Pseudoprime Statistics, Tables, and Data 
(Fermat, Miller-Rabin, Lucas, Fibonacci, Pell, Frobenius, Baillie-PSW)

by Dana Jacobsen, 31 March 2020

• PSP data from Jan Feitsma and William Galway, replicated to 1e12 with my code. SPSP-2's computed with my code from the Feitsma PSP-2 list. See the Galway link for full PSP/SPSP data to 2^64. Lucas pseudoprimes (LPSP) to 10¹⁵ were computed in 2019-2020 by Robert Baillie, Sam Wagstaff, and Andrew Fiori. SLPSP's verified identical to their data. All other data (ES Lucas, Perrin, Frobenius, etc.) generated with my code.

• Corrected: Robert Baillie informed me that my results to 10¹⁴ had missed 9 LPSPs and 1 SLPSP. The summary and data have been corrected, and many thanks for the correction and communication. This verification and new data were part of the 2019-2020 compputation by Robert Baillie, Sam Wagstaff, and Andrew Fiori of Lucas pseudoprimes to 10¹⁵. The missing data in my files were due to an issue with combining parallel result files.

• All code used from Math::Prime::Util (also known as "ntheory"; latest code in GitHub, GMP code also available). The module has a feature to efficiently loop over odd composites, using a segmented sieve internally. See the programs section at the end for examples. It is grossly less efficient than Jan Feitsma's more complicated methods for base-2 pseudoprimes. I am using my native 64-bit code, which is 3-6x faster than my GMP versions (though obviously limited to 2^64). Generation at ~2e13 for Lucas is about 1e10 range per hour per thread.

• The strong pseudoprimes to base 2 are a subset of the pseudoprimes to base 2. The strong Lucas pseudoprimes are a subset of the Lucas pseudoprimes. The extra strong Lucas pseudoprimes are a subset of the almost extra strong Lucas pseudoprimes. Because the LPSP and SLPSP tests use the Selfridge parameters while the AESLPSP and ESLPSP tests use the Baillie OEIS parameters, they have different pseudoprime sets. Below 2^64, the base-2 PSPs and SPSPs are a disjoint set from any of the listed Lucas PSPs. Above 2^64 we believe there is overlap but no examples have been found. The Frobenius(1,-1) pseudoprimes (A212424) are a subset of the odd Fibonacci pseudoprimes (A081264).

• We use the Selfridge parameters for the Lucas and strong Lucas pseudoprimes. The primary reference is Lucas Pseudoprimes by Baillie and Wagstaff (1980). When "Lucas pseudoprime" is seen in primality context, a variant of the Lucas-Selfridge test should be assumed unless otherwise specified. These are not the later Bruckman (1994) definition (OEIS A005845) which are not used since they correlate with the Fermat pseudoprimes, rather than anti-correlate like the Lucas-Selfridge pseudoprimes.

• The extra strong test is defined in Grantham 2000 theorem 2.3, as well as earlier preprints. It does not mention how to choose the parameter. The sequence above uses the parameter selection described in OEIS A217719 and Wikipedia:
  
        P = 3; Q = 1
      while (jacobi(P*P-4,n) != -1)
        P++;

  Thomas R. Nicely's code does what looks like the same test, but uses a different parameter selection (he stops when jacobi(D,n) != 0 rather than jacobi(D,n) == -1) which leads to overlap with the strong pseudoprime test, making it inappropriate for a combined (BPSW-style) test. Baillie and Wagstaff (1980) pages 1401-1406 spend some time justifying why one should use jacobi(D,n) = -1 when combining a Lucas and PSP test. The parameter choice shown earlier is the natural extension to the extra strong test and avoids all the problems Nicely discusses.
  
  MathWorld incorrectly defines the test as using Q = -1, contradicting their reference (Grantham 1998). Nicely discovered this error in 2005. Using Q = -1 results in a completely broken test. Note the comment later in their paragraph about using parameter (b,1) which is correct.
  
  Pari's documentation does not match the code. It indicates a strong Lucas test, but the implementation for years has been an almost extra strong test (see below). It indicates Q = -1 and Q = 1 in the same sentence. As expected, the code uses Q = 1. It indicates P is the smallest positive integer meeting the conditions, where in fact P starts at 3 and increments by 2.
  
  David Cleaver's mpz_prp software has code for the extra strong test. It does not define a parameter selection method. Note that previous to version 1.1 (3 July 2013) it had an error in the ES test. GMPY2 uses this code, so get at least version 2.0.1 if you use the ES test.

• The "almost extra strong" test is a generalization of the Lucas test used in Pari's is_pseudoprime. Use Crandall and Pomerance algorithm 3.6.7 to quickly get the V parameter, then perform an extra strong test only using V. This leads to some more pseudoprimes (since the U = 0 condition is ignored) but can be implemented such that it runs extremely fast, and still produces fewer pseudoprimes than the strong test.
  
  The version used in the table uses the Baillie OEIS parameter selection, so it is a superset of the OIES A217719 extra strong test. It takes about 2/3 the time of a strong Lucas-Selfridge test in my implementations. As n increases in size, the times converge however, so this is most useful for small n, depending on your big integer library.
  
  It is also possible to recover U_n using Crandall and Pomerance equation 3.13: U_n = D^-1 (2V_{n+1} - PV_n) allowing us to run the full extra-strong test at the cost of a single modular inversion. This computation is easy and fast in GMP, so we can get the full extra-strong test at essentially the same performance as the almost extra strong test.
  
  The Pari implementation has two differences. First, it does not include the required gcd check, probably because the overall implementation does trial division to 101 so it is unlikely to be necessary. More importantly, it uses a slightly different parameter selection: P += 2 in the method shown above, instead of P += 1. This leads to a different set of pseudoprimes, though the quantity produced seems to be similar. Here are the 7090 Pari LucasPsP pseudoprimes to 1e11.

• The BPSW test uses a combination of strong probable prime test with base 2 (the SPSP-2 column above) and a Lucas test. The reference by Baillie and Wagstaff indicate this should be a strong Lucas test with one of two specified parameter selections (the Selfridge method being the most commonly used). However, there are no BPSW pseudoprimes below 2^64 using any of the four Lucas tests above, and no known BPSW pseudoprime of any size. While many math packages say they uses a BPSW test, many seem to use a custom variant of the Lucas test and will produce a different set of pseudoprimes than the above. Note that the Bruckman (OEIS A005845) "Lucas pseudoprimes" should never be used here.

• There are additional, easy to apply, tests that can be performed to strengthen the Lucas tests. See Wikipedia's Lucas PSP page or section 6 of the Lucas Pseudoprimes paper.

• While the Fermat, Miller-Rabin, and Lucas tests are the most commonly used, there are other compositeness tests (with associated pseudoprimes).
  
  Of some interest are the Frobenius pseudoprimes, as they have much smaller error bounds: 1/7710 (cite) vs. 1/8 for the extra strong Lucas test (cite), 4/15 for strong Lucas test (cite), 1/4 for the Miller-Rabin test (cite). The computational cost is about 3x a Miller-Rabin test. Note that this bound applies to a specific strong QFT test with additional filters, while all current public implementations, including Perl/ntheory, use the "standard" Frobenius test. Read all of Crandall/Pomerance section 3.6 for more information. The examples shown of using (3,-5) or (P,2) do show it is a very effective test. Additional variants (EQFT, SQFT, MQFT) have been proposed -- I am only aware of one implementation of the variants (SQFT using RELIC for polynomial manipulation).
  
  If we create a Frobenius test with parameters (P,2) defining the monic quadratic polynomial x²-Px+2, then theorem 4.3 from Grantham 2000 shows that all pseudoprimes to the test will also be Fermat base-2 pseudoprimes, and theorem 4.9 that they are also Lucas (P,2) pseudoprimes. This allows us to verify there are no counterexamples to 2⁶⁴. P is selected as the first odd integer such that (D|n)=-1 where D = P*P-4Q. This is implemented as is_frobenius_pseudoprime($n) with no additional arguments.
  
  Khashin (2013) gives a different look at the Frobenius test. This has been implemented in my software and exhaustive testing has found no counterexamples to 2⁴³.
  
  Paul Underwood has created a Frobenius-style test: The Efficient Calculation of a Combined Fermat PRP and Lucas PRP Test (2014) (originally described in "Quadratic Compositeness Tests", preprint, 2012) that costs about 2.4x the cost of a Miller-Rabin test and has no known counterexamples (exhaustively verified to 2⁵⁰). This is basically an all-in-one BPSW type test, though using different bases for the two tests (this can be an advantage if combining with BPSW). Code in C, C+GMP, Perl, Javascript, and gwnum.
  
  The Fibonacci and Pell pseudoprimes have significant overlap with other tests, so there doesn't seem to be any compelling reason to use them.
  
  The Perrin pseudoprimes correlate strongly with the base 2 Fermat pseudoprimes and the Lucas test is computationally more efficient. The small number of pseudoprimes makes them interesting, but they are not usually used.
  
  Catalan pseudoprimes are quite rare (only 3 are known, with exhaustive testing to 10¹⁰), but are very expensive to calculate and are generally not used.

• Jan Feitsma's list of pages
• Jeff Gilchrist's Pseudoprime page
• ATH's list of Lucas pseudoprimes to 3*10^12
• The best known SPRP bases sets. This shows minimal sets required for deterministic testing of small numbers

Programs, wrap in perl -Mntheory=:all -E '...'

• PSP2: foroddcomposites { say if is_pseudoprime($_,2) } 1e10
• SPSP2: foroddcomposites { say if is_strong_pseudoprime($_,2) } 1e10
• LPSP: foroddcomposites { say if is_lucas_pseudoprime($_) } 1e10
• SLPSP: foroddcomposites { say if is_strong_lucas_pseudoprime($_) } 1e10
• AESLPSP: foroddcomposites { say if is_almost_extra_strong_lucas_pseudoprime($_) } 1e10
• Pari-like:foroddcomposites { say if is_almost_extra_strong_lucas_pseudoprime($_,2) } 1e10
• ESLPSP: foroddcomposites { say if is_extra_strong_lucas_pseudoprime($_) } 1e10
• Frobenius: foroddcomposites { say if is_frobenius_pseudoprime($_,1,-1) } 1e10
• Fibonacci: foroddcomposites { $e = (0,-1,1,1,-1)[$_%5]; say unless $e==0 || (lucas_sequence($_, 1, -1, $_+$e))[0] } 1e10
• Even Fibonacci: for (2..1e10) { my $n=$_<<1; $e = (0,-1,1,1,-1)[$n%5]; say $n unless !$e || (lucas_sequence($n,1,-1,$n+$e))[0] }
• Bruckman: foroddcomposites { say if (lucas_sequence($_,1,-1,$_))[1] == 1 } 1e10
• Pell: foroddcomposites { say unless ((lucas_sequence($_, 2, -1, $_))[0] - kronecker(2,$_)) % $_; } 1e10
• Perrin: forcomposites { say if is_perrin_pseudoprime($_) } 1e10
• Restricted Perrin: forcomposites { say if is_perrin_pseudoprime($_,1) } 1e10
• Catalan: forcomposites { say if is_catalan_pseudoprime($_) } 1e10

• Frobenius(P,2): foroddcomposites { say if is_frobenius_pseudoprime($_); } 1e10
• FK: foroddcomposites { say if is_frobenius_khashin_pseudoprime($_); } 1e10
• FU: foroddcomposites { say if is_frobenius_underwood_pseudoprime($_); } 1e10
• BPSW: foroddcomposites { say if is_bpsw_prime($_); } 1e10

Programs for older MPU:

• Frobenius: foroddcomposites { $k = kronecker(5,$_); ($U,$V)=lucas_sequence($_,1,-1,$_-$k); say if $U == 0 && $V == (($k==1)?2:$_-2) } 1e10
• Perrin: use Math::GMPz; { my $l = 2; my @A = map { Math::GMPz->new($_) } (3,0,2); sub perrin { my $n = shift; return 0 if $n <= 2; die if $n < $l; while ($n > $l) { @A = ($A[1],$A[2],$A[0]+$A[1]); $l++; } $A[-1] } } forcomposites { say unless perrin($_) % $_ } 1e10
• Catalan: use Math::GMPz; { my($c,$l)=(Math::GMPz->new(1),1); sub catalan { while ($_[0] > $l) { $l++; $c *= 4*$l-2; Math::GMPz::Rmpz_divexact_ui($c,$c,$l+1); } $c; } } foroddcomposites { $m = ($_-1)>>1; say if (catalan($m) % $_) == (($m&1) ? $_-2 : 2); } 1e10
• Catalan (alt): foroddcomposites { $m = ($_-1)>>1; say if (binomial($m<<1,$m) % $_) == (($m&1) ? $_-1 : 1); } 1e10