# Michael Monagan sent this bug fix of the isprime command to MUG:
# 'There is a bug in the Release 3 library routine isprime.
# On squares it goes into an infinite loop.  The bug is fixed
# in Release 4 and I think in one of the Release 3 patch libraries.
# Here is the source code for anyone who needs to install the fix -- just
# read in the code below.'

#-->	isprime(n)
#
#	Probabalistic primality test.	
#	Reference:  Knuth Vol 2, 2nd edition "Semi Numerical Algorithms".
#		Section 4.5.4, algorithm P
#
#	Version (May 1993): Returns false if n has a factor smaller
#	than 1000 or fails a strong (*) Fermat test with base 2 or
#	fails one Lucas test.
#
#	(*) we use a modified strong pseudo prime test, which is more
#	powerful than the Miller-Rabin algorithm described in Knuth.
#	For more details, see a future article on the Maple Technical
#	Newsletter.
#
#	Presently there are no composite numbers known to us that will
#	make isprime() return true.
#
#	The code for the Lucas test was written originally by Martin
#	Schoenert.
#
#
############################################################################
#
#			G.H. Gonnet (December 82)
#			M.B. Monagan (Updated 83, 84, 94)
#			G.H. Gonnet (October 88, extended test)
#			Gaston H. Gonnet (May 1993)
#                       $Notify: gonnet@inf.ethz.ch$
#
alias( TraceModQF = `isprime/TraceModQF` ):
alias( cyclotest = `isprime/cyclotest` ):
alias( sumxtor = `isprime/sumxtor` ):
alias ( w =`isprime/w` ):

w := [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]:

unprotect('isprime'):
isprime := proc(n)
local btor, nr, p, r;
option remember,system,
   
       `Copyright (c) 1993 Gaston Gonnet, Wissenschaftliches Rechnen, ETH Zurich. All rights reserved.`;

   if not type(n,integer) then
       if type(n,numeric) then
	     ERROR(`argument must be an integer`)
       else RETURN('isprime(n)')
            fi
       fi;
   if n < 2 then false
   elif has(w,n) then true

   #	Product of all primes in the range [1,100]
   elif igcd(2305567963945518424753102147331756070,n) <> 1 then false
   elif n < 10201 then true

   #	Product of all primes in the range [100,1000]
   elif igcd( 849696948923341811053233990918734996592606258664893273661154\
54263422038932707693909090694773095091375097869171186680288614993338250\
97682386722983737962963066757674131126736578936440788157186969893730633\
1306647862044862494925732402262739543736363903875260816675866125595683\
46306972204475122988482222285500626837863425199602259963013159456444700\
64720696621750477244528915927867113,n) <> 1 then
       false
   elif n < 1018081 then true
   else nr := igcd(210^5,n-1);
	 nr := igcd(nr^5,n-1);
	 r := iquo(n-1,nr);

	 # modified strong Fermat test
	 btor := modp( power(2,r), n );
	 if cyclotest(n,btor,2,r) = false or
	     irem(nr,3)=0 and cyclotest(n,btor,3,r) = false or
	     irem(nr,5)=0 and cyclotest(n,btor,5,r) = false or
	     irem(nr,7)=0 and cyclotest(n,btor,7,r) = false then RETURN(false)
	     fi;

	 # MBM (Nov/94) ensure that n is not a perfect square otherwise
	 # J(a/n)=J(a/m^2)=J(a/m)*J(a/m)=1 or 0, i.e. J(a/n) <> -1
	 if isqrt(n)^2=n then RETURN(false) fi;

	 # Lucas test
	 for p from 3 while numtheory[jacobi]( p^2-4, n ) <> -1 do od:
	 evalb( TraceModQF( p, n+1, n ) = [ 2, p ] );
        fi
end:
protect(isprime):


cyclotest := proc(n,btor,i,r)
local nr, k, t, s, x;
option 
   `Copyright (c) 1993 Gaston Gonnet, Wissenschaftliches Rechnen, ETH Zurich. All rights reserved.`;

nr := iquo(n-1,r);
for k from 0 while irem(nr,i,'t') = 0 do nr := t od;
x := modp(power(btor,nr), n);
if x=1 then RETURN(FAIL) fi;
to k do
   s := sumxtor(x,i,n);
   if s[1]=0 then RETURN(FAIL) elif s[2]=1 then RETURN(false) fi;
   x := s[2]
   od;
false
end:

sumxtor := proc(x,r,n)
local i1, i2, i3, r1, r2, r3, s, xj;
option 
   `Copyright (c) 1993 Gaston Gonnet, Wissenschaftliches Rechnen, ETH Zurich. All rights reserved.`;

if r < 1 then ERROR(`r is too small`)
elif r=1 then [1,x]
elif r < 6 then
    s := 1+x;
    xj := modp(x^2,n);
    to r-2 do s := s+xj; xj := modp(xj*x,n) od;
    [modp(s,n),xj]
else i1 := isqrt(r);
    i2 := iquo(r,i1,'i3');
    if i3=0 then
         r2 := procname(x,i2,n);
         r1 := procname(r2[2],i1,n);
         [modp(r2[1]*r1[1],n), r1[2]]
    else r3 := procname(x,i3,n);
         r2 := procname(x,i2,n);
         r1 := procname(r2[2],i1,n);
         [modp( r3[1] + r3[2] * modp(r2[1]*r1[1],n), n), modp(r3[2]*r1[2],n) ]
         fi
    fi;
end:

#############################################################################
##
##  Then isprime checks that n is a Lucas pseudoprime for p, choosen
##  so that the discriminant d=p^2/4-1 is an  quadratic nonresidue mod n.
##  I.e., isprime takes the root a = p/2+\sqrt{d} of x^2 - px + 1 in
##  the  ring Z_n[\sqrt{d}] and computes the  traces of a^n and a^{n+1}.
##  If n is a prime, this  ring is the field of  order n^2 and raising to
##  the nth power is conjugation, so trace(a^n)=p and trace(a^{n+1})=2.
##  However, these identities hold only for extremly few composite numbers.
##
##  Note that  this  test  for trace(a^n) = p  and  trace(a^{n+1}) = 2 is
##  usually formulated using the Lucas sequences  U_k = (a^k-b^k)/(a-b) and
##  V_k = (a^k+b^k)=trace(a^k), where one tests U_{n+1} = 0, V_{n+1} = 2.
##  However, the trace test is equivalent and requires fewer multiplications.
#############################################################################

TraceModQF := proc ( p, k, n )
   local  i, kk, trc, v;
   option 
       `Copyright (c) 1993 Gaston Gonnet, Wissenschaftl. Rechnen, ETH Zurich. All rights reserved.`;
   kk := k;
   for i while kk>1 do kk := iquo(kk+1,2,v[i]) od;
   trc := [ p, 2 ];
   for i from i-1 by -1 to 1 do
	if v[i]=1 then
      	     trc := modp( [trc[1]^2 - 2, trc[1]*trc[2] - p], n );
   	else trc := modp( [trc[1]*trc[2] - p, trc[2]^2 - 2], n );
   	     fi
	od;
   trc
end:

#savelib('isprime','w','cyclotest','sumxtor',\
   'TraceModQF','`isprime.m`'):