codepad

Language:

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

Private [?]

Run code