home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMBERTH.PAK / FACTORIN.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  14.2 KB  |  364 lines
  1.  
  2. (*:Version: Mathematica 2.0 *)
  3.  
  4. (*:Name: NumberTheory`FactorIntegerECM` *)
  5.  
  6. (*:Title: Integer Factoring Using Lenstra's Elliptic Curve Method *)
  7.  
  8. (*:Author: Ilan Vardi
  9.  
  10. *)
  11.  
  12. (*:Keywords: FactorInteger, Lenstra's Elliptic Curve Method
  13. *)
  14.  
  15. (*:Requirements: none. *)
  16.  
  17. (*:Warnings: The algorithm returns a single factor (not necessarily
  18.            a prime). To obtain a complete factorization, run
  19.            FactorInteger or FactorIntegerECM on the factor 
  20.            and cofactor. The algorithm is probabilistic, so there could
  21.            be a large variance in running times, even for similar
  22.            inputs. Note that SeedRandom[101] is used to generate 
  23.            pseudorandom numbers for the algorithm, so the program 
  24.            will always run the exactly the same on the same input.
  25.  
  26.            The input should under no circumstances be a prime number!           
  27.  
  28. *)
  29.  
  30. (*:Source: P.L. Montgomery, Speeding up the Pollard and Elliptic
  31.            Curve Methods of Factorization, Mathematics of 
  32.            Computation #48 (1987), pages 243-264
  33.  
  34. *)
  35.  
  36. (* :Limitations: It is assumed that the number was not factored
  37.            by FactorInteger, and so its smallest prime factor is
  38.            at least 10^12. The algorithm is optimized to find factors
  39.            of 17 digits or less, so it should factor most numbers of 40
  40.            digits or less (and such numbers will probably only have
  41.            two prime factors if they are hard to factor). 
  42.  
  43.            The program has some memory requirements. Stage one sets
  44.            up a table of 10,000 binary digits. Stage two requires a
  45.            table of up to 30,000 primes, and a table of 2400 numbers 
  46.            mod n.
  47. *)
  48.  
  49. (* :Discussion: The program is an implementations of EC as in Montgomery's 
  50.            paper above. The first stage uses the Elliptic curve 
  51.            parametrization of Peter Montgomery, which allows one to do
  52.            elliptic curve addition without modular inversion.
  53.            If {x,z} represents the rational number X = x/z of the
  54.            curve
  55.  
  56.            (1)   B Y^2 = X^3 + A X^2 + X
  57.  
  58.            and a = (A + 2)/4 (you pick a first, and can assume that 
  59.            a < Sqrt[n]).  You multiply a point with X = 2 by the LCM 
  60.            of numbers less than upperlimit = Prime[upperlimitp], where
  61.            upperlimit is best taken to be (according to Richard Crandall)
  62.  
  63.            Exp[Sqrt[Log[q] Log[Log[q]]/2]]
  64.  
  65.            for finding a factor < q. The default is therefore q=Sqrt[n].
  66.            At the end of the first stage, you generate an elliptic
  67.            curve in Weierstrass form with X coordinate corresponding
  68.            to the X coordinate generated by the first stage (you 
  69.            take Y = 1 in equation (1)). Let pt = {X,1}, then one
  70.            tries p*pt for 
  71.  
  72.            Prime[upperlimitp1] < p < Prime[upperlimitp2]
  73.  
  74.            This implementation uses Prime[upperlimitp2], about 25
  75.            times as large as Prime[upperlimitp1]. This was chosen 
  76.            so that stage one and stage two take about the same amount
  77.            of time.
  78.  
  79.            Instead of making (upperlimitp2-upperlimitp1) elliptic
  80.            curve additions, one can speed this up to about 
  81.            Sqrt[2/3 Prime[upperlimitp2]] additions as follows:
  82.  
  83.            If p = q w - r, w fixed and p*pt is the identity,
  84.            then (qw)*pt = r*pt, and the X coordinates corresponding
  85.            to these two are equal. One therefore picks w to be 
  86.            about Sqrt[Prime[upperlimitp1]] and constructs two tables
  87.            of X coordinates of (qw)*pt and r*pt (the qw table can
  88.            be omitted since the primes are taken in ascending order).
  89.            Letting w be divisible by 6, the remainders r are either
  90.            r = 1 or r = 5. Since pt and -pt have the same X coordinate,
  91.            one just makes a list of numbers of the form 6k+1 and
  92.            chooses p = q w + r or p = (q+1) w - (w - r) depending on
  93.            on Mod[r, 6] = 1 or 5. One then compares the X coordinates
  94.            of (qw)*pt and r*pt or (q+1)w *pt and (w-r)*pt. As in the 
  95.            Pollard rho method, you take cumulative produts of these
  96.            differences and test for GCD with n periodically. 
  97.  
  98.            Constructing the table of multiples of pt requires using
  99.            the addition rule for points on an elliptic curve. The
  100.            number of modular inversions can be reduced by taking a 
  101.            number of curves at a time and trading multiplications for
  102.            one modular inversion. 
  103.  
  104.            If n is small (n < 10^20), then the method will often 
  105.            succeed after one curve with second stage. In other words
  106.            the method is essentially Pollard p-1 with second stage.
  107.            For larger n one uses several curves at a time to reduce
  108.            GCDs and interpretation overhead in stage 2.
  109.  
  110. *)
  111.  
  112.  
  113. BeginPackage["NumberTheory`FactorIntegerECM`"]
  114.  
  115. FactorIntegerECM::usage = "FactorIntegerECM[n] finds a single factor
  116. of the composite number n using Lenstra's Elliptic Curve Method."
  117.  
  118. FactorSize::usage = "The option FactorIntegerECM[n, FactorSize -> q]
  119. optimizes the algorithm to find factors of size at most q. The 
  120. default is q = Sqrt[n]."
  121.  
  122. CurveNumber::usage = "The option FactorIntegerECM[n, CurveNumber -> b]
  123. specifies that the algorithm will use b curves at a time for each 
  124. iteration of the algorithm. This number must be a power of 2."
  125.  
  126. CurveCountLimit::usage ="The option FactorIntegerECM[n, CurveCountLimit->c]
  127. limits the total number of curves the algorithm will use. This allows the 
  128. user to abort the algorithm after a given number of iterations. The default is
  129. 10^4 curves. FactorIntegerECM[n, CurveCountLimit->c] returns n if no
  130. nontrivial factor has been found."
  131.  
  132.  
  133. Begin["`Private`"]  
  134.  
  135.  
  136. Off[PowerMod::ninv]         (* A factor is found by a failed division 
  137.                               (Mod n), so turn off the error message 
  138.                               for PowerMod[m,-1,n] when GCD[m,n] > 1 *) 
  139.  
  140. FactorIntegerECM[n_, opts___]:= 
  141. Block[{b1, b2, 
  142.        factorsize = FactorSize /. {opts} /. {FactorSize -> Sqrt[N[n]]},
  143.        upperlimitp1, upperlimitp2, sqrtn = Floor[Sqrt[N[n]]], 
  144.        digits, i, plen, 
  145.       curvenumber = CurveNumber /. {opts} /. 
  146.                      {CurveNumber -> If[n > 10^30, 8, If[n > 10^20, 4, 2]]},
  147.        a, x, x2, xw, w, digitsw, prod, j,
  148.        pt, ptw, pt6, inv3 = PowerMod[3, -1, n], gx, gxw, 
  149.        ppt, qptw, qqptw, qq, rr, q, r, ip, tablept, 
  150.        curvecountlimit = 
  151.        CurveCountLimit /. {opts} /. {CurveCountLimit -> 10^4},
  152.        curvecount = 0, stage   (* curvecount is the number of curves used 
  153.                                   so far, stage identifies stage of 
  154.                                   algorithm (1 or 2) *)
  155.        }, 
  156.        If[PrimeQ[n], Return[n]]; 
  157.        SeedRandom[101];   (* Same random seed *) 
  158.        curvenumber = 2^Round[N[Log[2, curvenumber]]]; 
  159.                      (* curvenumber must be a power of 2 *) 
  160.        b1 = Min[N[Exp[Sqrt[Log[factorsize]/2 Log[Log[factorsize]]]]], 
  161.                   13. 10^3];
  162.        upperlimitp1 = Round[N[b1/Log[b1]]];
  163.        upperlimitp2 = 21 upperlimitp1;
  164.        plen = upperlimitp2 - upperlimitp1; 
  165.        b1 = Prime[upperlimitp1]; 
  166.        b2 = Prime[upperlimitp2]; 
  167.        digits = Rest[IntegerDigits[LCM @@ 
  168.                     Range[Prime[upperlimitp1]], 2]];
  169.  
  170. (* w is about sqrt(6 b2) and is divisible by 6. *)
  171.  
  172.        w = 6 Ceiling[Sqrt[N[6 b2]]/6]; 
  173.        digitsw = Rest[IntegerDigits[w, 2]];
  174.        
  175.        factor = 1; 
  176.  
  177.        While[!(1 < factor < n) && (curvecount + 1) <= curvecountlimit,  
  178.  
  179. (* First stage *)
  180.           For[i = 1; x = a = {}, i <= curvenumber, i++,
  181.               curvecount++; 
  182.               sa = Random[Integer, {1, sqrtn}];
  183.               sx = mult[{2, 1}, digits, sa, n];
  184.               factor = GCD[sx[[2]], n];
  185.               If[1 < factor, Break[]];
  186.               AppendTo[a, sa]; AppendTo[x, sx];
  187.               ];
  188.           If[1 < factor, stage = 1; Continue[]];
  189.  
  190. (* Second stage *)  stage = 2; 
  191.  
  192.         (* Elliptic curve addition where a factor of n is detected
  193.            by an illegal inversion mod n. This is implemented by 
  194.            a Throw/Catch construction. *)
  195.  
  196.         Catch[x = Mod[(First /@ x) ListInverse[Last /@ x, n], n];
  197.               a = 4 a - 2;              
  198.               b = Mod[x (x (x + a) + 1), n];
  199.               t = Select[GCD[#, n]& /@ (b (a^2 - 4)), 1 < #&];
  200.               If[t != {}, factor = Min[t]; Throw[True]];
  201.               y = ListInverse[b, n];
  202.               x = Mod[(x + inv3 a) y, n];
  203.               g2 = Mod[(1 - inv3 a^2) y^2, n];
  204.               pt = PointECV[x, y, g2, n];
  205.  
  206. (*            pt corresponds to the x,z generated by stage 1 *)
  207.  
  208.               pt6 = 6 pt;
  209.               ptw = w pt; 
  210.               tablept = First /@ NestList[# + pt6&, pt, w/6];
  211.               ip = upperlimitp1;
  212.               q = Quotient[Prime[ip], w]; 
  213.               qptw = q ptw;
  214.               qqptw = {qptw, qptw + ptw};
  215.               r = Mod[Prime[ip], w];
  216.               For[i = 0, i < plen, i += upperlimitp1, 
  217.                   For[j = 1; prod = Table[1, {curvenumber}], 
  218.                       j <= upperlimitp1, j++, 
  219.                       r += (Prime[ip + 1] - Prime[ip]);
  220.                       ip++;
  221.                       While[r > w, 
  222.                             r -= w; q++; 
  223.                             qqptw = {qqptw[[2]], qqptw[[2]] + ptw}];
  224.                       If[q == 0, Continue[]];
  225.                       If[Mod[r, 6] == 1, 
  226.                          rr = r; qq = 1,
  227.                          rr = w - r; qq = 2];
  228.                       rr = (rr+5)/6;
  229.                       gx = tablept[[rr]];
  230.                       gxw = First[qqptw[[qq]]];
  231.                       prod = Mod[(gx - gxw) prod, n]; 
  232.                       ];
  233.                   factor = GCD[ListProduct[prod, n], n];  
  234.  
  235.           (*  The case factor = n is improbable so bactracking to
  236.               find a nontrivial factor has been omitted. *)
  237.  
  238.                   If[1 < factor < n, Break[]]]
  239.             ]];
  240.           If[factor == 1, n, factor]] 
  241.            
  242.               
  243. (* X coordinate of 2*P, P in Montgomery form *)
  244.  
  245. double[{x_, z_, a_, m_}]:= 
  246.        Block[{t = Mod[(x+z)^2, m], u = Mod[(x-z)^2, m], v},
  247.               v = t - u;
  248.               {Mod[t u, m], Mod[v (u + a v), m]}]
  249.  
  250. (* X coordinate of (2m+1)*P from m*P and (m+1)*P *)
  251.  
  252. doubleplusone[{x1_, z1_, xn_, zn_, xm_, zm_, m_}]:= 
  253.    Block[{u = Mod[(xm - zm) (xn + zn), m], 
  254.           v = Mod[(xm + zm) (xn - zn), m]},
  255.          {Mod[z1 (u + v)^2, m], Mod[x1 (u - v)^2, m]}]
  256.  
  257. (* X coordinate of {2m*P, (2m+1)*P} from {m*P, (m+1)*P} *)
  258.  
  259. even[{x1_, z1_, xn_, zn_, xm_, zm_, a_, m_}]:= 
  260. Join[{x1, z1}, double[{xn, zn, a, m}], 
  261.      doubleplusone[{x1, z1, xn, zn, xm, zm, m}], {a, m}]
  262.  
  263. (* X coordinate of {(2m+1)*P, (2m+2)*P} from {m*P, (m+1)*P} *)
  264.  
  265. odd[{x1_, z1_, xn_, zn_, xm_, zm_, a_, m_}]:= 
  266. Join[{x1, z1}, doubleplusone[{x1, z1, xn, zn, xm, zm, m}], 
  267.      double[{xm, zm, a, m}], {a, m}]
  268.  
  269. (* Compute k*P, where P has X coordinate x1:z1. Digits of k are 
  270.    computed only once by ECM algorithm *)
  271.  
  272. mult[{x1_, z1_}, digits_, a_, m_]:= 
  273.      Fold[If[#2 == 0, even[#1], odd[#1]]&,
  274.           Join[{x1, z1}, {x1,z1},double[{x1, z1, a, m}], {a, m}],
  275.                digits] [[{3,4}]]
  276.  
  277.  
  278. (* *******************************************************************
  279.  
  280.     Rules for adding points on the Elliptic curve 
  281.  
  282.          y^2 = x^3 + g2 x + g3 (mod m) 
  283.  
  284.     using lists of numbers to save GCDs. This implementation is 
  285.     not quite correct, since the test x1 == x2 will not detect if
  286.     equality holds in one component only (this is a low probability
  287.     event). This loss simplifies the program and saves time.
  288.  
  289.    ****************************************************************** *)
  290.  
  291. PointECV /: PointECV[x_,y_,g2_,m_] + IdentityEC:= PointECV[x,y,g2,m]
  292.  
  293. PointECV /:  PointECV[x1_,y1_,g2_,m_] + PointECV[x2_,y2_,g2_,m_]:=
  294.    Block[{slope, IdentityECQ = False},
  295.           If[x1 == x2,
  296.              If[Mod[y1, m] == Mod[-y2, m], IdentityECQ = True,
  297.                 slope = Mod[(3 x1^2 + g2) ListInverse[2 y1, m], m]],
  298.              slope = Mod[(y2 - y1) ListInverse[x2 - x1, m], m],
  299.               ]; 
  300.           If[SameQ[IdentityECQ, True], 
  301.              IdentityEC,
  302.              x3 = Mod[slope^2 - x1 - x2, m];
  303.              y3 = Mod[slope (x1 - x3) - y1, m];
  304.              PointECV[x3, y3, g2, m]
  305.                  ]
  306.             ]
  307.  
  308.  (* Multiply a point on the curve by an integer by repeated doubling, 
  309.     exactly like PowerMod[a,b,c].  *)
  310.  
  311. IdentityEC /: n_Integer * IdentityEC:= IdentityEC
  312.  
  313. PointECV /: 0 * PointECV[x_, y_, g2_, m_]:= IdentityEC
  314.  
  315. PointECV /: 1 * PointECV[x_, y_, g2_, m_]:= PointECV[x, y, g2, m]
  316.  
  317. PointECV /: 2 * PointECV[x_, y_, g2_, m_]:= 
  318.                 PointECV[x, y, g2, m] + PointECV[x, y, g2, m]
  319.  
  320. PointECV /: n_ * PointECV[x_, y_, g2_, m_]:= 
  321.                  (-n) * PointECV[x, -y, g2, m] /; n < 0
  322.  
  323. PointECV /: n_ * PointECV[x_, y_, g2_, m_]:= 
  324.   Fold[2 #1 + #2&, IdentityEC, 
  325.        IntegerDigits[n, 2] PointECV[x, y, g2, m]]  /;  n < 2^300 
  326.  
  327. (* For large n use base 64 to speed up the calculation by 25-> 50%.
  328.    If used repeatedly, precompute digits in main algorithm. *)
  329.  
  330. PointECV /: n_ * PointECV[x_, y_, g2_, m_]:= 
  331.   Fold[64 #1 + #2&, IdentityEC, 
  332.        NestList[# + PointECV[x, y, g2, m]&, 
  333.                 IdentityEC, 63] [[IntegerDigits[n, 64] + 1]]]
  334.  
  335.  
  336. (* Taking modular inverses of a list with a single inversion. *)
  337.  
  338. ListProduct[v_List, m_Integer]:= Fold[Mod[#1 #2, m]&, 1, v]
  339.  
  340. PrepareInverse[{a_Integer}, m_Integer]:= {1}
  341.  
  342. PrepareInverse[v_List, m_Integer]:=  
  343. Block[{u = Partition[v, Length[v]/2]},
  344.        Mod[Flatten[(PrepareInverse[#, m]& /@ u) 
  345.                    (ListProduct[#, m]& /@ Reverse[u])],
  346.                          m]]
  347.  
  348. ListInverse[v_List, m_Integer]:=   
  349.  Block[{pi = PrepareInverse[v, m], inv},
  350.         inv = PowerMod[First[pi] First[v], -1, m];
  351.         If[NumberQ[inv], Mod[inv pi, m], 
  352.            factor = Min[Select[GCD[#, m]& /@ v, 1 < #&]]; 
  353.            Throw[True]]]
  354.  
  355.  
  356. End[]   (* NumberTheory`FactorIntegerECM`Private` *)
  357.  
  358. Protect[FactorIntegerECM, FactorSize, CurveNumber, CurveCountLimit]
  359.  
  360. EndPackage[]   (* NumberTheory`FactorIntegerECM` *)
  361.  
  362.  
  363.  
  364.