GNU GMP mpz_probab_prime_p pseudoprimes

Thomas R. Nicely
Current e-mail address

Freeware copyright (C) 2018 Thomas R. Nicely. Released into the public domain by the author, who disclaims any legal liability arising from its use.

Last updated 0300 GMT 23 November 2018. Original posting December 2004.


The contents of this document are NOT to be interpreted as implying an error or bug in GCC, GNU C, or GMP (also known as GNU MP), the mpz_probab_prime_p function, or the Miller-Rabin algorithm. Instead, these results should be treated as information concerning the consequences and impact of the inherent mathematical limitations of the Miller-Rabin algorithm on the mpz_probab_prime_p function and its use. For detailed and in-depth discussions of pseudoprimes and the Miller-Rabin algorithm, see, for example, Prime numbers and computer methods for factorization, Hans Riesel (2nd ed., 1994, Birkhäuser, Boston) 84-140; also The new book of prime number records, Paulo Ribenboim (3rd ed., 1996, Springer-Verlag, New York) 19-178.

Note that the exponential operator is indicated herein by ** rather than ^.

General remarks

These results enumerate counterexamples (exceptions or false positives) returned by the GNU GMP mpz_probab_prime_p primality testing function (which is in fact a function for testing compositeness). Note that a return value of 1 from this function is in fact documented as indicating only that the multiple precision integer argument N (type mpz_t) is a "probable prime" (more accurately, a "strong probable prime"). However, many applications are written with the tacit assumption that a return value of 1 indicates that N is a prime "for all practical purposes", particularly if the number of bases (repetitions, Miller-Rabin tests) specified by the second argument REPS (a positive integer) is 5 or greater. In effect, we are addressing the reliability of that assumption.

Note also that if mpz_probab_prime_p returns a value of 2, it is guaranteeing that N is prime; in this case, trial divisors to the square root have been used, rather than the Miller-Rabin algorithm. However, this value is only returned for small primes, less than perhaps 1000000. Finally, a return value of zero indicates that the integer is unequivocally not a prime (N is composite or less than 2).

Let us define an integer N to be a member of (or to have the attribute) mpz_spspk if it is (a) composite, (b) the GMP function mpz_probab_prime_p(N, k)=1 (indicating a strong probable prime), and (c) mpz_probab_prime_p(N, k+1)=0 (guaranteeing that N is not prime). Such an integer will also be referred to as an mpz_spsp of order k, or simply an mpz_spspk---a strong pseudoprime to k bases for the function mpz_probab_prime_p---or generically as an mpz_spsp. It is characterized by the fact that the performance of k repetitions of Miller's test by mpz_probab_prime_p is insufficient to resolve N as a composite, whereas k+1 repetitions is sufficient. Note that this differs from my previous definition by the inclusion of condition (c), so that a given N has only one specific order (for a particular implementation of mpz_probab_prime_p).

Thus specifying k different bases for Miller's test on N will return k false positives, while specifying k+1 different bases will return a true negative (from the (k+1)st base). The actual values of the bases are chosen internally by the mpz_probab_prime_p algorithm; are local to the procedure and not visible to the user; and appear to vary with the value of N, and, to some extent, the version of GMP used. The internal consistency of the definition depends on the hypothesis that each particular mpz_probab_prime_p algorithm constructs a set of j+1 bases for a given N by appending a new and distinct (j+1)st base to the set produced by the choice for N and exactly j bases, for each j >= 1. No exceptions to this hypothesis have been observed.

If mpz_probab_prime_p(N, 1)=0, then N is not an mpz_spsp pseudoprime, but for classification purposes may be described as mpz_spsp0.

As first noted by Clifford Stern (27 July 2007), the specific mpz_spsp pseudoprimes depend to some extent on the specific version of GMP from which mpz_probab_prime_p is compiled. This reflects changes in the algorithm (principally in the internal choice of bases) and in the supporting routines (such as the pseudorandom number generator). Incidentally, the GMP version should be determined from the value of the library-defined constant gmp_version, rather than the header-defined constants __GNU_MP_VERSION, __GNU_MP_VERSION_MINOR, and __GNU_MP_VERSION_PATCHLEVEL; they are not always the same, due to installation anomalies.

GMP 4.0.1-4.1.4

The following results, obtained by the author, apply to versions 4.0.1, 4.1.2, and 4.1.4 of GMP (and probably to any intermediate versions as well). All known mpz_spsp's of order 10 or greater are presented, as well as the smallest known instance of each mpz_spsp of order less than 10. Note that only two 64-bit integers (209881968380913403 and 2370020490966486403) have been found which are of order ten or greater.

ORDER 12: The following integers are mpz_spsp12:

  N=43628926643584801291 (8 December 2004).
  N=99120207367868007331 (9 December 2004).
  N=2617057824672461696191 (15 December 2004).

ORDER 11: The following integers are mpz_spsp11:

  N=209881968380913403 (6 December 2004).
  N=19600159758889166731 (8 December 2004).
  N=455308758219013329403 (10 December 2004).
  N=1518132547139756696371 (14 December 2004).
  N=8031151951625760386288203 (11 December 2004).

ORDER 10: The following integers are mpz_spsp10:

  N=2370020490966486403 (6 December 2004).
  N=157351910012685225991 (9 December 2004).
  N=349228439450094648151 (10 December 2004).
  N=396594161005255290703 (10 December 2004).
  N=553755470998517690671 (10 December 2004).
  N=573990913701073033903 (10 December 2004).
  N=773464267476853150771 (12 December 2004).
  N=1240689310336284744781 (12 December 2004).
  N=1502965791107391010051 (14 December 2004).
  N=1574864695514005075351 (14 December 2004).
  N=1662627626869742333191 (14 December 2004).
  N=2794779165895069882231 (15 December 2004).
  N=3472182119018573666371 (15 December 2004).
  N=3513477856210217250391 (15 December 2004).
  N=3692761975706924194903 (15 December 2004).
  N=8000031293140601835421003 (11 December 2004).
  N=8063333906687749281446011 (12 December 2004).
  N=18066535549432303237551253 (11 December 2004).
  N=18115413527322036686653903 (12 December 2004).
  N=18123419880437997086265151 (12 December 2004).
  N=200000000135062271492802271468294969951 (16 December 2004).

ORDER 9: N=13856432344989391 (6 December 2004).
ORDER 8: N=1262316063544903 (6 December 2004).
ORDER 7: N=2237630988421 (6 December 2004).
ORDER 6: N=3398697686971 (6 December 2004).
ORDER 5: N=5563452248551 (6 December 2004).
ORDER 4: N=25366866661 is the first occurrence (7 December 2004).
ORDER 3: N=1152271 is the first occurrence (5 December 2004).
ORDER 2: N=1590751 is the first occurrence (5 December 2004).
ORDER 1: N=1039849 is the first occurrence (5 December 2004).

An exhaustive scan to 50e9 (GMP 4.0.1-4.1.4) has enumerated:

 543 mpz_spsp's of order 1.
  70 mpz_spsp's of order 2.
  19 mpz_spsp's of order 3.
   3 mpz_spsp's of order 4.
   0 mpz_spsp's of order 5 or greater.

Strong pseudoprimes to multiple bases

For purposes of comparison, the smallest N's which are strong pseudoprimes simultaneously to the first k prime bases are listed below (those marked with § are only first known occurrences---upper bounds). These were discovered by Carl Pomerance, John L. Selfridge, and Samuel S. Wagstaff, Jr. (The pseudoprimes to $25 \times 10**9$, Math. Comp. 35 (1980) 1003-1026) and by Gerhard Jaeschke (On strong pseudoprimes to several bases, Math. Comp. 61 (1993) 915-926), and are quoted from The new book of prime number records, Paulo Ribenboim (3rd ed., 1996, Springer-Verlag) 115. The last entry is Arnault's number ( Sur quelques tests probabilistes de primalité, François Arnault, thesis (1993) 28), a strong pseudoprime to all prime bases less than 200. Among these integers, only 1373653 (order 1) and Arnault's number (order 2) are mpz_spsp's in GMP 4.0.1-4.1.4. In GMP 4.2.1, 2152302898747, 3474749660383, and Arnault's number are of order one.

 k= 1, B=2: N=2047.
 k= 2, B=2,3: N=1373653.
 k= 3, B=2,3,5: N=25326001.
 k= 4, B=2,3,5,7: N=3215031751.
 k= 5, B=2,3,5,7,11: N=2152302898747.
 k= 6, B=2,3,5,7,11,13: N=3474749660383.
 k= 7, B=2,3,5,7,11,13,17: N=341550071728321.
 k= 8, B=2,3,5,7,11,13,17,19: N=341550071728321.
 k= 9, B=2,3,5,7,11,13,17,19,23: N=41234316135705689041.§
 k=10, B=2,3,5,7,11,13,17,19,23,29: N=1553360566073143205541002401.§
 k=11, B=2,3,5,7,11,13,17,19,23,29,31: N=56897193526942024370326972321.§
 k=46: N=803837457453639491257079614341942108138837688287558145837488917\
Jaeschke also notes that there are no strong pseudoprimes, below 10**12, common to the four prime bases 2, 13, 23, and 1662803 (cited from Prime numbers and computer methods for factorization, Hans Riesel (2nd ed., Birkhäuser Boston, 1994), p. 92).

GMP 4.2.1-4.3.1

The following results, obtained by the author (except as noted), apply to versions 4.2.1 and 4.3.1 of GMP (and probably to any intermediate versions as well).

All of the mpz_spsp's of GMP 4.0.1-4.1.4, enumerated above, are now more quickly resolved as composites. None is of order greater than 3; only 8063333906687749281446011 is of order 3, and only 349228439450094648151 is of order 2. Thirty are of order zero (not mpz_spsp pseudoprimes at all).

However, it appears (as expected) that the mpz_probab_prime_p of GMP 4.2.1-4.3.1 has simply substituted one set of mpz_spsp's for another, and that the two sets are of approximately the same size. To quickly test this conjecture, all composites of the form p(2p-1), where p and 2p-1 are prime and p < 10**8, were tested using separate codes compiled first with GMP 4.1.2 (under DJGPP 2.04 and GCC 4.1.2) and secondly with GMP 4.2.1 and GMP 4.3.1 (under MinGW with GCC 3.4.2). The results (counts of each order) were as follows.

MPZ_SPSP    GMP         GMP
 ORDER     4.1.2    4.2.1-4.3.1
   0      377353         377554  (not mpz_spsp pseudoprimes)
   1       35836          35772
   2        7996           7911
   3        1802           1782
   4         413            401
   5          89             86
   6          31             21
   7           4              3
   8           1              2
   9           1              0

An exhaustive scan to 50e9 (GMP 4.2.1-4.3.1) has revealed a total of:

 559 mpz_spsp's of order 1.
  68 mpz_spsp's of order 2.
  11 mpz_spsp's of order 3.
   0 mpz_spsp's of order 4.
   0 mpz_spsp's of order 5 or greater.
Clearly, there is no significant difference in the results from GMP 4.1.2 and GMP 4.2.1-4.3.1. There is a very slight advantage in favor of GMP 4.2.1-4.3.1, but this could well disappear (or be reversed) if the analysis were extended to a greater upper bound; my own expectation is that in each case the ratio #(4.3.1, k)/#(4.1.2, k) would approach unity asymptotically.

Up to 50e9 there are three more mpz_spsp pseudoprimes for GMP 4.2.1-4.3.1 than for GMP 4.0.1-4.1.4; however, GMP 4.2.1-4.3.1 produces fewer pseudoprimes above order one (and below 50e9) than GMP 4.0.1-4.1.4.

It appears that GMP 4.2.1-4.3.1 incorporates changes which alter the specific set of bases used by mpz_probab_prime_p(n, k) for each n and k. However, the global frequency of pseudoprimes is ultimately a consequence of the number of bases k, so the global frequency has not changed significantly. The details of mpz_probab_prime_p have changed (and probably will again, in some later version of GMP), but its ultimate nature has not: it is really a compositeness test, and only a probabilistic primality test. Unless, of course, you wish to increase the number of bases to k >= n/4, in which case Rabin's theorem guarantees that a return of 1 certifies primality; unfortunately, this version of the Miller-Rabin test is even more computationally expensive than using prime trial divisors to the square root of n.

As a further note, the following specific mpz_spsp pseudoprimes have been observed for GMP 4.2.1-4.3.1; only the smallest known occurrences are listed (for orders 1, 2, and 3, these are in fact the first occurrences).

Order            N                      Factorization
 13  189452997113368438678231 =  307776702427 * 615553404853 (MT, 2011)
 12    2867561004669023153611 =   75730588333 *  37865294167 (MT, 2011)
 11      35051193821423205091 =    4186358431 *   8372716861
 10       1858358298078962371 =     963939391 *   1927878781
  9        120175393101413203 =     245127919 *    490255837
  8            14910591535003 =       2730439 *      5460877
  7            47227901175703 =       4859419 *      9718837
  6            25768686677581 =       3589477 *      7178953
  5             2028576353203 =       1007119 *      2014237
  4              239626837621 =        346141 *       692281
  3                 465658903 =         15259 *        30517
  2                   1943521 =      61 * 151 * 211
  1                   1152271 =      43 * 127 * 211

The entries marked (MT, 2011) were discovered (using GMP 4.3.2, which I do not have installed) by Markus Tervooren and reported on 22 January 2011. I have verified their status under GMP 4.2.1 and 4.3.1; under GMP 4.1.2 they are order zero (not mpz_spsp pseudoprimes).

Among the above entries, only 1152271 (order three) is of order greater than one in GMP 4.0.1-4.1.4.

The integer 4972725568723375711 = 1576820467 * 3153640933 was found to be of order 5 for both GMP 4.0.1-4.1.4 and GMP 4.2.1-4.3.1.

GMP 5.0.1

On 15 February 2011, Andreas Höglund communicated to me the results of a similar analysis he has performed with version 5.0.1 of GMP. His results are available here.

GMP 6.1.2

On 22 November 2018, Paul Zimmermann communicated to me some findings of a similar analysis he has performed using version 6.1.2 of GMP. Among his results are the following:

Other primality tests

Finally, note that there is the possibility that a more effective "fast" primality test---O((ln N)**3), like the Miller-Rabin test, but with fewer pseudoprimes---can be constructed using one or more Miller-Rabin tests in conjunction with other compositeness tests of a different type. One prominent example is the Baillie-PSW algorithm (q.v.), which incorporates a compositeness test based on Lucas sequences. The author has verified (May, 2005) that, for N < 10**13, no counterexample exists for this algorithm, as implemented in, which includes source code, support routines, and a DOS/Wintel executable implementing the Baillie-PSW test. In January, 2007, with the aid of Nicely's code and William Galway's table of pseudoprimes, Martin Fuller determined that no Baillie-PSW pseudoprime (standard or strong) exists for N < 10**15. More recently, Jeff Gilchrist, using Nicely's code and a database of pseudoprimes prepared by Jan Feitsma, has verified (13 June 2009) that no Baillie-PSW pseudoprime (standard or strong) exists below 10**17. Furthermore, preliminary analysis by Gilchrist of Feitsma's database, further extended, indicates that no Baillie-PSW pseudoprime (standard or strong) exists below 2**64 (24 October 2009; double checking of this result is in progress), approximately 1.8446744e19. See The Baillie-PSW primality test for further details.