Fast Deterministic Prime Test for n Less Than Quintillion

The most common way to test whether a large number is prime is the Miller-Rabin test. If the test says a number is composite, it's definitely composite. Otherwise, the number is very likely, but not certain, to be prime. A pseudoprime is a composite number that slips past the Miller-Rabin test. (Actually, a strong pseudoprime. More on that below.)

Miller-Rabin Test

The Miller-Rabin test is actually a sequence of tests, one for each prime number. First, you run the test associated with 2, then the test associated with 3, then the one associated with 5, etc. If we knew the smallest numbers for which these tests fail, then for smaller numbers, we know for certain that they're prime if they pass. In other words, we can turn the Miller-Rabin test for probable primes into a test for provable primes.

Lower Bound on Failure

A recent result by Yupeng Jiang and Yingpu Deng finds the smallest number for which the Miller-Rabin test fails for the first nine primes. This number is:

N = 3,825,123,056,546,413,051

or more than 3.8 quintillion. So, if a number passes the first nine Miller-Rabin tests, and it's less than N, then it's prime. Not just a probable prime, but definitely prime. For a number n < N, this will be more efficient than running previously-known deterministic primality tests on n.

Python Implementation

Let's play with this in Python. The SymPy library implements the Miller-Rabin test in a function mr.
The following shows that N is composite and that it is a false positive for the first nine Miller-Rabin tests.

    from sympy.ntheory.primetest import mr

    N = 3825123056546413051
    assert(N == 149491*747451*34233211)
    ps = [2, 3, 5, 7, 11, 13, 17, 19, 23]
    print( mr(N, ps) )

This doesn't prove that N is the smallest number with these properties; we need the proof of Jiang and Deng for that. But assuming their result is right, here's an efficient deterministic primality test that works for all n less than N.

    def is_prime(n):
        N = 3825123056546413051
        assert(n < N)
        ps = [2, 3, 5, 7, 11, 13, 17, 19, 23]
        return mr(n, ps)

Jiang and Deng assert that N is also the smallest composite number to slip by the first 10 and 11 Miller-Rabin tests. We can show that N is indeed a strong pseudoprime for the 10th and 11th primes, but not for the 12th prime.

    print( mr(N, [29, 31]) )
    print( mr(N, [37]) )

This code prints True for the first test and False for the second. That is, N is a strong pseudoprime for bases 29 and 31, but not for 37.

Pseudoprimes and Strong Pseudoprimes

Fermat's little theorem says that if n is prime, then:

a^n-1 = 1 mod n

for all 0 < a < n. This gives a necessary but not sufficient test for primality. A (Fermat) pseudoprime for base a is a composite number n such that the above holds, an example of where the test is not sufficient.

The Miller-Rabin test refines Fermat's test by looking at additional necessary conditions for a number being prime. Often a composite number will fail one of these conditions, but not always. The composite numbers that slip by are called strong pseudoprimes or sometimes Miller-Rabin pseudoprimes.

Miller and Rabin's extra testing starts by factoring n-1 into 2 d where d is odd. If n is prime, then for all 0 < a < n either:

a^d = 1 mod n

or:

a^2kd = -1 mod n

for all k satisfying 0 ≤ k < s. If one of these two conditions holds for a particular a, then n passes the Miller-Rabin test for the base a.

It wouldn't be hard to write your own implementation of the Miller-Rabin test. You'd need a way to work with large integers and to compute modular exponents, both of which are included in Python without having to use SymPy.

Let us know your thoughts and questions in the comments.