fork download
copy 
  1. ; day 280
  2.  
  3. (import (rnrs hashtables (6)))
  4.  
  5. (define-syntax define-generator
  6.   (lambda (x)
  7.     (syntax-case x (lambda)
  8.       ((stx name (lambda formals e0 e1 ...))
  9.          (with-syntax ((yield (datum->syntax (syntax stx) 'yield)))
  10.            (syntax (define name
  11.              (lambda formals
  12.                (let ((resume #f) (return #f))
  13.                  (define yield
  14.                    (lambda args
  15.                      (call-with-current-continuation
  16.                       (lambda (cont)
  17.                         (set! resume cont)
  18.                         (apply return args)))))
  19.                  (lambda ()
  20.                    (call-with-current-continuation
  21.                     (lambda (cont)
  22.                       (set! return cont)
  23.                       (cond (resume (resume))
  24.                       (else (let () e0 e1 ...)
  25.                             (error 'name "unexpected return"))))))))))))
  26.         ((stx (name . formals) e0 e1 ...)
  27.           (syntax (stx name (lambda formals e0 e1 ...)))))))
  28.  
  29. (define-generator (primegen . args)
  30.   (define (oddify n) (if (odd? n) n (+ n 1)))
  31.   (let ((start (if (pair? args) (car args) 0)))
  32.     (when (<= start 2) (yield 2))
  33.     (when (<= start 3) (yield 3))
  34.     (let* ((pgen (primegen)) (p (and (pgen) (pgen))) (q (* p p))
  35.             (d (make-hashtable (lambda (x) x) =)))
  36.     (define (add m s)
  37.       (do () ((not (hashtable-contains? d m))) (set! m (+ m s)))
  38.       (hashtable-set! d m s))
  39.     (do () ((< start q))
  40.       (let ((x (* (quotient start p) p)))
  41.         (when (< x start) (set! x (+ x p)))
  42.         (when (even? x) (set! x (+ x p)))
  43.         (add x (+ p p)) (set! p (pgen)) (set! q (* p p))))
  44.     (do ((c (+ (oddify (max (- start 2) 3)) 2) (+ c 2))) (#f)
  45.       (cond ((hashtable-contains? d c)
  46.               (let ((s (hashtable-ref d c #f)))
  47.                 (hashtable-delete! d c)
  48.                 (add (+ c s) s)))
  49.             ((< c q) (yield c))
  50.             (else (add (+ c p p) (+ p p))
  51.                   (set! p (pgen))
  52.                   (set! q (* p p))))))))
  53.  
  54. (define (prime? n)
  55.   (define (expm b e m)
  56.     (define (times p q) (modulo (* p q) m))
  57.     (let loop ((b b) (e e) (x 1))
  58.       (if (zero? e) x
  59.         (loop (times b b) (quotient e 2) (if (odd? e) (times b x) x)))))
  60.   (define (inverse x m)
  61.     (let loop ((x x) (b m) (a 0) (u 1))
  62.       (if (zero? x)
  63.           (if (= b 1) (modulo a m) (error 'inverse "must be coprime"))
  64.           (let* ((q (quotient b x))) (loop (modulo b x) x u (- a (* u q)))))))
  65.   (define (jacobi a m)
  66.     (if (not (integer? a)) (error 'jacobi "must be integer")
  67.       (if (not (and (integer? m) (positive? m) (odd? m)))
  68.           (error 'jacobi "modulus must be odd positive integer")
  69.           (let loop1 ((a (modulo a m)) (m m) (t 1))
  70.             (if (zero? a) (if (= m 1) t 0)
  71.               (let ((z (if (member (modulo m 8) (list 3 5)) -1 1)))
  72.                 (let loop2 ((a a) (t t))
  73.                   (if (even? a) (loop2 (/ a 2) (* t z))
  74.                     (loop1 (modulo m a) a
  75.                            (if (and (= (modulo a 4) 3) (= (modulo m 4) 3)) (- t)
  76.                              t))))))))))
  77.   (define (isqrt n)
  78.     (if (not (and (positive? n) (integer? n)))
  79.         (error 'isqrt "must be positive integer")
  80.         (let loop ((x n))
  81.           (let ((y (quotient (+ x (quotient n x)) 2)))
  82.             (if (< y x) (loop y) x)))))
  83.   (define (square? n)
  84.     (let ((m (modulo n 128)))
  85.       (if (positive? (bitwise-and (* m #x8bc40d7d) (* m #xa1e2f5d1) #x14020a)) #f
  86.         (let ((large-mod (modulo n 3989930175))) ; (* 63 25 11 17 19 23 31)
  87.           (and (let ((m (modulo large-mod 63)))
  88.                  (zero? (bitwise-and (* m #x3d491df7) (* m #xc824a9f9) #x10f14008)))
  89.                (let ((m (modulo large-mod 25)))
  90.                  (zero? (bitwise-and (* m #x1929fc1b) (* m #x4c9ea3b2) #x51001005)))
  91.                (let ((m (* #xd10d829a (modulo large-mod 31))))
  92.                  (zero? (bitwise-and m (+ m #x672a5354) #x21025115)))
  93.                (let ((m (modulo large-mod 23)))
  94.                  (zero? (bitwise-and (* m #x7bd28629) (* m #xe7180889) #xf8300)))
  95.                (let ((m (modulo large-mod 19)))
  96.                  (zero? (bitwise-and (* m #x1b8bead3) (* m #x4d75a124) #x4280082b)))
  97.                (let ((m (modulo large-mod 17)))
  98.                  (zero? (bitwise-and (* m #x6736f323) (* m #x9b1d499) #xc0000300)))
  99.                (let ((m (modulo large-mod 11)))
  100.                  (zero? (bitwise-and (* m #xabf1a3a7) (* m #x2612bf93) #x45854000)))
  101.                (let ((root (isqrt n))) (if (= (* root root) n) root #f)))))))
  102.   (define (wheel? n limit)
  103.     (let ((ws (vector 1 2 2 4 2 4 2 4 6 2 6)))
  104.       (let loop ((f 2) (w 0))
  105.         (if (or (< n (* f f)) (< limit f)) #t
  106.           (if (zero? (modulo n f)) #f
  107.             (loop (+ f (vector-ref ws w)) (if (= w 10) 3 (+ w 1))))))))
  108.   (define (r-and-s n)
  109.     (let loop ((r 0) (s (- n 1)))
  110.       (if (odd? s) (values r s)
  111.         (loop (+ r 1) (/ s 2)))))
  112.   (define (strong-pseudoprime? n r s a)
  113.     (if (= (expm a s n) 1) #t
  114.       (let loop ((r r) (s s))
  115.         (if (zero? r) #f
  116.           (if (= (expm a s n) (- n 1)) #t
  117.             (loop (- r 1) (* s 2)))))))
  118.   (define (miller? n . xs)
  119.     (call-with-values
  120.       (lambda () (r-and-s n))
  121.       (lambda (r s)
  122.         (let loop ((xs xs))
  123.           (if (null? xs) #t
  124.             (if (not (strong-pseudoprime? n r s (car xs))) #f
  125.               (loop (cdr xs))))))))
  126.   (define (chain m f g x0 x1)
  127.     (let loop ((ms (digits m 2)) (u x0) (v x1))
  128.       (if (null? ms) (values u v)
  129.         (if (zero? (car ms)) (loop (cdr ms) (f u) (g u v))
  130.           (loop (cdr ms) (g u v) (f v))))))
  131.   (define (lucas? n)
  132.     (let loop ((a 11) (b 7))
  133.       (let ((d (- (* a a) (* 4 b))))
  134.         (if (square? d) (loop (+ a 2) (+ b 1))
  135.           (if (not (= (gcd n (* 2 a b d)) 1)) (loop (+ a 2) (+ b 2))
  136.             (let* ((x (modulo (- (* a a (inverse b n)) 2) n))
  137.                    (m (quotient (- n (jacobi d n)) 2))
  138.                    (f (lambda (u) (modulo (- (* u u) 2) n)))
  139.                    (g (lambda (u v) (modulo (- (* u v) x) n))))
  140.               (call-with-values
  141.                 (lambda () (chain m f g 2 x))
  142.                 (lambda (xm xm1) (zero? (modulo (- (* x xm) (* 2 xm1)) n))))))))))
  143.   (if (< n 2) #f (if (< n 1000000) (wheel? n 1000)
  144.     (and (wheel? n 1000) (miller? n 2 3) (lucas? n)))))
  145.  
  146. (define (sum-mod-prime-day limit)
  147.   (let ((pgen (primegen)))
  148.     (let loop ((sum (pgen)) (mod 1) (zs (list)))
  149.       (if (< limit mod) (reverse zs)
  150.         (if (prime? (modulo sum mod))
  151.             (loop (+ sum (pgen)) (+ mod 1) (cons mod zs))
  152.             (loop (+ sum (pgen)) (+ mod 1) zs))))))
  153.  
  154. (display (sum-mod-prime-day 365)) (newline)
  155. (display (length (sum-mod-prime-day 365))) (newline)
Success #stdin #stdout 0.29s 48312KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
(5 6 7 8 12 15 16 19 20 21 24 26 30 34 37 38 40 42 44 45 46 48 49 50 55 58 59 60 62 64 65 66 67 68 70 72 73 75 76 78 86 87 88 92 102 116 120 122 124 128 130 132 135 140 143 145 150 156 158 164 165 166 168 172 173 175 176 182 183 191 196 210 214 216 218 223 234 236 241 248 250 256 259 262 265 266 272 280 285 301 306 310 311 314 315 324 328 330 336 337 344 347 348 349 352 355 358 365)
108
https://ideone.com/TyDKiX
language:
Scheme (guile 2.2.4)
created:
6 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!