; pseudoprimes to bases 2 and 3

(define range
  (case-lambda
    ((stop) (range 0 stop (if (negative? stop) -1 1)))
    ((start stop) (range start stop (if (< start stop) 1 -1)))
    ((start stop step)
      (let ((le? (if (negative? step) >= <=)))
        (let loop ((x start) (xs (list)))
          (if (le? stop x) (reverse xs)
            (loop (+ x step) (cons x xs))))))
    (else (error 'range "too many arguments"))))

(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 (prime? n) ; baillie-wagstaff
  (define ps (list 2 3 5 7 11 13 17 19 23 29 31
    37 41 43 47 53 59 61 67 71 73 79 83 89 97))
  (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 (square? n)
    (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)))))
    (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 (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 (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)))))))))(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)))
        (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 (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))))))))))))))
  (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))))))

(define (psp23? n)
  (and (= (expm 2 (- n 1) n) 1)
       (= (expm 3 (- n 1) n) 1)
       (not (prime? n))))

(display (filter psp23? (range 2 10000)))