## Baillie-Wagstaff Pseudoprimality Tester

### March 1, 2013

There is a lot of code common to the two versions of the Lucas pseudoprimality test, so we will abstract much of the code to auxiliary functions. Here is the code to implement Selfridge’s selection of *P* and *Q*:

`(define (selfridge n)`

(let loop ((d-abs 5) (sign 1))

(let ((d (* d-abs sign)))

(cond ((< 1 (gcd d n)) (values d 0 0))

((= (jacobi d n) -1) (values d 1 (/ (- 1 d) 4)))

(else (loop (+ d-abs 2) (- sign)))))))

The biggest part of the code computes the Lucas chain. We compute and return a value `qkd`

which is used only in the strong version of the test so that we can share the code between the two tests. It’s imperative to compute the various intermediate values in the correct order; I made more than one mistake with out-of-order computations during the development of this code. For those not familiar with Scheme, the `let*`

version of `let`

computes its results left to right, with each computed value available to subsequent computations:

`(define (chain n u v u2 v2 d q m)`

(let loop ((u u) (v v) (u2 u2) (v2 v2) (q q) (qkd q) (m m))

(if (zero? m) (values u v qkd)

(let* ((u2 (modulo (* u2 v2) n))

(v2 (modulo (- (* v2 v2) (* 2 q)) n))

(q (modulo (* q q) n)))

(if (odd? m)

(let* ((t1 (* u2 v)) (t2 (* u v2))

(t3 (* v2 v)) (t4 (* u2 u d))

(u (+ t1 t2)) (v (+ t3 t4))

(qkd (modulo (* qkd q) n)))

(loop (modulo (quotient (if (odd? u) (+ u n) u) 2) n)

(modulo (quotient (if (odd? v) (+ v n) v) 2) n)

u2 v2 q qkd (quotient m 2)))

(loop u v u2 v2 q qkd (quotient m 2)))))))

With this, the standard version of the test is quite simple. The *m* value in the chain is (*n* + 1) / 2, and the only return value needed from *chain* is *u*:

`(define (standard-lucas-pseudoprime? n)`

(call-with-values

(lambda () (selfridge n))

(lambda (d p q)

(if (zero? p) (= n d)

(call-with-values

(lambda () (chain n 0 2 1 p d q (/ (+ n 1) 2)))

(lambda (u v qkd) (zero? u)))))))

For the strong version of the test, we need to split *n* + 1 into its constituent powers of two:

`(define (powers-of-two n)`

(let loop ((s 0) (n n))

(if (odd? n) (values s n)

(loop (+ s 1) (/ n 2)))))

Most of the strong test is reused from the standard test. The *m* value in the chain is ⌊*t* / 2⌋. The primary difference comes after the chain is computed:

`(define (strong-lucas-pseudoprime? n)`

(call-with-values

(lambda () (selfridge n))

(lambda (d p q)

(if (zero? p) (= n d)

(call-with-values

(lambda () (powers-of-two (+ n 1)))

(lambda (s t)

(call-with-values

(lambda () (chain n 1 p 1 p d q (quotient t 2)))

(lambda (u v qkd)

(if (or (zero? u) (zero? v)) #t

(let loop ((r 1) (v v) (qkd qkd))

(if (= r s) #f

(let* ((v (modulo (- (* v v) (* 2 qkd)) n))

(qkd (modulo (* qkd qkd) n)))

(if (zero? v) #t (loop (+ r 1) v qkd))))))))))))))

With all that, the *prime?* function is straight forward. We choose the primes less than a hundred for the trial division step of the algorithm, thougn many implementations prefer a limit of 1000; since *ps* is bound *inside* the `define`

but *outside* the `lambda`

, it is computed once, at compile time, and saved from one call of the function to the next. Here’s *prime?*:

`(define prime?`

(let ((ps (primes 100)))

(lambda (n)

(if (not (integer? n)) (error 'prime? "must be integer"))

(if (or (< n 2) (square? n)) #f

(let loop ((ps ps))

(if (pair? ps)

(if (zero? (modulo n (car ps))) (= n (car ps)) (loop (cdr ps)))

(and (strong-pseudoprime? n 2)

(strong-pseudoprime? n 3)

(strong-lucas-pseudoprime? n))))))))

We give several examples below. The first four are small calculations useful for testing; you can augment this implementation to print the results of intermediate calculations in strategic places, using them to debug your code. Then 323 and 5777 are, respectively, the first standard Lucas pseudoprime and the first strong Lucas pseudoprime, as you can see from the displayed results. The number that identified the error with the implementation in the previous exercise, 3825123056546413051, is tested next, so you can see that the function now works properly. The final test, 2^{89} − 1, is a moderate-sized prime:

`> (do ((ns (list 79 83 89 111 323 5777 3825123056546413051 (- (expt 2 89) 1)) (cdr ns))) ((null? ns))`

(let ((n (car ns)))

(display n) (display " ")

(display (strong-pseudoprime? n 2)) (display " ")

(display (strong-pseudoprime? n 3)) (display " ")

(display (standard-lucas-pseudoprime? n)) (display " ")

(display (strong-lucas-pseudoprime? n)) (display " ")

(display (prime? n)) (newline)))

83 #t #t #t #t #t

89 #t #t #t #t #t

111 #f #f #f #f #f

323 #f #f #t #f #f

5777 #f #f #t #t #f

3825123056546413051 #t #t #f #f #f

618970019642690137449562111 #t #t #t #t #t

We reused much code from several previous exercises, which you can see assembled at http://programmingpraxis.codepad.org/bHofkLUf. A Python version is also available.

Exhaustive checks by two independent sources (but both using Feitsma’s list of base-2 pseudoprimes) have been made to 2^64, which is over 10^20. This means for native 64-bit numbers, we can do a BPSW test and have a deterministic result. Noting also that we can do a single M-R test for numbers up to 291831, or two tests for up to 624,732,421 and have a deterministic result.

I have a BPSW test implemented in Perl in the Math::Prime::Util module, and one in C+GMP in the Math::Prime::Util::GMP module. I never got around to writing it in native C since I have deterministic M-R tests for all 64-bit values so it seemed of limited benefit (at best it would be a performance improvement for numbers larger than ~ 2^51).

[…] Pages: 1 2 […]

Thanks for the interesting article. When I tried it on a prime number I have tested recently, strong-lucas-pseudoprime? returned #f, so something is wrong (I suspect chain). A simple example is the prime 79:

(strong-lucas-pseudoprime? 79) => #f

Roland: Thanks for the compliment, and also for the bug report. The u/v chains were calculated correctly, but the error was the calculation of qkd; I multiplied qkd by q for each bit, but should have performed the multiplication only for the 1-bits. The text and the two versions of the program at codepad have now been fixed. The prime? function passes the test below for n = 1000000:

`(define (test n)`

(define (mesg . xs)

(display (car xs))

(do ((xs (cdr xs) (cdr xs))) ((null? xs) (newline))

(display " ") (display (car xs))))

(let loop ((i 2) (ps (primes n)))

(cond ((< n i) (void))

((and (prime? i) (null? ps))

(mesg "composite" i "reported as prime")

(loop (+ i 1) ps))

((and (prime? i) (< i (car ps)))

(mesg "composite" i "reported as prime")

(loop (+ i 1) ps))

((and (prime? i) (< (car ps) i))

(mesg "composite" i "reported as prime")

(loop (+ i 1) (cdr ps)))

((and (prime? i) (= (car ps) i))

(loop (+ i 1) (cdr ps)))

((and (not (prime? i)) (null? ps))

(loop (+ i 1) ps))

((and (not (prime? i)) (< i (car ps)))

(loop (+ i 1) ps))

((and (not (prime? i)) (< (car ps) i))

(loop (+ i 1) (cdr ps)))

((and (not (prime? i)) (= (car ps) i))

(mesg "prime" i "reported as composite")

(loop (+ i 1) ps)))))

Line 14 of the Python version should be “return n == n2 * n2”.

[…] we now have a function to compute Lucas sequences, we can take another look at testing primality using Lucas […]