Is there a PrimeQ whose accuracy guarantee you can adjust?

Question

Say I have a list of a million integers each with a million digits, and I want a crude sieve to see which have a chance at being prime.

Mathematica has a PrimeQ function, which appears to be slow because it performs a large number of trials of Miller-Rabin's primality test and the Lucas primality test before giving a Yes/No answer.

Now I could try looking at small factors first, but say this only removes a small fraction of my list. Is there a way to turn down the certainty of PrimeQ, or are there other, faster, tuneably-cruder ways to test for primality within Mathematica?

Answer 1

I'd like to expand on the answer given by @ChipHurst in the comments. My hope was that one of these approximate tests would be much faster than PrimeQ; however, that was not the case. Perhaps someone can code these methods more efficiently than shown here.

Fermat's Little Theorem: Prime p and any b such that GCD[b,p]=1 implies b^(p-1)=1, mod p. Use a prime base b with the Fermat test.

FermatPrimeQ[b_, n_] := Thread[PowerMod[b, n - 1, n] == 1]

Run Fermat's prime test with the first m primes as bases. If n passes for all bases, then return True. The larger m, the larger the probability that n is prime.

FermatLittleTheoremTest[n_?OddQ, m_] :=
   Block[{b = DeleteCases[Prime[Range[m]], n], k = 1, i = 1, len},
         len = Length[b];
         While[k <= len && (PowerMod[b[[i]], n - 1, n] == 1), k += 1; i += 1];
         len == k - 1]

Bressoud and Wagon, A Course in Computational Number Theory, define QuickPrimeQ[n], where primeProduct is the product of the first 100 primes. They combine a GCD test with a base b=2 Fermat test. Other tests also benefit from incorporating a GCD test.

primeProduct = Apply[Times, Prime[Range[2, 100]]];
QuickPrimeQ[n_] := GCD[primeProduct, n] == 1 && PowerMod[2, n - 1, n] == 1

A stronger test is described in Section 4.2 of Bressoud and Wagon. The following is an edit of the code given on page 116.

SpspSequence[b_, n_] :=
   Module[{s = IntegerExponent[n - 1, 2]},
      NestList[Mod[#^2, n] &, PowerMod[b, Quotient[n - 1, 2^s], n], s] /. n-1 -> -1]

StrongPseudoprimeTest[b_, 2] := True
StrongPseudoprimeTest[b_, n_?EvenQ] := False
StrongPseudoprimeTest[b_, n_?OddQ] :=
    (Union[#] == {1} || MemberQ[#, -1]) &[SpspSequence[b, n]]

They state on page 118 that "the probability of a composite integer passing the strong pseudoprime test using m random bases is at most 1/4^m", and "is usually much smaller than this".

StrongPseudoprimeProbability[n_?OddQ, m_] :=
   If[1 < GCD[n, primeProduct] < n, 0,
      If[VectorQ[RandomInteger[{2, n - 1}, m], 
         StrongPseudoprimeTest[#, n] &], 1 - 1/4^m, 0]]

Also on page 118 is MillerRabinPrimeQ[n,m], which applies m strong pseudoprime tests to n, as long as n passes the GCD test with the predefined primeProduct.

Options[MillerRabinPrimeQ] = {RandomBases -> True};

MillerRabinPrimeQ[n_Integer, m_Integer: 1, opts___Rule] :=
   Module[{rQ},
      If[1 < GCD[n, primeProduct] < n, False,
         rQ = RandomBases /. {opts} /. Options[MillerRabinPrimeQ];
         VectorQ[If[rQ, RandomInteger[{2, n - 1}, m], Prime[Range[m]]],
            StrongPseudoprimeTest[#, n] &]]]

The next test was mentioned by @ChipHurst. See Solovay-Strassen test, or page 193 of Bressoud and Wagon, or page 143 of Ribenboim's The New Book of Prime Number Records. The reference is: A Fast Monte-Carlo Test for Primality, SIAM Journal on Computing, v6, p84, 1977. Choose m>1 bases b such that 1<b<n and GCD[n,b]=1. Test each b for JacobiSymbol[b,n]=b^((n-1)/2), mod n. If a base b is found for which the congruence fails, then n is composite and probability 0 is returned. If no b fails, then n is prime with returned probability 1-1/2^k, where k is the number of bases tested.

SolovayStrassenPrimeProbability[n_?OddQ, m_] :=
   Block[{b = RandomInteger[{2, n - 1}, Max[m, 100]]},
      If[1 < GCD[n, primeProduct] < n, 0,
         If[VectorQ[b = Take[Pick[b, GCD[n, b], 1], m], 
           Divisible[JacobiSymbol[#, n] - PowerMod[#, Quotient[n - 1, 2], n], n] &],
           1 - 1/2^Length[b], 0]]]

Timings with random odd integers near 10^200 failed to produce results faster than PrimeQ. However, for tests of primes of order 10^2000, the recommended code is StrongPreudoprimeProbability (essentially the same as MillerRabinPrimeQ), which was faster than PrimeQ when m=1 or m=2. Many more tests are possible, your mileage may vary.