Thomas R. Nicely
Last updated 0700 GMT 13 January 2012. Originally posted 10 June 2005.
The Baillie-PSW (BPSW or BSW) primality test is actually a compositeness test, in the manner of Fermat's test and the Miller-Rabin test. It is named for Robert Baillie, Carl Pomerance, John L. Selfridge, and Samuel S. Wagstaff, Jr. The algorithm was apparently first conceived by Baillie (Baillie and Wagstaff, 1980), with refinements added by Selfridge (Pomerance, Selfridge, and Wagstaff, 1980). The procedure is as follows; note that some details may vary from one implementation to another (Martin, 2004).
As of this date, both the standard and strong tests are flawless; no counterexample (BPSW standard or strong pseudoprime) is known. With the aid of Richard G. E. Pinch's table of base-2 strong pseudoprimes, the author verified (May, 2005) that no BPSW pseudoprime (standard or strong) exists for N < 10^13. In January, 2007, with the aid of the author's code and William Galway's table of pseudoprimes, Martin Fuller determined that no BPSW pseudoprime (standard or strong) exists for N < 10^15. More recently, Jeff Gilchrist, using the author's code and a database of pseudoprimes prepared by Jan Feitsma, has verified (13 June 2009) that no BPSW pseudoprime (standard or strong) exists below 10^17. Furthermore, analysis (24 October 2009) by Gilchrist of an extension of Feitsma's database to 2^64 found no BPSW pseudoprime (standard or strong) below 2^64 (approximately 1.8446744e19). This lower bound of 2^64 for any BPSW pseudoprime has since been separately verified (11 July 2011) by Charles Greathouse; however, he also used Feitsma's database, so a completely independent check has not yet been carried out.
Note, however, that Carl Pomerance (1984) presented a heuristic argument that an infinite number of counterexamples exist to the standard test (and presumably to the strong test as well, based on similar reasoning), and even that (for sufficiently large x, dependent on µ) the number of BPSW pseudoprimes < x exceeds x^(1-µ), where µ is an arbitrarily small pre-assigned positive number. Nevertheless, not a single BPSW pseudoprime has yet been found. Consequently, the BPSW test carries an aura of dependability (justified or not) exceeding that of competing algorithms, such as multiple Miller-Rabin tests.
It is believed that some commercial mathematics software packages rely, in part or in whole, on the BPSW test for primality checking; see, for example, Ribenboim (1995/6, p. 142), also Martin (2004). In some instances, these packages appear to use the strong BPSW and/or Lucas test, or add additional Miller-Rabin tests.
BPSW requires O((log n)^3) bit operations, as do Fermat's test and the Miller-Rabin algorithm. However, a BPSW test typically requires roughly three to seven times as many bit operations as a single Miller-Rabin test. The strong version of BPSW differs only in replacing the standard Lucas-Selfridge test with the strong Lucas-Selfridge test. The strong Lucas-Selfridge test produces only roughly 30 % as many pseudoprimes as the standard version; for example, among the odd composites N < 10^6, there are 219 standard Lucas-Selfridge pseudoprimes, 58 strong Lucas-Selfridge pseudoprimes, and 46 base-2 strong pseudoprimes (these totals presume no screening with odd trial divisors). Since the strong Lucas-Selfridge test incurs roughly 10 % more running time, the strong tests appear to be more effective.
Selfridge specified the following parameters for the generation of the Lucas sequences: P = 1 and Q = (1 - D)/4, where D is the first integer in the sequence {5, -7, 9, -11, 13, -15, ...} for which GCD(D,N)=1 and the Jacobi symbol (D|N) = -1. Note, however, that if N is a perfect square, no such D will exist, and the search for D would continue all the way to ±sqrt(N), at which point GCD(D,N)=sqrt(N) would expose N as composite. Consequently, the algorithm also presumes the presence of a preliminary check for perfect squares (as well as even integers and integers < 3).
Additional conditions sometimes imposed (Riesel, 1994, p. 130), that N not divide Q and that D be square-free, were found to be irrelevant to this application of Lucas sequences.
There remains the theoretical concern that, for very large N, the required D would be so large in magnitude that the algorithm would become computationally useless. Given the exclusion of perfect square values of N, no such behavior has been observed in practice; indeed, the range of values observed for D is quite small---|D|max=47 for N < 10^6, |D|max=67 for 10^19 < N < 10^19 + 10^6. Baillie and Wagstaff (1980) presented analytical arguments to support this observation (Ribenboim, 1995/6, p. 142). However, the eventuality cannot be entirely excluded, and the author's code includes an overflow trap for D.
The author has translated the algorithm into GNU C, employing the GMP library, with a command-line executable compiled for the DOS/Wintel operating environment; see the accompanying zipfile, which includes the source codes and Wintel/DOS executables. To examine the complete library source modules or to recompile the codes, you will also need to download the support files trn.c and trn.h from trn.zip. The illustrative main routine computes pi(x) between specified bounds, and enumerates the corresponding pseudoprimes produced by various primality tests, including the standard and strong BPSW tests (none found to date), the standard, strong, and extra strong Lucas tests, and various others, including the Fermat and Miller-Rabin tests (base 2) as well as the GMP mpz_probab_prime_p function (see GNU GMP mpz_probab_prime_p pseudoprimes).
Code is also included in trn.c for the "extra strong" Lucas test, as developed by Zhaiyu Mo and James P. Jones (circa 1997), and described by Jon Grantham (1998). However, I have not employed the extra strong Lucas test in the BPSW test, as the Lucas sequence parameters are inconsistent with those of the Lucas-Selfridge tests; consequently, the extra strong Lucas pseudoprimes were not found to be disjoint from those of the Miller-Rabin test with base 2 (or any other single Miller-Rabin base employed). This behavior was observed for all bases employed for the extra strong Lucas test. For a more detailed discussion, see the documentation in the module iExtraStrongLucas in the support file trn.c included in trn.zip.
I would appreciate notification of any significant errors found in this document, the codes, or their results.