; generating large random primes
(define (drop-while pred? xs)
(let loop ((xs xs))
(if (or (null? xs) (not (pred? (car xs)))) xs
(loop (cdr xs)))))
(define (fortune xs)
(let loop ((n 1) (x #f) (xs xs))
(cond ((null? xs) x)
((< (rand) (/ n))
(loop (+ n 1) (car xs) (cdr xs)))
(else (loop (+ n 1) x (cdr xs))))))
(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 rand
(let ((seed 0.3141592654))
(lambda args
(set! seed
(if (pair? args)
(sin (car args))
(let ((x (* seed 147.0)))
(- x (floor x)))))
seed)))
(define (randint first past)
(inexact->exact (floor (+ (* (- past first) (rand)) first))))
(define (primes n) ; list of primes not exceeding n
(let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
(let loop ((i 0) (p 3) (ps (list 2)))
(cond ((< n (* p p))
(do ((i i (+ i 1)) (p p (+ p 2))
(ps ps (if (vector-ref bits i) (cons p ps) ps)))
((= i len) (reverse ps))))
((vector-ref bits i)
(do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
(vector-set! bits j #f)))
(else (loop (+ i 1) (+ p 2) ps))))))
(define prime? ; strong pseudoprime to prime bases less than 100
(let* ((ps (primes 100)) (p100 (apply * ps)))
(lambda (n)
(define (expm b e m)
(let loop ((b b) (e e) (x 1))
(if (zero? e) x
(loop (modulo (* b b) m) (quotient e 2)
(if (odd? e) (modulo (* b x) m) x)))))
(define (spsp? n a) ; #t if n is a strong pseudoprime base a
(do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
((odd? d) (if (= (expm a d n) 1) #t
(do ((r 0 (+ r 1)))
((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
(< r s)))))))
(if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f)
(do ((ps ps (cdr ps)))
((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))
(define (pocklington p k a)
(and (odd? p) (prime? p)
(< 1 k (* 2 (+ p 1)))
(positive? (modulo k p))
(let ((n (+ (* 2 k p) 1)))
(and (< 1 a n)
(= (expm a (/ (- n 1) 2) n) (- n 1))
(= (gcd (+ (expm a k n) 1) n) 1)
n)))) ; return n if prime, else #f
(define (ratchet p lo hi verbose?)
(define (pock p k a n)
(and (= (expm a (/ (- n 1) 2) n) (- n 1))
(= (gcd (+ (expm a k n) 1) n) 1)))
(let loop ((k (randint lo hi)))
(let ((n (+ (* 2 k p) 1)))
(cond ((pock p k 2 n) (when verbose? (display (list p k 2)) (newline)) n)
((pock p k 3 n) (when verbose? (display (list p k 3)) (newline)) n)
(else (loop (randint lo hi)))))))
(define (rand-prime lo hi . args)
(let ((verbose? (if (pair? args) (car args) #f)))
(if (< hi 10000) (fortune (drop-while (lambda (x) (< x lo)) (primes hi)))
(let loop ((p (fortune (cdr (primes 10000)))))
(if (<= lo p) p
(let ((k-lo (ceiling (/ (- lo 1) 2 p)))
(k-hi (floor (/ (- hi 1) 2 p))))
(loop (ratchet p (if (< (* 2 p) k-lo) 1 k-lo)
(min (* 2 p) k-hi) verbose?))))))))
(display (pocklington 2333 2001 2)) (newline)
(display (pocklington 2333 2017 2)) (newline)
(display (rand-prime #e1e49 #e1e50 #t)) (newline)