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