NTIC Primes – Probably

Section 12.2 Primes – Probably

Primality testing is full of little tidbits like those in the previous section, and tantalizingly devoid of easy methods that work for all special cases. Indeed, none of these paths lead us to reliable, reasonably fast discovery of large primes for cryptographic purposes, nor do other computationally infeasible methods like using Wilson's Theorem or other even stranger formulas (some of which appear later in this text).

🔗

Instead, what is typically done is to pick a number, and then use tests on it that do not guarantee primality!

🔗

Why would this work? The idea is that if a given number passes enough tests that do not guarantee primality but have a quite low false positive rate in practice, then the probability the number you have is composite is lower than the (very low) chance that your computer made an arithmetic error due to cosmic rays (though one still has to be careful of bugs like the one described in the discussion before Algorithm 12.1.7).

🔗

This is astonishing, but true. Then if you end up with a number that likely to be prime, you can always confirm its primality with one of the various slower tests I will not describe.

🔗

Subsection 12.2.1 Pseudoprimes

🔗

We start this discussion with our visual representation of powers (see Subsection 8.2.1).

🔗

Notice again here that Fermat's Little Theorem is visible in the second-to-last column. The graphic has been expanded, so that the last column is a slight restatement thereof, true for all

$a\text{:}$

$\begin{equation*} a^p\equiv a\text{ (mod }p)\text{.} \end{equation*}$

🔗

(See Exercise 12.7.3 and Exercise 9.6.3.) Go ahead and confirm it in the interactive version.

xxxxxxxxxx
 
import matplotlib.pyplot as plt
from matplotlib.ticker import IndexLocator, FuncFormatter
@interact
def power_table_plot(p=(11,prime_range(100))):
    mycmap = plt.get_cmap('gist_earth',p-1)
    myloc = IndexLocator(floor(p/5),.5)
    myform = FuncFormatter(lambda x,y: int(x+1))
    cbaropts = {  'ticks':myloc, 'drawedges':True, 'boundaries':srange(.5,p+.5,1)}
    P=matrix_plot(matrix(p-1,[mod(a,p)^b for a in range(1,p) for b in srange(p+1)]),cmap=mycmap, colorbar=True, colorbar_options=cbaropts, ticks=[myloc,myloc], tick_formatter=[None,myform])    
    show(P,figsize=6)

🔗

This is a useful criterion, as it works for all input, including multiples of the modulus. We can now use it to state a test for possible primality:

🔗

Fact 12.2.2.

If there is an $a$ such that $a^n\not\equiv a$ (mod $n$ ), then $n$ must be composite.

🔗

So if

$a^n\equiv a$ (mod

$n$ ) for a given

$n\text{,}$ it's at least possible that

$n$ is prime.

🔗

Definition 12.2.3.

If $a^n\equiv a$ (mod $n$ ), we say $n$ passes the base $a$ test.

🔗

It turns out that everyone from the ancient Chinese to Leibniz used this test for the base

$a=2$ to assert numbers are prime. And it doesn't do a bad job. As some former students pointed out, it's sort of like internet date matching for primes; it doesn't always work but can succeed reasonably often.

xxxxxxxxxx
 
@interact
def _(n=100):
    pretty_print(html("Here are the numbers through $%s$ that pass the base 2 test"%n))
    pretty_print(html("along with whether they are actually prime"))
    for i in [2..n]:
        if mod(2^i,i)==2:
            pretty_print(html(r"$2^{%s}\equiv 2\text{ (mod }%s)$ and the primality of $%s$ is %s"%(i,i,i,is_prime(i))))

🔗

We can change the numbers in the range of the preceding interact to check for more – say up to 1000, which allows exploring the following question.

🔗

Question 12.2.4.

Are there any numbers which satisfy the base $a$ test and are not prime?

🔗

To the surprise of many in the world of numbers, the answer is yes. The numbers

$n=341\text{,}$

$n=561\text{,}$ and

$n=645$ turn out to fall in that category.

xxxxxxxxxx
 
print("We factor 341 and get %s"%factor(341))
print("We factor 561 and get %s"%factor(561))
print("We factor 645 and get %s"%factor(645))

🔗

That's still not bad – out of one hundred seventy-one total such potential primes base

$2\text{,}$ only three of them actually are not prime, or about one and three quarters percent. That is unusual enough that we have a special name for composite numbers that pass one of the base

$a$ tests.

🔗

Definition 12.2.5. Pseudoprimes.

If $a^n\equiv a$ (mod $n$ ) but $n$ is not prime, we say it is a pseudoprime base $a\text{.}$

🔗

That is to say, if a number satisfies Fermat's Little Theorem, we think it is likely enough to be prime to call it a pseudoprime if it isn't. (Prime, that is.)

🔗

Perhaps unfortunately to cryptographers (though interestingly to pure mathematicians!), it turns out that there are infinitely many such pseudoprimes.

🔗

Fact 12.2.6.

If $n$ is a pseudoprime (base 2), then so is $2^n-1\text{.}$

🔗

We will get this result as a corollary of something stronger soon (see Corollary 12.4.3 and Theorem 12.4.2).

🔗

All the Fermat and Mersenne numbers pass the base

$2$ test, incidentally, though they are all quite large compared to a typical number you might try.

🔗

Subsection 12.2.2 Prime impostors, and how to avoid them

🔗

If we want to check things out more carefully, we can try to test for primality with a different base. In the next cell, we choose

$a=3\text{.}$

xxxxxxxxxx
 
for n in [341,561,645]:
    pretty_print(html(r"$3^{%s}\equiv %s\text{ (mod }%s)"%(n,mod(3,n)^n,n)))

🔗

As you can see, this exposes

$341$ and

$645$ as fakes. What about

$561\text{?}$ Let's try that one with base

$a=5$ as well.

xxxxxxxxxx
 
@interact
def _(p=(5,prime_range(50))):
    for pr in prime_range(next_prime(p)):
        pretty_print(html(r"$%s^{561}\equiv %s\text{ (mod }561)"%(pr,mod(pr,561)^561)))

🔗

Hmm, that's interesting. What if I change to a different prime base, like

$a=7$ or

$11\text{?}$ Try it above.

🔗

In the next cell, I get systematic. We should expect output if

$561$ doesn't pass the base

$a$ test for some

$a\text{.}$

xxxxxxxxxx
 
@interact
def _(p=(5,prime_range(1000))):
    pretty_print(html("The primes up to $%s$ for which $561$ fails the base $p$ test:"%p))
    for pr in prime_range(next_prime(p)):
        if mod(pr,561)^561!=pr:
            pretty_print(html(r"$%s^{561}\equiv %s\\text{ (mod }561)$"%(pr,mod(pr,561)^561)))

🔗

It appears that

$p^{561}\equiv p$ mod 561 for every prime

$p\text{!}$ Let's prove it.

🔗

Fact 12.2.7.

The number $561$ is a pseudoprime for every integer base $a\text{.}$

🔗

Proof.

We know that

\begin{equation*} 561=3\cdot 11\cdot 17\text{,} \end{equation*}

so by the Chinese Remainder Theorem

\begin{equation*} a^{561}\equiv a\text{ (mod }561) \end{equation*}

if and only if $a^{561}\equiv a$ holds for the prime power factors $3,11,17\text{;}$ so we will check them.

Remember, the exponents for these congruences live in the (mod $\phi(p)$) world, so we just need to check what $561$ is in each of those worlds. We get:

$561\equiv 1\text{ (mod }16=17-1)$ so $a^{561}\equiv a^1\text{ (mod }17)$
$561\equiv 1\text{ (mod }2=3-1)$ so $a^{561}\equiv a^1\text{ (mod }3)$
$561\equiv 1\text{ (mod }10=11-1)$ so $a^{561}\equiv a^1\text{ (mod }11)$

That is, for $p=3,11,17$ we see

\begin{equation*} a^{561}\equiv a^{1}\text{ (mod }p) \end{equation*}

Using Proposition 5.4.4, this congruence is always true!

By the way, we note that $a^{560}$ is not congruent to $1\text{,}$ which explains why we use $a^n\equiv a$ for these definitions.

🔗

Definition 12.2.8.

We call a number which is pseudoprime to every base $a\text{,}$ but is not a prime number a Carmichael number, in honor of the first person to actually produce such numbers, Robert Carmichael (in 1912).

🔗

So is

$561$ a Carmichael number? We saw the factorization above, but here it is again:

xxxxxxxxxx
 
factor(561)

🔗

The proof of Fact 12.2.7 suggests that to find a Carmichael number

$n\text{,}$ we might want to look at

$n$ which are a product of primes

$p_i$ such that

$n-1\equiv 1$ in the exponent world of

$p_i\text{.}$ It turns out that this is true, and one can prove something even more specific.

🔗

Proposition 12.2.9. Korselt's Theorem.

Carmichael numbers are precisely those composite $n$ for which $n$ is a product of at least two distinct primes $p_i$ (no squares)

$\begin{equation*} n=p_1p_2p_3\cdots p_k\text{ with }p_i\neq p_j \end{equation*}$

such that

$\begin{equation*} p_i-1\mid n-1 \end{equation*}$

for all the prime factors.

🔗

Proof.

Prime numbers satisfy almost all the conditions trivially. To show that $561$ is a Carmichael number we used this idea in the form $n\equiv 1$ (mod $\phi(p_i)$) for all three prime factors, and essentially the same argument applied to any number satisfying the hypotheses is a Carmichael number.

We will not prove the other half of this theorem (that all Carmichael numbers have this form). It is not hard, however, using a slight variant on the Euler $\phi$ function one can acquire from investigating $U_n$ for composite $n\text{.}$

🔗

Example 12.2.10.

Evaluate this Sage cell to see the previous result applied to identify another Carmichael number.

xxxxxxxxxx
 
n=29341
pretty_print(html("$%s$ is composite with factorization $%s$, but"%(n,factor(n))))
for fact,pow in factor(n):
    pretty_print(html(r"$%s^{%s}\equiv %s\text{ (mod }%s)$"%(fact,n,mod(fact,n)^n,n)))
pretty_print(html("and"))
for fact,pow in factor(n):
    pretty_print(html(r"$%s\equiv %s\text{ (mod }\phi(%s)=%s)$"%(n, mod(n,euler_phi(fact)), fact, euler_phi(fact))))