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