; 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)