ecmrpn2.1.gp

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
/*                          ecmrpn2.1.                             */
/*              ECM-Lenstra-Montgomery factorization               */
/*        for numbers of the form 'n = 2^(2m-1)-2^(m-1)+1'         */
/*              including right perfect numbers (RPN)              */
/*                                                                 */
/*  Language: Pari/GP .gp                                          */
/*  Compile (example) for ecmrpn2.1.gp:                            */
/*  gp2c-run -g ecmrpn2.1.gp                                       */
/*                                                                 */
/*  Included functions:                                            */
/*  ===================================                            */
/*  Redc(t, n, n2, ar, br, cr)                                     */
/*  pointAdd(px, pz, qx, qz, rx, rz, n, n2, ar, br, cr)            */
/*  pointDouble(px, pz, n, a24, n2, ar, br, cr)                    */
/*  ScalarMultiply(k, px, pz, n, a24, n2, ar, br, cr)              */
/*  ecm(m, B1, B2=100*B1, c=1, par=0, si=0, fa=[])                 */
/*  addhelp(ecm, "...")                                            */
/*  ecm_mt(m, B1, B2=100*B1, c=1, par=0, si=0, fa=[], thr=1)       */
/*  addhelp(ecm_mt, "...")                                         */
/*  ===================================                            */
/*                                                                 */
/*  For help type:                                                 */
/*  ?ecm                                                           */
/*  ?ecm_mt                                                        */
/*                                                                 */
/*  References:                                                    */
/*  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        */
/*                Free for non-commercial use                      */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
/* changes:                                                        */
/* v.1.2.:  better memory handling for stage 1.                    */
/* v.2.0.:  fast modular reduction for RPN                         */
/* v.2.1.:  'known factors' support                                */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/


Redc(t, n, n2, ar, br, cr) = {
  \\ this is not Montgomery's Redc core.
  \\ instead it is the modular reduction for RPN and similar numbers.
  my(a, b);
  while( t >= n2,       \\ t = bitand(t, ar) + (t >> cr) << br - t >> cr)
    a = bitand(t, ar);
    b = t >> cr;
    t = b << br;
    t-= b;
    t+= a;
  );
  if( t >= n, t-n, t)
};      \\ end Redc


pointAdd(px, pz, qx, qz, rx, rz, n, n2, ar, br, cr) = {
  my(t, u, v, upv, umv, x, z);
  t = px-pz;
  if(t < 0, t+= n);
  u = qx+qz;
  if(u > n, u-= n);
  u*= t;
  u = Redc(u, n, n2, ar, br, cr);
  t = px+pz;
  if(t > n, t-= n);
  v = qx-qz;
  if(v < 0, v+= n);
  v*= t;
  v = Redc(v, n, n2, ar, br, cr);
  upv = u+v;
  if(upv > n, upv-= n);
  umv = u-v;
  if(umv < 0, umv+= n);
  x = upv*upv;
  x = Redc(x, n, n2, ar, br, cr);
  x*= rz;
  x = Redc(x, n, n2, ar, br, cr);
  z = umv*umv;
  z = Redc(z, n, n2, ar, br, cr);
  z*= rx;
  z = Redc(z, n, n2, ar, br, cr);
  return([x, z])
};      \\ end pointAdd


pointDouble(px, pz, n, a24, n2, ar, br, cr) = {
  my(u2, v2, t, x, z);
  u2 = px+pz;
  if(u2 > n, u2-= n);
  u2*= u2;
  u2 = Redc(u2, n, n2, ar, br, cr);
  v2 = px-pz;
  if(v2 < 0, v2+= n);
  v2*= v2;
  v2 = Redc(v2, n, n2, ar, br, cr);
  t = u2-v2;
  if(t < 0, t+= n);
  x = u2*v2;
  x = Redc(x, n, n2, ar, br, cr);
  z = a24*t;
  z = Redc(z, n, n2, ar, br, cr);
  z+= v2;
  if(z > n, z-= n);
  z*= t;
  z = Redc(z, n, n2, ar, br, cr);
  return([x, z])
};      \\ end pointDouble


scalarMultiply(k, px, pz, n, a24, n2, ar, br, cr) = {
  my(lk, qx, qz, rx, rz, ik=100, ti=0, ts);
  lk=exponent(k);
  qx=px;
  qz=pz;
  [rx, rz] = pointDouble(px, pz, n, a24, n2, ar, br, cr);
  \\ bittest(k, i=lk) is always 1 - the neutral element of multiplication - so skip it :)
  forstep(i=lk-1, 0, -1,
    if(bittest(k, i),
      [qx, qz] = pointAdd(rx, rz, qx, qz, px, pz, n, n2, ar, br, cr);
      [rx, rz] = pointDouble(rx, rz, n, a24, n2, ar, br, cr),
      \\ else
      [rx, rz] = pointAdd(qx, qz, rx, rz, px, pz, n, n2, ar, br, cr);
      [qx, qz] = pointDouble(qx, qz, n, a24, n2, ar, br, cr)
    );
    if(default(debug) > 3,
      if(i%ik,,
        ti+=gettime();
        ts=ti\1000;
        printf("<%.2f%% %dsec.> ", (lk-i-1)/lk*100., ts)
      )
    )
  );
  return([qx, qz])
};      \\ end scalarMultiply


ecm(m, B1, B2=100*B1, c=1, par=0, si=0, fa=[]) = {
  my(ar, br, cr, n, n2, k, k2, b1, bl, sig=0, u, v, vmu, a, t, a24, px, pz, qx, qz, g=1);
  my(nb, dd, db, dv, beta, b, rx, rz, tx, tz, alpha, delta, step, f, trx, trz, ti, ts, ik=20);
  my(B1now, B1next, B1first, B1step = 100000);

  \\ build RPN out of m.
  n = 2^(2*m-1)-2^(m-1)+1;
  n2= 2*n;
  ar = 2^(2*m-1)-1;
  br = m-1;
  cr = 2*m-1;

  \\ stage 1 parameters.
  B1 = floor(B1);
  k = 1;
  b1 = sqrtint(B1);
  bl = log(B1);

  \\ stage 2 precomputation.
  B2 = floor(B2);
  nb = exponent(n)+1;      \\ #bits of 'n'.
   if(
    B2 < 3600, dd = 30; db = 4,       \\ elseif...
    B2 < 32400, dd = 90; db = 12,
    B2 < 176400, dd = 210; db = 24,
    B2 < 705600, dd = 420; db = 48,
    B2 < 6350400, dd = 1260; db = 144,
    B2 < 21344400, dd = 2310; db = 240,
    B2 < 42688800, dd = 4620; db = 480,
    B2 < 384199200, dd = 13860; db = 1440,
    B2 < 901800900, dd = 30030; db = 2880,
    B2 < 3607203600, dd = 60060; db = 5760,
    B2 < 16232416200, dd = 180180; db = 17280,
    dd = 510510; db = 46080
  );

  dv = vector(db);
  t=0;
  forstep( i = 1, dd/2, 2,
    if( gcd(i, dd) == 1,
      t++;
      dv[t] = i

    )
  );

  \\ compute a B1-powersmooth integer 'k'
  forprime(p = 2, b1, k*= p^(bl\log(p)));

  \\ for memory issues we set B1first = min(B1, B1step).
  B1first = min(B1, B1step);
  k*= factorback(primes([b1+1, B1first]));
  gettime();
  if(si, c = 1);


  \\ curve loop
  for(i = 1, c,
    while( sig==0,          \\ iferr then choose another sigma
      if(si,
        sig = si; setrand(si),                       \\else
        sig = random([6,2^32-1])        \\ 5 < sigma < 2^32
      );

      \\ 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,                   \\ parameterization used in Alpertron's ECM. sigma 2: ...
        u = 2*sig/(3*sig^2-1);
        a = (-3*u^4-6*u^2+1)/(4*u^3);
        px = 3*u^2+1;
        pz = 4*u
      );

      if( par == 3,                   \\ no parameterization: set a24 = sigma. sigma 3: ...
        a = 4*sig-2;
        px = 2;
        pz = 1
      );

      if( par < 1 || par > 3,         \\ 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
      );

      if(default(debug)>2, printf("sigma=%d:%d D=%d\n", par, sig, dd));

      \\ Normalizing step done once at the beginning of stage 1.
      \\ As 'px', 'pz' and 'a24' may be rationals at this point, an 'impossible inverse' error
      \\ can occur by mod 'n', usually if 'n' has small factors. Pari/GP's 'iferr' handles
      \\ this and returns a factor before stage 1 begins. This is marked as:
      \\ stage 0 factor found ?

      iferr(px%= n, E, sig = 0; si = 0; g = lcm(gcd(lift(component(E,2)),n),g));  \\ inverse mod handling.
      iferr(pz%= n, E, sig = 0; si = 0; g = lcm(gcd(lift(component(E,2)),n),g));
      a24 = a+2;
      a24/= 4;
      iferr(a24%= n, E, sig = 0; si = 0; g = lcm(gcd(lift(component(E,2)),n),g));
      for(i = 1, #fa, g\= gcd(g, fa[i]));     \\ in case of known factors, remove them.
      if(g > 1 && g < n, return([g, sig, i, 0]))
    );

    \\ stage 1 - standard double-add ladder

    k2=k;
    B1next = B1first;
    B1now = B1next;

    if( B1now == B1,
      [qx, qz] = scalarMultiply(k, px, pz, n, a24, n2, ar, br, cr),     \\ else

      \\  for memory issues we split B1 into B1steps.
      while ( B1now < B1,
        [px, pz] = scalarMultiply(k2, px, pz, n, a24, n2, ar, br, cr);
        B1now = B1next;
        B1next = min(B1, B1next + B1step);
        k2 = factorback(primes([B1now+1, B1next]));
      );
      qx = px;
      qz = pz
    );

    g = gcd(n, qz);         \\ factor found stage 1 ?
    if(default(debug) > 3, print);
    if(default(debug) > 1, print("c",i,"s1 <",gettime()," ms.>"));
    for(i = 1, #fa, g\= gcd(g, fa[i]));     \\ in case of known factors, remove them.
    if(g > 1 && g < n, return([g, sig, i, 1]));

    \\ stage 2 - advanced standard continuation

      g = 1;
      beta = vector(db);
      delta = vector(db);
      b = (B1\dd+1)*dd;     \\ 'b' divisble by 'dd' and next greater 'B1'.

    for(i = 1, db,      \\ create table delta[i], beta[i].
       delta[i] = scalarMultiply(dv[i], qx, qz, n, a24, n2, ar, br, cr);
       beta[i] = delta[i][1]*delta[i][2];
       beta[i] = Redc(beta[i], n, n2, ar, br, cr)
    );

    step = scalarMultiply(dd, qx, qz, n, a24, n2, ar, br, cr);
    [rx, rz] = scalarMultiply(b, qx, qz, n, a24, n2, ar, br, cr);
    [tx, tz] = scalarMultiply(b-dd, qx, qz, n, a24, n2, ar, br, cr);
    if(default(debug) > 3, print);
    if(default(debug) > 1, print("c",i,"s2pre <",gettime()," ms.>"));
    \\ step 2 main loop
    forstep(r = b, B2, dd,
      alpha = rx*rz;
      alpha = Redc(alpha, n, n2, ar, br, cr);
      for(i = 1, db,
        \\ do a pseudo primality test only if n is => 1024 bits.
        if(nb < 1024 || ispseudoprime(r-dv[i]) || ispseudoprime(r+dv[i]),
          t = rx-delta[i][1];
          if(t < 0, t+= n);
          f = rz+delta[i][2];
          if(f > n, f-= n);
          f*= t;
          f = Redc(f, n, n2, ar, br, cr);
          f-= alpha;
          f+= beta[i];
          if(f > n, f-= n, f < 0, f+= n);
          g*= f;
          g = Redc(g, n, n2, ar, br, cr)
        )
      );
      trx = rx;
      trz = rz;
      [rx, rz] = pointAdd(rx, rz, step[1], step[2], tx, tz, n, n2, ar, br, cr);
      tx = trx;
      tz = trz;
      if(default(debug) > 3,
        if(r%(dd*ik),,
          ti+=gettime();
          ts=ti\1000;
          printf("<%.2f%% %dsec.> ", (r-b)/B2*100., ts)
        )
      )
    );

    g = gcd(n, g);          \\ factor found stage 2 ?
    if(default(debug) > 3, print);
    if(default(debug) > 1, print("c",i,"s2 <",gettime()," ms.>"));
    for(i = 1, #fa, g\= gcd(g, fa[i]));     \\ in case of known factors, remove them.
    if(g > 1 && g < n, return([g, sig, i, 2]));

    sig = 0        \\ Stage 2 done. Reset sig for more curves.
  );
};     \\ end ecm


addhelp(ecm, "ecm(m, B1, {B2=100*B1}, {c=1}, {par=0}, {si=0}, {fa=[]}) elliptic curve method for numbers of the form \n\t'n = 2^(2*m-1)-2^(m-1)+1' including right perfect numbers (RPN) (V2.1). \n\nInputs are: \n'm' exponent of the number to be factored 'n = 2^(2*m-1)-2^(m-1)+1'. 'n' must be composite. \n'B1' the bound for stage 1. 'B1' <= 10^10 must be even. \nIf specified: \n'B2' {default=100*B1} the bound for stage 2. 'B2' <= 10^14. \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\t3: no parameterization: set a24 = sigma. \n'si' {default=0} sets the sigma-value for one curve, if omitted a pseudo-random value 5 < 'si' < 2^32 is choosen.\n'fa' {default=[]} vector of one or more known factors e.g.: '[<factor1>, <factor2>, <...>]'. \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 entries is returned. \nFrom left to right these entries are: \n\t- a non trivial factor of 'n', \n\t- the generating sigma value, \n\t- the curve and \n\t- the stage where the factor was found. \n\nThe debug level for increased output can be set with '\\g n' command while 'n' ranges from 0 to 20. \n\t'n' >= 2, timings for stage 1 and 2. \n\t'n' >= 3, sigma and stage 2 D-values. \n\t'n' >= 4, more detailed timings (for curves running several hours). \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 4Gbytes in this case. ");


ecm_mt(m, B1, B2=100*B1, c=1, par=0, si=0, fa=[], thr=1) = {
  my(i=0, R, S);
  forstep(b = 0, c-1, thr,
    S = vector(min(thr, c-b), s, random([6,2^32-1]));
    R = parapply(s -> ecm(m, B1, B2, 1, par, s, fa), S);
    for(t=1, #R,
      i++;
      if(R[t],
        R[t][3] = i;
        return(R[t])
      )
    );
    if(default(debug) > 0, printf("<c%d/%d> ", i, c))
  )
};      \\ end ecm_mt


addhelp(ecm_mt, "ecm_mt(m, B1, {B2=100*B1}, {c=1}, {par=0}, {si=0}, {fa=[]}, {thr=1}) framework for ecm to run multi-threaded (if available). \n\n Inputs are the same as for ecm (see ?ecm) except: \n'thr' {default=1} sets the number of curves to run in parallel \n\nFor memory management it is recommended to increase the 'threadsizemax' before running ecm_mt(). \nTo do this, enter for example: 'default(threadsizemax, 1024*1048576)' which sets the max size to 1Gbyte for each thread in this case.");


export(ecm);        \\ export to the parallel world.
export(scalarMultiply);
export(pointDouble);
export(pointAdd);
export(Redc);


\\ end ecmrpn2.1.gp