A better way to make the Fermat test more accurate is to realize that if an odd number n is prime, then the number 1 has just two square roots modulo n: 1 and -1. So the square root of an-1, a(n-1)/2 (since n will be odd), is either 1 or -1. (We actually could calculate which it should be using the Jacobi symbol, see the glossary page on Euler PRP's, but we wish to develop a stronger test here.) If (n-1)/2 is even, we can easily take another square root... Let's make this into an algorithm:
Write n-1 = 2sd where d is odd and s is non-negative: n is a strong probable-prime base a (an a-SPRP) if either ad ≡ 1 (mod n) or (ad)2r ≡ -1 (mod n) for some non-negative r less than s.
Again all integers n > 1 which fail this test are composite; integers that pass it might be prime. The smallest odd composite SPRP's are the following.
A test based on these results is quite fast, especially when combined with trial division by the first few primes. If you have trouble programming these results Riesel [Riesel94, p100] has PASCAL code for a SPRP test, Bressoud has pseudocode [Bressoud89, p77], and Langlois offers C-Code. See the glossary page "Strong PRP" for more information.
It has been proven ([Monier80] and [Rabin80]) that the strong probable primality test is wrong no more than 1/4th of the time (3 out of 4 numbers which pass it will be prime). Jon Grantham's "Frobenius pseudoprimes" can be used to create a test (see [Grantham98]) that takes three times as long as the SPRP test, but is far more than three times as strong (the error rate is less than 1/7710).
Individually these tests are still weak (and again there are infinitely many a-SPRP's for every base a>1 [PSW80]), but we can combine these individual tests to make powerful tests for small integers n>1 (these tests prove primality):
To make a quick primality test from these results, start by dividing by the first few primes (say those below 257); then perform strong primality tests base 2, 3, ... until one of the criteria above is met. For example, if n < 25,326,001 we need only check bases 2, 3 and 5. This is much faster than trial division because someone else has already done much of the work, but will only work for small numbers (n < 1016 with the data above).
Note that these results can be strengthened by not treating them as separate tests, but rather realizing we are finding square root of -1. For example, n = 46,856,248,255,981 is a 2 and 7 pseudoprime, but
2(n-1)/4 ≡ 34456063004337 (mod n), and
7(n-1)/4 ≡ 21307242304265 (mod n).
The square of both of these is -1. If n were prime, then it would have only two square root and the above would be equal or negatives of each other; yet gcd(n,34456063004337-21307242304265) = 4840261 and we have factored n.
Finally, there is a fair amount more that could (and should) be said. We could discuss Euler pseudoprimes and their relationship with SPRP's. Or we could switch
to the "plus side" and discuss Lucas pseudoprimes, or Fibonacci pseudoprimes, or the important combined tests... but that would take a chapter of a book--and it has already been well written by Ribenboim [Ribenboim95]. Let us end this section with one last result:
Miller's Test [Miller76]: If the extended Riemann hypothesis is true, then if n is an a-SPRP for all integers a with 1 < a < 2(log n)2, then n is prime.
The extended Riemann hypothesis is far too complicated for us to explain here--but should it be proven, then we would have a very simple primality test. Until it is proven, we can at least expect that if n is composite, we should be able to find an a that shows it is composite (a witness) without searching "too long." Most surveys cover Miller's test (often with the constant 70 from [Osterle1979] as Miller's article just said O((log n)2)); the improvable constant 2 is due to Bach [Bach85], see also [CP2001, pp. 129-130]. Note that heuristically Bach and Huelsbergen [BH1993] argue that we should be able to replace the bound in Miller's test with a bound near:
(log 2)-1 log n log log n.
Note that there is no finite set of bases that will work in Miller's test. In fact, if for n composite we let W(n) denote the least witness for n (the least a which shows n is composite), then there are infinitely many composite n with
W(n) > (log n)1/(3 log log log n) [AGP94]