fork(1) download
copy 
  1. ; pseudoprimes to bases 2 and 3
  2.  
  3. (define range
  4.   (case-lambda
  5.     ((stop) (range 0 stop (if (negative? stop) -1 1)))
  6.     ((start stop) (range start stop (if (< start stop) 1 -1)))
  7.     ((start stop step)
  8.       (let ((le? (if (negative? step) >= <=)))
  9.         (let loop ((x start) (xs (list)))
  10.           (if (le? stop x) (reverse xs)
  11.             (loop (+ x step) (cons x xs))))))
  12.     (else (error 'range "too many arguments"))))
  13.  
  14. (define (expm b e m)
  15.   (define (m* x y) (modulo (* x y) m))
  16.   (cond ((zero? e) 1)
  17.         ((even? e) (expm (m* b b) (/ e 2) m))
  18.         (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
  19.  
  20. (define (prime? n) ; baillie-wagstaff
  21.   (define ps (list 2 3 5 7 11 13 17 19 23 29 31
  22.     37 41 43 47 53 59 61 67 71 73 79 83 89 97))
  23.   (define (expm b e m)
  24.     (define (m* x y) (modulo (* x y) m))
  25.     (cond ((zero? e) 1)
  26.           ((even? e) (expm (m* b b) (/ e 2) m))
  27.           (else (m* b (expm (m* b b) (/ (- e 1) 2) m))))) 
  28.   (define (square? n)
  29.     (define (isqrt n)
  30.       (if (not (and (positive? n) (integer? n)))
  31.           (error 'isqrt "must be positive integer")
  32.           (let loop ((x n))
  33.             (let ((y (quotient (+ x (quotient n x)) 2)))
  34.               (if (< y x) (loop y) x)))))
  35.     (let ((m (modulo n 128)))
  36.       (if (positive? (bitwise-and (* m #x8bc40d7d) (* m #xa1e2f5d1) #x14020a)) #f
  37.         (let ((large-mod (modulo n 3989930175))) ; (* 63 25 11 17 19 23 31)
  38.           (and (let ((m (modulo large-mod 63)))
  39.                  (zero? (bitwise-and (* m #x3d491df7) (* m #xc824a9f9) #x10f14008)))
  40.                (let ((m (modulo large-mod 25)))
  41.                  (zero? (bitwise-and (* m #x1929fc1b) (* m #x4c9ea3b2) #x51001005)))
  42.                (let ((m (* #xd10d829a (modulo large-mod 31))))
  43.                  (zero? (bitwise-and m (+ m #x672a5354) #x21025115)))
  44.                (let ((m (modulo large-mod 23)))
  45.                  (zero? (bitwise-and (* m #x7bd28629) (* m #xe7180889) #xf8300)))
  46.                (let ((m (modulo large-mod 19)))
  47.                  (zero? (bitwise-and (* m #x1b8bead3) (* m #x4d75a124) #x4280082b)))
  48.                (let ((m (modulo large-mod 17)))
  49.                  (zero? (bitwise-and (* m #x6736f323) (* m #x9b1d499) #xc0000300)))
  50.                (let ((m (modulo large-mod 11)))
  51.                  (zero? (bitwise-and (* m #xabf1a3a7) (* m #x2612bf93) #x45854000)))
  52.                (let ((root (isqrt n))) (if (= (* root root) n) root #f)))))))
  53.   (define (jacobi a m)
  54.     (if (not (integer? a)) (error 'jacobi "must be integer")
  55.       (if (not (and (integer? m) (positive? m) (odd? m)))
  56.           (error 'jacobi "modulus must be odd positive integer")
  57.           (let loop1 ((a (modulo a m)) (m m) (t 1))
  58.             (if (zero? a) (if (= m 1) t 0)
  59.               (let ((z (if (member (modulo m 8) (list 3 5)) -1 1)))
  60.                 (let loop2 ((a a) (t t))
  61.                   (if (even? a) (loop2 (/ a 2) (* t z))
  62.                     (loop1 (modulo m a) a
  63.                            (if (and (= (modulo a 4) 3)
  64.                                     (= (modulo m 4) 3))
  65.                                (- t) t))))))))))
  66.   (define (strong-pseudoprime? n a)
  67.     (let loop ((r 0) (s (- n 1)))
  68.       (if (even? s) (loop (+ r 1) (/ s 2))
  69.         (if (= (expm a s n) 1) #t
  70.           (let loop ((r r) (s s))
  71.             (cond ((zero? r) #f)
  72.                   ((= (expm a s n) (- n 1)) #t)
  73.                   (else (loop (- r 1) (* s 2)))))))))(define (selfridge n)
  74.     (let loop ((d-abs 5) (sign 1))
  75.       (let ((d (* d-abs sign)))
  76.         (cond ((< 1 (gcd d n)) (values d 0 0))
  77.               ((= (jacobi d n) -1) (values d 1 (/ (- 1 d) 4)))
  78.               (else (loop (+ d-abs 2) (- sign)))))))
  79.   (define (lucas p q m n) ; right-to-left
  80.     (define (even e o) (if (even? n) e o))
  81.     (define (mod n) (if (zero? m) n (modulo n m)))
  82.     (let ((d (- (* p p) (* 4 q))))
  83.       (let loop ((un 1) (vn p) (qn q) (n (quotient n 2))
  84.                  (u (even 0 1)) (v (even 2 p)) (k (even 1 q)))
  85.         (if (zero? n) (values u v k)
  86.           (let ((u2 (mod (* un vn))) (v2 (mod (- (* vn vn) (* 2 qn))))
  87.                 (q2 (mod (* qn qn))) (n2 (quotient n 2)))
  88.             (if (even? n) (loop u2 v2 q2 n2 u v k)
  89.               (let* ((uu (+ (* u v2) (* u2 v)))
  90.                      (vv (+ (* v v2) (* d u u2)))
  91.                      (uu (if (and (positive? m) (odd? uu)) (+ uu m) uu))
  92.                      (vv (if (and (positive? m) (odd? vv)) (+ vv m) vv))
  93.                      (uu (mod (/ uu 2))) (vv (mod (/ vv 2))))
  94.                 (loop u2 v2 q2 n2 uu vv (* k q2)))))))))
  95.   (define (powers-of-two n)
  96.     (let loop ((s 0) (n n))
  97.       (if (odd? n) (values s n)
  98.         (loop (+ s 1) (/ n 2)))))
  99.   (define (strong-lucas-pseudoprime? n)
  100.     ; assumes odd positive integer not a square
  101.     (call-with-values
  102.       (lambda () (selfridge n))
  103.       (lambda (d p q)
  104.         (if (zero? p) (= n d)
  105.           (call-with-values
  106.             (lambda () (powers-of-two (+ n 1)))
  107.             (lambda (s t)
  108.               (call-with-values
  109.                 (lambda () (lucas p q n t))
  110.                 (lambda (u v k)
  111.                   (if (or (zero? u) (zero? v)) #t
  112.                     (let loop ((r 1) (v v) (k k))
  113.                       (if (= r s) #f
  114.                         (let* ((v (modulo (- (* v v) (* 2 k)) n))
  115.                                (k (modulo (* k k) n)))
  116.                           (if (zero? v) #t (loop (+ r 1) v k))))))))))))))
  117.   (if (not (integer? n)) (error 'prime? "must be integer"))
  118.   (if (or (< n 2) (square? n)) #f
  119.     (let loop ((ps ps))
  120.       (if (pair? ps)
  121.           (if (zero? (modulo n (car ps))) (= n (car ps)) (loop (cdr ps)))
  122.           (and (strong-pseudoprime? n 2)
  123.                (strong-pseudoprime? n 3)
  124.                (strong-lucas-pseudoprime? n))))))
  125.  
  126. (define (psp23? n)
  127.   (and (= (expm 2 (- n 1) n) 1)
  128.        (= (expm 3 (- n 1) n) 1)
  129.        (not (prime? n))))
  130.  
  131. (display (filter psp23? (range 2 10000)))
Success #stdin #stdout 0.86s 12224KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
(919 1105 1729 1999 2039 2465 2701 2711 2821 3319 3511 4591 4831 5471 5711 5839 6079 6199 6601 6679 6791 6959 6991 7039 7411 8311 8911 9151 9319 9391 9431 9439 9551 9631 9791 9839 9871)
https://ideone.com/chL78D
language:
Scheme (chicken 4.13)
created:
4 years ago
visibility:
secret

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!