*(Continues this post.)*

So we have this integer $n$, and we would like to know whether it is prime or not. We have used the Rabin-Miller compositeness test on it a couple times and it has always turned negative, so the probability that $n$ is not prime is excruciatingly small. Still, we have not *proved* that $n$ is prime; in order to do that, we need a real *primality* test. Basically, we are looking for a converse to Fermat’s theorem. Such a converse was found by Pocklington and improved by Lehmer. Essentially, it is based on the fact that an integer $a$ of multiplicative order $n-1$ modulo $n$ (*i.e.*, such that $a^{n-1} \equiv 1 \pmod n$ and $a^e \not \equiv 1 \pmod n$ for all positive $e < n-1$) exists if and only if $n$ is prime. More precisely:

Let $p$ be a prime divisor of $n-1$, and $a_p$ be such that $a_p^{n-1} \equiv 1 \pmod n$ and $\gcd(a_p^{(n-1)/p}-1,n) = 1$. Then for any divisor $d$ of $n$ we have $d \equiv 1 \pmod {p^{\alpha_p}}$, where $p^{\alpha_p}$ is the largest power of $p$ which divides $n-1$.

*Proof:* It is clearly sufficient to prove the result only when $d$ is a prime divisor of $n$ (since any divisor of $n$ is a product of prime divisors). Let then $d$ be a prime divisor of $n$. Since $a_p^{n-1} \equiv 1 \pmod n$, $a_p$ is in the group of units modulo $n$ and so is coprime to $n$. Thus it is also coprime to $d$, and we have $a_p^{d-1} \equiv 1 \pmod d$. On the other hand, since $\gcd(a_p^{(n-1)/p}-1,n) = 1$ we see that $a_p^{(n-1)/p} \not \equiv 1 \pmod d$ (otherwise $d$ divides $a_p^{(n-1)/p}-1$ and $n$, contradicting the $\gcd$). Letting $e$ be the multiplicative order of $a_p$ modulo $d$ we obtain that $e | d-1$ (since $a_p^{d-1} \equiv 1 \pmod d$), $e \not | (n-1)/p$ (since $a_p^{(n-1)/p} \not \equiv 1 \pmod d$), and $e | n-1$ (since $a_p^{n-1} \equiv 1 \pmod n$ and $d$ divides $n$ we also have $a_p^{n-1} \equiv 1 \pmod d$).

So $e|n-1$ and $e \not | (n-1)/p$. This means that if we remove only one factor $p$ from $n-1$, $e$ no longer divides it. In turn this means that all factors $p$ of $n-1$ are in $e$ also, and so that $p^{\alpha_p} | e$. Since also $e | d-1$ this means that $p^{\alpha_p} | d-1$, and so $d \equiv 1 \pmod{p^{\alpha_p}}$, as desired. $\Box$

This gives us a first primality test: if for all prime divisors $p$ of $n-1$ we can find an integer $a_p$ satisfying the conditions above, then $n$ is prime. Indeed, for any divisor $d$ of $n$ we then have $d \equiv 1 \pmod {p^{\alpha_p}}$ for all prime divisors $p$ of $n-1$, and so $d \equiv 1 \pmod{n-1}$, which means that either $d = 1$ or $d \ge n$. We remark also that if the same integer serves as $a_p$ for all $p$, then its multiplicative order modulo $n$ is $n-1$ (and as we remarked above, such an element exists if and only if $n$ is prime).

This test has two drawbacks: the first and more serious is that it requires that we know all the prime factors of $n-1$, the second is that we need a way to find a suitable $a_p$ for all $p$.

The first problem can be mitigated as follows: instead of testing for all prime divisors of $n-1$ we only factor it up to $\sqrt{n}$. More precisely, we write $n-1 = F\cdot U$, where $F > \sqrt{n}$ and $\gcd(F,U) = 1$. Then we only need the prime factors of $F$: the same argument as above yields that if $d$ is a divisor of $n$ we have $d \equiv 1 \pmod F$, and since $F > \sqrt{n}$ we must have $d = 1$ or $d > \sqrt{n}$, so $n$ must be prime.

Finally, we need to find an integer $a_p$ satisfying the conditions above. Remember that we only do this when we are already certain that $n$ is prime and want to obtain a rigorous proof. So we can assume that $n$ is prime, and the group of units modulo $n$ is cyclic of order $n-1$. If we take $a_p$ to be a generator of this group, we obtain of course $a_p^{n-1} \equiv 1 \pmod n$. On the other hand, we have $a_p^{(n-1)/p} \not \equiv 1 \pmod n$, and so $a_p^{(n-1)/p}-1$ is not a multiple of $n$, and so is coprime with $n$ since $n$ is prime.

There are $\varphi(n-1)$ generators in the group of units modulo $n$, so if we pick an element at random it will work with probability at least $\varphi(n-1)/(n-1)$. Since, in the words of Hardy and Wright, $\varphi(n)$ is “always ‘nearly $n$’,” this doesn’t seem too bad. We have in fact

\[ \varphi(n-1) > C\frac{n}{\ln \ln n} \]

for some constant $C$, and so we have a good chance of finding a suitable $a_p$ after about $\frac{10}{C}\ln \ln n$ tries.

This algorithm can be implemented as follows (we do not use the factorisation $n-1 = F\cdot U$ in order to keep the code simple):

def testpocklington(n,p): t = range(2,n-1); while len(t) > 0: a = t[ZZ.random_element(len(t))] if a**(n-1) % n == 1 and gcd(a**((n-1)//p)-1, n) == 1: return True t.remove(a) return False def pocklington(n): for p in prime_divisors(n-1): if not testpocklington(n,p): return False return True

And so we can use this in conjunction with the Rabin-Miller compositeness test for a full primality test:

def my_is_prime(n): if rabinmiller(n,2): return False return pocklington(n)

which works as intended:

sage: for i in range(10,10000): ....: if is_prime(i) and not my_is_prime(i) or not is_prime(i) and my_is_prime(i): ....: print(i) ....: sage: