Is a probabilistic primality test implemented in Mathematica?

Question

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? — Henrik Schumacher, Jun 13 '18 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. — Daniel Lichtblau, Jun 13 '18 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. — Валерий Заподовников, Apr 15 '23 at 09:50

score 1 · Answer 1 · edited Sep 09 '18 at 14:53

1

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.

edited Sep 09 '18 at 14:53

Johu

4,918
16
43

answered Sep 09 '18 at 14:00

PaulH

11
2

2

Would you please provide a source for these claims? – Henrik Schumacher Sep 09 '18 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. – PaulH Sep 09 '18 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. – PaulH Sep 10 '18 at 06: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. – PaulH Sep 11 '18 at 13:20
Mathematica now uses 2020's Baillie–PSW, see https://arxiv.org/abs/2006.14425. No base 3 test. – Валерий Заподовников Apr 15 '23 at 09:47

score 0 · Answer 2 · answered Dec 12 '21 at 21:20

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
)*

Is a probabilistic primality test implemented in Mathematica?

2 Answers2