(*:Mathematica Version: 2.0 *)

(*:Package Version: 1.0 *)

(* :Copyright: Copyright 1992-1999,  Wolfram Research, Inc. *)

(*:Name: NumberTheory`FactorIntegerECM` *)

(*:Title: Integer Factoring Using Lenstra's Elliptic Curve Method *)

(*:Author: Ilan Vardi
*)

(*:Summary: 
This package provides a function for finding a single factor of a
composite number using Lenstra's elliptic curve method.
*)

(*:Keywords: FactorInteger, Lenstra's Elliptic Curve Method
*)

(*:Requirements: None. *)

(*:Warnings: The algorithm returns a single factor (not necessarily
           a prime). To obtain a complete factorization, run
           FactorInteger or FactorIntegerECM on the factor 
           and cofactor. The algorithm is probabilistic, so there could
           be a large variance in running times, even for similar
           inputs. Note that SeedRandom[101] is used to generate 
           pseudorandom numbers for the algorithm, so the program 
           will always run exactly the same on the same input.

           The input should under no circumstances be a prime number!           

*)

(*:Source: P.L. Montgomery, Speeding up the Pollard and Elliptic
           Curve Methods of Factorization, Mathematics of 
           Computation #48 (1987), pages 243-264.

*)

(* :Limitations: It is assumed that the number was not factored
           by FactorInteger, and so its smallest prime factor is
           at least 10^12. The algorithm is optimized to find factors
           of 17 digits or less, so it should factor most numbers of 40
           digits or less (and such numbers will probably only have
           two prime factors if they are hard to factor). 

           The program has some memory requirements. Stage one sets
           up a table of 10,000 binary digits. Stage two requires a
           table of up to 30,000 primes, and a table of 2400 numbers 
           mod n.
*)

(* :Discussion: The program is an implementations of EC as in Montgomery's 
           paper above. The first stage uses the elliptic curve 
           parametrization of Peter Montgomery, which allows one to do
           elliptic curve addition without modular inversion.
           If {x,z} represents the rational number X = x/z of the
           curve

           (1)   B Y^2 = X^3 + A X^2 + X

           and a = (A + 2)/4 (you pick a first, and can assume that 
           a < Sqrt[n]).  You multiply a point with X = 2 by the LCM 
           of numbers less than upperlimit = Prime[upperlimitp], where
           upperlimit is best taken to be (according to Richard Crandall)

           Exp[Sqrt[Log[q] Log[Log[q]]/2]]

           for finding a factor < q. The default is therefore q=Sqrt[n].
           At the end of the first stage, you generate an elliptic
           curve in Weierstrass form with the X coordinate corresponding
           to the X coordinate generated by the first stage (you 
           take Y = 1 in equation (1)). Let pt = {X,1}, then one
           tries p*pt for 

           Prime[upperlimitp1] < p < Prime[upperlimitp2]

           This implementation uses Prime[upperlimitp2], about 25
           times as large as Prime[upperlimitp1]. This was chosen 
           so that stage one and stage two take about the same amount
           of time.

           Instead of making (upperlimitp2-upperlimitp1) elliptic
           curve additions, one can speed this up to about 
           Sqrt[2/3 Prime[upperlimitp2]] additions as follows:

           If p = q w - r, w fixed and p*pt is the identity,
           then (qw)*pt = r*pt, and the X coordinates corresponding
           to these two are equal. One therefore picks w to be 
           about Sqrt[Prime[upperlimitp1]] and constructs two tables
           of X coordinates of (qw)*pt and r*pt (the qw table can
           be omitted since the primes are taken in ascending order).
           Letting w be divisible by 6, the remainders r are either
           r = 1 or r = 5. Since pt and -pt have the same X coordinate,
           one just makes a list of numbers of the form 6k+1 and
           chooses p = q w + r or p = (q+1) w - (w - r) depending on
           on Mod[r, 6] = 1 or 5. One then compares the X coordinates
           of (qw)*pt and r*pt or (q+1)w *pt and (w-r)*pt. As in the 
           Pollard rho method, you take cumulative products of these
           differences and test for GCD with n periodically. 

           Constructing the table of multiples of pt requires using
           the addition rule for points on an elliptic curve. The
           number of modular inversions can be reduced by taking a 
           number of curves at a time and trading multiplications for
           one modular inversion. 

           If n is small (n < 10^20), then the method will often 
           succeed after one curve with second stage. In other words
           the method is essentially Pollard p-1 with second stage.
           For larger n one uses several curves at a time to reduce
           GCDs and interpretation overhead in stage 2.

*)


BeginPackage["NumberTheory`FactorIntegerECM`"]

FactorIntegerECM::usage = "FactorIntegerECM[n] finds a single factor
of the composite number n using Lenstra's Elliptic Curve Method."

FactorSize::usage = "FactorSize is an option for FactorIntegerECM.
FactorSize -> q optimizes the algorithm to find factors of size at most q. The
default is Automatic; q is chosen to be Sqrt[n], where n is the integer to be
factored."

CurveNumber::usage = "CurveNumber is an option for FactorIntegerECM.
CurveNumber -> b specifies that the algorithm will use b curves at a time 
for each iteration of the algorithm. This number must be a power of 2.
The default is Automatic; b is chosen to be either 2, 4, or 8 (increasing
with n, the integer to be factored)."

CurveCountLimit::usage ="CurveCountLimit is an option for FactorIntegerECM.
CurveCountLimit->c limits the total number of curves the algorithm will use
to c. This allows the user to abort the algorithm after a given number of
iterations. The default is 10^4 curves.  FactorIntegerECM[n, CurveCountLimit->c]
returns n if no nontrivial factor has been found."


Unprotect[FactorIntegerECM, FactorSize, CurveNumber, CurveCountLimit]
Begin["`Private`"]  

Options[FactorIntegerECM] =
{FactorSize -> Automatic, CurveNumber -> Automatic, CurveCountLimit -> 10^4} 

FactorIntegerECM[n_, opts___]:= 
Block[{ninv, factorsize, curvenumber, curvecountlimit,
       b1, upperlimitp1, upperlimitp2, plen, b2, digits, w, digitsw, factor,
       curvecount = 0, (* curvecount is the number of curves used so far *)
       stage,   (* stage identifies stage of algorithm (1 or 2) *)
       sqrtn, i, x, a, sa, sx,
       inv3, b, t, y, g2, pt, pt6, ptw, tablept, ip, q, qptw, qqptw, r,
       j, prod, rr, qq, gx, gxw
       },
       (* Added the test for EvenQ because the stage 1 loop was infinite
		for n = 4, 8, 10, and 12.  This test also protects against
		Min::nord errors;  for n = 4 and 6, there were
	        Min::nord errors because Log[Log[factorsize]] is negative
		for factorsize < E.  (ECM 961016)
	*)
       If[EvenQ[n], Return[2]];	
       (* Added the test for the following composites, because the
		algorithm does not terminate for these values.  (ECM 961016) *)
       If[MemberQ[{1, 9, 25, 27, 49, 81, 125}, n],
		Return[FactorInteger[n][[1, 1]]]
       ];
       If[PrimeQ[n], Return[n]]; 
     (* A factor is found by a failed division 
     		(Mod n), so turn off the error message 
   		for PowerMod[m,-1,n] when GCD[m,n] > 1 *) 
       ninv = (Head[PowerMod::ninv] === $Off);
       Off[PowerMod::ninv]; 
       {factorsize, curvenumber, curvecountlimit} =
		{FactorSize, CurveNumber, CurveCountLimit} /. {opts} /.
			Options[FactorIntegerECM];
       If[factorsize === Automatic,
		factorsize = Sqrt[N[n]]];
       If[curvenumber === Automatic,
		curvenumber = If[n > 10^30, 8, If[n > 10^20, 4, 2]]];
       SeedRandom[101];   (* Same random seed *) 
       curvenumber = 2^Round[N[Log[2, curvenumber]]]; 
                     (* curvenumber must be a power of 2 *) 
       (* Note: the following gives a negative value for b1
	  for factorsize < E *)
       b1 = Min[N[Exp[Sqrt[Log[factorsize]/2 Log[Log[factorsize]]]]], 
                  13. 10^3];
       upperlimitp1 = Round[N[b1/Log[b1]]];
       upperlimitp2 = 21 upperlimitp1;
       plen = upperlimitp2 - upperlimitp1; 
       b1 = Prime[upperlimitp1]; 
       b2 = Prime[upperlimitp2]; 
       digits = Rest[IntegerDigits[LCM @@ 
                    Range[Prime[upperlimitp1]], 2]];

       (* w is about sqrt(6 b2) and is divisible by 6. *)

       w = 6 Ceiling[Sqrt[N[6 b2]]/6]; 
       digitsw = Rest[IntegerDigits[w, 2]];
       sqrtn = Floor[Sqrt[N[n]]];
       inv3 = PowerMod[3, -1, n];
       
       factor = 1; 

       While[!(1 < factor < n) && (curvecount + 1) <= curvecountlimit,  

        (* First stage *)
        For[i = 1; x = a = {}, i <= curvenumber, i++,
              curvecount++; 
              sa = Random[Integer, {1, sqrtn}];
              sx = mult[{2, 1}, digits, sa, n];
              factor = GCD[sx[[2]], n];
              If[1 < factor, Break[]];
              AppendTo[a, sa]; AppendTo[x, sx];
        ]; (* end For *)
        If[1 < factor, stage = 1; Continue[]];

        (* Second stage *)  stage = 2; 

        (* Elliptic curve addition where a factor of n is detected
           by an illegal inversion mod n. This is implemented by 
           a Throw/Catch construction. *)

        (* Note: ListInverse sets factor *)
        Catch[x = Mod[(First /@ x) ListInverse[Last /@ x, n], n];
              a = 4 a - 2;              
              b = Mod[x (x (x + a) + 1), n];
              t = Select[GCD[#, n]& /@ (b (a^2 - 4)), 1 < #&];
              If[t != {}, factor = Min[t]; Throw[True]];
              (* Note: ListInverse sets factor *)
              y = ListInverse[b, n];
              x = Mod[(x + inv3 a) y, n];
              g2 = Mod[(1 - inv3 a^2) y^2, n];
              pt = PointECV[x, y, g2, n];

	     (* pt corresponds to the x,z generated by stage 1 *)

              pt6 = 6 pt;
              ptw = w pt; 
              tablept = First /@ NestList[# + pt6&, pt, w/6];
              ip = upperlimitp1;
              q = Quotient[Prime[ip], w]; 
              qptw = q ptw;
              qqptw = {qptw, qptw + ptw};
              r = Mod[Prime[ip], w];
              For[i = 0, i < plen, i += upperlimitp1, 
                  For[j = 1; prod = Table[1, {curvenumber}], 
                      j <= upperlimitp1, j++, 
                      r += (Prime[ip + 1] - Prime[ip]);
                      ip++;
                      While[r > w, 
                            r -= w; q++; 
                            qqptw = {qqptw[[2]], qqptw[[2]] + ptw}];
                      If[q == 0, Continue[]];
                      If[Mod[r, 6] == 1, 
                         rr = r; qq = 1,
                         rr = w - r; qq = 2];
                      rr = (rr+5)/6;
                      gx = tablept[[rr]];
                      gxw = First[qqptw[[qq]]];
                      prod = Mod[(gx - gxw) prod, n]; 
                  ]; (* end For *)
                  factor = GCD[ListProduct[prod, n], n];  

                  (*  The case factor = n is improbable so bactracking to
                      find a nontrivial factor has been omitted. *)

                  If[1 < factor < n, Break[]]
	      ] (* end For *)
        ] (* end Catch *)

       ]; (* end While *)
       If[!ninv, On[PowerMod::ninv]];
       If[factor == 1, n, factor]
] (* end FactorIntegerECM *) 
           
              
(* X coordinate of 2*P, P in Montgomery form *)

double[{x_, z_, a_, m_}]:= 
       Block[{t = Mod[(x+z)^2, m], u = Mod[(x-z)^2, m], v},
              v = t - u;
              {Mod[t u, m], Mod[v (u + a v), m]}]

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

doubleplusone[{x1_, z1_, xn_, zn_, xm_, zm_, m_}]:= 
   Block[{u = Mod[(xm - zm) (xn + zn), m], 
          v = Mod[(xm + zm) (xn - zn), m]},
         {Mod[z1 (u + v)^2, m], Mod[x1 (u - v)^2, m]}]

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

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

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

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

(* Compute k*P, where P has X coordinate x1:z1. Digits of k are 
   computed only once by ECM algorithm *)

mult[{x1_, z1_}, digits_, a_, m_]:= 
     Fold[If[#2 == 0, even[#1], odd[#1]]&,
          Join[{x1, z1}, {x1,z1},double[{x1, z1, a, m}], {a, m}],
               digits] [[{3,4}]]


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

    Rules for adding points on the Elliptic curve 

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

    using lists of numbers to save GCDs. This implementation is 
    not quite correct, since the test x1 == x2 will not detect if
    equality holds in one component only (this is a low probability
    event). This loss simplifies the program and saves time.

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

PointECV /: PointECV[x_,y_,g2_,m_] + IdentityEC:= PointECV[x,y,g2,m]

PointECV /:  PointECV[x1_,y1_,g2_,m_] + PointECV[x2_,y2_,g2_,m_]:=
   Block[{slope, IdentityECQ = False},
          If[x1 == x2,
             (* Note: ListInverse sets factor *)
             If[Mod[y1, m] == Mod[-y2, m], IdentityECQ = True,
                slope = Mod[(3 x1^2 + g2) ListInverse[2 y1, m], m]
	     ],
                slope = Mod[(y2 - y1) ListInverse[x2 - x1, m], m]
          ]; 
          If[SameQ[IdentityECQ, True], 
             IdentityEC,
             x3 = Mod[slope^2 - x1 - x2, m];
             y3 = Mod[slope (x1 - x3) - y1, m];
             PointECV[x3, y3, g2, m]
          ]
   ]

 (* Multiply a point on the curve by an integer by repeated doubling, 
    exactly like PowerMod[a,b,c].  *)

IdentityEC /: n_Integer * IdentityEC:= IdentityEC

PointECV /: 0 * PointECV[x_, y_, g2_, m_]:= IdentityEC

PointECV /: 1 * PointECV[x_, y_, g2_, m_]:= PointECV[x, y, g2, m]

PointECV /: 2 * PointECV[x_, y_, g2_, m_]:= 
                PointECV[x, y, g2, m] + PointECV[x, y, g2, m]

PointECV /: n_ * PointECV[x_, y_, g2_, m_]:= 
                 (-n) * PointECV[x, -y, g2, m] /; n < 0

PointECV /: n_ * PointECV[x_, y_, g2_, m_]:= 
  Fold[2 #1 + #2&, IdentityEC, 
       IntegerDigits[n, 2] PointECV[x, y, g2, m]]  /;  n < 2^300 

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

PointECV /: n_Integer * PointECV[x_, y_, g2_, m_]:= 
  Fold[64 #1 + #2&, IdentityEC, 
       NestList[# + PointECV[x, y, g2, m]&, 
                IdentityEC, 63] [[IntegerDigits[n, 64] + 1]]]


(* Taking modular inverses of a list with a single inversion. *)

ListProduct[v_List, m_Integer]:= Fold[Mod[#1 #2, m]&, 1, v]

PrepareInverse[{a_Integer}, m_Integer]:= {1}

PrepareInverse[v_List, m_Integer]:=  
Block[{u = Partition[v, Length[v]/2]},
       Mod[Flatten[(PrepareInverse[#, m]& /@ u) 
                   (ListProduct[#, m]& /@ Reverse[u])],
                         m]]

ListInverse[v_List, m_Integer]:=   
 Block[{pi = PrepareInverse[v, m], inv},
        inv = PowerMod[First[pi] First[v], -1, m];
        If[NumberQ[inv],
	   Mod[inv pi, m], 
           factor = Min[Select[GCD[#, m]& /@ v, 1 < #&]]; 
           Throw[True]
	]
 ]


End[]   (* NumberTheory`FactorIntegerECM`Private` *)

Protect[FactorIntegerECM, FactorSize, CurveNumber, CurveCountLimit]

EndPackage[]   (* NumberTheory`FactorIntegerECM` *)



