[kujswynr] Revisiting Fermat 20

We repeat testing the primality of the twentieth Fermat number on a slightly faster computer.

? f=2^2^20+1;
? ispseudoprime(f)
time = 9h, 20min, 50,056 ms.
0

This computation can be sped up considerably by adding a flag:

? ispseudoprime(f,1)
time = 4h, 46min, 24,280 ms.
0

Without the flag, Pari/GP does the BPSW primality test, which begins with a Miller-Rabin test with base 2.  But Fermat numbers are special: they always pass Miller-Rabin base 2, regardless of whether they are prime or composite.  Interestingly, base-2 Miller-Rabin passes (false positive) extremely quickly.

? lift(Mod(2,f)^(f-1))
time = 9,917 ms.
1

This is because modular squarings number #22 through #1048576 are all just computing 1*1=1.

Without the flag, after base-2 Miller-Rabin passes, BPSW moves on to the Lucas probable prime test (so named because it uses Lucas sequences, not be confused with the totally different Lucas primality test which relies on the factorization of p-1), which is more complicated than Miller-Rabin so takes 9 hours.

With the flag 1, ispseudoprime instead does 1 Miller-Rabin trial with a random base (not 2, probably).  Unless we are unlucky, it will report that F20 is definitely composite.

Performing the Pepin test, which isn't probabilistic, is pretty easy to do in Pari/GP.  It does take some effort to look up the formula and interpret the result.

? r=Mod(3,f)^((f-1)/2);
time = 3h, 48min, 14,148 ms.
? r+1==0
0

We're not sure why the Pepin test was faster than Miller-Rabin to a random base, as they both seem to be the same amount of work.

Incidentally, the Pepin test is identical to the (second) Lucas primality test mentioned above which relies on the factorization of p-1.  p-1 is easy to factor for Fermat numbers.

We can verify the final residue with the remainders published in "The Twentieth Fermat Number is Composite" by Jeff Young and Duncan A. Buell.  The values published in their paper are in octal, a fact which they confusingly do not explicitly state.

? printf("%o",lift(r)%(2^36))
175517362761
? printf("%o",lift(r)%(2^36-1))
411337412531
? printf("%o",lift(r)%(2^35-1))
161572365764

The values are in octal probably following a tradition established by Selfridge and Hurwitz for F14.  However, Selfridge and Hurwitz do explicitly state in their paper that their values are in octal.

For reference, here is the final F20 residue modulo several bases in decimal:

? lift(r)%(2^36)
16865158641
? lift(r)%(2^36-1)
35626292569
? lift(r)%(2^35-1)
15265819636