; lucas pseudoprimes
(define-syntax fold-of
(syntax-rules (range in is)
((_ "z" f b e) (set! b (f b e)))
((_ "z" f b e (v range fst pst stp) c ...)
(let* ((x fst) (p pst) (s stp)
(le? (if (positive? s) <= >=)))
(do ((v x (+ v s))) ((le? p v) b)
(fold-of "z" f b e c ...))))
((_ "z" f b e (v range fst pst) c ...)
(let* ((x fst) (p pst) (s (if (< x p) 1 -1)))
(fold-of "z" f b e (v range x p s) c ...)))
((_ "z" f b e (v range pst) c ...)
(fold-of "z" f b e (v range 0 pst) c ...))
((_ "z" f b e (x in xs) c ...)
(do ((t xs (cdr t))) ((null? t) b)
(let ((x (car t)))
(fold-of "z" f b e c ...))))
((_ "z" f b e (x is y) c ...)
(let ((x y)) (fold-of "z" f b e c ...)))
((_ "z" f b e p? c ...)
(if p? (fold-of "z" f b e c ...)))
((_ f i e c ...)
(let ((b i)) (fold-of "z" f b e c ...)))))
(define-syntax list-of (syntax-rules ()
((_ arg ...) (reverse (fold-of
(lambda (d a) (cons a d)) '() arg ...)))))
(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 ((n2 (isqrt n)))
(= (* n2 n2) n)))
(define (expm b e m)
(define (m* x y) (modulo (* x y) m))
(cond ((zero? e) 1)
((even? e) (expm (m* b b) (/ e 2) m))
(else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
(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 (primes . args) ; (primes [lo] hi) inclusive at both ends
(let* ((lo (if (null? (cdr args)) 0 (car args)))
(hi (if (null? (cdr args)) (car args) (cadr args))))
(cond ((and (<= lo 100000) (<= hi 1000000)) ; simple sieve
(let* ((max-index (quotient (- hi 3) 2))
(v (make-vector (+ max-index 1) #t)))
(let loop ((i 0) (ps (list 2)))
(let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
(cond ((< hi (* p p))
(let loop ((j i) (ps ps))
(cond ((< max-index j)
(let loop ((ps (reverse ps)))
(if (<= lo (car ps)) ps
(loop (cdr ps)))))
((vector-ref v j)
(loop (+ j 1) (cons (+ j j 3) ps)))
(else (loop (+ j 1) ps)))))
((vector-ref v i)
(let loop ((j startj))
(when (<= j max-index)
(vector-set! v j #f) (loop (+ j p))))
(loop (+ i 1) (cons p ps)))
(else (loop (+ i 1) ps)))))))
((< lo (sqrt hi))
(let* ((r (inexact->exact (ceiling (sqrt hi))))
(r (if (even? r) r (+ r 1))))
(append (primes lo r) (primes r hi))))
(else ; segmented sieve
(let* ((lo (if (even? lo) lo (- lo 1)))
(b 50000) (bs (make-vector b #t))
(r (inexact->exact (ceiling (sqrt hi))))
(ps (cdr (primes r)))
(qs (map (lambda (p)
(modulo (* -1/2 (+ lo 1 p)) p)) ps))
(zs (list)) (z (lambda (p) (set! zs (cons p zs)))))
(do ((t lo (+ t b b))
(qs qs (map (lambda (p q) (modulo (- q b) p))
ps qs)))
((<= hi t)
(let loop ((zs zs))
(if (<= (car zs) hi) (reverse zs)
(loop (cdr zs)))))
(do ((i 0 (+ i 1))) ((= i b)) (vector-set! bs i #t))
(do ((ps ps (cdr ps)) (qs qs (cdr qs))) ((null? qs))
(do ((j (car qs) (+ j (car ps)))) ((<= b j))
(vector-set! bs j #f)))
(do ((j 0 (+ j 1))) ((= j b))
(if (vector-ref bs j) (z (+ t j j 1))))))))))
(define (strong-pseudoprime? n a)
(let loop ((r 0) (s (- n 1)))
(if (even? s) (loop (+ r 1) (/ s 2))
(if (= (expm a s n) 1) #t
(let loop ((r r) (s s))
(cond ((zero? r) #f)
((= (expm a s n) (- n 1)) #t)
(else (loop (- r 1) (* s 2)))))))))
(let ((primes487000 (primes 487000)))
(list-of n ; strong pseudoprimes to base 2 A001262
(n range 3 487000)
(strong-pseudoprime? n 2)
(not (member n primes487000))))
; (2047 3277 4033 4681 8321 15841 29341 42799 49141 52633 65281
; 74665 80581 85489 88357 90751 104653 130561 196093 220729
; 233017 252601 253241 256999 271951 280601 314821 357761
; 390937 458989 476971 486737)
(let ((primes227000 (primes 227000)))
(list-of n ; strong pseudoprimes to base 3 A020229
(n range 3 227000 2)
(strong-pseudoprime? n 3)
(not (member n primes227000))))
; (121 703 1891 3281 8401 8911 10585 12403 16531 18721 19345
; 23521 31621 44287 47197 55969 63139 74593 79003 82513 87913
; 88573 97567 105163 111361 112141 148417 152551 182527 188191
; 211411 218791 221761 226801)
(define (selfridge n)
(let loop ((d-abs 5) (sign 1))
(let ((d (* d-abs sign)))
(cond ((< 1 (gcd d n)) (values d 0 0))
((= (jacobi d n) -1) (values d 1 (/ (- 1 d) 4)))
(else (loop (+ d-abs 2) (- sign)))))))
(define (lucas p q m n) ; right-to-left
(define (even e o) (if (even? n) e o))
(define (mod n) (if (zero? m) n (modulo n m)))
(let ((d (- (* p p) (* 4 q))))
(let loop ((un 1) (vn p) (qn q) (n (quotient n 2))
(u (even 0 1)) (v (even 2 p)) (k (even 1 q)))
; (display un) (display " ") (display vn) (display " ")
; (display qn) (display " ") (display n) (display " ")
; (display u) (display " ") (display v) (display " ")
; (display k) (newline)
(if (zero? n) (values u v k)
(let ((u2 (mod (* un vn))) (v2 (mod (- (* vn vn) (* 2 qn))))
(q2 (mod (* qn qn))) (n2 (quotient n 2)))
(if (even? n) (loop u2 v2 q2 n2 u v k)
(let* ((uu (+ (* u v2) (* u2 v)))
(vv (+ (* v v2) (* d u u2)))
(uu (if (and (positive? m) (odd? uu)) (+ uu m) uu))
(vv (if (and (positive? m) (odd? vv)) (+ vv m) vv))
(uu (mod (/ uu 2))) (vv (mod (/ vv 2))))
(loop u2 v2 q2 n2 uu vv (* k q2)))))))))
(define (standard-lucas-pseudoprime? n)
; assumes odd positive integer not a square
(call-with-values
(lambda () (selfridge n))
(lambda (d p q)
(if (zero? p) (= n d)
(call-with-values
(lambda () (lucas p q n (+ n 1)))
(lambda (u v qkd) (zero? u)))))))
(define (powers-of-two n)
(let loop ((s 0) (n n))
(if (odd? n) (values s n)
(loop (+ s 1) (/ n 2)))))
(define (strong-lucas-pseudoprime? n)
; assumes odd positive integer not a square
(call-with-values
(lambda () (selfridge n))
(lambda (d p q)
(if (zero? p) (= n d)
(call-with-values
(lambda () (powers-of-two (+ n 1)))
(lambda (s t)
(call-with-values
(lambda () (lucas p q n t))
(lambda (u v k)
(if (or (zero? u) (zero? v)) #t
(let loop ((r 1) (v v) (k k))
(if (= r s) #f
(let* ((v (modulo (- (* v v) (* 2 k)) n))
(k (modulo (* k k) n)))
(if (zero? v) #t (loop (+ r 1) v k))))))))))))))
(let ((primes48000 (primes 48000)))
(list-of n ; standard lucas pseudoprimes A217120
(n range 3 48000 2)
(not (square? n))
(standard-lucas-pseudoprime? n)
(not (member n primes48000))))
; (323 377 1159 1829 3827 5459 5777 9071 9179 10877 11419 11663
; 13919 14839 16109 16211 18407 18971 19043 22499 23407 24569
; 25199 25877 26069 27323 32759 34943 35207 39059 39203 39689
; 40309 44099 46979 47879)
(let ((primes325000 (primes 325000)))
(list-of n ; strong lucas pseudoprimes A217255
(n range 3 325000 2)
(not (square? n))
(strong-lucas-pseudoprime? n)
(not (member n primes325000))))
; (5459 5777 10877 16109 18971 22499 24569 40309 58519 75077
; 97439 113573 115639 130139 155819 158399 161027 162133
; 176399 176471 189419 192509 197801 224369 230691 243629
; 253259 268349 288919 313499 324899)
(define prime?
(let ((ps (primes 100)))
(lambda (n)
(if (not (integer? n)) (error 'prime? "must be integer"))
(if (or (< n 2) (square? n)) #f
(let loop ((ps ps))
(if (pair? ps)
(if (zero? (modulo n (car ps))) (= n (car ps)) (loop (cdr ps)))
(and (strong-pseudoprime? n 2)
(strong-pseudoprime? n 3)
(strong-lucas-pseudoprime? n))))))))
(do ((xs (list #e1e3 #e1e6 #e1e9 #e1e12) (cdr xs))) ((null? xs))
(let loop ((x (car xs)) (ps (primes (car xs) (+ (car xs) 100000))))
(when (pair? ps)
(when (and (= x (car ps)) (not (prime? x)))
(display "prime reported as composite ") (display x) (newline))
(when (and (prime? x) (not (= x (car ps))))
(display "composite reported as prime ") (display x) (newline))
(loop (+ x 1) (if (= x (car ps)) (cdr ps) ps)))))
(with-input-from-file "23primes" (lambda ()
(do ((x (read) (read))) ((eof-object? x))
(when (prime? x) (display x) (newline)))))