# Is there a PrimeQ whose accuracy guarantee you can adjust?

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?

-
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 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 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 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 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 at 16:52