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.

3What 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 – 2014-04-23T06:48:32.3732I 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. – 2014-04-23T07:59:59.847Yes. This is a question specifically about tuning the efficiency of PrimeQ, or whether there are alternatives in Mathematica that allow for this. The context is just for illustration (to avoid people answering things like "try lots of small divisors!"). – JeremyKun – 2014-04-23T13:47:22.847

2In 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 – 2014-04-23T14:00:38.270

I recognize that the title was potentially misleading, and I've updated it to reflect the content better. – JeremyKun – 2014-04-23T15:49:37.130

5(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 – 2014-04-23T16:32:46.4276Here'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 – 2014-04-23T16:52:36.853http://mathworld.wolfram.com/Pseudoprime.html – Enrique Pérez Herrero – 2015-07-05T15:16:10.830