(define (expm b e m)
(define (times x y) (modulo (* x y) m))
(let loop ((b b) (e e) (r 1))
(if (zero? e) r
(loop (times b b) (quotient e 2)
(if (odd? e) (times b r) r)))))
(define (iroot k n)
(let loop ((u n))
(let* ((s u)
(t (+ (* (- k 1) s)
(quotient n (expt s (- k 1)))))
(u (quotient t k)))
(if (< u s) (loop u) s))))
(define (digits n . args)
(let ((b (if (null? args) 10 (car args))))
(let loop ((n n) (d '()))
(if (zero? n) d
(loop (quotient n b)
(cons (modulo n b) d))))))
(define square?
(let ((q11 (make-vector 11 #f))
(q63 (make-vector 63 #f))
(q64 (make-vector 64 #f))
(q65 (make-vector 65 #f)))
(do ((k 0 (+ k 1))) ((< 5 k))
(vector-set! q11 (modulo (* k k) 11) #t))
(do ((k 0 (+ k 1))) ((< 31 k))
(vector-set! q63 (modulo (* k k) 63) #t))
(do ((k 0 (+ k 1))) ((< 31 k))
(vector-set! q64 (modulo (* k k) 64) #t))
(do ((k 0 (+ k 1))) ((< 32 k))
(vector-set! q65 (modulo (* k k) 65) #t))
(lambda (n)
(if (not (integer? n)) (error 'square? "must be integer")
(if (< n 1) #f
(if (not (vector-ref q64 (modulo n 64))) #f
(let ((r (modulo n 45045)))
(if (not (vector-ref q63 (modulo r 63))) #f
(if (not (vector-ref q65 (modulo r 65))) #f
(if (not (vector-ref q11 (modulo r 11))) #f
(let ((q (iroot 2 n)))
(if (= (* q q) n) q #f))))))))))))
(define (spsp? n a)
(do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
((odd? d)
(let ((t (expm a d n)))
(if (or (= t 1) (= t (- n 1))) #t
(do ((s (- s 1) (- s 1))
(t (expm t 2 n) (expm t 2 n)))
((or (zero? s) (= t (- n 1)))
(positive? s))))))))
(define (pseudoprime? n)
(and (spsp? n 2) (spsp? n 3) (spsp? n 5)))
(define (td n limit)
(let twos ((n n) (fs (list)))
(cond ((= n 1) fs)
((even? n) (twos (/ n 2) (cons 2 fs)))
(else (let odds ((n n) (f 3) (fs fs))
(cond ((< limit f) fs)
((< n (* f f)) (reverse (cons n fs)))
((zero? (modulo n f))
(odds (/ n f) f (cons f fs)))
(else (odds n (+ f 2) fs))))))))
(define (rho n limit)
(let loop ((t 2) (h 2) (d 1) (c 1) (limit limit))
(define (f x) (modulo (+ (* x x) c) n))
(cond ((zero? limit) #f)
((= d 1) (let ((t (f t)) (h (f (f h))))
(loop t h (gcd (- t h) n) c (- limit 1))))
((= d n) (loop 2 2 1 (+ c 1) (- limit 1)))
((pocklington? d) d)
(else (rho d (- limit 1))))))
(define (pocklington? n)
(if (not (pseudoprime? n)) #f
(if (< n #e1e12) (= (car (td n #e1e6)) n)
(let* ((n-1 (- n 1))
(root2 (iroot 2 n)) (root3 (iroot 3 n))
(fs (td n-1 1000)) (f (apply * fs)))
(let f-loop ((f f) (fs fs))
(if (< f root3)
(let ((p (rho (/ n-1 f) #e1e6)))
(if p (f-loop (* f p) (cons p fs)) 'unknown))
(let a-loop ((a 2))
(if (not (= (expm a n-1 n) 1)) #f
(let loop ((fs fs))
(let ((g (gcd (- (expm a (/ n-1 (car fs)) n) 1) n)))
(cond ((< 1 g n) #f)
((= g n) (a-loop (+ a 1)))
((pair? (cdr fs)) (loop (cdr fs)))
((< root2 f) #t)
(else (let* ((ds (digits n f))
(c2 (car ds)) (c1 (cadr ds))
(s (- (* c1 c1) (* 4 c2))))
(not (square? s)))))))))))))))
(define (mersenne n) (- (expt 2 n) 1))
(display (pocklington? (mersenne 89))) (newline)
(display (pocklington? (mersenne 521))) (newline)
(display (pseudoprime? 3825123056546413051)) (newline)
(display (pocklington? 3825123056546413051)) (newline)