
(*:Version: Mathematica 2.0 *)

(*:Name: NumberTheory`FactorIntegerECM` *)

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

(*:Author: Ilan Vardi

*)

(*: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 the 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 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 produts 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 = "The option FactorIntegerECM[n, FactorSize -> q]
optimizes the algorithm to find factors of size at most q. The 
default is q = Sqrt[n]."

CurveNumber::usage = "The option FactorIntegerECM[n, 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."

CurveCountLimit::usage ="The option FactorIntegerECM[n, CurveCountLimit->c]
limits the total number of curves the algorithm will use. 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."


Begin["`Private`"]  


Off[PowerMod::ninv]         (* 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 *) 

FactorIntegerECM[n_, opts___]:= 
Block[{b1, b2, 
       factorsize = FactorSize /. {opts} /. {FactorSize -> Sqrt[N[n]]},
       upperlimitp1, upperlimitp2, sqrtn = Floor[Sqrt[N[n]]], 
       digits, i, plen, 
      curvenumber = CurveNumber /. {opts} /. 
                     {CurveNumber -> If[n > 10^30, 8, If[n > 10^20, 4, 2]]},
       a, x, x2, xw, w, digitsw, prod, j,
       pt, ptw, pt6, inv3 = PowerMod[3, -1, n], gx, gxw, 
       ppt, qptw, qqptw, qq, rr, q, r, ip, tablept, 
       curvecountlimit = 
       CurveCountLimit /. {opts} /. {CurveCountLimit -> 10^4},
       curvecount = 0, stage   (* curvecount is the number of curves used 
                                  so far, stage identifies stage of 
                                  algorithm (1 or 2) *)
       }, 
       If[PrimeQ[n], Return[n]]; 
       SeedRandom[101];   (* Same random seed *) 
       curvenumber = 2^Round[N[Log[2, curvenumber]]]; 
                     (* curvenumber must be a power of 2 *) 
       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]];
       
       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];
              ];
          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. *)

        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]];
              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]; 
                      ];
                  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[]]]
            ]];
          If[factor == 1, n, factor]] 
           
              
(* 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,
             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_ * 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` *)



