Primality testing, part 1: Compositeness tests, Fermat and Rabin-Miller

How to check that a given integer is prime? Of course we can do the naive method of trial division up to $\lceil \sqrt{n} \rceil$, but this becomes impossible in practice if $n$ is at all large, so we need something better. (Note that trial division is actually used in practice, up to a point, in order to quickly conclude that a given integer is composite if it has some small prime factors.)

A common strategy is to first perform a compositeness test, which can either say with certainty that the given integer is composite, or say that it is prime with high but not 100% probability. Then, once we are, in the words of Cohen, “morally certain” that our integer is prime, we do a real primality test in order to prove it. The reason for this strategy is that in general compositeness tests are easy but primality tests are hard, so it makes sense to first acquire a high degree of certainty that a number is prime before doing a primality test. In this post we will see two compositeness tests, the Fermat test, which is based on Fermat’s little theorem, and the Rabin-Miller test, and experiment with them using Sage (for a change). In a next post we will see primality tests.

The Fermat test

Recall Fermat’s little theorem: if $p$ is a prime and $a$ is any integer in $\{1,\dots,p-1\}$, we have
\[ a^{p-1} \equiv 1 \pmod p. \]
This gives a first compositeness test: to test a given integer $n$, we pick an $a$ at random in $\{1,\dots,n-1\}$, and if $a^{n-1} \not \equiv 1 \pmod n$ we know for certain that $n$ is composite. If, however, we do have $a^{n-1} \equiv 1 \pmod n$, we cannot say that $n$ is prime. For example $15$ is not prime but $4^{14} \equiv 1 \pmod {15}$.

We say that $n$ is a Fermat pseudoprime in the base $a$ if we have $a^{n-1} \equiv 1 \pmod n$. From the above it is clear that a prime $p$ is a Fermat pseudoprime in all bases not divisible by $p$, and so if for an integer $n$ there is any base not divisible by $n$ and in which $n$ is not a Fermat pseudoprime, we can say with certainty that $n$ is composite. If on the other hand $n$ is a Fermat pseudoprime in any such base, we cannot say that $n$ is prime. Of course it would be foolish to test if $n$ is a Fermat pseudoprime in all the bases not divisible by $n$, because that would amount to testing all integers in $\{1,\dots,n-1\}$, and we might as well do trial division in the first place. However, this test turns out to be surprisingly accurate even if we test only a small number of bases. We can implement this test as follows:

def testfermat(n,a):
    return a**(n-1) % n != 1

def fermat(n,k):
    t = range(2,n-1)
    for i in range(k):
        a = t[ZZ.random_element(len(t))]
        if testfermat(n,a):
            return True
        t.remove(a)
    return False

The function fermat(n,k) tests whether $n$ is composite by checking whether it is a Fermat pseudoprime in $k$ bases chosen at random, and returns False if and only if it is a Fermat pseudoprime in all of them. Note that since this is a compositeness test, we return True if $n$ is composite, and False if its compositeness could not be ascertained. Also we only test for bases in $\{2,\dots,n-2\}$, since $n$ is always a Fermat pseudoprime in the base $1$ and almost always in the base $n-1$. If we run this test on all integers between $10$ and $9999$ with $2$ bases each time:

sage: for i in range(10,10000):
....:     if not fermat(i,2) and not is_prime(i):
....:         print(i)
....:         
65
91
703
1729
1921
2821
5611
8911

only $8$ numbers are not detected as composite even though they are actually composite (of course results will vary slightly with each run, since the bases are chosen at random). Using $3$ bases instead, we have on average $5$ misdetected numbers. With $5$ bases we have $2$, and with $10$ bases we almost never have any misdetection.

The Rabin-Miller test

One way to improve the Fermat test in order to eliminate most misdetections is the Rabin-Miller test. This test was first proposed by Miller in 1976, and improved to its current form by Rabin in 1980. Recall that if $p$ is prime, we have by Fermat’s theorem
\[ a^{p-1} \equiv 1 \pmod p, \]
for all integers $a$ relatively prime with $p$. Considering only odd primes $p$, the idea of the algorithm is to factor $p-1$ as $p-1 = 2^eq$, with $q$ odd. Then one of the following must be true: either $a^q \equiv 1 \pmod p$, or there exists $i \in \{0,\dots,e-1\}$ such that $a^{2^iq} \equiv -1 \pmod p$.

Proof: We have
\[\begin{align*}
a^{p-1} &= a^{2^eq} \\
&= a^{2\times 2^{e-1}q} \\
&= (a^{2^{e-1}q})^2
\end{align*} \]
and so by Fermat’s theorem we have
\[ (a^{2^{e-1}q})^2 \equiv 1 \pmod p, \]
whence
\[ a^{2^{e-1}q} \equiv \pm 1 \pmod p. \]
If $a^{2^{e-1}q} \equiv -1 \pmod p$ we are done, otherwise we have by the same argument
\[ a^{2^{e-2}q} \equiv \pm 1 \pmod p, \]
and so on. Assuming that we never get a result of $-1$ we arrive at
\[ a^{2^0q} = a^q \equiv 1 \pmod p, \]
as desired. $\Box$

An integer $n$ is said to be a strong pseudoprime in the base $a$ if one of the two conditions above hold. That is, if letting $n-1 = 2^eq$ with $q$ odd, we have either $a^q \equiv 1 \pmod n$ or $a^{2^iq} \equiv -1 \pmod n$ for some integer $i \in \{0,\dots,e-1\}$. As for the Fermat test, we have seen that a prime $p$ is a strong pseudoprime in any base not divisible by $p$, and so if for any $n$ there exists a base not divisible by $n$ and in which $n$ is not a strong pseudoprime, we know for certain that $n$ is composite. Again, if $n$ is a strong pseudoprime in a given base $a$ this does not imply that $n$ is prime. In fact there is a less than 25% probability that $n$ is composite if $a$ is chosen at random. This means that if we try, say, $25$ values of $a$, the probability of error becomes less than $1/4^{25} = 1/2^{50}$, which is in general lower than the probability of a hardware error.

We can implement this test as follows:

def testrabinmiller(n,a):
    e = 0
    q = n-1
    while q % 2 == 0:
        e += 1
        q = q/2
    if a**q % n == 1:
        return False
    for i in range(e):
        if a**(q*2**i) % n == n-1:
            return False
    return True

def rabinmiller(n,k):
    t = range(2, n-1)
    for i in range(k):
        a = t[ZZ.random_element(len(t))]
        if testrabinmiller(n,a):
            return True
        t.remove(a)
    return False

And we now have:

sage: for i in range(10,10000):                                                                    
....:     if not rabinmiller(i,2) and not is_prime(i):
....:         print(i)
....:         
3367
8321

only two misdetected numbers, instead of eight.

Leave a Reply

Your email address will not be published. Required fields are marked *