/* for gmp-ecm version 7.x: for parameter sigma = 0:s */
/* also for gmp-ecm version 6.x: for sigma = s */
FindGroupOrder(p,s)={
  my(K,v,u,x,b,a,A,E);
  K = Mod(1,p);
  v = K*(4*s);
  u = K*(s^2-5);
  x = u^3;
  b = 4*x*v;
  a = (v-u)^3*(3*u+v);
  A = a/b-2;
  x = x/v^3;
  b = x^3 + A*x^2 + x;
  E = ellinit([0,b*A,0,b^2,0],K);
  return(ellcard(E));
}

FindGroupOrderA(p,A)={
  my(K, d, a, b, E);
  K = Mod(1,p);
  d = K*((A+2)/4);
  a = K*(4*d-2);
  b = K*(16*d+2);
  E = ellinit([0,a/b,0,1/b^2,0],K);
  return(ellcard(E));
}

/* for parameter sigma = 1:s */
FindGroupOrderParam1(p,s)={
  return(FindGroupOrderA(p, 4*s^2/2^64-2));
}

/* for parameter sigma = 2:s */
FindGroupOrderParam2(p,s)={
  my(K,E,P,x,y,x3,A);
  K = Mod(1,p);
  E = ellinit([0,36],K);
  [x,y] = ellmul(E, [-3,3], s);
  x3 = (3*x+y+6)/(2*(y-3));
  A = -(3*x3^4+6*x3^2-1)/(4*x3^3);
  return(FindGroupOrderA(p, A));
}

/* for parameter sigma = 3:s */
FindGroupOrderParam3(p,s)={
   return(FindGroupOrderA(p, 4*s/2^32-2));
}

FindGroupOrderParam(p, sigma, param) = {
  if (param == 0, return(FindGroupOrder(p, sigma)));
  if (param == 1, return(FindGroupOrderParam1(p, sigma)));
  if (param == 2, return(FindGroupOrderParam2(p, sigma)));
  if (param == 3, return(FindGroupOrderParam3(p, sigma)));
  print("Invalid parametrization: ", param);
}

/*
# check if a prime p is found with bounds B1 and B2,
# for parameter 'param' and sigma in [sigma_min,sigma_max-1]
# check_found (31622776601683800097, 11000, 1873422, 0, 1000)
# check_found (31622776601683800097, 11000, 1873422, 1, 1000)
# check_found (31622776601683800097, 11000, 1873422, 2, 1000)
# check_found (31622776601683800097, 11000, 1873422, 3, 1000)
*/
check_found(p, B1, B2, param, sigma_max) = {
  my(e2=0,e3=0,tries=0,found=0,sigma,f);
  for(sigma=0,sigma_max-1,
    iferr(f = factor(FindGroupOrderParam(p, sigma, param)),
          E, next(), 1);
    f = factor(FindGroupOrderParam(p, sigma, param));
    tries += 1;
    if (f[1,1] != 2,
      print(" * Error 1,1 != 2");
      print("factors = ",f);
      return();
    );
    e2 += f[1,2];
    if (f[2,1] == 3,
      e3 += f[2,2];
    );
    ms=matsize(f)[1];
    if (f[ms-1,1] <= B1 && f[ms,1] <= B2,
      found += 1;
    );
  );
  printf("tries=%d, found=%d, %0.8f %0.8f %0.8f \n",tries,found,1.0*e2/tries,1.0*e3/tries,2.0^(e2/tries)*3.0^(e3/tries));
}

/* check all parametrizations 0, 1, 2, 3 */
check_found_all(p, B1, B2, sigma_max) = {
  for (param=0,3,
    check_found(p,B1,B2,param,sigma_max);
  );
}

/*
sample run:

check_found_all(31622776601683800097, 11000, 1873422, 1000)
*/