PROGRAMMING  TIPS

What Is a Prime Number?

The notion of primality is, in many respects, elusive. In this column we describe several ways in which Mathematica can help us answer the age-old question, “Is n prime?”

by Stan Wagon

Although it is one of the most elementary concepts of mathematics, the notion of primality still holds many secrets. In this article we will delve into several algorithms that determine primality. It is important that such an algorithm be fast, and for that one must give up on the notion of absolute correctness. Still, there is an algorithm that is very fast and is perfect for all integers up to (10^16), and possibly much farther. This algorithm is used by the built-in PrimeQ command and we will discuss it in detail here. Along the way we will learn about efficient approaches to linear recurrences, as well as ways to deal with the situation where absolute certainty regarding a number's primality is required.

.  Probable Primes

Prime numbers have many beautiful properties. One of the most famous is that, if p is an odd prime, then (2^(p - 1)) ≡  1 (mod p).This is known as Fermat's little theorem and is in fact more general: if gcd(a, p) = 1, then (a^(p - 1)) ≡  1 (mod p). The proof is simple and can be found in any number theory textbook. Because modular powers can be computed very quickly (one uses the base-2 digits of p – 1 to break the problem down into steps; by reducing mod p at each step the numbers stay small), this is a fine initial test for primality.

p = Prime[10000]; PowerMod[2, p - 1, p]

A composite, odd n for which (2^n)–1 ≡  1 (mod n) 
      is called a 2-pseudoprime. The following shows that there are none among the first 100 odd integers beyond one million.

Select[Range[(10^6) + 1, (10^6) + 200, 2], PowerMod[2, # - 1, #] == 1 &]
         == Select[Range[(10^6) + 1, (10^6) + 200, 2], PrimeQ]

Of course, 2-pseudoprimes do exist. There are 23 of them under 10000.

\(Select[Range[3, 10000, 2], 
	
      \(! PrimeQ[#]\)\  && \ PowerMod[2, \ # - 1, #] == 1\ &]\)

Similarly one can define b-pseudoprimes to be odd composite n for which (b^n) ≡ 1 (mod n).

\(Select[Range[3, 10000, 2], 
	
      \(! PrimeQ[#]\)\  && \ PowerMod[13, \ # - 1, #] == 1\ &]\)

Pseudoprime tests are very fast, and so they are perfect for an initial investigation into the primality question for a large integer n. But the failure rate is too great and so more sophisticated methods must be used for something as important as PrimeQ. Note that 561 (= 3·11·17) is both a 2-psp and a 13-psp. In fact, 561 is a b-pseudoprime for every b relatively prime to 561. Such numbers are called Carmichael numbers (it is now known that there are infinitely many of them) show that primality cannot be verified by just trying lots of b-pseudoprime tests.

This sort of test has nothing to do with probability theory. If a number passes several pseudoprime tests then one might say that it is probably prime, but the procedure is totally deterministic. Sometimes randomness is introduced into the choice of b, but then different runs might yield different results. To ensure repeatability, the built-in prime tester makes no use of random numbers.

. Strong Pseudoprimes

A very pretty and easy enhancement of the pseudoprime test comes by considering square roots. This leads to the notion of strong pseudoprime. Suppose (b^(n - 1)) ≡ 1 (mod n) where n is odd. Write n  -  1 as (2^s) m and observe that (b^(n - 1)) is the last term of the sequence (b^m), (b^(2  m)), . . . , (b^((2^s)m)), where each term of the sequence is obtained by squaring its predecessor. But if n is indeed prime then not only is the last term equal to 1 (modulo n), but the last term that is not 1 must be -1. This is because the only square roots of +1 modulo a prime p are ±1 (proof: \(TraditionalForm
      \`x\^2\  - \ 1\  = \ \((x\  - \ 1)\) \((x\  + \ 1)\)\), and if the right side is 0, then p divides it, in which case p divides one of the two factors so x ≡ ±1). So if p is prime the b-sequence just given either consists entirely of 1s or looks like { . . . , –1, +1, . . . , +1}. Call n a b-strong pseudoprime (b-spsp) if its b-sequence looks like this but n is composite. Here is code to find such numbers. First we isolate an odd-part routine that returns {m,  s} as just defined.

OddPart[n_] := 
      Module[{m = n, s = 0}, While[EvenQ[m], m /= 2; s++]; {m, s}]

OddPart /@ {100, 101, 102}

Now we can look at some b-sequences. In the code that follows NestList is used to iterate the mod-n squaring function the number of times specified by oddPart. For legibility we replace n  -  1 by -1.

\(spspSequence[b_, n_] := 
      Module[{m, s}, 
		{m, s} = OddPart[n - 1]; 
		
        NestList[Mod[#\^2, n]&, PowerMod[b, m, n], s] /. n - 1 -> \(-1\)]
        \)

When we look at the 1000th or 10000th prime we see the proper behavior.

spspSequence[2, 7919]

spspSequence[2, 104729]

The smallest 2-pseudoprime, 341, is quickly unmasked by this test, since the 2-sequence has a 1 preceded by a 32.

spspSequence[2, 341]

Here is complete code to recognize b-strong pseudoprimes. We treat n = 1 as a special case. We need only check that –1 is a member of the b-sequence with its last entry dropped, for its presence anywhere guarantees that the sequence has the form typical of primes.

\(StrongPseudoprime[b_, n_] := 
      Module[{bSeq, m, s}, 
		{m, s} = OddPart[n - 1]; 
		
          bSeq = NestList[PowerMod[#, 2, n]\ &, PowerMod[b, m, n], s]; 
		
          Union[bSeq] == {1} ∨ MemberQ[Drop[bSeq, \(-1\)], n - 1]] /; 
        n > 1
\)StrongPseudoprime[_, 1] = False;

We can now select all the 2-spsps below 10,000, where the term 2-spsp is used for an odd integer that passes the 2-strong pseudoprime test but is not in fact prime.

Select[Range[10000], 
      ! (PrimeQ[#]) ∧ 
          StrongPseudoPrime[2, #] &]

We see that this test too is fooled by some rather small integers. But in this situation the combination of several strong pseudoprime tests leads to a very strong test. For example, the first integer that is a b-strong pseudoprime for b = 2, 3, 5, and 7 is 3,215,031,751 [PSW]

Daniel Bleichenbacher [Bleichenbacher, 1996] has made a detailed investigation of integers that are b-spsps for several b. His largest example is 68528663395046912244223605902738356719751082784386681071, which is a b-spsp for all b ≤ 100. Here is the proof, where the psp is given in factored form.

Table[
      StrongPseudoPrime[i, 
        18215745452589259639*4337082250616490391*867416450123298079], {i, 
        101}]

In such a case we say that 101 is the first witness to the compositeness of n. Unlike the case for pseudoprimes, all composites have a first witness that is relatively small. Well, this is under the assumption of the extended Riemann hypothesis, but the evidence for the truth of this statement is good. More precisely, under ERH every composite integer n is not a b-strong pseudoprime for some b less than 2(log ((n))^2). In short, ERH => “primes are in 𝒫”, where 𝒫 is the class of problems solvable in polynomial time. See [Bach and Shallit, 1996].

Mathematica's PrimeQ function starts off (after some checking for small divisors) by using a 2-spsp test. If the input fails the test, it is proved composite. If it passes, then a 3-spsp test is used, and then one other test, to be described in section 4.

3. Perrin's Pretty Primality Test

As a warmup to the more important Lucas test, we discuss a simple and very pretty test. Consider the recurrence defined by the code that follows; it was in fact studied by Lucas, but it is now known as the Perrin recurrence, after R. Perrin who examined it in 1899. The code that follows serves as the mathematical definition, but also works as a recursive definition in Mathematica. In such cases it is essential to cache the values as they are computed, and this is done by the (Subscript[a, n]) = (Subscript[a, n - 2]) + (Subscript[a, n - 3]) phrase. While subscripted variables do not work in all situations, they can be used here, as opposed to the more familiar a[n_] setup.

\(a\_0 = 3; a\_1 = 0; a\_2 = 2; 

    a\_n_ := \(a\_n = a\_\(n - 2\) + a\_\(n - 3\)\)\){(Subscript[a, 0]),  (Subscript[a, 1]),  (Subscript[a, 2]),  (Subscript[a, 3]),  (Subscript[a, 4]),  (Subscript[a, 5])}

Note that to clear the definition one would use Clear[Subscript]. In order to allow a to work directly as a function too, it is convenient to add the following two lines.

a[n_] := (Subscript[a, n])Attributes[a]  =  Listable;

Here are the first 20 values.

a[Range[0, 20]]

Start from the left and count 0, 1, 2, 3, . . .  . You will notice that at the primes, and only there, the a-value is divisible by the index (with the exception of (Subscript[a, 1])). We can check it.

\(pr\  = \ Prime[Range[20]]; 
a[pr]/pr\)

For a proof, we will need some elementary facts about linear recurrences. If p(x) is the characteristic polynomial of the recurrence — p(x) = (x^3)– x – 1 in this case — then (Subscript[a, n])is equal to \(TraditionalForm
      \`\(c\_1\) α\^n\  + \ \(c\_2\) β\^n\  + \ 
        \(c\_3\) γ\^n\), where α, β, and γ are the three roots of the polynomial and the coefficients (Subscript[c, i]) are determined by the initial conditions. The roots of a cubic are easy to find and we can then solve for the coefficients. We suppress the large symbolic output of these two steps.

\({α, β, γ} = x /. Solve[x\^3 - x - 1 == 0, x]; 

    soln = Solve[
        Thread[
			\ \ {3, 0, 2} == 
            c\_1\ α\^{0, 1, 2} + c\_2\ β\^{0, 1, 2} + 
              c\_3\ γ\^{0, 1, 2}], {c\_1, c\_2, c\_3}]; \)

Checking the numerical values shows that that the coefficients are surely {1, 1, 1}.

soln //N //Chop

Well, it is easy to see that setting each (Subscript[c, i])= 1 does indeed work. Setting n = 0 yields (Subscript[c, 1])+ (Subscript[c, 2]) + (Subscript[c, 3]) = 3, which is satisfied. Setting n = 1 yields  0  =  (Subscript[a, 1]) =  (Subscript[c, 1])α + (Subscript[c, 2]) β + (Subscript[c, 3])γ = α + β + γ, which holds because 0 is the coefficient of (x^2) in the polynomial. And it is easy to see, by expanding (α + β + ((γ))^2), that  (α^2)+ (β^2)+ (γ^2) = 2. So we have that (Subscript[a, n])  =  (α^n)  +  (β^n)  +  (γ^n). In fact, for n ≥ 10, it is true that (Subscript[a, n]) is simply the nearest integer to (α^n), where α is the real root, which can be simplified to the following form.

\(α = 
      \(\((\@27 - \@23)\)\^\(1/3\) + 
          \((\@27 + \@23)\)\^\(1/3\)\)\/\(2\^\(1/3\)\ \@3\); 
    
{a[Range[9, 20]]\ , Round[α\^Range[9, 20]]}\)

Now, if p is prime, then  (α^p)  +  (β^p)  +  (γ^p) and (α + β + ((γ))^p) are in the subring of algebraic integers of the number field determined by the roots of p(x). Therefore the binomial theorem and elementary facts about divisibility of (({{p}, {i}})) by p tell us that  (α^p)  +  (β^p)  +  (γ^p) ≡ (α + β + ((γ))^p) (mod p) in this ring. But α  +  β  +  γ = 0, so (α^p)  +  (β^p)  +  (γ^p) has the form pρ where ρ is in the aforementioned ring.  But (α^p)  +  (β^p)  +  (γ^p) is an integer — it equals (Subscript[a, p])— and ρ is an integer linear combination of elements of an integral basis (see [Bach and Shallit, 1996, §6.7]) so ρ must in fact be an integer too. Therefore (Subscript[a, p])is an integer divisible by p.

Whenever one has a simple formula that holds for primes, it is natural to ask if it holds for any composites. Call n a Perrin pseudoprime if n is composite and n divides (Subscript[a, n]). They exist, but they are not easy to find. Indeed, none was known until 1982, when Adams and Shanks found the first one. Here are the results of a quick search up to 100.

\(Select[Range[2, 100], 
      \(! PrimeQ[#]\) && Mod[a\_#, #] == 0\ &]\)

The recursive definition makes even a modest search up to 1000 difficult because of recursion limits. One could program the definition in a Do-loop. But it is much more efficient to take an approach based on matrices. First we digress to discuss how to raise matrices to very large powers, and that requires a digression on the raising of integers to very large powers modulo n. Recall that modular powers for integers are very fast when PowerMod is used.

PowerMod[1007, 34000000000, 100001]

The algorithm underlying this is simple. As an example, suppose the exponent is 26. Write 26 in base 2 and view the goal, (b^n), as being (b^11010). Now start with 1 and scan the base-2 digits from left to right. When a 0 is read, square the current value modulo m. When a 1 is read, square and then multiply by b, again reducing modulo m. The sequence will look like the following, where the exponent is in base 2: (b^0), (b^1), (b^10), (b^11), (b^110), (b^1100), (b^1101), (b^11010).

The fact that a modular reduction takes place after each digit is scanned means that the numbers never get large. If we wished to implement this ourselves, here is a fast way to do it. Fold is used to scan the base-2 digits and do the right thing.

\(PowerModManual[b_, n_, m_] := 
      Fold[Mod[#1\^2*If[#2 == 1, b, 1], m]\ &, 
		1, IntegerDigits[n, 2]]
        \)

{PowerModManual[123, 456, 789], PowerMod[123, 456, 789]}

Now it is simple to generalize the idea to modular powers of integer matrices. We start with mat, as opposed to the identity matrix.

\(MatrixPowerMod[mat_, \ n_, \ m_]\  := \ 
      Fold[
\ \ Mod[If[#2\  == \ 1, #1.#1.mat, #1.#1], m]\ &, 
\ \ mat, 
        Rest[IntegerDigits[n, 2]]]\)

Here is an example that uses the matrix  (({{0, 1}, {1, 1}})), which defines the Fibonacci recurrence, in that (({{0, 1}, {1, 1}}))(({{a}, {b}})) = (({{b}, {a + b}})). Using ∞ as the modulus means that we get the unreduced numbers. Taking the 100th power should yield the 100th Fibonacci number.

MatrixPowerMod[{{0, 1}, {1, 1}}, 100, ∞]

Fibonacci[100]

Now we can formulate the Perrin recurrence in terms of multiplication by (({{0, 1, 0}, {0, 0, 1}, {1, 1, 0}})). Here, for example, is how to get (Subscript[a, 7919]) (mod 7919). Because 7919 is prime, we expect 0.

Mod[First[MatrixPowerMod[(({{0, 1, 0}, {0, 0, 1}, {1, 1, 0}})),7919,7919]].{3, 0, 2},7919]

Of course, we can formulate this method for any linear recurrence. The LinearRecurrence function that follows accepts the coefficients of the recurrence ({1, 1, 0} for Perrin's rule), the initial condition ({3, 0, 2} for Perrin), the desired index n, and a modulus m (which defaults to ∞, giving the unreduced numbers). It forms the matrix, takes the mod-m power efficiently, and picks out the correct element of the matrix product with the initial vector. We also arrange the function to handle a list in its third argument.

LinearRecurrence[coeffs_, init_, n_Integer, m_:∞] := 
  Mod[First[MatrixPowerMod[
    Append[RotateRight /@    
      Drop[IdentityMatrix[Length[coeffs]], -1],
    coeffs], n, m]] . init, m] /; n ≥ Length[coeffs]

LinearRecurrence[coeffs_, init_, n_Integer, m_:∞] := Mod[
init[[n+1]],m] /; n < Length[coeffs]		

LinearRecurrence[coeffs_, init_, n_List, m_:∞] := 
  LinearRecurrence[coeffs, init, #, m]& /@ n

As a check we look at the millionth prime.

LinearRecurrence[{1, 1, 0}, {3, 0, 2}, Prime[(10^6)], Prime[(10^6)]]

And now we can search for Perrin pseudorpimes.

Select[Range[2, 500], 
      ! PrimeQ[#] && 
          LinearRecurrence[{1, 1, 0}, {3, 0, 2}, #, #] == 0 &]

The first one is quite large. Since I know where it is, I can focus the search very well!

\(Select[Range[271200, 271500], 
      \(! PrimeQ[#]\) && 
		
          LinearRecurrence[{1, 1, 0}, {3, 0, 2}, #, #] == 0\ &]\)

It is striking, especially when compared to other pseudoprime tests, how large the first failure of the Perrin test is. We can see that 271441, which is (521^2), is not a 2-pseudoprime:

PowerMod[2, 271440, 271441]

So we might be tempted into thinking that Perrin + 2-psp has the potential of being a good primality test. But in fact there are many Perrin pseudoprimes that are even 2-strong pseudoprimes. The first such is 27,664,033.

\(n = 3037*9109; 
LinearRecurrence[{1, 1, 0}, {3, 0, 2}, n, n]\)StrongPseudoprime[2, n]

True

One can go farther with Perrin's idea. For example, one can run the recurrence backwards as well as forwards and discover that (Subscript[a, -p]) ≡ 1 (mod p)when p is prime. This and other variations, as well as references to papers with more detailed investigations, can be found in chapter 18 of [Wagon, 1998]. Such ideas lead to a Perrin test that is even stronger, in the sense that the first failure is quite large. But let us leave Perrin and move on to Lucas.

. The Lucas Primality Test

Mathematica's PrimeQ uses the 2- and 3-strong pseudoprime tests together with a Lucas test that we will now describe. No composite is known that passes these three tests, and all cases up to (10^16) have been checked.

A Lucas sequence is the sequence of integers obtained from the Fibonacci-lik recurrence (Subscript[V, 0])  =  2, (Subscript[V, 1])  =  P, (Subscript[V, n])  =  P (Subscript[V, n - 1]) -  Q (Subscript[V, n - 2]), where Δ, defined to be (P^2) - 4  Q, is congruent to 0 or 1 (mod 4) and is not a perfect square. Since the special case Q = 1 will be our primary interest, we make that the default. As with the Perrin sequence, the V-sequence is related to the roots of a polynomial: if α and β are the roots of the quadratic (x^2) - P x + Q , then (Subscript[V, n])= (α^n) + (β^n).

\(\(Lucas[P_, Q_:  1]\)[0]\  := \ 2; 
\(Lucas[P_, Q_:  1]\)[1]\  := \ P; 
    
\(Lucas[P_, Q_:  1]\)[n_]\  := \ 
      \(\(Lucas[P, Q]\)[n]\  = 
		\ 
        P\ \ \(Lucas[P, Q]\)[n - 1]\  - \ Q\ \(Lucas[P, Q]\)[n - 2]\)\)

Because of recursion limits, the preceding definition will not work well except for quite small values of n. We will use the matrix method for serious work. We next summarize the facts that lead to an efficient prime test using the V-sequence.

THEOREM 1. Suppose p is prime, P is relatively prime to p,  Δ is not a square modulo p (where Δ = (P^2) - 4 ), and {(Subscript[V, n])}
is the Lucas sequence that starts with (2, P).

(a).  (Subscript[V, p + 1]) ≡ 2 (mod p)

(b).  (Subscript[V, 2  i])  =  (Subscript[V, i]^2)  -  2

(c).  (Subscript[V, 2  i + 1]) = (Subscript[V, i]) (Subscript[V, i + 1]) - P

(d).  2  (Subscript[V, i + 1])  =  P (Subscript[V, i]) + Δ (Sqrt[((Subscript[V, i]^2)  -  4)/Δ])

We will not prove these identities; the reader can find proofs in [Bressoud, 1989, chap. 12]. These identities are, in essence, equivalent to the following result, a quadratic field analog of Fermat's little theorem that puts the primality test we are about to describe on a basis very similar to that of ordinary pseudoprimes (see [Bleichenbacher, 1996]).

THEOREM 2. Let p, P, and Δ be as in Theorem 1. Let α be the larger root of (x^2) - P x + 1  =  0; α = (^1)\(/\_2\)(P + (Sqrt[Δ])). Then (α^n)≡  \(TraditionalForm\`α\&_\ \ \((mod\ p)\)\).

Our first task is to define a routine to get the least P that works. All the work that follows will be restricted to odd inputs. Further, if P = ±1 or ±2, then (Subscript[V, n]) is a periodic sequence, so we start our search for P at 3.

LucasParameter[n_?OddQ] := 
      Module[{P = 3}, 
        While[JacobiSymbol[(P^2) - 4, n] == 1 || GCD[n, P] ≠ 1, 
          P++]; P]

We can use this routine to illustrate Theorem 2. In the code that follows we use PowerMod to rationalize the 2 in the denominator. But because (Sqrt[k]) has the internal form Power[k, Rational[1, 2]], we must hide the square root from the rule that transforms rationals. Then we add the square root back at the end.

\(n = 61; P = LucasParameter[n]; Δ = P\^2 - 4; 

    α = 
      x /. Last[
          Solve[x\^2 - P\ x + 1 == 0, x] /. P\^2 - 4 -> Δ]
          \)

ExpandAll[{α, (α^n) /. (Sqrt[21]) :> -(Sqrt[21])}]
             /. (Sqrt[21]) :> hide /. 
        Rational[m_, 2] :> PowerMod[2,  -1,  n] /. 
      hide -> (Sqrt[21])

The reader can check some small odd composite values of n to see that the equation of Theorem 2 fails (to do this, add "/. m_Integer :> Mod[m, n]" to the end of the code). The Lucas test that follows is essentially a computationally efficient was of checking the congruence of Theorem 2.

The value of P that LucasParameter gets will typically be quite small.

Union[LucasParameter /@ Range[3, 1000, 2]]

Now we can verify the Lucas condition (Theorem 1(a)) for primality in the case of 7.

\(n\  = \ 7; 

    ls = \(\(Lucas[LucasParameter[n]]\)\ [#]&\)\ /@\ Range[0, n + 1]; 

    DisplayForm[GridBox[{Range[0, \ n + 1], ls, Mod[ls, n]}]]\)

Code for GridBox below

\(n\  = \ 7; 	
    ls = \(\(Lucas[LucasParameter[n]]\)\ [#]&\)\ /@\ Range[0, n + 1]; 

    DisplayForm[
      StyleBox[GridBox[

          			{Prepend[Range[0, \ n + 1], "n"], 
				
            Prepend[ls, "V"], 
	\ \ \ \ 	
            Prepend[	Mod[ls, n], \ "V (mod 7)"]}, 
	GridFrame -> 1, \ 
          RowLines -> 1, ColumnAlignments -> {Center, Right}, \ \ 
          ColumnLines\  -> \ {1, 0}], \ FontSize -> 9, \ 
        Background -> GrayLevel[0.85]]]\)

And here is the last term of the V-sequence for the first 25 odd primes.

Mod[Lucas[LucasParameter[#]] [# + 1], #] & /@ 
      Prime[Range[2, 26]]

But much more is true, and some ideas of the theory of strong pseudoprimes can be used [Baillie and Wagstaff, 1980]. Theorem 1(b) means that we can get (Subscript[V, n + 1]) by first computing (Subscript[V, s]) where s is the odd part of n+1, and then repeatedly squaring and subtracting 2 mod-n. Moreover, if n is really a prime and (V^2) -  2 ≡  2 (mod n), then (V^2)  ≡  4, and this means that V must be one of ±2 (mod n); this is analogous to the situation in the strong pseudoprime test, where the fact that the modular square root of 1 is one of ±1 is used to great advantage. Further, the V-values that occur right after those just mentioned are important too: still assuming p is prime, if (Subscript[V, 2  i]) ≡  ±2, then \(TraditionalForm
      \`\(V\_\(2  i + 1\)\  ≡ \ \(±P\)\ \)\); this follows from Theorem 1(d). And, finally, these next V-values are easily computed from Theorem 1(c) — (Subscript[V, 2  i + 1]) = (Subscript[V, i]) (Subscript[V, i + 1]) - P — so we can easily check this condition on the next V-value.

Now we can form a strong Lucas pseudoprime test as follows:

     1. Given n, let P be the Lucas parameter as defined earlier.
     2. Write n + 1 as s (2^r), where s is odd.
     3. Get (Subscript[V, s])and (Subscript[V, s + 1]) modulo n by the matrix method.
     4. Form the two V-sequences indexed by {s, 2s, 4s, 8s, . . . , n+1} and {s+1, 2s+1, 4s+1, 8s+1, . . . , n+2} modulo n.
     5. Check that (Subscript[v, n + 1]) = 2, and that either the first sequence consists of all 2s or the first 2 in the sequence is preceded by a -2.
     6. Check that the first 2 or -2 in the first sequence is matched by a P or -P in the next sequence.

It would be more efficient to not compute the whole V-sequence, but to check the condition term-by-term. But, in part because we wish to look at the whole sequences here, we simply follow that approach. Exercise: Program the Lucas test more efficiently, checking each term of the sequence rather than precomputing the entire sequence. Here is code that implements these ideas. The lspspSequence function uses LinearRecurrence to get started and then NestList to continue, and returns two lists of V-values, the first indexed by {s, 2s, 4s, 8s, . . . , n+1} and the second indexed by {s+1, 2s+1, 4s+1, 8s+1, . . . , n+2}.

\(lspspSequence[P_, n_?OddQ] := 
      Module[{oddData = OddPart[n + 1]}, 
		
        Transpose[
          NestList[
              Mod[{#[[1]]\^2 - 
                      2, #[[1]]\ #
                          [[2]] - P}, 
                  n]&, 
				\ 
              LinearRecurrence[{\(-1\), P}, {2, P}, 
                oddData[[1]] + {0, 1}, 
                n], oddData[[2]]] /. {
              n - 1 -> \(-1\), n - 2 -> \(-2\)}]]\)

Here is an example. The parameter P is 5 and a -2 shows up with the corresponding value in the sequence of next terms being the expected 7919  -  5, or 7914.

\(P\  = \ LucasParameter[7919]; 
    
{P, lspspSequence[LucasParameter[7919], 7919]}\)

Now we can write some code to check the strong Lucas pseudoprime conditions in general. We will use the term strong Lucas pseudoprime do describe a composite integer that passes the test as just described, but the reader should be aware that the terminology in the literature is not standard. Others have used the term QF-psp and strong QF-psp to emphasize that this idea is a pseudoprime construction in a certain quadrtic number field.

\(StrongLucasPSP[n_] := 
      Module[{P = LucasParameter[n], lSeq, pos}, 
		
        lSeq = lspspSequence[P, n]; 
		
        lSeq[[1, \(-1\)]] == 2 ∧ 

          		\((lSeq[[{1, 2}, 1]] == {
                  2, P} || 
					\ 
              \((pos = 
                  Position[lSeq[[1]], 
                    \(-2\)]; 
						
                pos ≠ {} ∧ 
                  lSeq[[2, 
                      pos[[1, 1
                        ]]]] == n - P)
                \))\)]\)

We first check that a dozen random primes under 10 million pass the test.

12 == 
      Length[Select[Table[Prime[Random[Integer, {2, (10^7)}]], {12}], 
          StrongLucasPSP]]

And now we search for strong Lucas pseudoprimes according to this definition. We avoid a Select[Range[]] construction since that requires forming a huge range; instead we use a Do-loop. We add a Print statement to monitor progress since this is a time-consuming computation, and such a trace eases the task of stopping and restarting if desired.

\(LucasPSPs = {}; 

    Do[If[Mod[i + 1, 2000] == 0, Print[{i, LucasPSPs}]]; 
	\ \ \ \ \ 
      If[GCD[i, 6] == 1 ∧ \(! PrimeQ[i]\) ∧ 
          StrongLucasPSP[i], 
		\ \ \ \ \ \ AppendTo[LucasPSPs, i]], {i, 
        5, 100000, 2}]; \)

LucasPSPs

So there are 12 strong Lucas psps under 100,000. But they all would have their compositeness detected by a 2-spsp test! Indeed, even the simple 2-psp test suffices.

(1 == PowerMod[2, # - 1, #] &) /@ LucasPSPs

Now, what the built-in PrimeQ function does on input n is, in essence, the following:

     1. If n has a small divisor, return False.
     2. If n is proved to be composite by the 2-strong pseudoprime test, return False.
     3. If n is proved to be composite by the 3-strong pseudoprime test, return False.
     4. If n is proved to be composite by the strong Lucas pseudoprime test, return False.
     5. Return True.

So, if PrimeQ declares that n is prime, is n definitely prime? Good question! The common belief is that such a combination of simple tests will not truly characterize primality (note that if such a finite combination of tests did characterize primality, that would mean that the primes are in the complexity class 𝒫). On the other hand, Daniel Bleichenbacher [Bleichenbacher, 1996] has computed a list of all composite integers under (10^16) that are both 2- and 3-strong pseudoprimes; there are 52,593 of them and he checked that PrimeQ is correct (i.e., False) for all of them. In other words, under (10^16) there is no counterexample to the assertion that primality is determined by the conjunction of the 2-spsp test, the 3-spsp test, and the strong Lucas pseudoprime test. Mathematica's PrimeQ is actually a little different than the algorithm described here in that a slightly different method is used for the choice of Lucas parameter. But I checked that all the 2- and 3-spsps under (10^16) fail the Lucas test as presented here too. So, for either method of choosing the Lucas parameter, we can be sure that the prime-query algorithm based on 2-spsp, 3-spsp, and the Lucas test tells the truth up to (10^16).

. Certifiable Primality

References

[Bach and Shallit, 1996] Eric Bach and Jeffrey Shallit, Algorithmic Number Theory, Volume 1, Efficient Algorithms, MIT Press, Cambridge.

[Baillie and Wagstaff, 1980] R. Baillie and S. S. Wagstaff, Lucas pseudoprimes, Math. Comp. 35, 1391–1417.

[Bleichenbacher, 1996]  D. Bleichenbacher, Efficiency and Security of Cryptosystems based on Number Theory, Ph. D dissertation, Swiss Federal Institute of Technology, ETH Diss. No. 11404, Zürich.

[Bressoud, 1989]  D. M. Bressoud, Factorization and Primality Testing, Springer, New York.

[Schlafly and Wagon, 1994]  A. Schlafly and S. Wagon, Carmichael's conjecture is valid below (10^(10, 000, 000)), Mathematics of Computation 63, 415–419.

[Wagon, 1998] Stan Wagon, Mathematica in Action, Springer/TELOS, New York.