Baillie-Wagstaff Pseudoprimality Tester

March 1, 2013

We examined the Baillie-Wagstaff probabilistic primality checker based on Lucas pseudoprimes in a previous exercise, where we got it wrong; among other things it incorrectly reported that 3825123056546413051 = 149491 × 747451 × 34233211 was prime. We’ll fix that today. Robert Baillie and Samuel Wagstaff Jr. describe the test in their article “Lucas Pseudoprimes” in Mathematics of Computation, Volume 35, Number 152, pages 1391-1417, October 1980; see also Thomas Nicely’s web page devoted to the test. This is the same algorithm used in the PrimeQ function in Mathematica. If at any point the math gets daunting, feel free to skip ahead to the code, which is short and simple.

Lucas numbers are defined by a recurrence formula similar to Fibonacci numbers, where Ln = Ln−1 + Ln−2 with L1 = 1 and L2 = 3; the Lucas numbers are 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, … (Sloane’s A000204). Lucas numbers have the rather startling property that, if n is prime, Ln ≡ 1 (mod n). But the converse is not true, and composite numbers n such that Ln ≡ 1 (mod n) are known as Lucas pseudoprimes; the first few Lucas pseudoprimes are 705, 2465, 2737, 3745, 4171, … (Sloane’s A005845).

Lucas numbers are a special case of Lucas sequences. If P and Q are integers such that the discriminant D = P2 − 4Q, then the roots of x2Px + Q = 0 are a = (P + √D) / 2 and b = (P − √D) / 2 (that’s the quadratic equation you learned in middle school). There are two Lucas sequences Un(P,Q) = (anbn) / (ab) and Vn(P,Q) = an + an for n ≥ 0, which can be computed by the recurrence equations Um(P,Q) = PUm−1(P,Q) − QUm−2(P,Q) and Vm(P,Q) = PVm−1(P,Q) − QVm−2(P,Q). The Lucas numbers are given by the sequence V(1, -1) and the Fibonacci numbers are given by the sequence U(1, -1).

It is possible to compute the nth element of a Lucas sequence in logarithmic time rather than linear time using a double-and-add algorithm similar to the peasant algorithm for multiplication; such a computation is known as a Lucas chain. The rules are U2n = Un Vn, U2n+1 = Un+1 VnQn, V2n = Vn2 − 2 Qn, and V2n+1 = Vn+1 VnP Qn. For our purposes, we will be making all computations mod n.

Thus we have a method to determine if a number n is a Lucas pseudoprime. Choose a suitable P and Q, use a Lucas chain to compute Un(P,Q), and declare n either prime or a Lucas pseudoprime if Un(P,Q) ≡ 0 (mod n). There are several methods to compute suitable P and Q, of which the most commonly-used is due to John Selfridge: choose D as the smallest number in the sequence 5, -7, 9, -11, 13, -15, 17, … for which the Jacobi symbol (D / n) = -1, then set P = 1 and Q = (1 – D) / 4.

There is a strong version of the Lucas pseudoprimality test that bears the same relationship to the standard version of the Lucas pseudoprimality test described above as Gary Miller’s strong Fermat pseudoprimality test bears to the standard Fermat pseudoprimality test. Represent n as d · 2s + 1, then compute Ud(P,Q) and Vd(P,Q). Then n is either prime or a Lucas pseudoprime if Ud(P,Q) ≡ 0 (mod n) or if Vd 2r(P,Q) ≡ 0 (mod n) for any r with 0 ≤ r < s.

Finally, the Baillie-Wagstaff test of the primality of n consists of the following steps: 1) if n is not an integer, report a domain error; 2) if n < 2, report n is not prime; 3) if n is a square, report n is not prime; 4) if n is divisible by any primes less than some convenient limit (e.g. 100 or 1000), report n is not prime; 5) if n is not a strong Fermat pseudoprime to base 2, report n is not prime; 6) (optional, added by Mathematica) if n is not a strong Fermat pseudoprime to base 3, report n is not prime; 7) if n is not a Lucas pseudoprime (standard or strong), report n is not prime; 8) otherwise, report n is prime.

The Baillie-Wagstaff pseudoprimality test, whether or not the base-3 test is included and whether the standard or strong Lucas test is chosen, is extremely strong: there are no known errors. Exhaustive checks have been performed to 1017, and the test has been used for years in several high-quality primality proving programs that would have found an error if one occurred. There is a $620 reward for a counter-example to the test ($500 from Selfridge, $100 from Wagstaff, and $20 from Carl Pomerance), plus the certain fame that will accrue to the finder.

Your task is to write a Baillie-Wagstaff primality checker as described above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.


Pages: 1 2

6 Responses to “Baillie-Wagstaff Pseudoprimality Tester”

  1. danaj said

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

  2. Roland said

    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

  3. programmingpraxis said

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

  4. Dang Griffith said

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

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

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: