; day 280 (import (rnrs hashtables (6))) (define-syntax define-generator (lambda (x) (syntax-case x (lambda) ((stx name (lambda formals e0 e1 ...)) (with-syntax ((yield (datum->syntax (syntax stx) 'yield))) (syntax (define name (lambda formals (let ((resume #f) (return #f)) (define yield (lambda args (call-with-current-continuation (lambda (cont) (set! resume cont) (apply return args))))) (lambda () (call-with-current-continuation (lambda (cont) (set! return cont) (cond (resume (resume)) (else (let () e0 e1 ...) (error 'name "unexpected return")))))))))))) ((stx (name . formals) e0 e1 ...) (syntax (stx name (lambda formals e0 e1 ...))))))) (define-generator (primegen . args) (define (oddify n) (if (odd? n) n (+ n 1))) (let ((start (if (pair? args) (car args) 0))) (when (<= start 2) (yield 2)) (when (<= start 3) (yield 3)) (let* ((pgen (primegen)) (p (and (pgen) (pgen))) (q (* p p)) (d (make-hashtable (lambda (x) x) =))) (define (add m s) (do () ((not (hashtable-contains? d m))) (set! m (+ m s))) (hashtable-set! d m s)) (do () ((< start q)) (let ((x (* (quotient start p) p))) (when (< x start) (set! x (+ x p))) (when (even? x) (set! x (+ x p))) (add x (+ p p)) (set! p (pgen)) (set! q (* p p)))) (do ((c (+ (oddify (max (- start 2) 3)) 2) (+ c 2))) (#f) (cond ((hashtable-contains? d c) (let ((s (hashtable-ref d c #f))) (hashtable-delete! d c) (add (+ c s) s))) ((< c q) (yield c)) (else (add (+ c p p) (+ p p)) (set! p (pgen)) (set! q (* p p)))))))) (define (prime? n) (define (expm b e m) (define (times p q) (modulo (* p q) m)) (let loop ((b b) (e e) (x 1)) (if (zero? e) x (loop (times b b) (quotient e 2) (if (odd? e) (times b x) x))))) (define (inverse x m) (let loop ((x x) (b m) (a 0) (u 1)) (if (zero? x) (if (= b 1) (modulo a m) (error 'inverse "must be coprime")) (let* ((q (quotient b x))) (loop (modulo b x) x u (- a (* u q))))))) (define (jacobi a m) (if (not (integer? a)) (error 'jacobi "must be integer") (if (not (and (integer? m) (positive? m) (odd? m))) (error 'jacobi "modulus must be odd positive integer") (let loop1 ((a (modulo a m)) (m m) (t 1)) (if (zero? a) (if (= m 1) t 0) (let ((z (if (member (modulo m 8) (list 3 5)) -1 1))) (let loop2 ((a a) (t t)) (if (even? a) (loop2 (/ a 2) (* t z)) (loop1 (modulo m a) a (if (and (= (modulo a 4) 3) (= (modulo m 4) 3)) (- t) t)))))))))) (define (isqrt n) (if (not (and (positive? n) (integer? n))) (error 'isqrt "must be positive integer") (let loop ((x n)) (let ((y (quotient (+ x (quotient n x)) 2))) (if (< y x) (loop y) x))))) (define (square? n) (let ((m (modulo n 128))) (if (positive? (bitwise-and (* m #x8bc40d7d) (* m #xa1e2f5d1) #x14020a)) #f (let ((large-mod (modulo n 3989930175))) ; (* 63 25 11 17 19 23 31) (and (let ((m (modulo large-mod 63))) (zero? (bitwise-and (* m #x3d491df7) (* m #xc824a9f9) #x10f14008))) (let ((m (modulo large-mod 25))) (zero? (bitwise-and (* m #x1929fc1b) (* m #x4c9ea3b2) #x51001005))) (let ((m (* #xd10d829a (modulo large-mod 31)))) (zero? (bitwise-and m (+ m #x672a5354) #x21025115))) (let ((m (modulo large-mod 23))) (zero? (bitwise-and (* m #x7bd28629) (* m #xe7180889) #xf8300))) (let ((m (modulo large-mod 19))) (zero? (bitwise-and (* m #x1b8bead3) (* m #x4d75a124) #x4280082b))) (let ((m (modulo large-mod 17))) (zero? (bitwise-and (* m #x6736f323) (* m #x9b1d499) #xc0000300))) (let ((m (modulo large-mod 11))) (zero? (bitwise-and (* m #xabf1a3a7) (* m #x2612bf93) #x45854000))) (let ((root (isqrt n))) (if (= (* root root) n) root #f))))))) (define (wheel? n limit) (let ((ws (vector 1 2 2 4 2 4 2 4 6 2 6))) (let loop ((f 2) (w 0)) (if (or (< n (* f f)) (< limit f)) #t (if (zero? (modulo n f)) #f (loop (+ f (vector-ref ws w)) (if (= w 10) 3 (+ w 1)))))))) (define (r-and-s n) (let loop ((r 0) (s (- n 1))) (if (odd? s) (values r s) (loop (+ r 1) (/ s 2))))) (define (strong-pseudoprime? n r s a) (if (= (expm a s n) 1) #t (let loop ((r r) (s s)) (if (zero? r) #f (if (= (expm a s n) (- n 1)) #t (loop (- r 1) (* s 2))))))) (define (miller? n . xs) (call-with-values (lambda () (r-and-s n)) (lambda (r s) (let loop ((xs xs)) (if (null? xs) #t (if (not (strong-pseudoprime? n r s (car xs))) #f (loop (cdr xs)))))))) (define (chain m f g x0 x1) (let loop ((ms (digits m 2)) (u x0) (v x1)) (if (null? ms) (values u v) (if (zero? (car ms)) (loop (cdr ms) (f u) (g u v)) (loop (cdr ms) (g u v) (f v)))))) (define (lucas? n) (let loop ((a 11) (b 7)) (let ((d (- (* a a) (* 4 b)))) (if (square? d) (loop (+ a 2) (+ b 1)) (if (not (= (gcd n (* 2 a b d)) 1)) (loop (+ a 2) (+ b 2)) (let* ((x (modulo (- (* a a (inverse b n)) 2) n)) (m (quotient (- n (jacobi d n)) 2)) (f (lambda (u) (modulo (- (* u u) 2) n))) (g (lambda (u v) (modulo (- (* u v) x) n)))) (call-with-values (lambda () (chain m f g 2 x)) (lambda (xm xm1) (zero? (modulo (- (* x xm) (* 2 xm1)) n)))))))))) (if (< n 2) #f (if (< n 1000000) (wheel? n 1000) (and (wheel? n 1000) (miller? n 2 3) (lucas? n))))) (define (sum-mod-prime-day limit) (let ((pgen (primegen))) (let loop ((sum (pgen)) (mod 1) (zs (list))) (if (< limit mod) (reverse zs) (if (prime? (modulo sum mod)) (loop (+ sum (pgen)) (+ mod 1) (cons mod zs)) (loop (+ sum (pgen)) (+ mod 1) zs)))))) (display (sum-mod-prime-day 365)) (newline) (display (length (sum-mod-prime-day 365))) (newline)