MAT 338 Day 20 2011

2449 days ago by kcrisman

Today and next time we will finish up the cryptographic side of things - via the more theoretical issues of finding primes and factoring big numbers.

As we have seen, it is not terribly hard to find lots of small primes easily - using the Sieve of Eratosthenes, or by making numbers coprime to known primes and then factoring them.  The problem is that almost every effort to find lots of big ones has been stymied.  Primes simply do not follow nice enough rules to find them.  This is despite the fact that they seem to follow very nice rules on average, which we will explore next quad.  

One example of an attempt to do this we mentioned last time.  Fermat thought that $2^{2^n}+1$ would always be prime, but as noted, in 1732 Euler proved that this is not true for $n=5$:

for n in [0..7]: html("If $n=%s$, then $2^{2^n}+1=2^{2^%s}+1=%s$ factors as $%s$"%(n,n,2^(2^n)+1,factor(2^(2^n)+1))) 
       
If n=0, then 2^{2^n}+1=2^{2^0}+1=3 factors as 3

If n=1, then 2^{2^n}+1=2^{2^1}+1=5 factors as 5

If n=2, then 2^{2^n}+1=2^{2^2}+1=17 factors as 17

If n=3, then 2^{2^n}+1=2^{2^3}+1=257 factors as 257

If n=4, then 2^{2^n}+1=2^{2^4}+1=65537 factors as 65537

If n=5, then 2^{2^n}+1=2^{2^5}+1=4294967297 factors as 641 * 6700417

If n=6, then 2^{2^n}+1=2^{2^6}+1=18446744073709551617 factors as 274177 * 67280421310721

If n=7, then 2^{2^n}+1=2^{2^7}+1=340282366920938463463374607431768211457 factors as 59649589127497217 * 5704689200685129054721
If n=0, then 2^{2^n}+1=2^{2^0}+1=3 factors as 3

If n=1, then 2^{2^n}+1=2^{2^1}+1=5 factors as 5

If n=2, then 2^{2^n}+1=2^{2^2}+1=17 factors as 17

If n=3, then 2^{2^n}+1=2^{2^3}+1=257 factors as 257

If n=4, then 2^{2^n}+1=2^{2^4}+1=65537 factors as 65537

If n=5, then 2^{2^n}+1=2^{2^5}+1=4294967297 factors as 641 * 6700417

If n=6, then 2^{2^n}+1=2^{2^6}+1=18446744073709551617 factors as 274177 * 67280421310721

If n=7, then 2^{2^n}+1=2^{2^7}+1=340282366920938463463374607431768211457 factors as 59649589127497217 * 5704689200685129054721

Nobody knows if there are any more primes in this sequence; the information above suggest that the prime factors are quite large, though.  There is a special test called Pépin's test that tests Fermat numbers for primality.  It is equivalent to checking whether $3$ is a primitive root of $2^{2^n}+1$, but proving it is just a little beyond us right now.

However, we can at least prove what seems obvious above - namely, that lots of primes come out of Fermat numbers.  First, we need a lemma, whose proof is just multiplication.

Lemma:  $$\text{If }\ell=jk,\text{ then }2^\ell-1=2^{jk}-1=\left(2^j+1\right)\left(\left(2^j\right)^{k-1}-\left(2^j\right)^{k-2}-\left(2^j\right)^{k-3}-\cdots-\left(2^j\right)-1\right)$$

Proposition: $$F_n=2^{2^n}+1\text{ and }F_m=2^{2^m}+1\text{ are coprime if }m\neq n$$

  • First, notice that the numbers are very closely related; if $n<m$, then $F_n-1$ divides $F_m-1$.  In fact, $$2^{2^m}=\left(2^{2^n}\right)^{2^{m-n}}\; .$$
  • Because of this, using the lemma with $j=2^n$ and $k=2^{m-n}$, we get$$2^{2^m}-1=\left(2^{2^n}+1\right)\left(\left(2^{2^n}\right)^{2^{m-n}-1}-\left(2^{2^n}\right)^{2^{m-n}-2}-\cdots -\left(2^{2^n}\right)^1-1 \right)$$ so $$F_n=2^{2^n}+1\mid 2^{2^m}-1=F_m-2$$
  • So any number $d$ that divides $F_n$ also divides $F_m-2$, hence a common divisor of $F_n$ and $F_m$ could only be two or one.
  • But both numbers are odd, so the gcd is 1.

Another early attempt at this was an idea of Marin Mersenne.  Mersenne was a Minim monk who not only acted as a clearinghouse for scientific knowledge in early 17th century France (particularly between Pascal, Fermat, Descartes, and their friends) but also wrote major theological and music theoretical treatises of his own.  He suggested that one try searching for primes of the form $2^p-1$, where $p$ is itself prime.  

for p in prime_range(100): html("If $p=%s$, then $2^p-1=2^{%s}-1=%s$ factors as $%s$"%(p,p,2^p-1,factor(2^p-1))) 
       
If p=2, then 2^p-1=2^{2}-1=3 factors as 3

If p=3, then 2^p-1=2^{3}-1=7 factors as 7

If p=5, then 2^p-1=2^{5}-1=31 factors as 31

If p=7, then 2^p-1=2^{7}-1=127 factors as 127

If p=11, then 2^p-1=2^{11}-1=2047 factors as 23 * 89

If p=13, then 2^p-1=2^{13}-1=8191 factors as 8191

If p=17, then 2^p-1=2^{17}-1=131071 factors as 131071

If p=19, then 2^p-1=2^{19}-1=524287 factors as 524287

If p=23, then 2^p-1=2^{23}-1=8388607 factors as 47 * 178481

If p=29, then 2^p-1=2^{29}-1=536870911 factors as 233 * 1103 * 2089

If p=31, then 2^p-1=2^{31}-1=2147483647 factors as 2147483647

If p=37, then 2^p-1=2^{37}-1=137438953471 factors as 223 * 616318177

If p=41, then 2^p-1=2^{41}-1=2199023255551 factors as 13367 * 164511353

If p=43, then 2^p-1=2^{43}-1=8796093022207 factors as 431 * 9719 * 2099863

If p=47, then 2^p-1=2^{47}-1=140737488355327 factors as 2351 * 4513 * 13264529

If p=53, then 2^p-1=2^{53}-1=9007199254740991 factors as 6361 * 69431 * 20394401

If p=59, then 2^p-1=2^{59}-1=576460752303423487 factors as 179951 * 3203431780337

If p=61, then 2^p-1=2^{61}-1=2305843009213693951 factors as 2305843009213693951

If p=67, then 2^p-1=2^{67}-1=147573952589676412927 factors as 193707721 * 761838257287

If p=71, then 2^p-1=2^{71}-1=2361183241434822606847 factors as 228479 * 48544121 * 212885833

If p=73, then 2^p-1=2^{73}-1=9444732965739290427391 factors as 439 * 2298041 * 9361973132609

If p=79, then 2^p-1=2^{79}-1=604462909807314587353087 factors as 2687 * 202029703 * 1113491139767

If p=83, then 2^p-1=2^{83}-1=9671406556917033397649407 factors as 167 * 57912614113275649087721

If p=89, then 2^p-1=2^{89}-1=618970019642690137449562111 factors as 618970019642690137449562111

If p=97, then 2^p-1=2^{97}-1=158456325028528675187087900671 factors as 11447 * 13842607235828485645766393
If p=2, then 2^p-1=2^{2}-1=3 factors as 3

If p=3, then 2^p-1=2^{3}-1=7 factors as 7

If p=5, then 2^p-1=2^{5}-1=31 factors as 31

If p=7, then 2^p-1=2^{7}-1=127 factors as 127

If p=11, then 2^p-1=2^{11}-1=2047 factors as 23 * 89

If p=13, then 2^p-1=2^{13}-1=8191 factors as 8191

If p=17, then 2^p-1=2^{17}-1=131071 factors as 131071

If p=19, then 2^p-1=2^{19}-1=524287 factors as 524287

If p=23, then 2^p-1=2^{23}-1=8388607 factors as 47 * 178481

If p=29, then 2^p-1=2^{29}-1=536870911 factors as 233 * 1103 * 2089

If p=31, then 2^p-1=2^{31}-1=2147483647 factors as 2147483647

If p=37, then 2^p-1=2^{37}-1=137438953471 factors as 223 * 616318177

If p=41, then 2^p-1=2^{41}-1=2199023255551 factors as 13367 * 164511353

If p=43, then 2^p-1=2^{43}-1=8796093022207 factors as 431 * 9719 * 2099863

If p=47, then 2^p-1=2^{47}-1=140737488355327 factors as 2351 * 4513 * 13264529

If p=53, then 2^p-1=2^{53}-1=9007199254740991 factors as 6361 * 69431 * 20394401

If p=59, then 2^p-1=2^{59}-1=576460752303423487 factors as 179951 * 3203431780337

If p=61, then 2^p-1=2^{61}-1=2305843009213693951 factors as 2305843009213693951

If p=67, then 2^p-1=2^{67}-1=147573952589676412927 factors as 193707721 * 761838257287

If p=71, then 2^p-1=2^{71}-1=2361183241434822606847 factors as 228479 * 48544121 * 212885833

If p=73, then 2^p-1=2^{73}-1=9444732965739290427391 factors as 439 * 2298041 * 9361973132609

If p=79, then 2^p-1=2^{79}-1=604462909807314587353087 factors as 2687 * 202029703 * 1113491139767

If p=83, then 2^p-1=2^{83}-1=9671406556917033397649407 factors as 167 * 57912614113275649087721

If p=89, then 2^p-1=2^{89}-1=618970019642690137449562111 factors as 618970019642690137449562111

If p=97, then 2^p-1=2^{97}-1=158456325028528675187087900671 factors as 11447 * 13842607235828485645766393

Certainly this isn't always giving us primes, but it's not bad.  You can help the world search for more such Mersenne primes if you leave your personal computer on and connected to the Internet, via the Great Internet Mersenne Prime Search (GIMPS).  Random computers in labs at the University of Central Missouri and UCLA have found some of the largest known primes this way.  Currently, there hasn't been a new one found in over a year, but expect that to change soon.  The largest known such primes are VERY large - well over ten million digits! - with a monetary prize having been awarded to GIMPS for finding these huge ones; they shared it with many of the people who made it possible.

The reason this is conceivable is because of a special test, known as the Lucas-Lehmer test, which applies just to numbers $2^p-1$.  One start with $x_0=4$, and then keeps feeding $x$ into the loop $$x_{n+1}=\text{residue of }x_n^2-2\text{ modulo }2^p-1\; ;$$ if the number you get after doing this $p-2$ times is still divisible by $2^p-1$ (i.e., zero), then $2^p-1$ is in fact prime.  

With $p=5$ and $2^p-1=31$, we would start with $x_0=4$; doing it $5-2=3$ times gives

  1. $4^2-2=14$ modulo $31$ is $14$
  2. $14^2-2=194$ modulo $31$ is $8$
  3. $8^2-2$ modulo $31$ is $0$

And of course $31$ is indeed prime.

@interact def _(p=(71,prime_range(100))): test = 4 num = 2^p-1 for i in range(p-2): test=(test^2-2)%num html("The test says "+str(bool(test==0))) html("And in fact $2^{%s}-1=%s$ primality is "%(p,num)+str(is_prime(num))) 
       

Click to the left again to hide and once more to show the dynamic interactive window

 

Again, proving this is slightly beyond our capabilities right now.  Once again, we can prove the lesser result that Mersenne numbers are coprime.

Proposition:

$$\text{Mersenne numbers }2^p-1\text{ and }2^q-1\text{ with coprime exponents are themselves coprime.}$$

  • By way of contradiction, let $d>1$ be the gcd of the two numbers $2^p-1$ and $2^q-1$.
  • Let's investigate the order of $2\neq 1$ in $U_{d}$.  (Why is 2 even in this group?)
  • By definition of divisibility, $$2^p\equiv 1\text{ mod }(d)\text{ and }2^q\equiv 1\text{ mod }(d)$$
  • We know that $2^k\equiv 1$ means that $k$ is a multiple of $|2|$ by our group theory.  (For instance, use Lagrange's theorem on the group generated by $2$.)
  • So $p$ and $q$ are both multiples of $|2|$, which means that $|2|=1$.  This is a contradiction, so our assumption that $d>1$ was wrong.

Primality testing is full of little tidbits like this, 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 I hope to share later in the course).  

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 you use enough tests that do not guarantee primality but have a quite low false positive rate in practice, then the number you have is more likely to be a prime than the chance that your computer made an arithmetic error.   This is astonishing, but true; and if you end up with a number that likely to be prime, you can always check it with one of the various slower tests I will not describe.

 

@interact def _(p=(11,prime_range(100))): P=matrix_plot(matrix(p,[mod(a,p)^b for a in [1..p] for b in srange(1,p+1)]),cmap='gist_earth') P.show() 
       

Click to the left again to hide and once more to show the dynamic interactive window

Notice again here that Fermat's Little Theorem is visible in the second-to-last column.  

The last column is a slight restatement of Fermat's Little Theorem.  It turns out that of course if $a\equiv 0$ mod ($p$) then $a^{p-1}$ can't be $1$, but no matter what $a$ is, it is always true that $$a^p\equiv a\text{ mod }(p)\, .$$  

What makes this useful is that we can now use it to state a test for possible primality:

  • Fact: 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$, it's at least possible that $n$ is prime.
  • Definition: If $a^n\equiv a$ mod ($n$), we say $n$ passes the base $a$ test.
  • Definition: If $a^n\equiv a$ mod ($n$) but $n$ is not prime, we say it is a pseudoprime base $a$.

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.

It turns out that everyone from the ancient Chinese (e.g. the times of Master Sun from the CRT) to Leibniz used this test for the base $a=2$ to assert numbers are prime.  And it doesn't do a bad job:

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

Click to the left again to hide and once more to show the dynamic interactive window

We can change the numbers in the range above to check for more - say up to 1000.  Are there any that satisfy this test and are not prime?

To the surprise of many in the world of numbers, $n=341$, $n=561$, and $n=645$ turn out to fall in that category!  

html("We factor $341$ and get $%s$"%factor(341)) html("We factor $561$ and get $%s$"%factor(561)) html("We factor $645$ and get $%s$"%factor(645)) 
       
We factor 341 and get 11 * 31
We factor 561 and get 3 * 11 * 17
We factor 645 and get 3 * 5 * 43
We factor 341 and get 11 * 31
We factor 561 and get 3 * 11 * 17
We factor 645 and get 3 * 5 * 43

That's still not bad - out of 171 total such potential pseudoprimes base 2, only 3 of them actually are not prime, or about one and three quarters percent!  

Interestingly, it turns out that there are infinitely many such pseudoprimes:

Fact:

$$\text{If }n\text{ is a pseudoprime (base 2), then so is }2^n-1\; .$$  We will get this as a corollary of something stronger below.  All the Fermat and Mersenne numbers pass this test, incidentally, though they are all quite large compared to a typical number you might try.

If we want to check things out more carefully, of course, we can try to test for primality with a different base:

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

3^{561}\equiv 3\text{ mod }(561)

3^{645}\equiv 108\text{ mod }(645)
3^{341}\equiv 168\text{ mod }(341)

3^{561}\equiv 3\text{ mod }(561)

3^{645}\equiv 108\text{ mod }(645)

As you can see, this exposes 341 and 645 as fakes.  What about 561?  Let's try with base 5 as well.

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

Click to the left again to hide and once more to show the dynamic interactive window

Hmm, that's interesting.  What if I add some primes, like 7 or 11?  Try it above.

In the next cell, I get systematic.  We should expect output if $561$ doesn't pass the test.

@interact def _(p=(5,prime_range(1000))): 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: html("$%s^{561}\equiv %s\\text{ mod }(561)"%(pr,mod(pr,561)^561)+"</p>") 
       

Click to the left again to hide and once more to show the dynamic interactive window

It appears that $p^{561}\equiv p$ mod 561 for every prime $p$!  This can in fact be proved.

  • By the CRT, $$a^{561}\equiv a\text{ mod }(561)$$ precisely if this congruence holds for all the prime powers factors of 561.
  • We know that $$561=3\cdot 11\cdot 17\; ,$$ so these are the ones to check.  
  • The exponents for these congruences live in the mod ($\phi(p)$) world, so we just check what $561$ is in each of those worlds.
    • $$561\equiv 1\text{ mod }(16=17-1)$$
    • $$561\equiv 1\text{ mod }(2=3-1)$$
    • $$561\equiv 1\text{ mod }(10=11-1)$$
  • That is, for $p=3,11,17$ we have that $$a^{561}\equiv a^{1}\text{ mod }(p)$$.
  • Thus, by the CRT, this formula is always true!

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

factor(561) 
       
3 * 11 * 17
3 * 11 * 17

This example suggests that to find a Carmichael number, 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$.  It turns out that this is true - in fact, we can prove something even more accurate:

Proposition (Korselt's Theorem):

Carmichael numbers are precisely those composite $n$ for which $n$ is a product of all distinct primes $p_i$ (no squares) $$n=p_1p_2p_3\cdots p_k\text{ with }p_i\neq p_j$$ such that $$p_i-1\mid n-1$$ for all the prime factors.

Certainly prime numbers satisfy this trivially!  In this case, 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 that is the entire proof that such numbers are Carmichael numbers.  

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$.  

Here's another example:

n=29341 html("$%s$ is composite with factorization $%s$, but"%(n,factor(n))+"</p>") for fact,pow in factor(n): html("$%s^{561}\equiv %s\\text{ mod }(561)"%(fact,mod(fact,n)^n)) print 'but' for fact,pow in factor(n): html("$%s\equiv %s\\text{ mod }(\phi(%s)=%s)"%(n,mod(n,euler_phi(fact)),fact,euler_phi(fact))+"</p>") 
       
29341 is composite with factorization 13 * 37 * 61, but

13^{561}\equiv 13\text{ mod }(561) 37^{561}\equiv 37\text{ mod }(561) 61^{561}\equiv 61\text{ mod }(561) but 29341\equiv 1\text{ mod }(\phi(13)=12)

29341\equiv 1\text{ mod }(\phi(37)=36)

29341\equiv 1\text{ mod }(\phi(61)=60)

29341 is composite with factorization 13 * 37 * 61, but

13^{561}\equiv 13\text{ mod }(561) 37^{561}\equiv 37\text{ mod }(561) 61^{561}\equiv 61\text{ mod }(561) but 29341\equiv 1\text{ mod }(\phi(13)=12)

29341\equiv 1\text{ mod }(\phi(37)=36)

29341\equiv 1\text{ mod }(\phi(61)=60)

Two good exercises the book gives us are to check that $1729$ and $2821$ are also Carmichael numbers, and to find one of the form $7\cdot 23 \cdot p$ for a prime $p$.

In fact, as was proved about 20 years ago, there are (sadly?) infinitely many Carmichael numbers.  So now what?

To answer this, we turn to another result visible in our modular power graphic.

@interact def _(p=(11,prime_range(100))): P=matrix_plot(matrix(p-1,[mod(a,p)^b for a in [1..p-1] for b in srange(1,p)]),cmap='gist_earth') P.show() 
       

Click to the left again to hide and once more to show the dynamic interactive window

Again, Fermat's Little Theorem is the right-hand column.  What's that pattern in the middle?

Theorem:

$$a^{(p-1)/2}\equiv \pm 1\text{ mod }(p)\text{ for any prime modulus }p$$  

  • Since $a^{p-1}\equiv 1$,
  • Since $a^{(p-1)/2}$ is a solution to $x^2\equiv 1$,
  • Since $x^2\equiv 1$ factors as $p\mid x^2-1=(x+1)(x-1)$, 
  • Since, given that $p$ is prime, that means $p\mid x-1$ or $p\mid x+1$,
  • Then $x\equiv \pm 1$, so $a^{(p-1)/2}\equiv 1$ mod ($p$) !

What is the use for us of this theorem?  If we are testing $n$ for primality but $$a^{(n-1)/2}\not\equiv \pm 1\text{ mod }(n)\; ,$$ then that number is definitely not prime. 

Let's try this on our pesky Carmichael number, once again starting with base $a=2$.

mod(2,561)^((561-1)/2) 
       
1
1

Not again!  Try another base - maybe $a=3$?

mod(3,561)^((561-1)/2) 
       
441
441

Phew!  So this does help us test at least a little better.

A slightly stronger variant of this test is called Miller's test base $a$ for primality.  

  • Start by taking $n-1$, dividing by two, and checking the power $$a^{(n-1)/2}\text{ mod }(n)\; .$$  If $n$ actually is prime, then this will be $\pm 1$.  If the result is $-1$ we say it passes for sure.
  • If it hasn't passed yet, but the power is $+1$, continue dividing the power by two and then taking $a$ to that power; if $n$ is actually prime, then the result will be $\pm 1$, and again we say $n$ passes if the result is $-1$.
  • Continue - divide by two, check the same thing.
  • If we have divided $n-1$ by all possible powers of two, then we say $n$ passes as long as $a$ to that power is $\pm 1$; certainly this will be true if $n$ is prime and we haven't already had $n$ pass yet. 

Here is another number pseudoprime base 2 - but it does not pass this test, which is good since it's composite.

n=1387 html("We know $%s$ is composite because it factors as $%s$"%(n,factor(n))) html("Let's check $2^{(%s-1)/2}$ modulo $%s$: it's $%s$"%(n,n,mod(2,n)^((n-1)/2))) 
       
We know 1387 is composite because it factors as 19 * 73
Let's check 2^{(1387-1)/2} modulo 1387: it's 512
We know 1387 is composite because it factors as 19 * 73
Let's check 2^{(1387-1)/2} modulo 1387: it's 512

Looking good... But let's try another such pseudoprime, just to be sure.

n=2047 html("We know $%s$ is composite because it factors as $%s$"%(n,factor(n))) html("Let's check $2^{(%s-1)/2}$ modulo $%s$: it's $%s$"%(n,n,mod(2,n)^((n-1)/2))) 
       
We know 2047 is composite because it factors as 23 * 89
Let's check 2^{(2047-1)/2} modulo 2047: it's 1
We know 2047 is composite because it factors as 23 * 89
Let's check 2^{(2047-1)/2} modulo 2047: it's 1

As we can see, this shows that $n=2047$ passes the first part of Miller's test base 2, and that there is no further to go because $(2047-1)/2=1023$ is odd.  So, as far as we know thus far, $2047$ is prime.  

This is a little frustrating!  Composite numbers can pass this primality test too.  So we come up with another name.

Definition: We call a composite number $n$ that passes Miller's test base $a$ a strong pseudoprime base $a$.  

The bad news is that strong pseudoprimes exist.  In fact, we can prove a theorem about it:

Theorem: If $n$ is a pseudoprime base $2$, then $2^n-1$ is a strong pseudoprime base $2$.

  • Let $n$ be composite and odd, but $$2^{n}\equiv 2\text{ mod }(n)$$
  • Then we can cancel 2, and get $$2^{n-1}\equiv 1\text{ mod }(n)$$
  • So $$2^{n-1}-1=nk\text{ for some (odd) integer }k$$
  • So $$\left(2^{n}-1\right)-1=2^n-2 =2\left(2^{n-1}-1\right)=2nk$$
  • Now apply Miller's test to $2^n-1$.
    • $$2^{((2^{n}-1)-1)/2}=2^{nk}=\left(2^n\right)^k\equiv 1^k\equiv 1\text{ mod }(2^n-1)$$
    • It passes Miller's test.
  • All that remains is show $2^n-1$ is composite if $n$ is composite; this is the lemma we used in proving that Fermat numbers were coprime, so we are done.

This has two corollaries:

  • Since $(\pm 1)^2=1$, we have proved that if $n$ is a pseudoprime base $a$, so is $2^n-1$
  • There are infinitely many strong pseudoprimes base 2.  

For instance, we now know that $2^{341}-1$ must fall in that category, and since the second number below is odd, this confirms it:

n=2^341-1 print mod(2,n)^((n-1)/2),(n-1)/2 
       
1
223974474217780421055744228056844427812164549723464953489998910096379187\
1180160945380877493271607115775
1 2239744742177804210557442280568444278121645497234649534899989100963791871180160945380877493271607115775

BUT there are not strong Carmichael numbers!  In fact: 

Theorem:

If $n$ is an odd composite positive integer, then $n$ passes Miller's test for at most $(n-1)/4$ bases $a$ between 1 and $n-1$.

The proof of this is doable for us.  It counts numbers of solutions of $x^\ell-1$ modulo various prime powers and combines them with the CRT to give a good counting argument.  It is long enough to omit here, though.

Needless to say, no one could compute that many bases to prove primality!  But Rabin used this fact to suggest a test for a probable prime with probability of failure less than $(\frac{1}{4})^k$ for any desired $k$; simply run Miller's test for $k$ different bases less than $n-1$, and if a number passes all of them, that is the probability.  

For 100 primes, this is the probability that would come out.

(1./4)^100 
       
6.22301527786114e-61
6.22301527786114e-61

So if you run the test for 100 primes, you are in pretty decent shape.  

You can also always use some slow test to prove primality.  That is what is called a certificate of primality, and although you may not believe it, programs that reliably generate reasonably large (100-200 digits, right now) primes and can verify it are hot items on the virtual shelves of those who care about such things.  

Finally, let's see this in action.  Remember that we wanted keys larger than 1024 bits for at least a semblance of security in RSA?  Here we go with a start:

p=next_probable_prime(randrange(2^1024)) q=next_probable_prime(randrange(2^1024)) n=p*q print p,q,n 
       
133731079371886953587741222083781152978109465241569454709483753758437097\
157524648618405567429736328245111949726102618480846163359131881238397522\
465408746143098984759036965241016744842409092448605381821281981700532335\
533319937788464622880315050991622802890119477612137433256064692848313835\
437231389442427180867
100073944325582449386487435903499304985342069793973786815972644980328009\
677997238912232151257808561069433669958157190579068464095068937262429989\
860294296973577466087788054412554377311666557138304811293677718085830459\
922920018283609804729362028789383130538222444584313373678299500264538587\
742341729178750665511
133829965916622625486679322125979389452067414613190039529322107447545642\
350459504726220211312426068980551491058852922450324907547090473709717609\
185596762507056630748423257113594044460822524064092037021851441328955661\
740812482520295240308652546291452722247397877365374892115617739289746389\
468864768240925264837310849099788951435508569839862479594743331833972537\
914175034166484597996164875824655956692777703145717762123873427062634254\
158069621281408713500893090348125935774600220579410033860478284982298426\
539491984227874696828379529917239847091827529716770691000886696600853237\
71243791647307400860537576842131815978037
133731079371886953587741222083781152978109465241569454709483753758437097157524648618405567429736328245111949726102618480846163359131881238397522465408746143098984759036965241016744842409092448605381821281981700532335533319937788464622880315050991622802890119477612137433256064692848313835437231389442427180867 100073944325582449386487435903499304985342069793973786815972644980328009677997238912232151257808561069433669958157190579068464095068937262429989860294296973577466087788054412554377311666557138304811293677718085830459922920018283609804729362028789383130538222444584313373678299500264538587742341729178750665511 13382996591662262548667932212597938945206741461319003952932210744754564235045950472622021131242606898055149105885292245032490754709047370971760918559676250705663074842325711359404446082252406409203702185144132895566174081248252029524030865254629145272224739787736537489211561773928974638946886476824092526483731084909978895143550856983986247959474333183397253791417503416648459799616487582465595669277770314571776212387342706263425415806962128140871350089309034812593577460022057941003386047828498229842653949198422787469682837952991723984709182752971677069100088669660085323771243791647307400860537576842131815978037

The $p$ and $q$ we get above are just probable primes.  Verifying them could take a little longer!

time is_prime(p) 
       
True
Time: CPU 7.00 s, Wall: 7.15 s
True
Time: CPU 7.00 s, Wall: 7.15 s
time is_prime(q) 
       
True
Time: CPU 6.83 s, Wall: 6.84 s
True
Time: CPU 6.83 s, Wall: 6.84 s