# Is a probabilistic primality test implemented in Mathematica?

Imagine that I need to verify that a specific integer, like $n = 10^{1347}+3049$ is prime.

ProvablePrimeQ takes a lot of time, and, even if it will eventually work for this integer, it will work forever for somewhat larger integers.

PrimeQ works in 0.156 seconds, but, according to specification, it uses "the multiple Rabin-Miller test in bases 2 and 3 combined with a Lucas pseudoprime test". As far as I understand, its answer for every specific number is deterministic - either correct or false (no probabilities involved). It may happen that it returns false answer for $n=10^{1347}+3049$, and, in this case, it will still return false answer even if I call PrimeQ[n] a million times.

I do not understand why use Rabin-Miller test in bases 2 and 3 only? The theory says that if we select bases AT RANDOM, we get the correct answer with probability $3/4$. We can then repeat the test $N$ times to reduce the probability of error down to $(1/4)^N$. This way I can control the error probability.

The question is: Is there a randomised primality test implemented in Mathematica, where I can control the error probability as a parameter? Or where the error probability is given in the specification and I can reduce it by repeating the test?

• Where in the documentation did you find the claim, that PrimeQ might return incorrect answers? Commented Jun 13, 2018 at 15:38
• There is no randomized primality test built in as a system function. The reason to use bases 2 and 3, then go to the Lucas test, is that there is no compelling reason to believe alternatives offer any improvement. (1) A randomized set of Miller-Rabin tests would mean any incorrect outcome would be effectively unrepeatable. (2) A bunch of runs with probability 1-(1/4)^n incorrect is not really a great substitute for a faster tandem that has no known failure cases. Commented Jun 13, 2018 at 15:50
• "As far as I understand, its answer for every specific number is deterministic - either correct or false (no probabilities involved)." That is not enough. A primality test is deterministic if it outputs True when the number is a prime and False when the input is composite with probability 1. Deterministic algorithm is ECPP and APR-CL, e.g. Commented Apr 15, 2023 at 9:50

PrimeQ is using Miller-Rabin tests with bases 2 and 3 to ensure that n is not a square, which is a necessary test, as otherwise the choice for Q in the Lucas test will not succeed. Other similar implementations use only the Miller-Rabin test with base 2 and check if n is square afterwards. See the wiki page on Baillie-PSW primality test.

• Would you please provide a source for these claims? Commented Sep 9, 2018 at 14:06
• The wiki page on Baillie-PSW does indeed not talk about implementation. Also the Mathematica documentation does not clarify this. I have found the remark on the Mathematica implementation in Python source code I found some time ago on the web. I do not know the source at the moment (but I will dig further). I would like to know why testing with Miller-Rabin bases 2 and 3 ensures, that n is squarefree. Commented Sep 9, 2018 at 14:31
• I just checked if any square would pass Fermat tests base 2 and 3, up to about 5*10^18. So, testing with Fermat bases 2 and 3 (or even better Miller-Rabin tests bases 2 and 3) prevents squares to pass. Commented Sep 10, 2018 at 6:08
• @Henrik Schumacher: I posted a question on link and got the answer that it is very unlikely that, if a number passes the Miller-Rabin tests with bases 2 and 3, the number is a square, but that you cannot rule it out. Commented Sep 11, 2018 at 13:20
• Mathematica now uses 2020's Baillie–PSW, see arxiv.org/abs/2006.14425. No base 3 test. Commented Apr 15, 2023 at 9:47

Not actually and answer to the question, but may be of interest. Code below is a slight modification of code in this book.

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

StrongPseudoprimeQ[_, 1] = False;
StrongPseudoprimeQ[b_, n_] := With[{bSeq = SpspSequence[b,n]},
Union[bSeq] == {1} || MemberQ[bSeq, -1]
]


StrongPseudoprimeQ above applies the Miller-Rabin test using base n. You can use that to apply Miller-Rabin to any bases you want. For example:

n = 1373653;
StrongPseudoprimeQ[2, n]
StrongPseudoprimeQ[3, n]
StrongPseudoprimeQ[5, n]
(*
True
True
False
)*


Related content is here, here and here.