\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\            ECM-Lenstra-Montgomery Implementation v.1.1          \\
\\  Language: Pari/GP .gp                                          \\
\\  Compile (example) for ecmm.gp:                                 \\
\\  gp2c-run -g ecmm.gp                                            \\
\\                                                                 \\
\\  Included functions:                                            \\
\\  ===================================                            \\
\\  pointAdd(px, pz, qx, qz, rx, rz, n)                            \\
\\  pointDouble(px, pz, n, a24)                                    \\
\\  scalarMultiply(k, px, pz, n, a24)                              \\
\\  ecm(n, B1, B2=100*B1, c=1, par=0, si=0)                        \\            
\\  addhelp(ecm, "...")                                            \\
\\  factorecm(n, dig=10, par=0)                                    \\
\\  addhelp(factorecm, "...")                                      \\
\\  ===================================                            \\
\\                                                                 \\ 
\\  For help type:                                                 \\
\\  ?ecm                                                           \\ 
\\  ?factorecm                                                     \\ 
\\                                                                 \\
\\  Reference:                                                     \\
\\  https://www.rieselprime.de/ziki/Elliptic_curve_method          \\
\\  https://rosettacode.org/wiki/Fermat_numbers#Go                 \\ 
\\  https://www.alpertron.com.ar/ECM.HTM                           \\
\\  https://members.loria.fr/PZimmermann/records/ecm/params.html   \\
\\                                                                 \\
\\  Martin Hopf 2022     mailto : athlontiger@googlemail.com       \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

  
pointAdd(px, pz, qx, qz, rx, rz, n) = { local(t, u, v, upv, umv, x, z);
  t = px-pz;
  u = qx+qz;
  u*= t;
  t = px+pz;
  v = qx-qz;
  v*= t;
  upv = u+v;
  umv = u-v;
  x = upv*upv;
  x*= rz;
  x%= n;
  z = umv*umv;
  z*= rx;
  z%= n;
  return([x, z])
};      \\ end pointAdd
  

pointDouble(px, pz, n, a24) = { local(u2, v2, t, x, z);
  u2 = px+pz;
  u2*= u2;
  v2 = px-pz;
  v2*= v2;
  t = u2-v2;
  x = u2*v2;
  x%= n;
  z = a24*t;
  z+= v2;
  z*= t;
  z%= n;
  return([x, z])
};      \\ end pointDouble


scalarMultiply(k, px, pz, n, a24) = { local(sk, lk, qx, qz, rx, rz);
  sk=binary(k);
  lk=#sk;
  qx=px;
  qz=pz;
  [rx, rz] = pointDouble(px, pz, n, a24);
  for(i=2, lk,      \\ sk[i=1] is always 1 so don't use it :) 
    if(sk[i],
      [qx, qz] = pointAdd(rx, rz, qx, qz, px, pz, n);
      [rx, rz] = pointDouble(rx, rz, n, a24),           
      \\ else
      [rx, rz] = pointAdd(qx, qz, rx, rz, px, pz, n);
      [qx, qz] = pointDouble(qx, qz, n, a24)
    )
  );
  return([qx, qz])
};      \\ end scalarMultiply


ecm(n, B1, B2=100*B1, c=1, par=0, si=0) = {
  local(k, b1, bl, u, v, vmu, a, t, a24, px, pz, qx, qz, g);
  local(dd, beta, s, d, d2, b, tx, tz, step, alpha, limit, f, trx, trz);

  k = 1;
  b1 = sqrtint(B1); 
  bl = log(B1);
  
  \\ compute a B1-powersmooth integer 'k'
  forprime(p = 2, b1, k*= p^floor(bl/log(p)));
  k*= factorback(primes([b1+1,B1]));
  gettime();
  if(si,c=1); 
  for(i = 1, c,
    if(si,sig=si,sig=random([6,2^32-1]));      \\ 5 > sigma > 2^32
    if(default(debug)>2, print("sigma=",sig));
    
    \\ generate a new curve in Montgomery form 
      
    if( par == 1,                   \\ 64bit parameterization equals GMP-ECM -sigma 1: ...
      a = 4*sig^2/2^64-2;
      px = 2; 
      pz = 1
    );

    if( par == 2,                   \\ Alpertron's parameterization -sigma 2: ...
      u = 2*sig/(3*sig^2-1);
      pz = 4*u;
      a = (-3*u^4-6*u^2+1)/(4*u^3);
      a24 = (a+2)/4;
      px = 3*u^2+1
    );

    if( par < 1 || par > 2,       \\ default Suyama's parameterization equals GMP-ECM -sigma: 0: ...
      par = 0;     
      u = sig*sig; u-= 5; v = 4*sig; vmu = v-u; a = vmu*vmu; a*= vmu; t = 3*u;        
      t+= v; a*= t; t = 4*u; t*= u; t*= u; t*= v; a/= t; a-= 2; 
      px = u*u; px*= u; t = v*v; t*= v; px/= t; pz = 1
    );
      
    \\ stage 1

    a24 = a+2;
    a24/= 4;

    [qx, qz] = scalarMultiply(k, px, pz, n, a24);

    g = gcd(n, qz);
    if(default(debug) > 1, print("c",i,"s1 <",gettime()," ms.>"));
    if(g > 1 && g < n, return([g, sig, i, 1]));

    \\ stage 2 
    dd = sqrtint(B2);
    beta = vector(dd+1);
    s = vector(2*dd+2);
    
    [s[1], s[2]] = pointDouble(qx, qz, n, a24);
    [s[3], s[4]] = pointDouble(s[1], s[2], n, a24);
    beta[1] = s[1]*s[2];
    beta[1]%= n;
    beta[2] = s[3]*s[4];
    beta[2]%= n;
    for(d = 3, dd,
      d2 = 2*d;  
      [s[d2-1], s[d2]] = pointAdd(s[d2-3], s[d2-2], s[1], s[2], s[d2-5], s[d2-4], n);
      beta[d] = s[d2-1]*s[d2];
      beta[d]%= n
    );    
    g = 1;
    b = B1-1;
    [rx, rz] = scalarMultiply(b, qx, qz, n, a24);
    t = 2*dd;
    t = b-t;
    [tx, tz] = scalarMultiply(t, qx, qz, n, a24);
    step = 2*dd;
    forstep(r = B1-1, B2, step,
      alpha = rx*rz;
      alpha%= n;
      limit = r + step;
      forprime(q = r+1, limit,
        d = (q-r)/2;
        t = rx-s[2*d-1];
        f = rz+s[2*d];
        f*= t;
        f-= alpha;
        f+= beta[d];
        g*= f;  
        g%= n
      );
      trx = rx;
      trz = rz;
      [rx, rz] = pointAdd(rx, rz, s[2*dd-1], s[2*dd], tx, tz, n);
      tx = trx;    
      tz = trz
    );
    g = gcd(n, g);
    if(default(debug) > 1, print("c",i,"s2 <",gettime()," ms.>"));
    if(g > 1 && g < n, return([g, sig, i, 2]))
  );
};     \\ end ecm


addhelp(ecm, "ecm(n, B1, {B2=100*B1}, {c=1}, {par=0}, {si=0}) performs the elliptic curve method by Lenstra and Montgomery. \n\nInputs are: \n'n' the number to be factored. \n'B1' the bound for stage 1. \nIf specified: \n'B2' {default=100*B1} indicates the bound for stage 2. \n'c' {default=1} indicates the number of curves to run. \n'par' {default=0} sets the parameterization for the curve. Values are: \n\t0: Suyama's parameterization equals GMP-ECM -sigma 0:... . \n\t1: parameterization for 64-bit processors equals GMP-ECM -sigma 1:... . \n\t2: parameterization used in Alpertron's ECM. \n'si' {default=0} sets the sigma-value for one curve, otherwise a pseudo-random value between 5 and 2^32 is choosen. \n\nNote on pseudo-random numbers: if a new Pari/GP-session is started 'random()' will always create the \nsame sequence of pseudo-random numbers. However if you use 'setrand(n)' at the begin of a session \nwith 'n' is a wildly typed number < 2^64, a different sequence of pseudo-random numbers will be created. \n\nIf a factor (not necessarily prime) was found within the specified boundaries, a vector with 4 elements is returned. \nFrom left to right these elements are: \n\t- a non trivial factor of 'n', \n\t- the corresponding sigma value, \n\t- the curve and \n\t- the stage where the factor was found. \n\nThe debug level can be set with '\\g n' command while 'n' ranges from 0 to 20. \n\tIf debug level is set >= 2, timings for stage 1 and 2 are shown. \n\tIf debug level is set >= 3, additionally each sigma-value is shown. \n\nFor memory management it is recommended to increase the 'parisizemax' before running ecm(). \nTo do this, enter for example: 'default(parisizemax, 4096*1048576)' which sets the max size to 4Gbyte in this case. \n\nAvoid 'n' having small factors, as they often result in early 'impossible inverse mod' errors.");  


factorecm(n, dig=9, par=0) = {
  local(B1, B2, c, M, i, t, d, cd, R, v, bdig, s=10^6, p=2);
  i=2;d = 0; 
  bdig = dig*log(10)/log(2);
  printf("\n");

  \\ trial-divide to find small factors <= s.
  while( p <= s, 
    if( n%p, 
      p = nextprime(p+1),       \\ else
      printf("P%d=%d found by trial division.\n", floor(log(p)/log(10))+1, p);
      n/=p;
    )
  );

  \\ Parameter-Matrix [bits,B1,curves] for bits[30-200] step 10 bits. See: https://members.loria.fr/PZimmermann/records/ecm/params.html
  M=[30, 1358, 2; 40, 1630, 10; 50, 12322, 9; 60, 21906, 21; 70, 32918, 66; 80, 76620, 119; 90, 183850, 219; 100, 445658, 339; 110, 1305196, 439; 120, 3071166, 649; 130, 4572524, 1507; 140, 9267682, 2399; 150, 22025674, 3159; 160, 35158748 ,6076; 170, 47862548, 14038; 180, 153319098, 12017; 190, 410593604, 13174;  200, 491130496, 29584];

  if(bdig < M[1,1], 
    i = 1,          \\ else
    if(bdig > M[#M[,1],1],  
      i = #M[,1],   \\ else
      while(bdig > M[i,1], i++)  
    )    
  );

  printf("\n");
  while(!ispseudoprime(n) && i <= #M[,1],
    c = M[i,3]-d;
    B1 = M[i,2];
    B2 = M[i,1]*B1;
    cd = floor(log(n)/log(10))+1;
    t = M[i,1]*log(2)/log(10);
    if(c, 
      if(c==1,v="",v="s");
      printf("          %d curve%s on C%d with B1=%d and B2=%d (t=%.2f).\n", c, v, cd, B1, B2, t);
    );
    if(c && R=ecm(n, B1, B2, c, par),
        if(ispseudoprime(R[1]),v="P",v="C");
        printf("\n%s%d=%d found @ curve %d stage %d sigma %d:%d\n\n", v, floor(log(R[1])/log(10))+1, R[1], R[3], R[4], par, R[2]);
        n/= R[1];
        d+= R[3],     \\ else
        i++;
        d=0
      
    ) 
  );
  printf("P%d factoring complete!\n\n",floor(log(n)/log(10))+1)
};    \\ end factorecm


addhelp(factorecm, "factorecm(n, {dig=9}, {par=0}) factoring routine for ecm() function. \n\nInputs are: \n'n' the number to be factored. \n If specified: \n'dig' {default=9} the (minimum) number of digits that an expected factor of 'n' has. \n'par' {default=0} parameterization value for which ecm() initializes the curves. \n\nFor more details see '?ecm'.");
 
\\ end ecmm.gp