Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
What you are asking is a matter of on-going research since the time of the first digital computers. It is an inherently difficult problem. Even when restricted to context of Mathematica, the question is too broad. I'm voting to close it. – m_goldberg Apr 23 '14 at 6:48
I am not sure if I agree that the question is too broad, given that it takes PrimeQ as its performance reference and is asking whether it is possible to perform a test that is both cheaper and weaker. Arguably the answer to the question is "no" (unless you code it yourself), but the question itself is not ill-posed. I certainly don't agree with the downvote. – Oleksandr R. Apr 23 '14 at 7:59
In other words, I want to know if Mathematica has a function which I imagine would be ProbablyPrimeQ[n, certainty] which never errs when saying No, but may err when saying Yes with probability scaling down with certainty. I have not been able to find such a function, so I ask here. – JeremyKun Apr 23 '14 at 14:00
(1) Mathematica does not have any parameters to fiddle with for PrimeQ. Best I can recommend is a quick sieve against small primes, then use a few M-R tests and skip Lucas, as I think that's the slower one. (2) I think this one should be reopened in case anyone wants to provide code, along the above lines or otherwise. – Daniel Lichtblau Apr 23 '14 at 16:32
Here's something that has an adjustable parameter: ProbablyPrimeQ[p_, m_:100] := OddQ[p] && VectorQ[RandomInteger[{1, p}, m], Divisible[JacobiSymbol[#, p] - PowerMod[#, (p-1)/2, p], p]&] – Chip Hurst Apr 23 '14 at 16:52

1 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] :=
      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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.