MATHEMATICA IN ACTION for issue 8-4
 of Mathematica in Education and Research

Certified Primes

An elegant recursive method for getting large integers that are certifiably prime outperforms traditional ways of generating mere “probable primes”.

by  Stan Wagon

Suppose you want a random 100-digit prime. The obvious way to proceed is to use PrimeQ, as in the following code, which generates a random d-digit prime by choosing the first prime after a random d-digit integer.

PrimeAfter[n_] := 
        2 /; -2 ≤ n ≤ 1;
\(\(PrimeAfter[n_] := 
        Module[{j = n + 1 + Mod[n, 2]}, 
          While[¬ PrimeQ[j], j += 2]; 
          j];\)
		\)
\(\(RandomPrime[d_] := 
        Module[{n}, 
\ \ While[
            n = PrimeAfter[Random[Integer, 10\^\(d - {1, 0}\)]]; 
            Log[10. , n] > d]; n];\)\)

RandomPrime[5]

rp = RandomPrime[50]

The problem with this approach is that it is conceivable (though no examples are known; see [Wagon 1999]) that PrimeQ would return True on a composite integer n. Thus it is of interest to find a method that produces integers that are certifiably prime. The primary motivation is esthetic, since it is annoying to have to say “very probably prime” when one wishes to say “prime”; in practical applications to cryptography one tends to not worry about the chance of getting a nonprime because this probability can be made very small with a few more pseudoprime tests. One could use the routines in the standard package PrimeQ` to validate the output of RandomPrime. The ProvablePrimeQ function in that package uses elliptic curve methods and works reasonably well, but it is slow and its theoretical underpinnings are quite advanced.

Needs["NumberTheory`PrimeQ`"]

ProvablePrimeQ[rp] // Timing

The method that we will describe in this article (due to U. Maurer [Maurer 1995]) produces certified primes very quickly and is both simple and elegant.

.  The Theory

Maurer’s work relies on the following result, discovered by H. C. Pocklingon in 1914. The main idea is similar to the idea of Pratt certificates discussed in an earlier column of mine [Wagon 1997], but with the important new point that the witness, b, need only work for a partial factorization of n. We give a simplified version of Pocklington’s theorem where the partial factorization consists of a single prime q; this is adequate to the task at hand. For the general version, which handles the case that the factored part of n - 1 has several primes (as opposed to just one), see [Bressoud and Wagon, 2000].

POCKLINGTON’S THEOREM.  Suppose n - 1 has the prime factor q, where q > (Sqrt[n]) and the following hold:

1.  (2^(n - 1)) ≡ 1 (mod n)

	2.  gcd((2^((n - 1)/q))-1, n)=1.

Then n is prime.

PROOF.  Suppose n is not prime and p is its smallest prime divisor; then p < (Sqrt[n]) and p is odd (because, by (1), n is odd). Let d be the mod-p order of 2. Then, by standard results, d divides p - 1 (recall that Fermat’s little theorem tells us that (2^(p - 1)) ≡ 1 (mod p)). Further, by (1), d divides n - 1. But d does not divide (n - 1)/q, because doing so would violate (2) (for then (2^((n - 1)/q))≡1 (mod p) and so p would divide both n and (2^((n - 1)/q)) - 1). Therefore q divides d. But d divides p - 1, so q divides p - 1, which means that  p - 1 ≥ q ≥ (Sqrt[n]), a contradiction.   □

Of course, any integer b can serve as a witness in place of 2; but we may as well restrict to the choice of 2 because, if n is truly prime, the chance that 2 is not a witness is very small (about 1/(Sqrt[n])).

And now here is Maurer’s method of generating certified primes. Suppose we want p, a d-digit prime that is certifiably prime. Select a prime q having ⌈(d/2)⌉ + 1 many digits (q should be a certified prime: this can be done recursively). Then choose a random integer r < q (choosing repeatedly) until p = r q + 1 can be declared probably prime by PrimeQ (or your favorite pseudoprime test). Because q is prime and greater than (Sqrt[p]), we have enough of a factorization of p - 1 to apply Pocklington’s Theorem and check that 2 is a witness to p’s primality. Below (10^16) we can simply use PrimeQ, which is known to tell the truth in that domain.

.  The Implementation

Because PrimeQ is known to be 100% correct below (10^16), we can consider primes below that bound to be self-certifying. However, experiments show that it is faster to lower the crossover point to (10^7). The QuickCompositeQ function performs a fast check for compositeness with a 2-pseudoprime test, preceded by a check for small divisors. When looking for the random multiplier r, we take some care (by looking at the mod-3 class of q) to put it in the right mod 6 class so that 1 + 
          q r ≡ ±1  (mod 6), and so has a chance of being prime. We do this by choosing r and then subtracting Mod[r, 6] to get a multiple of 6; then we add the small correction (0 or 4 or 2) based on Mod[q, 3]. But if this fails to yield a prime value of 1 + r q there is really no need to pick another large random number r; it is sufficient to add 6 to r’s current value and try again (and again and again). The use of the decimal point in the bounds in the first Random speeds things up by avoiding rational arithmetic; and the If statement that follows avoids the mod-3 reduction of q whenever the random choice from {0, 1} is 0. Mod[q,3,2] uses the new modular offset feature in version 4; in this case the modulus is 3, so the offset of 2 gives a residue in {2, 3, 4}.

\(primeProd = Times@@Prime[Range[3, 150]]; 


    QuickCompositeQ[n_] := 
      1 < GCD[n, primeProd] < n || PowerMod[2, n - 1, n] ≠ 1; 


    Options[CertifiedPrime] = {ShowCertificate -> False}; 


    CertifiedPrime[d_, opts___] := 
      \(If[\(ShowCertificate /. {opts}\) /. 
            Options[CertifiedPrime], Flatten, First]\)[
        PrimeAndCertificate[d]]; 


    PrimeAndCertificate[d_] := 
      Module[
{p, r, q, 
          q1 = PrimeAndCertificate[Ceiling[d\/2] + 1]}, 
        
q = q1[[1]]; 
        
r = 
          Random[Integer, {10. \^\(d - 1\), 10. \^d}\/q]; 
        
r = r - Mod[r, 6] + adjust; 
		
        While[r = r + 6; p = 1 + q\ r; 
\ \ 
          QuickCompositeQ[p] || GCD[p, PowerMod[2, r, p] - 1] ≠ 1]; 
        
{p, q1}]; 



    PrimeAndCertificate[d_ /; d < 8] := {q\_0 = RandomPrime[d]; 
        adjust = 1 - PowerMod[q\_0, \(-1\), 6]; q\_0}\)

CertifiedPrime[80]

cert = CertifiedPrime[80, ShowCertificate -> True]

And we can easily check that all the proper conditions hold throughout the certificate.

\(Table[{p, q} = 
        cert[[{i, 
            i + 1}]]; 
	\ \ \ \ \ \ \ \ \ \ \ {q > \@p, 
        PowerMod[2, p - 1, p], 
        GCD[p, PowerMod[2, \((p - 1)\)/q, \ p] - 1]}, \ 
	{i, \ 
        Length[cert] - 1}]\)

We can easily get quite large primes.

CertifiedPrime[300]

.  Timing Surprises

Now for a surprise: Maurer’s method is faster than standard methods of generating probable primes, while at the same time having the great advantage that there is no doubt about the status of the output!

We will take several steps to make a fair timing comparison. The randomness in the process makes a true comparison difficult, but the results in this section will provide solid evidence that Maurer’s algorithm is just as fast as, and often faster than, standard methods for getting a random probable prime. We will optimize the CertifiedPrime code in several ways:

• There is no need to keep track of the certificate (the smaller primes q). For security reasons this is sound, since knowledge of the factors of p - 1 can help one factor a number of the form p (Subscript[p, 1]). The prime returned will be a certified prime, but neither the user nor anyone else will be able to check its certificate since it is too difficult to find the large prime factor of p - 1. 

• We replace the recursive method with an iterative method: first we generate the list of digit-sizes that will show up; then we reverse this list and work upward, using Fold with #1 representing q, the certified prime in hand, and #2 representing d, the number of digits.

• We will choose r to be in the right congruence class modulo 30, so that 1 + q r is congruent to one of {1, 7, 11, 13, 17, 19, 23, 29} modulo 30. In fact, if we note the mod-30 class of the first prime generated ((Subscript[q, 0])), then we can arrange that each subsequent q has the same mod-30 class. This is handled with the adjust variable, which is used to adjust the random r upward after it is reduced downward to be 0 (mod 30).

• Quotient is faster than Ceiling.

• The gcd part of the quick test for compositeness is done in two pieces, the first using \(TraditionalForm\`7·11·13·17·19\), thus allowing it to more often work with machine integers.

To be scrupulously fair we should avoid PrimeQ in the generation of the smallest random primes. But this is done only once for each call, and the timing issue on this point is insignificant.

p7To150 = 
        Times @@ Prime[Range[7, 150]];

\(\(CertifiedPrime[d_]\  := \ 
        Module[{q0, adjust, r, 
            p, 
\ \ \ \ \ \ \ \ sizes = 
              Rest[Reverse[
                  FixedPointList[If[# > 10, Quotient[# + 1, 2] + 1, #] &, 
                    d]]]}, 

 (*\ 
            Get\ q0, \ 
            the\ smallest\ of\ the\ certified\ primes\ *) 
q0 = 
            Random[Integer, \ {10. \^\(First[sizes] - 1\), \ 
                10. \^First[sizes]}]; 
q0 += \ 
            1 - Mod[q0, 2]; 
While[\(! PrimeQ[q0]\), 
            q0 += 2]; 

If[d > 10, 
            adjust = 
              1 + PowerMod[\(-q0\), \(-1\), 
                  30030]]; 
Fold[\((\ \ r = \((\(\((# + 
                            Mod[\(-#\), 30030] + adjust)\) &\)\ [
                      Random[Integer, {10. \^\(#2 - 1\), 10. \^#2}\/\(\(\ \)\(#1\)\)]])\); 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ While[
                  p = 1 + #1\ *r; 
\((1 < GCD[p, p7To150] < 
                        p)\)\  ∨ 
\ \((PowerMod[2, 
                          p - 1, p] ≠ 
                        1)\) ∨ \ 
\((GCD[p, 
                          PowerMod[2, r, p]\  - \ 1] ≠ 1)\), 
                  r += 30030]; 
p)\) &\ , q0, 
            Rest[sizes]]];\)\)

CertifiedPrime[100] // Timing

These optimizations yield a nicely fast routine. Here is the average time needed for a 100-digit prime, based on 1000 trials.

trials = 1000;
Plus @@ Table[Timing[CertifiedPrime[100]][[1]], {trials}]/trials

If we wish to compare this with standard methods of generating random probable primes (sometimes known as “industrial grade primes”), we have several choices (see [Wagon 1997, 1999] for the definitions).

1. Use the built-in PrimeQ, which calls n prime if it passes the 2- and 3-strong pseudoprime tests and a Lucas pseudoprime test.

2. Use a method that uses only a single 2-strong psp test.

3. Use 2- and 3- strong psp tests.

4. Use 2-, 3-, and 5-strong psp tests.

Method 2 is surely the fastest and will be adequately reliable for large integers, where the probability of running into a 2-strong psp is very small. Methods 3 and 4 add some extra security, while method 1 is probably the best since not a single counterexample is currently known.

Coding the strong pseudoprime test gives us an opportunity to show off two new and useful version 4 functions: NestWhile and IntegerExponent. To perform a b-strong pseudoprime test on n one first forms s = (b^m) (mod n) where m is the odd part of n - 1. Then one squares repeatedly (mod n) until the exponent reaches n - 1. If n - 1 shows up before the penultimate step, n passes the test; if it does not, n fails and is definitely composite. Thus we wish to iterate the modular squaring operation, exiting the Nest sequence as soon as a condition is met. There was no simple way to do this in version 3. But now we can use:

n - 1 == 
        NestWhile[PowerMod[#, 2, n]&, s, # ≠ n - 1 &, 1, r - 1]

The third argument tells us under what condition to continue the iteration; the fourth argument gives the maximum number of iterations. If the result of this is n - 1, the number passes; otherwise it fails.

So here is code that returns True if n passes the b-strong pseudoprime test, followed by the PrimeAfterSPSP and RandomPrimeSPSP functions that make use of the StrongPseudoprimeTest code. We split the gcd step into two as before.

\(\(Attributes[StrongPseudoprimeTest] = 
        Listable;\)
\)
\(\(StrongPseudoprimeTest[b_Integer, n_]\  := \ 
        Module[
{r = IntegerExponent[n - 1, 2], s}, 
	
          s = PowerMod[b, \(n - 1\)\/2\^r, n]; 
	
          s == 1 ∨ 
            s == n - 1 ∨ 
	\ \ \ \ \ n - 1 == 
              NestWhile[PowerMod[#, 2, n] &, s, # ≠ n - 1\  &, 1, 
                r - 1]];\)\)

Here are some examples; 2047 is the first composite that passes the 2-strong pseudoprime test; 1373653 is the first composite that passes both the 2- and 3-strong psp tests.

StrongPseudoprimeTest[
      2, {101, 541,  7919,  2047,  1373653}]

But 2047 is not a 3-strong pseudoprime.

StrongPseudoprimeTest[
      3, {101, 541,  7919,  2047,  1373653}]

Now we modify the PrimeAfter and RandomPrime code so that the 2-strong pseudoprime test is used.

PrimeAfterSPSP2[n_] := 
        Module[{j = n + 1 + Mod[n, 2]}, 
While[
            1 < GCD[j, 15015] < j ∨ 
              1 < GCD[j, p7To150] < 
                j  ∨  
¬ StrongPseudoprimeTest[2, j]  , 
j += 
              2]; 
j];

RandomPrimeSPSP2[d_]  :=  
        Module[{n}, 
          While[n = 
              PrimeAfterSPSP2[
                Random[Integer, 
                  (10^(d - {1, 0}))]];  
  Log[10. , 
                n] > d];  n];

Plus @@ Table[Timing[RandomPrimeSPSP2[100]][[1]], {trials}]/trials

The timing experiment shows that the generation of certified primes takes essentially the same amount of time as the generation of probably primes. And when more sophisticated probable prime generators are used, such as PrimeQ, then the certified-prime algorithm takes less time!

References

[Bressoud and Wagon 2000]  D. Bressoud and S. Wagon, A Course in Computational Number Theory, Springer-Verlag, New York.
[Maurer 1995] U. Maurer, Fast generation of prime numbers and secure public-key cryptographic parameters, Journal of Cryptology 9, 123–155.
[Wagon 1997] S. Wagon, What is a prime number?, Mathematica in Education and Research 6:4, 54–61.
[Wagon 1999] S. Wagon, Mathematica in Action, Springer/TELOS, New York.

About the Author

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition, and is especially appreciative of the new ways Mathematica allows us to look at mathematical objects. He is author or coauthor of several books, including the just-released second edition of Mathematica in Action, Animating Calculus, VisualDSolve, and Which Way Did the Bicycle Go?. He teaches a high-altitude course, Rocky Mountain Mathematica, every summer in the mountains of Colorado.

Stan Wagon
Department of Mathematics and Computer Science
Macalester College
St. Paul, MN 55105
wagon@macalester.edu