/* Written  4:11 pm  Sep  7, 1989 by ilan@Gang-of-Four.Stanford.EDU in uxe.cso.uiuc.edu:sci.math.symbolic */
/* ---------- "Deterministic version of PrimeQ" ---------- */
I have just finished implementing the Atkin-Morain deterministic
primality test. I ended up including a number of other functions
relating to prime numbers. Enjoy!

-Ilan Vardi

---------------------------------------------------------------------------



(*          Program by Ilan Vardi, September 4, 1989. 

            Copywright Wolfram Research, Inc.         
            P.O. Box 6059, Champaign, IL 61820, 217-398-0700   *)


BeginPackage["PrimeQNew`"] 

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

Options[PrimeQNew] := {Probabilistic -> False,
                       Certificate -> False,
                       SmallPrime -> 10^8
                      }

PrimeQCertificate::usage = "PrimeQCertificate[n] prints out a certificate
that n is prime using Atkin-Morain if n > SmallPrime and Pratt's
certificate if n < SmallPrime, or that n is composite."


Options[PrimeQCertificate] := {SmallPrime -> 10^8}  

PrimeQCertificateCheck::usage = "PrimeQCertificateCheck[n] checks that
PrimeQCertificate[n] is correct."

Options[PrimeQCertificateCheck] := {SmallPrime -> 10^8}  

PointEC::usage = "PointEC[x,y,g2,g3,p] is a point on the elliptic curve
y^2 = x^3 + g2 x + g3  (mod p). You can add two points or multiply a point 
by an integer."

PointECQ::usage ="PointEC[x,y,g2,g3,p] is True if y^2 = x^3 + g2 x + g3 
(mod p), False otherwise."

ClassNumber::usage = "ClassNumber[d] gives the number of inequivalent
quadratic forms of discriminant d, where d < 0 is a square free integer
of the form 4 n + 1."

ClassList::usage = "ClassList[d] gives a list of inequivalent quadratic
forms of discriminant d, where d < 0 is a square free integer of the
form 4 n + 1, and  a quadratic form a x^2 + b x y + c y^2 is represented
as { a, b, c}."

SquareFreeQ::usage = "SquareFreeQ[n] is True if n is a square free
integer, False otherwise."

SqrtMod::usage = "SqrtMod[ a, p] is the square root of a (mod p), where
p is a prime number and JacobiSymbol[ a, p] = 1."

PrimeCount::usage ="PrimeCount[x] gives the number of prime numbers
less than or equal to x."

PrimitiveRoot::usage ="PrimitiveRoot[p] gives the smallest generator
a of the multiplicative group (mod p), where p is a prime number."


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])."


PrimeQNew[p_Integer, opts___] :=
          If[(Probabilistic /. {opts} /. Options[PrimeQNew]),
             PrimeQ[p], 

(*  
    PrimeQNew[p] has to wait until PrimeQCertificate has finished to
    return its answer. Therefore a certificate has been generated.
    Since PrimeQCertificate caches its value, run PrimeQCertificate[p]
    to get the certificate. 

    PrimeQNew uses the old (probabilistic) 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 PrimeQNew will not 
    return a value (will go in an infinite loop). Therefore, PrimeQNew cannot
    be used to find counterexamples to the old PrimeQ.
*)

            If[ Not[(Certificate /. {opts} /. Options[PrimeQNew])],
                Block[{},
                      PrimeQCertificate[ p, opts];
                      Return[PrimeQ[p]]
                      ],  

(*  Or else get a certificate along with the answer   *)

       {PrimeQ[p], PrimeQCertificate[ p, opts]}
               ]
              ]

PrimeQCertificate[ 2, opts___] := 2

PrimeQCertificate[p_Integer, opts___] := PrimeQCertificate[p, opts] =
    If[Not[PrimeQ[p]], 
(*
     Certificate of compositeness. Apply[PowerMod, #] to 
     PrimeQCertificate[p] to prove that p is composite.

     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. 
*)
       Block[{ a = 1, q = (p - 1)/2,  r, s = -1, exp = 2},
              While[ EvenQ[q], q /= 2]; 
              While[ Mod[ s + 1, n] == 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; 
v                                   r = PowerMod[ s, 2, p] 
                                 ]
                              ]
                       ]
                     ]  ;
              Return[ { s, exp, p}]
          ],              
(*  
      Certificate of Primality  

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

      Pratt's certificate gives a proof that a is a generator 
      (primitive root) of the multiplicative group (mod p) which
      proves that p is a prime.  You need to factor p-1, which is 
      feasible for small p. 
     
      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.
     
   *)
      If[PrimeQ[p] && p  > 2 &&
         p <= (SmallPrime /. {opts} /. Options[PrimeQCertificate]), 
               Block[{ a = 2, flist}, 
                     flist = First[Transpose[FactorInteger[ p  - 1]]];
                     While[ PowerMod[ a, p - 1, p] != 1 ||  
                            Select[ Map[ PowerMod[ a, (p-1)/#, p] &, flist], 
                                    # == 1 &
                                  ] != {}, 
                            a++
                           ];
                     Return[{p, a, 
                             Map[ PrimeQCertificate[ #, opts] &,
                                 flist 
                                ]
                            }
                           ]
                   ], 
(*   
      For large primes p > SmallPrime use the Atkin-Morain test:

      The Atkin-Morain deterministic primality test. At present my version
      runs only up to sixty digit primes. The main problem is that polynomial 
      arithmetic (mod p) had to be written in top level for large primes
      and so finding roots of polynomials (mod p) is too slow above that.
      
      This can be avoided by precomputing Hilbert Polynomials
      for d > -3000 computing the roots in terms of radicals over the complex
      numbers and storing them in a table. Expressing these radicals (mod p) 
      will be very fast.

      The other slow part of the algorithm is addition on an elliptic
      curve, but this can't be speeded up significantly.

      The Atkin-Morain test returns a recursive certificate for the 
      probable 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).

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

           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 Goldwasser--Kilian (see Lenstra, Elliptic Curve algorithms, 
       World Congress 1986).

 (iii) Each q has its recursive certificate following it. So if 
       the first 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]] 

NextPrime == CertificateM /CertificateK && NextPrime > 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 F. Morain, INRIA technical report #911.

A heuristic running time of this algorithm is (Log[n])^6, see
Lenstra & Lenstra, ``Algorithms in Number Theory'', Handbook
of Theoretical Computer Science, North Holland, to appear.

*) 
 Block[{d,qr,t,qplus,qminus,m,q,moverq,hp,hr,ell,
        g2,g3,xec,yec2,p1,p2,point,c,primetime = 50},
        d = -7;
        While[True,
              While[JacobiSymbol[d,p] != 1 || Not[SquareFreeQ[d]], d -= 4];
                    qr= QuadraticRepresentation[d,p];
                    If[qr == False, d -= 4; Continue[]];
                    t = Abs[2 qr[[1]]]; s = Abs[2 qr[[2]]];
                    gcdp = GCD[t + 2, s]; gcdm = GCD[t - 2, s];
                    qplus =  LargeFactor[p + 1 + t, p ];
                    If[IntegerQ[qplus], m = p + 1 + t; q = qplus,
                       qminus = LargeFactor[p + 1 - t,p];
                       If[IntegerQ[qminus], m = p + 1 - t; q = qminus,
                          d -= 4; Continue[]
                         ]
                      ];
                 Break[]
               ]; 
           hr = RootMod[HilbertPolynomial[d], p];
         ell = Mod[ hr PowerMod[1728 - hr,-1,p], p]; 
         g2 = Mod[3 ell,p]; g3= Mod[2 ell,p];
          SeedRandom[101];                   (* Same random factor  *)
          While[True,
                xec = Random[Integer,{1,p}]; 
                yec2 = Mod[ PowerMod[xec,3,p] + Mod[g2 xec,p] + g3,p];
                If[JacobiSymbol[yec2,p] == -1, 
                   Continue[],
                   point = PointEC[xec, SqrtMod[yec2,p],g2,g3,p];  
                   moverq = m/q;
                   p1 = moverq point;  
                   If[SameQ[p1, IdentityEC], 
                      Continue[], 
                      p2 = q p1; 
                      If[SameQ[p2, IdentityEC],
                         Break[],
    (* 
       Apparently the choice of c below can be made deterministic
       for class number 1  
    *)
                         c = Random[Integer,{1,p}]; 
                         While[JacobiSymbol[c,p] == 1,
                               c = Random[Integer,{1,p}]];
                         g2 = Mod[g2 c^2, p]; g3 = Mod[g3 c^3, p];
                         Continue[]
                        ]
                      ]
                 ]
             ]; 
       Return[Prepend[PrimeQCertificate[q, opts], 
      {CertificatePrime -> p, 
       CertificatePoint -> point,  
       CertificateK -> m / q, 
       CertificateM -> m,  
       NextPrime -> q, 
       Discriminant -> 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[ p_Integer, opts___] :=
      If[ Not[PrimeQ[p]], 
          Apply[ PowerMod, PrimeQCertificate[ p, opts]], 
     
(*   
     Checking the Pratt certificate amounts to recomputing it. So
     this is how it's checked here.
*)
    If[ p <= (SmallPrime /. {opts} /. Options[PrimeQCertificate]),
        PrimeQCertificate[ p, opts], 

(*  
    Checking the Atkin-Morain certificate, see above.  
    Note that the certificate is short, but checking it will be slow
    (up to 60% of original running time). 
*)

     Block[{pqc = PrimeQCertificate[ p, opts]},   
     {Map[ { CertificatePrime /. #,
            CertificateK CertificatePoint /. #, 
            CertificateM CertificatePoint /. #,
            PointECQ[(CertificatePoint /. #)] &&
            (NextPrime /. #) > 
            N[( (CertificatePrime /. #)^(1/4) + 1)^2] &&
            (NextPrime /. #) == (CertificateM /CertificateK /. #),
            NextPrime /. #
           } &, 
          Drop[ pqc,  {-3,-1}]
         ],
       PrimeQCertificate[ pqc[[-3]], opts]
      } ] 
     ]]

ModularInvariantj[z_, prec_Integer, range_Integer] := 
  Block[{q = N[Exp[ 2 Pi I z], prec],e4,eta,i},
        e4 = 1 + 
             240 Table[DivisorSigma[3,i], {i,range}] . 
                 NestList[ q # &, q, range - 1];
        eta = 1 + Apply[Plus,
                        Table[ (-1)^i q^(i (3 i - 1)/ 2) (1 + q^i), 
                        {i, Floor[N[(1 + Sqrt[1 + 24 range]) / 6]]}
                             ]
                       ];
        Return[(e4/eta^8)^3 / q]
       ]
   
HilbertPolynomial[d_Integer] := HilbertPolynomial[d] =
      Block[{hprec = 300, hrange = 50, j},
    (*  hprec = Round[10 N[Log[10,d]]]; hrange = Round[5 N[Log[10,d]]];*)
      Round[
            N[
              CoefficientList[
                          Expand[Apply[Times,
                                       x - 
                                   Map[ModularInvariantj[#, hprec, hrange]&,
                                        Map[ 
                                        (- #[[2]] + Sqrt[d]) / (2 #[[1]]) &, 
                                                 ClassList[d]
                                                 ]
                                             ]
                                          ]
                                    ],
                              {x}],
                hprec
              ]
             ] .
        Table[x^j,{j,0,Length[ClassList[d]]}]
       ]
                        
(* 
   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.
*)


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 = Mod[f,p], s, t, i },
  SeedRandom[101];                   (* Same random factor  *)
  While[True, 
       s = Apply[
                 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], Break[],
         s =  PolyPowerMod[s, (p-1)/2, {fp, p}] - 1 ;
         t = GCDMod[fp, s, p]; 
         If[ 0< Exponent[t,x], Break[]]
        ]
       ];
    Return[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 = Mod[f, p], gp = Mod[g, p], q,r,monic},
                     monic = Coefficient[gp, x, Exponent[gp,x]];
                     If[monic != 1, 
                        gp = Mod[PowerMod[monic, -1, p] gp, p]];
                     q = Mod[PolynomialQuotient[fp, gp, x], p];
                     r = Mod[PolynomialRemainder[fp, gp, x], p];
                     If[Exponent[r, x] <= 0,
                        If[r == 0, Return[gp], Return[1]],
                     Return[GCDMod[q,r,p]]
                     ]
             ]  /; Exponent[g,x] >= 0

(*  The following is supposed to be faster since you don[t take 
    a GCD like in the above, but I haven't been able to get rid
    of the bugs.     

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


*)

(* Repeated squaring method. This takes up a lot of time and is a 
   bottle neck for the program at present. For example if p = 3 2^324 -1
   = 100 digit number it takes 100 seconds to raise a polynomial
   to power (p-1)/2 (mod degree 4 polynomial)    
 
   Digits only exists in Version 1.2 *)

PolynomialRemainderMod[f_, g_, x_, p_Integer] := 
    Mod[PolynomialRemainder[Mod[f, p] , Mod[g, p], x], p];
   


PolyPowerMod[f_,n_Integer,{g_,p_Integer}] :=
 Last[
  Accumulate[
      PolynomialRemainderMod[#1 #2, g, x, p] &,
      First[Transpose[
         Select[
           Transpose[
              {NestList[
                 PolynomialRemainderMod[# #, g, x, p] &,
                 Mod[f,p],
                 Length[Digits[n,2]]-1],
               Reverse[Digits[n,2]]}
                    ],
          #[[2]] == 1 &
            ]
          ]
        ]
   ]
  ]   /; n>0



(* 
   This tries to find a nontrivial large prime factor of n.  
   It is unobvious how to choose PrimeTable correctly to get
   best performance. According to Morain, this part of the program
   ends up being the bottleneck.

  *)
 
PrimeTable := Block[{i},Table[Prime[i], {i,500}]];   (* or 5000 *)

LargeFactor[m_Integer,n_Integer] :=
         Block[{ q, t},
                   t = FactorInteger[ 2 m, FactorComplete -> False];
                   q = If[ IntegerQ[Last[t]], Last[t], Last[t] [[1]]]; 
                If[PrimeQ[q] &&  N[(n^(1/4)+ 1)^2] < q < m, 
                   Return[q],
                   q = m;
                   Scan[ While[ Mod[q, #] == 0, q /= #] &, PrimeTable];
                   If[PrimeQ[q] &&  N[(n^(1/4)+ 1)^2] < q < m, 
                      Return[q]
                     ]
                   ]
                ]


(*  This version seems to create core dumps if time is too small *)

(*

LargeFactor[n_Integer, m_Integer, time_Integer] :=   
              Block[{fi},
               fi = Check[TimeConstrained[FactorInteger[n],time], 0]//Timing;
                    If[ SameQ[fi[[2]], 0], 
                       Return[Null],
                      Scan[If[N[(m^(1/4) + 1)^2] < #[[1]] < n, 
                           Return[#[[1]] ] ]&,
                             fi[[2]] ]    
                      ] 
                   ]
*)



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


QuadraticRepresentation[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, 
               Return[{x p - b y/2, y/2}]
               ];   
            Return[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[]
                    ]
                  ]
               ];
          Return[g[[2]]]
        ]   


(* Shanks modular Sqrt algorithm, Knuth, Vol. 2, page 626. Note that
   finding a non residue below takes < 2 (Log [n])^2 time by GRH, and
   this avoids using Random[Integer,{1,p}]. *)

SqrtMod[d_Integer,p_Integer] := 
       Block[{e,q,z,n,v,w,k,t,a,count=0},
             e = 1; q = (p-1)/2;
             While[EvenQ[q], e++; q /= 2];
             If[ e == 1, 
                 z = -1,    
                 n = 2;                                
                 While[JacobiSymbol[n,p] != -1, n++];  
                 z = PowerMod[n,q,p]
               ];
               v = PowerMod[d,(q+1)/2,p]; 
               w = Mod[PowerMod[d,-1, p] Mod[ v v ,p], p]; 
               While[w != 1, count++; 
                     k = 1; t = Mod[w w,p];
(*  The following loop is very slow for primes of the form 
    h 2^n + 1, large  *)
                     While[t != 1, k++; t = Mod[t t,p]]; 
                     a = PowerMod[z,2^(e-k-1),p];
                     e = k;
                     v = Mod[v a, p];
                     z = Mod[a a,p];
                     w = Mod[w z,p];
                    ]; 
                Return[v]
              ]             /; JacobiSymbol[d,p] == 1


(*    You have to compute the class group of Q[Sqrt[d]]. Since d  
      will be small (< 2000), you don't need to be very efficient   *)

ClassList[d_Integer] := 
     Block[{a,b,c,list},
           a = 1; list = {};
           While[a <= N[Sqrt[-d/3]],
                 b = 1 - a;
                 While[b <= a, 
                       c= (b^2 - d)/(4 a);
                       If[Mod[c,1] == 0 && c >= a && GCD[a,b,c] == 1 &&
                            Not[a == c && b < 0],
                          AppendTo[list,{a,b,c}]
                         ];
                        b++
                       ];
                   a++
                  ];
             Return[list]
            ]     /;           d < 0 && Mod[d,4] == 1  && SquareFreeQ[d] 

ClassNumber[d_Integer] := Length[ClassList[d]];

SquareFreeQ[d_Integer] := Select[FactorInteger[d], #[[2]] > 1 &] == {};



(* *******************************************************************

    Elliptic curves, y^2 = x^3 + g2 x + g3 (mod p). Mostly how
    to add two points on an elliptic curve.

   ****************************************************************** *)

PointECQ[x_,y_,g2_,g3_,p_] := SameQ[ 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[{s,IdentityECQ,},
                   IdentityECQ = False;
                   If[ Mod[x1,p] == Mod[x2,p],
                     If[ Mod[y1,p] == Mod[-y2,p], 
                         IdentityECQ = True,
       s = Mod[
              Mod[3 Mod[x1 x1,p] + g2,p]  * PowerMod[2 y1, -1, p], p]
              ] ,
           s = Mod[Mod[y2 - y1,p]  * PowerMod[x2-x1, -1,p],p],
              ];
                  If[ SameQ[IdentityECQ , True], 
                  Return[IdentityEC],
                  x3=  Mod[ Mod[s s, p] - x1 - x2,p];
                  y3=  Mod[Mod[-y1,p] + Mod [s (x1-x3),p],p];
                  Return[PointEC[x3,y3,g2,g3,p]]
                 ]
            ]
             
 (* Multiply a point on the curve by an integer by repeated doubling, 
    exactly like PowerMod[a,b,c] 

************************************************************************

Note: This routine can be written in essentially one line
      as below. Unfortunately, this seems to run slower
      (around twice as slow for 100 digit numbers). It's
      unclear why that is (perhaps because you are storing 
      your partial results in the NestList). It would be nice 
      to fix it up.

 PointEC /:  0 * PointEC[x_,y_,g2_,g3_,p_] := IdentityEC
 PointEC /:  1 * PointEC[x_,y_,g2_,g3_,p_] := PointEC[x,y,g2,g3,p]
 PointEC /:  n_Integer * PointEC[x_,y_,g2_,g3_,p_] :=
               Block[{dgt = Digits[n,2]}, 
               Return[
   Reverse[dgt] . NestList[ # + # &, PointEC[x,y,g2,g3,p],Length[dgt]-1]
                    ] ]   /; n>1

**************************************************************************

*)



 PointEC /:  0 * PointEC[x_,y_,g2_,g3_,p_] := IdentityEC
 PointEC /:  1 * PointEC[x_,y_,g2_,g3_,p_] := PointEC[x,y,g2,g3,p]
 PointEC /:  n_Integer * PointEC[x_,y_,g2_,g3_,p_] :=
               Block[{dgt,w,sm,i},
                      dgt = Digits[n,2]; 
                      w = PointEC[x,y,g2,g3,p];
                    sm = IdentityEC;
                      i = Length[dgt];
                      While[ i >= 1,
                             If[ dgt[[i]] == 1,
                                 sm += w
                               ];
                             w = w+w;
                             i--
                           ];
                      Return[sm]
                     ]    /; n>1

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

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

(*  ********************************************************  *)  



(*
      PrimeCount[x] gives the number of primes less than or equal x
      (pi(x) in number theory terms). The prime number theorem states 
      that this is aymptotically  LogIntegrate[x]. The Riemann Hypothesis 
      asserts that the absolute error is  < constant * Sqrt[x] Log[x] 
      (good luck!).

      This algorithm is valid only up to the largest n for which
      Prime[n] is defined.

*)  


PrimeCount[x_] := PrimeCount[ Floor[x]] /; Not[IntegerQ[x]]

PrimeCount[n_Integer] := 0    /; n <= 1

PrimeCount[n_Integer] :=   
              Block[{ t, diff, m},
                    t = Round[ N[Log[n]]];
                    m = Round[ N[LogIntegral[n]]]; 
                    diff = n - Prime[m]; 
                    While[ Abs[diff] > t^2 / (1.16)^2 ,
                           m += Round[ diff /  t ];
                           diff = n - Prime[m]  ; 
                         ];
                    If[diff < 0, 
                       While[ n < Prime[m], m--], 
                       While[ n >= Prime[m], m++]; m--
                      ];
                    Return[m]
               ]                     /; n > 1



(* As above, PrimitiveRoot finds the smallest primitive root (mod p) *)

 
PrimitiveRoot[p_Integer] :=  
               Block[{ a = 2, flist}, 
                     flist = First[Transpose[FactorInteger[ p  - 1]]];
                     While[ PowerMod[ a, p - 1, p] != 1 ||  
                            Select[ Map[ PowerMod[ a, (p-1)/#, p] &, flist], 
                                    # == 1 &
                                  ] != {}, 
                            a++
                           ];
                      Return[a]
                     ]  /; PrimeQ[p] && p > 2

EndPackage[]
/* End of text from uxe.cso.uiuc.edu:sci.math.symbolic */
