Programming Praxis

We begin with the functions that implement the four conditions. Two tests are sufficient for the four conditions, because Condition 1 and Condition 3 are essentially the same test, and likewise Condition 2 and Condition 4. For Condition 3 and Condition 4, the additional r1 or r2 acts as just another “prime” in the test of Condition 1 and Condition 2. Both functions are straight forward, but the Lucas test is ugly because of the multiple return values:

(define (fermat-test? n ps)
  (let loop ((ps ps) (a 1))
    (if (null? ps) #t
      (if (and (= (expm a (- n 1) n) 1)
               (= (gcd (expt a (/ (- n 1) (car ps))) n) 1))
          (loop (cdr ps) 1)
          (loop ps (+ a 1))))))

(define (lucas-test? n ps)
  (call-with-values
    (lambda () (selfridge n))
    (lambda (d pp qq)
      (let loop ((ps ps) (p pp) (q qq))
        (if (null? ps) #t
          (call-with-values
            (lambda () (lucas p q n (+ n 1)))
            (lambda (un+1 vn+1 qn+1)
              (call-with-values
                (lambda () (lucas p q n (/ (+ n 1) (car ps))))
                (lambda (un+1/p vn+1/p qn+1/p)
                  (if (and (zero? un+1) (= (gcd un+1/p n) 1))
                      (loop (cdr ps) pp qq)
                      (loop ps (+ p 2) (+ p q 1))))))))))))

After a primality test to eliminate composites and removing factors of two, the outer loop runs a 2,3,5,7-wheel to find factors less than b1, stopping if it finds enough. If not, the inner loop runs Pollard rho tests up to b2 iterations on both the n−1 and n+1 cofactors, again stopping if it finds enough. Finally, the function reports failure if not enough factors are found, or checks the appropriate set of conditions:

(define (bls-prime? n . args)
  (define (last-pair xs) (if (null? (cdr xs)) xs (last-pair (cdr xs))))
  (define (cycle xs) (set-cdr! (last-pair xs) xs) xs)
  (define wheel (list 2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8
    4 2 4 2 4 8 6 4 6 2 4 6 2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10))
  (define (rho n c b)
    (define (f y) (modulo (+ (* y y) c) n))
    (define (g p x y) (modulo (* p (abs (- x y))) n))
    (let stage1 ((x 2) (y (+ 4 c)) (z (+ 4 c)) (j 1) (q 2) (p 1))
      (if (= j b) #f ; timeout
        (if (= x y) (rho n (+ c 1) (- b j)) ; cycle
          (if (= j q) (let ((t (f y))) (stage1 y (f y) z (+ j 1) (* q 2) (g p y t)))
            (if (positive? (modulo j 100)) (stage1 x (f y) z (+ j 1) q (g p x y))
              (let ((d (gcd p n)))
                (if (= d 1) (stage1 x (f y) y (+ j 1) q (g p x y))
                  (if (and (< 1 d n) (bls-prime? d)) d ; factor
                    (let stage2 ((k 1) (z (f z)))
                      (if (= k 100) (rho n (+ c 1) (- b j))
                        (let ((d (gcd (- z x) n)))
                          (if (= d 1) (stage2 (+ k 1) (f z))
                            (if (= d n) (rho n (+ c 1) (- b j))
                              (if (bls-prime? d) d ; factor
                                (rho n (+ c 1) (- b j)))))))))))))))))
  (define (remove-twos n)
    (let loop ((f 1) (r n))
      (if (odd? r) (values f (list 2) r)
        (loop (* f 2) (/ r 2)))))
  (define (join x xs) (if (member x xs) xs (cons x xs)))
  (define (enough? n p f1 f2)
    (< n (* (max (+ (* p f1) 1) (- (* p f2) 1)) (+ (* p p f1 f2 1/2) 1))))
  (if (not (prime? n)) #f ; sanity check
    (let-values (((b1) (if (pair? args) (car args) 20000000))
                 ((b2) (if (and (pair? args) (pair? (cdr args)))
                           (cadr args) 10000000))
                 ((f1 f1s r1) (remove-twos (- n 1)))
                 ((f2 f2s r2) (remove-twos (+ n 1))))
      (let loop ((p 3) (ws (cons 2 (cons 2 (cons 4 (cycle wheel))))))
        (cond ((and (< p b1) (not (enough? n p f1 f2)))
              (case (modulo (+ n 1) p)
              ((0) (while (zero? (modulo r2 p))
                (set! f2 (* f2 p)) (set! f2s (join p f2s)) (set! r2 (/ r2 p))))
              ((2) (while (zero? (modulo r1 p))
                (set! f1 (* f1 p)) (set! f1s (join p f1s)) (set! r1 (/ r1 p)))))
              (when (< r1 (* p p))
                (set! f1 (- n 1)) (set! f1s (join r1 f1s)) (set! r1 1))
              (when (< r2 (* p p))
                (set! f2 (+ n 1)) (set! f2s (join r2 f2s)) (set! r2 1))
              (loop (+ p (car ws)) (cdr ws)))
        (else (set! p (min p b1))
              (when verbose?
                (display (list 'wheel n p f1 f1s r1 f2 f2s r2)) (newline))
              (let loop ((more1? #t) (more2? #t) (more3? #t))
                (cond ((and (not (enough? n p f1 f2)) (or more1? more2?))
                        (if more1?
                          (let ((f (rho r1 1 b2)))
                              (if (not f) (loop #f #t #t)
                                (begin (set! f1 (* f1 f)) (set! f1s (join f f1s))
                                         (set! r1 (/ r1 f)) (loop #t #t#t))))
                              (let ((f (rho r2 1 b2)))
                                (if (not f) (loop #f #f #t)
                                  (begin (set! f2 (* f2 f)) (set! f2s (join f f2s))
                                         (set! r2 (/ r2 f)) (loop #f #t #t))))))
                      (more3?
                        (when verbose?
                          (display (list 'rho n p f1 f1s r1 f2 f2s r2)) (newline))
                        (loop #f #f #f))
                      ((not (enough? n p f1 f2)) #f) ; failure to prove
                      ((< r1 f1)
                        (when verbose? (display "fermat only") (newline))
                        (fermat-test? n f1s)) ; condition 1 only
                      ((< r2 f2)
                        (when verbose? (display "lucas only") (newline))
                        (lucas-test? n f2s)) ; condition 2 only
                      (else
                        (when verbose? (display "fermat and lucas") (newline))
                        (and (fermat-test? n (cons r1 f1s)) ; 1 and 3
                             (lucas-test? n (cons r2 f2s)))))))))))) ; 2&4

The function takes n and, optionally, the two bounds, so you can adjust the limits depending on the size of n; the default limits allow most primes to sixty digits to be proved quickly, and sometimes much larger primes. Here’s an example, which first performs trial division then finds one factor of n−1 and two factors of n+1 using Pollard rho to prove a 75-digit prime in under half a minute:

> (time (bls-prime? 851486813098510619326357783811917873528153342227994200571828
674136369388533))
(wheel 8514868130985106193263577838119178735281533422279942005718286741363693885
33 20000000 27473484 (34171 67 3 2) 30993040893485173534101382402461874639858320
926024314956626129912623 331524754 (31859 43 11 2) 25683958824338969853405208583
96950574400104500878369403968636001571)
(wheel 426612541 53 60 (5 3 2) 7110209 2 (2) 213306271)
(rho 426612541 53 60 (5 3 2) 7110209 2 (2) 213306271)
fermat and lucas
(wheel 1106008481 29 160 (5 2) 6912553 6 (3 2) 184334747)
(rho 1106008481 29 160 (5 2) 6912553 6 (3 2) 184334747)
fermat and lucas
(wheel 712747153 17 1872 (13 3 2) 380741 14 (7 2) 50910511)
(rho 712747153 17 1872 (13 3 2) 380741 14 (7 2) 50910511)
fermat and lucas
(rho 851486813098510619326357783811917873528153342227994200571828674136369388533
20000000 11720532819362844 (426612541 34171 67 3 2) 726491556503144935303282197
75568844892110944591345979575003 261342420969838665149595122 (712747153 11060084
81 31859 43 11 2) 3258127057745363430685369369279126297395107564947)
fermat and lucas
(time (bls-prime? 85148681309851061932635778381191787352815334222799420057182867
4136369388533))
    908 collections
    27862 ms elapsed cpu time, including 77 ms collecting
    28062 ms elapsed real time, including 75 ms collecting
    7645064288 bytes allocated, including 7647832128 bytes reclaimed
#t

We stole many functions from the Standard Prelude and previous exercises. You can see and run the entire program at http://programmingpraxis.codepad.org/cKNnptBX.

Pages: 1 2