## Combined N±1 Primality Prover

### March 7, 2014

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 *r*1 or *r*2 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 *b*1, stopping if it finds enough. If not, the inner `loop`

runs Pollard rho tests up to *b*2 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