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