In his e-mail message of 08-Sep-92 Andries E. Brouwer writes:

gap> IsPrime(274177*67280421310721);
true

This is  obviously a  bug.  Let me discuss  it a  little  bit.  Note that
274177 * 67280421310721 is 2^64 + 1.

First note that the problem is not that GAP cannot *factorize* Fermat
non-primes. GAP's function 'Factors' doesn't even *try* to factor this
integer, because it believes it is a prime. That is to say that the
problem is in 'IsPrime' not in 'Factors'. If 'IsPrime' is changed so
that it returns the correct answer, 'Factors' has no problems at all to
split this integer.

Now let me describe how 'IsPrime( N )' works. First it checks that N
does not have any prime factors less than 1000. 2^64+1 obviously passes
this test.

Next 'IsPrime' tests whether N is a (strong) pseudoprime to base 2.  That
is, 'IsPrime' tests whether  2^(N-1) = 1 mod  N.   In  this case we  have
2^128 = 1 mod (2^64+1), so 2^64+1 passes this test too.

Finally 'IsPrime' performs another similar test, however this  time using
GF(N^2),  instead  of  GF(N)  as above.  To do  this  'IsPrime' selects a
polynomial X^2 - P X + 1, and works in (Z_N)[X]/(X^2 - P X + 1).  If N is
a prime and X^2 - P  X  + 1 is irreducible,  this is the  field  with N^2
elements.  The multiplicative group of this field has order N^2-1, so any
element of this field has an order that divides N^2-1.  So 'IsPrime' test
whether X^(N^2-1) is trivial in the ring (Z_N)[X]/(X^2  - P X + 1).   The
hope is  of  course that  X  will have an order that is not  a divisor of
N^2-1 if N is not a prime.

Now I have to say how X^2 - P X + 1, i.e.,  the parameter P, is selected.
'IsPrime' takes the smallest  P such that Jacobi( P^2-4,  N ) = -1.  If N
is a  prime  this  implies  that X^2 - P X +  1 is  irreducible modulo N.
Because Jacobi( -3, 2^64+1  )  =  -1, 'IsPrime' takes P = 1 in this case.
But X^2 - X + 1 is the sixth cylotomic polynomial, thus the element X has
order  6 in Z[X]/(X^2 - X  + 1), and  thus an  order dividing  6  in  any
(Z_N)[X]/(X^2 - X + 1).  So the whole test  comes down to testing whether
6 divides N^2-1, which is equivalent to gcd( N, 6 ) = 1.

The corresponding problem for the (strong) pseudoprime test would be to
test if N is a (strong) pseudoprime w.r.t. the base -1.

The fix is very simple of course. Instead of selecting the first
positive p such that Jacobi( p, N ) = -1, 'IsPrime' should take the first
such p that is greater than 1.

If somebody wants to fix the problem right away, here is the necessary
change.

--- integer.g   1992/03/19 18:16:39
+++ integer.g   1992/09/10 12:38:07
@@ -344,7 +347,7 @@
     fi;

     # find a quadratic nonresidue $d = p^2/4-1$ mod $n$
-    p := 1;  while Jacobi( p^2-4, n ) <> -1  do p := p+1;  od;
+    p := 2;  while Jacobi( p^2-4, n ) <> -1  do p := p+1;  od;

# for a prime $n$ the trace of $(p/2+\sqrt{d})^n$ must be $p$
# and the trace of $(p/2+\sqrt{d})^{n+1}$ must be 2

Martin.

PS: Now comes the strange part. When I wrote 'IsPrimeInt' I knew about
this problem, that certain polynomial would lead to very weak tests, and
should be avoided. How could I fail to see that x^2 - x + 1 is such a
polynomial? Well this is more than 2 years ago, so my memory is a bit
hazy. Still it leaves me puzzled.

--
Martin Sch"onert, Martin.Schoenert@Math.RWTH-Aachen.DE, +49 241 804551
Lehrstuhl D f"ur Mathematik, Templergraben 64, RWTH, D 51 Aachen, Germany