
(*:Version: Mathematica 2.0 *)

(*:Name: NumberTheory`PrimeQ` *)

(* :Title:  PrimeQ *)

(* :Author: Ilan Vardi *)

(*:Keywords: 
    Primality Proving, Certificate of Primality, Pratt's Certificate,
    Elliptic Curves, Complex Multiplication, Goldwasser-Kilian,
    Atkin-Morain
*)


(* :Source: Francois Morain, ``Implementation of the Atkin-Goldwasser-Kilian
            Primality Testing Algorithm,'' INRIA Research Report #911,
            October 1988.

            A.O.L. Atkin and F. Morain, ``Elliptic Curves and 
            Primality Proving,'' to appear, Mathematics of Computation.

            D.A. Cox, ``Primes of the Form x^2 + n y^2'', Wiley 1989.

            Stan Wagon, ``Mathematica in Action,'' Freeman 1991.

            D. Bressoud, ``Factorization and Primality Testing,''
            Springer Verlag 1989.
*)

(* :Warning: 

  
    ProvablePrimeQ[p] has to wait until PrimeQCertificate has
    finished to return its answer, i.e., until a certificate has been
    generated. Since PrimeQCertificate caches its value, running 
    PrimeQCertificate[p] after ProvablePrimeQ[p] will not take 
    any extra time. (There have been problems with caching of 
    certificates.)

    ProvablePrimeQ uses PrimeQ and only proves that the old PrimeQ 
    has returned a correct answer. If PrimeQ is incorrect, i.e., 
    it claims that a composite number is prime, then ProvablePrimeQ will 
    either not return a value (will go in an 
    infinite loop) or else it will find a nontrivial divisor of a number 
    that is claimed to be prime by PrimeQ. Thus ProvablePrimeQ can
    be used to find counterexamples to PrimeQ.
*)


(* :Limitation:  The Elliptic Curves method will not work for very 
                 small primes, so SmallPrime indicating when this
                 method takes over from Pratt's method has a minimum
                 value of 1000.

                 It looks like this method will return an answer in 
                 reasonable time for inputs of up to 200 digits, i.e.,
                 this takes about 6 hours on a DEC workstation.
                 
*)  

(* :Discussion: 

      This package does primality proving. In fact, the package will
      produce a certificate of primality, i.e., a short set of data
      that is a witness to the primality of a prime. This certificate
      can easily be checked in the sense that the theory involved in 
      implementing a check is much easier than the theory that produced
      the certificate. This is analogous to the certificate of 
      compositeness given by exhibiting a nontrivial factorization.
      The package also includes a certificate of compositeness, though 
      it is more subtle than just giving a factor.
    
*********************************************************************
               Pratt's certificate of Primality
*********************************************************************

      For small prime p <= SmallPrime use Pratt's certificate of primality:

      Pratt's certificate gives a proof that a is a generator 
      (primitive root) of the multiplicative group (mod p) which,
      along with the fact that a has order p-1, proves that p is a
      prime.  You need to factor p-1, which is feasible for p with 
      20 digits or less. 
     
      The certificate shows that PowerMod[a, p-1, p] = 1 and 
      PowerMod[ a, (p-1)/q, p] != 1 for all primes q dividing p-1
      (=> a is a primitive root). These primes q are the children of 
      {p,a} in the tree and are proved prime recursively. You have
      to check that all prime factors of p - 1 are included.

      This part of the program uses some suggestions of Ferris Wheeler
      as in the PrimitiveRoot[]  program in NumberTheoryFunctions.m
      package.

      See Stan Wagon's book for details on this test.

***************************************************************
    Atkin-Goldwasser-Kilian-Morain certificate of primality
***************************************************************

      For large primes p > SmallPrime use the Atkin-Morain test:

      The Atkin-Morain test returns a recursive certificate for the 
      prime p. This consists of a list of:

 (i)  A point (CertificatePoint) on an Elliptic Curve: PointEC[x,y,g2,g3,p],
      in other words, a solution of the equation

          y^2 = x^3 + g2 x + g3  (mod p)

      for some numbers g2 and g3.

 (ii) A prime q with q > (p^(1/4) + 1)^2, such that for some
      other number k (CertificateK), and m = k q (CertificateM)
      such that k is not equal to 1 and such that 

           m PointEC[x,y,g2,g3,p] 

      is the identity on the curve, but 

           k PointEC[x,y,g2,g3,p]

       is not the identity. This guarantees primality of p
       by a theorem of Goldwasser and Kilian (see Lenstra, 
       Elliptic Curve algorithms, World Congress 1986). So one needs 
       only to show that q is also prime in order to show that p is 
       prime, and q is shown prime using the same method.

 (iii) Each q has its recursive certificate following it. So if 
       the smallest q is known to be prime, all the numbers are
       certified prime up the chain.

   Thus to check the certificate at level q, do

      CertificateK CertificatePoint /. 
        Select[PrimeQCertificate[p, opts], #[[1]] == q &] [[2]] 

      CertificateM CertificatePoint /. 
        Select[PrimeQCertificate[p, opts], #[[1]] == q &] [[2]] 

      CertificateNextPrime == CertificateM /CertificateK &&
        CertificateNextPrime > ((N[q]^(1/4) + 1)^2) /.
           Select[PrimeQCertificate[p, opts], #[[1]] == q &] [[2]]  

   The first should not be the identity, the second should be the
   identity and the third should be True.

   This algorithm is described in Francois Morain's INRIA technical
   report #911.  A more recent version is in the Atkin-Morain paper.

     
***************************************************************
              Certificate of Compositeness 
***************************************************************
   
   The certificate of compositeness shows that either 
   PowerMod[a, p - 1, p] != 1, or else a != 1 or -1 and
   PowerMod[a, 2, p] == 1, which proves that p is not a prime. 

   Applying PowerMod to PrimeQCertificate[p] proves p composite
   by checking this condition for the value of a given by the 
   certificate.

*)


BeginPackage["NumberTheory`PrimeQ`", "NumberTheory`NumberTheoryFunctions`"]  


ProvablePrimeQ::usage = "ProvablePrimeQ[n] determines if n
is prime or not and can provide a certificate (proof) of this."

Certificate::usage = "Certificate is an option to ProvablePrimeQ,
which prints a certificate for the primality or compositeness of n if
and only if Certificate -> True."

SmallPrime::usage = 
"SmallPrime is an option for ProvablePrimeQ, PrimeQCertificate, and
PrimeQCertificateCheck, which gives a lower bound for using the Atkin-Morain
primality test."

PrimeQCertificate::usage = "PrimeQCertificate[n] prints out a certificate
that n is prime or that n is composite."

TrialDivisionLimit::usage =
"TrialDivisionLimit is an option to ProvablePrimeQ and PrimeQCertificate that
specifies how many primes should be used in the trial division part of 
PrimeQCertificate[n]."

PollardPTest::usage = 
"PollardPTest is an option to ProvablePrimeQ and PrimeQCertificate which 
determines whether the Pollard p-1 factoring algorithm in the factorization
part of PrimeQCertificate[n] is used."

PollardRhoTest::usage =
"PollardRhoTest is an option to ProvablePrimeQ and PrimeQCertificate which 
determines whether the Pollard rho factoring algorithm in the factorization
part of PrimeQCertificate[n] is used."

PrimeQCertificateCheck::usage = "PrimeQCertificateCheck[cert, n]
checks that the certificate cert is a correct certificate for the
primality or compositeness of n. It returns True if the certificate
is correct. If it does not return True then it returns a counterexample
to PrimeQ[], i.e., a nontrivial factor of a number declared prime
by PrimeQ."

ModularInvariantj::usage = "ModularInvariantj[z] gives the j invariant
of the complex number z, Im[z] > 0."

HilbertPolynomial::usage = "HilbertPolynomial[d] gives the irreducible
polynomial of ModularInvariantj[(-1 + Sqrt[d])/2], where d < 0 is a 
square free integer of the form 4 n + 1. The splitting field is
the Hilbert Class Field (maximal unramified Abelian extension) of
Q(Sqrt[d])."

x::usage= "Indeterminate in the Hilbert polynomial."

PointEC::usage = "PointEC[a, b, g2, g3, p] gives a point
on the elliptic curve y^2 = x^3 + g2 x + g3 (mod p)."

PointECQ::usage = "PointECQ[PointEC[a, b, g2, g3, p]] returns True if
b^2 = a^3 + g2 a + g3 (mod p) so that the point is on the elliptic
curves that it claims it is."

CertificatePrime::usage = "CertificatePrime is used to certify primality."

CertificatePoint::usage = "CertificatePoint is used to certify primality."

CertificateK::usage = "CertificateK is used to certify primality."

CertificateM::usage = "CertificateM is used to certify primality."

CertificateNextPrime::usage = "CertificateNextPrime is used to certify primality."

CertificateDiscriminant::usage = "CertificateDiscriminant is used to certify primality."

PrimeQMessages::usage = "PrimeQMessages -> True is an option to
PrimeQCertificate which prints out the progress of the algorithm."


Begin["`Private`"]

Options[PrimeQCertificate] =
{SmallPrime -> 10^10, Certificate -> False, PrimeQMessages->False,
 PollardPTest -> Automatic, PollardRhoTest -> Automatic,
 TrialDivisionLimit -> Automatic}
 
Options[ProvablePrimeQ] = Options[PrimeQCertificate]

ProvablePrimeQ[p_Integer, opts___]:=
          If[!TrueQ[(Certificate /. {opts} /. Options[ProvablePrimeQ])],
               Block[{},
                 PrimeQCertificate[p, opts];
                 Return[PrimeQ[p]]
	       ],
	      (*  Or else get a certificate along with the answer.   *)
       	      {PrimeQ[p], PrimeQCertificate[ p, opts]}
          ]
        

PrimeQCertificate::false =
	"Warning: PrimeQCertificate has detected a counterexample to PrimeQ."
PrimeQCertificate::smprm = "Warning: SmallPrime set to 1000."
PrimeQCertificate::trdiv = "Warning: TrialDivisionLimit set to default ``."
PrimeQCertificate::pollp = "Warning: PollardPTest set to default ``."
PrimeQCertificate::pollr = "Warning: PollardRhoTest set to default ``."

PrimeQCertificate[p_Integer, opts___]:= PrimeQCertificate[p, opts]=
   Block[{a, q, r, s, exp, smallprime, flist,
          d = -7, qr,t,qplus,qminus,m,moverq,hr,ell,
          sqrtp = Round[Sqrt[N[p]]], 
          classtablelen = Length[ClassTable], tableq = True, 
          trialvector, triallen, tdr, qplustr, pminustr, tdmod,
          trialdivisionlimit, pollardptest, pollardrhotest,
          g2,g3,xec,yec2,p1,p2,point,c,di = 1, hp, curve, messages},


    Which[p == 2, 2, 

(* p Composite *)
          Not[PrimeQ[p]], 
(*
             Certificate of compositeness  

*)
   a = 1; q = (p - 1)/2;  s = -1;  exp = 2; 
    While[EvenQ[q], q /= 2]; 
          While[Mod[s + 1, p] == 0, 
                    a++; 
                     If[PowerMod[a, p - 1, p] != 1, 
                          exp = p - 1;
                          s = a;  
                          Break[], 
                            r = PowerMod[a, q, p];
                            If[r != 1, 
                            While[r != 1, 
                                  s = r; 
                                  r = PowerMod[s, 2, p] 
                                 ]
                              ]
                       ]
                     ];
              {s, exp, p},

(* Prime *) 

     PrimeQ[p], 

     smallprime = SmallPrime /. {opts} /. Options[PrimeQCertificate];
 
   (* restrict SmallPrime to >= 1000 *)

     If[!NumberQ[smallprime] || smallprime < 1000,
        Message[PrimeQCertificate::smprm];
        smallprime=1000
        ];

     (* Pratt's certificate *)

If[p <= smallprime,
   a = 2; flist = First /@ FactorInteger[p-1];
   While[JacobiSymbol[a, p] == 1 ||
         Scan[If[PowerMod[a, (p-1)/#, p] == 1,
                 Return[True]]&, flist],
                        a++]; 
                If[PowerMod[a, p-1, p] == 1, 
         {p, a, PrimeQCertificate[#, opts] & /@ flist}
                   ],
 
       (* Elliptic curves certificate  *)

(* First build a vector of p + 1 mod small primes in order to 
   facilitate trial division. *)

messages = PrimeQMessages /. {opts} /. Options[PrimeQCertificate]; 

If[messages, Print[p]; Print["Create trial division vector"]];

trialdivisionlimit = TrialDivisionLimit /. {opts} /. Options[PrimeQCertificate]; 
If[trialdivisionlimit =!= Automatic &&
   (!NumberQ[trialdivisionlimit] || Head[trialdivisionlimit] =!= Integer ||
    trialdivisionlimit <= 0),
	trialdivisionlimit = Which[p > 10^150, 20000,
				   p > 10^60, 10000, True, 5000];
    	Message[PrimeQCertificate::trdiv, trialdivisionlimit]];

If[trialdivisionlimit === Automatic, 
	trialdivisionlimit = Which[p > 10^150, 20000, 
                                   p > 10^60, 10000, True, 5000]]; 

trialdivisionlimit = Min[trialdivisionlimit, 
                         2 Sqrt[N[p]]/Log[N[p]]];

trialvector = Mod[p+1, Prime[Range[trialdivisionlimit]]]; 

triallen = Length[trialvector];

(* What kind of factoring tests should be used in LargeFactor. *)

pollardptest = PollardPTest /. {opts} /. Options[PrimeQCertificate]; 
If[!MemberQ[{Automatic, True, False}, pollardptest],
	pollardptest = (p < 10^150);
	Message[PrimeQCertificate::pollp, pollardptest]];
If[pollardptest === Automatic, pollardptest = p < 10^150];

pollardrhotest = PollardRhoTest /. {opts} /. Options[PrimeQCertificate]; 
If[!MemberQ[{Automatic, True, False}, pollardrhotest],
	pollardrhotest = (p < 10^60);
	Message[PrimeQCertificate::pollr, pollardrhotest]];
If[pollardrhotest === Automatic, pollardrhotest = p < 10^60];

If[messages, Print[{"Trial division vector done", triallen}]];  

    While[True, 
         While[JacobiSymbol[d,p] != 1, 
                If[di < classtablelen, di++; d = 1 - 4 ClassTable[[di]],
                   If[di == classtablelen || tableq, 
                      d = -7; di++, tableq = False;
                      {d, di} += {-4, 1}]];
                If[di > classtablelen, 
                   While[MemberQ[ClassTable, (1 - d)/4] || 
                         !SquareFreeQ[d], {d, di} += {-4, 1}]]
                ];
          If[di > classtablelen && ClassNumber[d] > 30, 
             {d, di} += {-4, 1};
             While[MemberQ[ClassTable, (1 - d)/4] && 
                   !SquareFreeQ[d], {d, di} += {-4, 1}];
             Continue[]];
           If[messages, Print[{di, d}]];
           qr= QuadraticRepresentation4[d, p];
           If[qr === False, 
              If[di < classtablelen, di++; d = 1 - 4 ClassTable[[di]],
                 If[di == classtablelen || tableq, 
                    d = -7; di++; tableq = False,
                    {d, di} += {-4, 1}]];
              If[di > classtablelen, 
                 While[MemberQ[ClassTable, (1 - d)/4] || 
                       !SquareFreeQ[d], {d, di} += {-4, 1}]];
                 Continue[]];                      
            t = Abs[2 qr[[1]]]; s = Abs[2 qr[[2]]];
            qplustr = p + 1 + t;
            qminustr = p + 1 - t;

(* Trial division of p + 1 + t, p + 1 - t, is done in parallel
   by comparing t mod small primes with the stored values 
   of p + 1 mod small primes in trialvector. This allows one 
   to do both trial divisions at once and by modular operations
   with a number of size sqrt[p]. See Atkin-Morain for details. *)

  If[messages, Print["trying Trial Division"]];  


tdmod = Mod[t, Prime[Range[triallen]]];

Scan[While[Mod[qminustr, #] == 0, qminustr = Quotient[qminustr, #]]&, 
     Prime[Flatten[Position[tdmod - trialvector, 0]]]];

Scan[While[Mod[qplustr, #] == 0, qplustr = Quotient[qplustr, #]]&,   
     Prime[Flatten[Position[
           MapIndexed[If[#1 == 0 || {#1} == Prime[#2],0,1]&,
                         tdmod + trialvector], 0]]]];
           qplus =  
           LargeFactor[qplustr, p+1+t, p, pollardptest, pollardrhotest];
            If[IntegerQ[qplus], m = p + 1 + t; q = qplus,
               qminus = 
               LargeFactor[qminustr, p+1-t, p, pollardptest, pollardrhotest];
               If[IntegerQ[qminus], m = p + 1 - t; q = qminus,
                  If[di < classtablelen, di++; d = 1 - 4 ClassTable[[di]],
                     If[di == classtablelen || tableq, 
                        d = -7; di++; tableq = False,
                      {d, di} += {-4, 1}]];
                      If[di > classtablelen, 
                         While[MemberQ[ClassTable, (1 - d)/4] ||
 
                               !SquareFreeQ[d], {d, di} += {-4, 1}]];
                   Continue[]]]; 
                Break[]
          ]; 
          If[messages,  Print[{di, d}]]; 
          hp = If[di > Length[HilbertPolynomialTable]
                       || di > classtablelen,
                  HilbertPolynomial[d], 
                  HilbertPolynomialTable[[di]]];
          SeedRandom[101];                   (* Same random factor  *)
          hr = RootMod[hp, p]; 
          ell = Mod[hr PowerMod[1728 - hr,-1,p], p]; 
          g2 = Mod[3 ell,p]; g3= Mod[2 ell,p]; 
          xec = 2; 
          While[True, 
                yec2 = Mod[PowerMod[xec,3,p] + g2 xec + g3, p];
                While[JacobiSymbol[yec2,p] == -1,
                      xec++;
                      yec2 = Mod[PowerMod[xec,3,p] + g2 xec + g3, p];
                      ];
                point = PointEC[xec, SqrtMod[yec2,p],g2,g3,p];  
                moverq = m/q; 

(* Point will be expressed in projective coordinates so as 
   to avoid doing GCD operations. The X coordinate will be 
   written as a ratio x:z. The identity element occurs when
   z = 0. See Bressoud's book, p. 213.
      One has to check that no nontrivial factors of p 
   occur, i.e., that 1 < GCD[z, p] < p does not occur. *)

                p1 = mult[{xec, 1}, moverq, g2, g3, p]; 
                If[p1[[2]] == 0, xec++; Continue[]]; 
                p2 = multcheckprime[p1, q, g2, g3, p]; 
                If[GCD[p2[[1]], p] > 1, 
                   Message[PrimeQCertificate::false];
                   Return[{p, GCD[p, p2[[1]]], g2, g3}]];

(* Check that second to last operation did not give a nontrivial
   factor of p. If it does then PrimeQ[p] is wrong. 
   p2[[1]] : p2[[2]] gives X coordinate of q * pt which should
   be the identity (p2[[2]] == 0), p2[[3]] : p2[[4]] gives the
   X coordinate of (q + 1) * pt. If a nontrivial divisor of
   p occured in the computation of q * pt, then it would remain
   in the the Z coordinate of (q+1) * pt (i.e., GCD[p2[[4]], p]>1).
*)

                If[p2[[3]] == 0, Break[]]; 
                c = 2; While[JacobiSymbol[c,p] == 1, c++];
                g2 = Mod[g2 c^2, p];
                g3 = Mod[g3 c^3, p]; 
                curve = 2;
                xec = 2
                ];  
      Prepend[PrimeQCertificate[q, opts], 
      {CertificatePrime -> p, 
       CertificatePoint -> point,  
       CertificateK -> m / q, 
       CertificateM -> m,  
       CertificateNextPrime -> q, 
       CertificateDiscriminant -> d}]
      ]]]
             (*  End of PrimeQCertificate[p]   *)





(*  *************************************************************
               Checking the certificate 
    ************************************************************

*)

(*  
    If p is composite check that PowerMod[a, p - 1, p] != 1 or
    that ( PowerMod[ a, 2, p] == 1 for a != -1, 1 ), for some a. 
*)

PrimeQCertificateCheck::badarg := "Certificate does not correspond
to number being tested."

PrimeQCertificateCheck[2, 2] := True

PrimeQCertificateCheck[_, 2] := False

PrimeQCertificateCheck[list_List, p_] :=
  Block[{am, pratt}, 
        If[Not[PrimeQ[p]], 

(* Make sure that certificate corresponds to number. *)

  If[!VectorQ[list] || p != list[[3]], 
     Message[PrimeQCertificateCheck::badarg]; Return[]]; 

Return[
         Length[list] == 3 && 

 (* Either a^(p-1) != 1 (mod p)  *)

         ((Rest[list] == {p-1, p} && Apply[PowerMod, list] != 1) ||

 (* or there is an element a != 1 or -1 such that a^2 = 1 (mod p). *)

         (Rest[list] == {2, p} && (First[list] != 1 && First[list] != p-1)&&
          Apply[PowerMod, list] == 1))]]; 
     
    (* Checking for primality: pratt is the Pratt portion of the test, i.e.,
       the last 3 entries in the certificate and am the rest of the 
       certificate, i.e., Atkin-Morain portion. *)

(* Check that certificate matches a certificate of primality, not
   compositeness. *)

If[VectorQ[list], 
   Message[PrimeQCertificateCheck::badarg]; Return[]]; 

      pratt = Take[list, -3];  am = Drop[list, -3]; 


(* Make sure that certificate corresponds to number *)

If[Not[p == First[pratt] || (p == (CertificatePrime /. First[am]))], 
   Message[PrimeQCertificateCheck::badarg]; Return[]]; 



(* Testing the Pratt certificate is equivalent to recomputing it.
   
   A nontrivial failed inversion during elliptic curve addition 
   means that a counterexample to PrimeQ has been found, i.e., 
   PrimeQ[p] is true, but a nontrivial factor of p is found. 
   This test is implemented using a Throw/Catch construction. *)

 Catch[answer = 
       (p == First[pratt] &&                  (* If only Pratt's test *)
        pratt == PrimeQCertificate[First[pratt], 
                                   SmallPrime-> p + 1]) ||
  (     (* Elliptic curve test *)
(((CertificateNextPrime /. Last[am]) == First[pratt])) &&
(pratt == PrimeQCertificate[First[pratt], SmallPrime->First[pratt+1]]) &&
   (p == First[pratt] || 
   (And @@ Map[PointECQ[(CertificatePoint /. #)] &&
              (!((CertificateK CertificatePoint /. #) === IdentityEC)) &&
              (CertificateM CertificatePoint /. #) === IdentityEC &&
              (CertificateM /. #) > (CertificateNextPrime /. #) > 
               ((N[CertificatePrime /. #]^(1/4) + 1)^2) &&
    (CertificateK /. #)  (CertificateNextPrime /. #) == (CertificateM /. #)&,
             am]))
   && (Table[CertificateNextPrime, {Length[am] -1}] /. Drop[am, -1]) ==
      (Table[CertificatePrime, {Length[am] -1}] /. Rest[am]))];
    If[Length[PrimeQCounterexample] == 3, 
       PrimeCounterexample, 
       answer]
   ]


   (* Computing the j invariant using the Dedekind eta function. 
      See Atkin--Morain, Section 3.5. *)

   
ModularEta[z_, prec_Integer, range_Integer]:= 
Block[{q = Exp[2 N[Pi, prec] I N[z, prec]], i}, 
       Exp[N[Pi, prec] I N[z, prec]/12] * 
       (1 + 
       Plus @@ Table[(-1)^i (q^(i (3 i -1)/2) + q^(i (3 i + 1)/2)), 
                     {i, range}])
      ]

ModularInvariantj[z_, prec_Integer, range_Integer] := 
Block[{f2 = 
       2^12 (ModularEta[2 z, prec, range] /
              ModularEta[z, prec, range])^24},
      (f2 + 16)^3/f2]   


(* Computing the Hilbert Polynomial. The specifications for
   the precision in computing the j-invariant and the number
   of terms to use in the series for eta(z) are given in 
   Atkin--Morain, Section 8.1. *)

HilbertPolynomial[d_Integer]:= 
      Block[{hprec, hrange, s, classlist = ClassList[d]},  
            hprec = Ceiling[10 + 2 Sqrt[N[-d]] * 
                            (Plus @@ N[1/(First /@ classlist)])]; 
            hrange = 3; 
      Expand[Fold[#1 x + #2&, 1, 
           Rest[Reverse[Round[N[CoefficientList[Expand[Apply[Times,
           x - Map[ModularInvariantj[
               (- #[[2]] + Sqrt[N[d, hprec]]) / (2 #[[1]]) , 
                    hprec, Round[hrange Sqrt[N[#[[1]]]]]]&, 
           classlist]]], {x}], hprec]]]]]]]
                        
(* 
   Find a root of the monic reducible polynomial f[x] (mod p). 
   This doesn't need Berlekamp's algorithm since Class Field Theory
   guarantees that the Hilbert Polynomial splits completely (mod p). 
   However this seems to be unnacceptable for larger degrees, e.g., 11.
*)

(* PolynomialMod does the wrong thing, so use PolyMod instead. *)

PolyMod[f_, p_]:= PolynomialMod[Expand[f], p]

RootMod[f_, p_Integer]:= Mod[-Coefficient[f,x,0], p]  /; Exponent[f,x] == 1

RootMod[f_,p_Integer]:= 
        Block[{cl,b,c},
              cl = Map[Mod[#, p] &, CoefficientList[f,x]];
              b = cl[[2]]; c = cl[[1]];
              Return[Mod[(-b + SqrtMod[b b - 4 c, p]) PowerMod[2, -1,p],p]]
             ]   /; Exponent[f,x] == 2

RootMod[f_, p_Integer] := 
  Block[{exp = Exponent[f, x], fp = PolyMod[f,p], s, t, i, q},
  While[True, 
        s = Plus @@ Table[Random[Integer,{0,p}] x^i, {i,0,exp-2}] + 
            x^(exp - 1); 
        t = GCDMod[s, fp, p]; 
       If[0 < Exponent[t,x] < exp, Break[],
         t =  GCDMod[fp, PolyPowerMod[s, (p-1)/2, {fp, p}] - 1, p];
            If[ 0 < Exponent[t,x] < exp, Break[]]
        ]]; 
    q = PolyMod[PolynomialQuotient[fp, t, x], p];
    t = If[Exponent[q, x] < Exponent[t, x], q, t];
    RootMod[t, p]
   ]       /; Exponent[f,x] > 2


(* You need a different GCD for polynomials (mod p). E.g., 
   GCD[x^3 - 1, x-2] = 1, but PolynomialRemainder[x^3 -1, x - 2, x] = 7.
   In other words, using the built-in GCD would not allow you to find
   the root 2 of x^3 - 1 == 0 (mod 7).      *)
   
GCDMod[f_,g_,p_Integer] := 
    GCDMod[g,f,p] /; Exponent[f,x] < Exponent[g,x] 

GCDMod[f_, g_, p_Integer] := f  /;  Exponent[g,x] == -Infinity

GCDMod[f_, g_, p_Integer] := 
          Block[{fp = PolyMod[f, p], gp = PolyMod[g, p], 
                 q,r,monic},
                     monic = Coefficient[gp, x, Exponent[gp,x]];
                     If[monic != 1, 
                        gp = PolyMod[PowerMod[monic, -1, p] gp, p]];
                     q = PolyMod[PolynomialQuotient[fp, gp, x], p];
                     r = PolyMod[PolynomialRemainder[fp, gp, x], p];
                     If[Exponent[r, x] <= 0,
                        If[r == 0, Return[gp], Return[1]],
                     Return[GCDMod[gp,r,p]]
                     ]
             ]  /; Exponent[g,x] >= 0


PolynomialRemainderMod[f_, g_, x_, p_Integer] := 
    PolyMod[PolynomialRemainder[PolyMod[f, p], 
                  PolyMod[g, p], x], p];
   
PolyPowerMod[f_,n_Integer,{g_,p_Integer}] :=
Fold[PolynomialRemainderMod[#1^2 #2, g, x, p] &,
     1, 
     Block[{prm = PolynomialRemainderMod[f, g, x, p]},
            If[# == 0, 1, prm]& /@ IntegerDigits[n, 2]]] /; n > 0



(* 
   This tries to find a nontrivial large prime factor of n.  
   (This part of the program is the bottleneck since large numbers require 
   nontrivial factoring routines). The problem in this implementation
   is the fact that trial division is so slow. This part of the 
   program should be speeded up a lot when the complete version of 
   FactorInteger is implemented in the kernel.

   The larger the number, the simpler the factoring routine used. 
   For very large numbers, one uses only trial division. 
*)
 

LargeFactor[qq_, m_Integer,n_Integer, pptest_, prtest_] :=
         Block[{q = qq, pp, pr, nn = (N[n]^(1/4)+1)^2}, 
                If[PrimeQ[q] &&  nn < q < m, Return[q]];
                If[pptest, 
                 If[messages, Print["trying PollardP"]];   
                   pp = GCD[q, PowerMod[2, LCM @@ Range[3000], q]-1];
                   q = Quotient[q, pp];
                   If[PrimeQ[q] &&  nn < q < m, 
                      Return[q]]]; 
                 If[prtest, 
                 If[messages,  Print["trying PollardRho"]];  
                   If[!(N[10^5 nn] < q < m), Return[]];
                 pr = PollardRhoPrimeQ[q, 2, Min[2000, Round[N[n]^(1/4)]]];
                   If[pr != q, q = Quotient[q, pr]];
                   If[PrimeQ[q] &&  nn < q < m, 
                      Return[q]]];
              ]


(* Finding a solution of  4p = x^2 + |d| y^2  if it exists.  *)


QuadraticRepresentation4[d_Integer,p_Integer] :=
      Block[{b,prec,u,v,w,x,y,lr, inv},
            b = SqrtMod[d,p];
            If[EvenQ[b], b = p - b];
            prec = Max[Round[N[4/3 Log[10,p]]], 50];   
            v = HighPoint[ N[ (-b +  Sqrt[d])/(2 p), prec]];
            x= v[[2]]; y = v[[1]];
            If[(x p - b y/2)^2 - y^2 d/4 == p, 
               {x p - b y/2, y/2}, 
               False]
           ]    /;     d < 0 && Mod[d,4] == 1   

HighPoint[z_Complex] :=  
   Block[{g = IdentityMatrix[2], abs, t=z, round}, 
         While[True,
               round = Round[Re[t]]; 
               If[round != 0,
                      t = t - round;  
                      g = {{1,-round},{0,1}}. g,
                      abs = 1/(Re[t]^2 + Im[t]^2);
                  If[abs > 1, 
                     t =  - Conjugate[t] abs; 
                     g = {{0,-1},{1,0}} . g,
                     Break[]
                    ]
                  ]
               ];
          g[[2]]
        ]   

(* This version of PollardRho doesn't do intermediate GCD's,,
   since you want to find all the small factors of a large 
   number. *)

PollardRhoPrimeQ[n_, c_, max_]:= 
  Block[{x1 = 2, x2 = 2^2 + c, i = 0, j, gcd, 
         range = 1, prod = 1, factor = 1},
         While[i <= max, 
               Do[x2 = Mod[x2^2 + c, n];
                  prod = Mod[prod (x1 - x2), n], {range}];
               i += range; 
               x1 = x2; range *= 2;
               x2 = Nest[Mod[#^2 + c, n]&, x2, range]];
           factor = GCD[prod, n]; 
        If[factor > 1, factor, n]]


(* ClassTable stores the fundamental discriminants d > -24000 with
   class number 20 or less. They are ordered by class number and by 
   size. Since d is negative of the form 4k + 1, d is stored as 
   (1 - d)/4 in order to save space. This large list of 1186 numbers
   is required so that weaker factoring algorithms can be used on 
   large numbers (so that the algorithms searches more numbers). 
   If this is too large, you can reduce it by taking the first 489
   elements (class number <= 12). 
*)

ClassTable = 
  {2, 3, 5, 11, 17, 41, 4, 9, 13, 23, 29, 31, 47, 59, 67, 101, 107, 6, 8, 15, 
   21, 27, 35, 53, 71, 77, 83, 95, 125, 137, 161, 221, 227, 10, 14, 39, 49, 
   51, 55, 65, 73, 81, 89, 109, 121, 139, 149, 157, 167, 179, 181, 191, 199, 
   239, 251, 257, 307, 311, 347, 353, 359, 377, 389, 12, 20, 26, 32, 33, 45, 
   57, 87, 111, 131, 143, 155, 171, 173, 185, 197, 237, 263, 281, 431, 437, 
   467, 551, 587, 671, 22, 62, 85, 103, 113, 129, 177, 193, 209, 211, 265, 
   275, 287, 301, 305, 317, 329, 337, 341, 391, 401, 461, 479, 491, 557, 571, 
   611, 629, 641, 697, 731, 809, 857, 881, 941, 18, 38, 56, 63, 116, 117, 
   122, 147, 203, 207, 215, 291, 293, 371, 381, 407, 447, 497, 503, 521, 545, 
   563, 617, 677, 755, 767, 797, 977, 1151, 1277, 1481, 24, 28, 46, 74, 75, 
   93, 99, 145, 146, 163, 229, 235, 245, 247, 249, 261, 283, 289, 299, 325, 
   335, 361, 409, 413, 415, 433, 443, 449, 451, 485, 487, 499, 509, 515, 517, 
   535, 541, 577, 581, 599, 605, 613, 647, 653, 667, 679, 689, 707, 737, 749, 
   751, 811, 829, 839, 851, 877, 899, 947, 971, 991, 1031, 1049, 1067, 1081, 
   1097, 1187, 1211, 1217, 1271, 1367, 1397, 1427, 1487, 1577, 50, 92, 105, 
   123, 141, 206, 272, 297, 323, 356, 395, 501, 701, 791, 815, 827, 887, 911, 
   1007, 1061, 1091, 1121, 1181, 1247, 1361, 1511, 1607, 1691, 1721, 1931, 
   2141, 2201, 2267, 2657, 30, 36, 40, 76, 80, 104, 153, 159, 175, 195, 201, 
   213, 231, 279, 411, 427, 445, 455, 459, 473, 481, 591, 661, 725, 773, 785, 
   787, 823, 909, 917, 921, 953, 965, 1021, 1057, 1109, 1145, 1157, 1229, 
   1283, 1291, 1379, 1403, 1417, 1451, 1529, 1565, 1601, 1667, 1781, 1841, 
   1847, 1859, 1871, 1907, 2057, 2237, 2327, 2537, 2621, 3461, 42, 68, 165, 
   242, 321, 326, 327, 365, 383, 425, 507, 567, 635, 683, 713, 743, 801, 837, 
   875, 935, 983, 1013, 1295, 1421, 1541, 1637, 1757, 1877, 1901, 1967, 2111, 
   2321, 2351, 2411, 2447, 2747, 3251, 3317, 3527, 3671, 3917, 58, 64, 82, 
   136, 164, 172, 183, 189, 253, 267, 309, 314, 339, 343, 355, 373, 379, 387, 
   441, 452, 489, 523, 539, 559, 589, 597, 623, 639, 649, 657, 659, 681, 699, 
   717, 739, 757, 761, 779, 847, 863, 865, 905, 931, 937, 959, 989, 1009, 
   1037, 1039, 1073, 1103, 1117, 1133, 1139, 1147, 1175, 1189, 1193, 1199, 
   1207, 1227, 1237, 1241, 1259, 1289, 1315, 1325, 1327, 1343, 1349, 1381, 
   1399, 1439, 1441, 1453, 1459, 1547, 1559, 1567, 1571, 1621, 1651, 1661, 
   1679, 1697, 1711, 1733, 1739, 1741, 1747, 1777, 1823, 1889, 1921, 1973, 
   2033, 2039, 2081, 2087, 2099, 2197, 2207, 2251, 2285, 2339, 2381, 2417, 
   2461, 2501, 2651, 2677, 2687, 2699, 2729, 2789, 2837, 2927, 2951, 3077, 
   3161, 3611, 3791, 3821, 4001, 4451, (* last classno = 12 *)
   48, 66, 152, 158, 182, 255, 363, 375, 
   417, 477, 533, 536, 593, 665, 741, 771, 923, 1001, 1127, 1161, 1337, 1355, 
   1445, 1655, 1811, 1991, 2387, 2435, 2867, 2897, 2957, 2981, 3011, 3587, 
   3947, 4241, 5141, 54, 72, 98, 112, 128, 134, 176, 202, 225, 303, 351, 382, 
   463, 471, 531, 537, 543, 584, 607, 627, 643, 695, 733, 807, 895, 927, 929, 
   967, 1047, 1079, 1111, 1165, 1201, 1257, 1273, 1313, 1317, 1431, 1493, 
   1625, 1631, 1745, 1767, 1775, 1787, 1979, 2009, 2047, 2153, 2225, 2279, 
   2309, 2357, 2531, 2579, 2591, 2603, 2807, 3037, 3167, 3197, 3257, 3359, 
   3371, 3401, 3551, 4217, 4547, 4637, 4661, 5057, 5387, 5771, 60, 110, 188, 
   243, 315, 332, 357, 392, 405, 561, 662, 675, 711, 833, 893, 951, 1025, 
   1055, 1251, 1307, 1331, 1391, 1457, 1497, 1517, 1523, 1553, 1643, 1805, 
   1865, 1887, 2117, 2177, 2195, 2261, 2477, 2561, 2567, 2615, 2663, 2681, 
   2771, 2993, 3041, 3191, 3287, 3491, 3581, 3707, 3713, 3797, 3911, 3977, 
   4151, 4211, 4367, 4481, 4511, 4631, 4847, 4967, 5177, 5501, 100, 102, 118, 
   140, 166, 200, 224, 226, 236, 254, 256, 262, 285, 290, 345, 399, 496, 505, 
   524, 549, 553, 595, 685, 703, 766, 793, 841, 845, 859, 861, 883, 901, 949, 
   955, 985, 1011, 1043, 1045, 1063, 1095, 1107, 1129, 1171, 1172, 1225, 
   1235, 1243, 1279, 1297, 1299, 1341, 1351, 1471, 1477, 1499, 1525, 1531, 
   1549, 1579, 1589, 1599, 1649, 1687, 1693, 1731, 1751, 1759, 1763, 1799, 
   1829, 1837, 1849, 1873, 1895, 1927, 1937, 1939, 1949, 1955, 1961, 1981, 
   1999, 2011, 2021, 2071, 2075, 2129, 2137, 2159, 2161, 2171, 2179, 2209, 
   2215, 2305, 2371, 2377, 2399, 2441, 2459, 2467, 2489, 2549, 2551, 2557, 
   2597, 2641, 2647, 2659, 2701, 2711, 2741, 2767, 2777, 2795, 2801, 2881, 
   2891, 2909, 2929, 2999, 3007, 3065, 3097, 3131, 3149, 3187, 3209, 3215, 
   3281, 3299, 3331, 3341, 3377, 3449, 3455, 3457, 3539, 3593, 3601, 3637, 
   3677, 3691, 3749, 3767, 3847, 3851, 3887, 3929, 4007, 4049, 4087, 4133, 
   4139, 4181, 4307, 4331, 4337, 4357, 4379, 4601, 4679, 4721, 4727, 4787, 
   4799, 4987, 4997, 5039, 5099, 5351, 5429, 5459, 5561, 5711, 5849, 5897, 
   96, 248, 273, 393, 416, 446, 633, 831, 987, 1085, 1112, 1137, 1163, 1371, 
   1551, 1595, 1613, 1707, 1727, 1971, 2135, 2183, 2471, 2813, 2861, 3227, 
   3407, 3521, 3695, 3737, 4175, 4457, 4577, 4991, 5267, 5891, 84, 130, 132, 
   170, 284, 302, 346, 422, 423, 482, 512, 513, 542, 573, 579, 687, 715, 759, 
   777, 886, 913, 1075, 1077, 1205, 1221, 1329, 1385, 1475, 1557, 1583, 1597, 
   1685, 1709, 1831, 1835, 1893, 1929, 1943, 1957, 2051, 2083, 2101, 2127, 
   2147, 2281, 2303, 2391, 2407, 2421, 2587, 2631, 2855, 2911, 2921, 2963, 
   3017, 3047, 3059, 3071, 3163, 3181, 3203, 3307, 3329, 3347, 3437, 3487, 
   3497, 3541, 3557, 3629, 3667, 3779, 3811, 4031, 4043, 4097, 4157, 4259, 
   4283, 4351, 4409, 4571, 4757, 4781, 4913, 5009, 5207, 5261, 5417, 5477, 
   5567, 5611, 5627, 5737, 5837, 5867, 5921, 5981, 78, 90, 230, 266, 386, 
   458, 525, 585, 615, 836, 866, 867, 902, 1005, 1035, 1082, 1265, 1287, 
   1382, 1415, 1701, 2105, 2231, 2243, 2405, 2723, 2825, 3773, 3833, 4091, 
   4187, 4253, 4325, 4385, 4421, 4877, 5297, 5303, 5321, 5801, 114, 154, 316, 
   374, 424, 435, 472, 562, 602, 621, 622, 651, 686, 746, 747, 752, 775, 776, 
   783, 789, 805, 842, 849, 879, 892, 914, 939, 995, 1003, 1004, 1089, 1093, 
   1099, 1101, 1135, 1183, 1191, 1214, 1223, 1255, 1333, 1433, 1489, 1513, 
   1535, 1585, 1617, 1689, 1783, 1791, 1793, 1867, 1903, 1909, 1913, 1917, 
   1963, 1985, 2017, 2063, 2089, 2123, 2189, 2249, 2263, 2287, 2299, 2333, 
   2335, 2341, 2361, 2393, 2423, 2439, 2449, 2495, 2507, 2521, 2539, 2543, 
   2573, 2575, 2627, 2629, 2705, 2839, 2845, 2849, 2857, 2885, 2915, 2939, 
   2971, 2987, 2989, 3005, 3035, 3079, 3083, 3089, 3091, 3117, 3125, 3147, 
   3151, 3221, 3233, 3239, 3241, 3289, 3311, 3389, 3413, 3451, 3577, 3583, 
   3617, 3623, 3665, 3689, 3809, 3839, 3901, 3923, 3941, 3971, 3989, 4037, 
   4099, 4109, 4121, 4127, 4171, 4177, 4229, 4231, 4267, 4297, 4391, 4411, 
   4441, 4477, 4517, 4541, 4549, 4589, 4591, 4771, 4861, 4889, 4981, 5021, 
   5051, 5147, 5171, 5189, 5221, 5273, 5309, 5327, 5347, 5399, 5431, 5441, 
   5471, 5597, 5617, 5639, 5651, 5681, 5861, 5987}


(* Table of Hilbert polynomials corresponding to fundamental 
   discriminants of class number <= 3. *)

HilbertPolynomialTable = 
  {3375 + x, 32768 + x, 884736 + x, 884736000 + x, 147197952000 + x, 
   262537412640768000 + x, -121287375 + 191025*x + x^2, 
   -134217728000 + 117964800*x + x^2, 6262062317568 + 5541101568*x + x^2, 
   -3845689020776448 + 10359073013760*x + x^2, 
   130231327260672000 + 427864611225600*x + x^2, 
   148809594175488000000 + 1354146840576000*x + x^2, 
   -3845689020776448000000 + 4545336381788160000*x + x^2, 
   11946621170462723407872000 + 823177419449425920000*x + x^2, 
   531429662672621376897024000000 + 19683091854079488000000*x + x^2, 
   -108844203402491055833088000000 + 2452811389229331391979520000*x + x^2, 
   155041756222618916546936832000000 + 15611455512523783919812608000*x + x^2, 
   12771880859375 - 5151296875*x + 3491750*x^2 + x^3, 
   1566028350940383 - 58682638134*x + 39491307*x^2 + x^3, 
   374643194001883136 - 140811576541184*x + 30197678080*x^2 + x^3, 
   549755813888000000000 - 41490055168000000*x + 2691907584000*x^2 + x^3, 
   337618789203968000000000 - 6764523159552000000*x + 129783279616000*x^2 + 
    x^3, 67408489017571610198016 - 53041786755137667072*x + 
    12183160834031616*x^2 + x^3, 
   5310823021408898698117644288 + 277390576406111100862464*x + 
    65873587288630099968*x^2 + x^3, 
   201371843156955365376000000000 + 90839236535446929408000000*x + 
    89611323386832801792000*x^2 + x^3, 
   8987619631060626702336000000000 - 5083646425734146162688000000*x + 
    805016812009981390848000*x^2 + x^3, 
   56176242840389398230218488594563072 + 368729929041040103875232661504*x + 
    6647404730173793386463232*x^2 + x^3, 
   15443600047689011948024601807415148544 - 
    121567791009880876719538528321536*x + 364395404104624239018246144*x^2 + 
    x^3, 4671133182399954782798673154437441310949376 - 
    6063717825494266394722392560011051008*x + 
    3005101108071026200706725969920*x^2 + x^3, 
   83303937570678403968635240448000000000 - 
    139712328431787827943469744128000000*x + 
    81297395539631654721637478400000*x^2 + x^3, 
   308052554652302847380880841299197952000000000 - 
    6300378505047247876499651797450752000000*x + 
    39545575162726134099492467011584000*x^2 + x^3, 
   167990285381627318187575520800123387904000000000 - 
    151960111125245282033875619529124478976000000*x + 
    34903934341011819039224295011933392896000*x^2 + x^3, 
   149161274746524841328545894969274007552000000000 + 
    39181594208014819617565811575376314368000000*x + 
    123072080721198402394477590506838687744000*x^2 + x^3}




(* Elliptic Curves *)


PointEC /: PointECQ[PointEC[pt___]] := PointECQ @@ PointEC[pt]


PointECQ[x_,y_,g2_,g3_,p_] := Mod[y^2 - x^3 - g2 x - g3, p] === 0

PointEC /:  PointEC[x1_,y1_,g2_,g3_,p_] == PointEC[x2_,y2_,g2_,g3_,p_]:=
              (Mod[x1,p] == Mod[x2,p]) && (Mod[y1,p] == Mod[y2,p]) 


(* 
    Addition on the Elliptic curve y^2 = x^3 + g2 x + g3 (mod p).
    This is the part of the algorithm that gets used the most.
*)

PointEC /:  PointEC[x_Integer,y_Integer,g2_,g3_,p_] + IdentityEC := 
                      PointEC[x,y,g2,g3,p]

PointEC /: PointEC[x1_,y1_,g2_,g3_,p_] + PointEC[x2_,y2_,g2_,g3_,p_] :=
Block[{slope, IdentityECQ = False, gcd},
       If[Mod[x1 - x2, p] == 0,
          If[Mod[y1 + y2, p] == 0,
             IdentityECQ = True,
       slope = Mod[(3 x1^2 + g2) PowerMod[2 y1, -1, p], p]],
       slope = Mod[(y2 - y1) PowerMod[x2 - x1, -1, p], p],
              ]; 
 (* A failed inversion implies that a counterexample to PrimeQ. *) 

       If[!IntegerQ[slope], 
          gcd = GCD[y1, p]; 
          If[1 < gcd < p,   
             Message[PrimeQCertificate::false];
             PrimeQCounterexample = {False, p, gcd, False}; 
             Throw[True]];
          gcd = GCD[x2 - x1, p]; 
          If[1 < gcd < p,   
             Message[PrimeQCertificate::false];
             PrimeQCounterexample = {False, p, gcd, False}; 
             Throw[True]]
           ];

(*   Otherwise continue with addition. *)

             If[SameQ[IdentityECQ, True], 
                     IdentityEC,
                     x3 =  Mod[slope^2 - x1 - x2, p];
                     y3 =  Mod[slope (x1-x3) - y1, p];
                     PointEC[x3,y3,g2,g3,p]
                 ]
            ]
             
 (* Multiply a point on the curve by an integer by repeated doubling, 
    exactly like PowerMod[a,b,c].  *)

IdentityEC /:  n_Integer * IdentityEC := IdentityEC

PointEC /:  0 * PointEC[x_,y_,g2_,g3_,p_]:= IdentityEC

PointEC /:  1 * PointEC[x_,y_,g2_,g3_,p_]:= PointEC[x,y,g2,g3,p]

PointEC /: 2 * PointEC[x_, y_, g2_, g3_, p_]:= 
                PointEC[x, y, g2, g3, p] + PointEC[x, y, g2, g3, p]

PointEC /: n_ * PointEC[x_, y_, g2_, g3_, p_]:= 
                 (-n) * PointEC[x, -y, g2, g3, p] /; n < 0

PointEC /: n_ * PointEC[x_, y_, g2_, g3_, p_]:= 
  Fold[2 #1 + #2&, IdentityEC, 
       IntegerDigits[n, 2] PointEC[x, y, g2, g3, p]]  /; 2 < n  <= 2^100 

(* For large n use base 16 to speed up the calculation by 25-> 50%.
   If used repeatedly, precompute digits in main algorithm. *)


PointEC /: n_ * PointEC[x_, y_, g2_, g3_, p_]:= 
  Fold[16 #1 + #2&, IdentityEC, 
       NestList[# + PointEC[x, y, g2, g3, p]&, 
                IdentityEC, 15] [[IntegerDigits[n, 16] + 1]]] /; n > 2^100



(* Elliptic curve multiplication in homogeneous form (no GCD's),
   as in Bressoud's book, p. 213. *)

(* X coordinate of 2*P. *)

double[{x_, z_, g2_, g3_, m_}]:= 
Block[{x2 = Mod[x^2, m], z2 = Mod[z^2, m], bz3},
       bz3 = g3 Mod[z2 z, m]; 
       {Mod[(x2 - Mod[g2 z2, m])^2 - 8 x bz3, m],
        Mod[4 z Mod[Mod[x2 x, m] + g2 Mod[x z2, m] + bz3, m], m]}
      ]    


(* X coordinate of (2i+1)*P from i*P and (i+1)*P. *)

doubleplusone[{x1_, z1_, xi_, zi_, xi1_, zi1_, g2_, g3_, m_}]:= 
Block[{zizi1 = Mod[zi zi1, m], xi1zi = Mod[xi1 zi, m],
       xizi1 = Mod[xi zi1, m]},
      {Mod[z1 Mod[PowerMod[xi xi1 - g2 zizi1, 2, m] -
                  4 g3 Mod[zizi1 (xizi1 + xi1zi), m], m], m],
       Mod[x1 PowerMod[xi1zi - xizi1, 2, m], m]}
     ] 
   

(* X coordinate of {2i*P, (2i+1)*P} from {i*P, (i+1)*P}. *)

even[{x1_, z1_, xn_, zn_, xm_, zm_, g2_, g3_, m_}]:= 
Join[{x1, z1}, double[{xn, zn, g2, g3, m}], 
     doubleplusone[{x1, z1, xn, zn, xm, zm, g2, g3, m}], {g2, g3, m}]

(* X coordinate of {(2i+1)*P, (2i+2)*P} from {i*P, (i+1)*P} *)

odd[{x1_, z1_, xn_, zn_, xm_, zm_, g2_, g3_, m_}]:= 
Join[{x1, z1}, doubleplusone[{x1, z1, xn, zn, xm, zm, g2, g3, m}], 
     double[{xm, zm, g2, g3,  m}], {g2, g3, m}]

(* Compute n*P, where P has X coordinate x1:z1. *)

mult[{x1_, z1_}, n_, g2_, g3_,  m_]:= 
     Fold[If[#2 == 0, even[#1], odd[#1]]&,
          Join[{x1, z1}, {x1,z1},double[{x1, z1, g2, g3, m}], {g2, g3, m}],
               Rest[IntegerDigits[n, 2]]] [[{3,4}]]

(* This is the same as mult but multiplies together all the 
   intermediate zi's in order to see whether a spurious 
   factor of p exists. *) 

multcheckprime[{x1_, z1_}, n_, g2_, g3_,  m_]:= 
  Fold[Join[{Mod[#1[[1]] #1[[5]] #1[[7]], m]}, 
             If[#2 == 0, even[Rest[#1]], odd[Rest[#1]]]]&,
       Join[{1}, {x1, z1}, {x1,z1}, double[{x1, z1, g2, g3, m}], {g2, g3, m}],
            Rest[IntegerDigits[n, 2]]] [[{1, 4, 5}]]

Protect[ProvablePrimeQ, PrimeQCertificateCheck, (* PrimeQCertificate, *)
        SmallPrime, Certificate, TrialDivisionLimit, PollardPTest,
        PollardRhoTest, ModularInvariantj, HilbertPolynomial, fact,
        CertificatePrime, CertificatePoint, CertificateK, CertificateM,
        CertificateNextPrime, CertificateDiscriminant, PointEC, PointECQ]


End[]


EndPackage[]
