// Returns either 2*P or GCD(N,x1-x2) =/= 1
function double(P,a,N)  
   x,y,z := Explode(P);
   if z eq 0 then  // point at infinity
      return P;
   end if;
   g,_,y2inv := XGCD(N,Integers()!(2*y));
   if g ne 1 then
      return g;
   end if;
   xx := ((x^2-a)^2 - 8*x)*y2inv^2;
   yy := ((3*x^2 + a)*(x - xx) - 2*y^2)*y2inv;
   return [xx,yy,1];
end function;

// Returns P + Q or GCD(N,x1-x2) =/= 1
function add(P,Q,a,N)
   if P eq Q then
      return double(P,a,N);
   end if;
   x1,y1,z1 := Explode(P);
   x2,y2,z2 := Explode(Q);
   if z1 eq 0 then
      return Q;
   elif z2 eq 0 then
      return P;
   end if;
   if x1 eq x2 and y1 eq -y2 then
      return [0,1,0];
   end if;
   g,_,inv := XGCD(N,Integers()!(x1-x2));
   if g ne 1 then
      return g;
   end if;
   lambda := (y1-y2)*inv;
   nu := y1 - lambda*x1;
   x3 := lambda^2 -x1-x2;
   y3 :=  -lambda*x3-nu;
   return [x3,y3,1];
end function;

// Try to compute R=m*[0,1,1] on y^2=x^3+ax+1; returns
// either R or GCD(N,some denominator) =/= 1 if not possible.
function multiply(m,a,N)
   // Points are represented as triples [x,y,z] with z either 0 or 1.
   P := [IntegerRing(N)|0,1,1];   
   R := [IntegerRing(N)|0,1,0];
   while m ne 0 do
      if IsOdd(m) then
         R := add(R,P,a,N);
         if Type(R) eq RngIntElt then
            return R;
         end if;
      end if;         
      m := Floor(m/2);
      P := double(P,a,N);
      if Type(P) eq RngIntElt then
         return P;
      end if;
   end while;
   return R;
end function;

intrinsic ECM1(N::RngIntElt, m::RngIntElt, a::RngIntElt) -> RngIntElt
{Try to find a B-power smooth factor of N using Lenstra's ECM with given a
and m=lcm(1,...,B). Returns N on failure.}
   printf "Trying a = %o: \t", a;
   if GCD(4*a^3 + 27, N) ne 1 then
      print "Split using discriminant.";
      return GCD(4*a^3 + 27, N);
   end if;
   R := multiply(m,a,N);
   if Type(R) eq RngIntElt then
      printf "Failed to compute mP. ";
      if R lt N then
         print "Split using denominator.";
         return R;
      end if;
         print "Denominator gives no factor.";
   end if;
   print "Computed mP (no factor found).";
   return N;
end intrinsic;


intrinsic ECM(N::RngIntElt, B::RngIntElt, maxtries::RngIntElt) -> RngIntElt
{Try to find a B-power smooth factor of N using Lenstra's ECM. 
Returns N on failure.}
   m := LCM([1..B]);
   for i in [1..maxtries] do
      a := Random(N);
      M := ECM1(N,m,a);
      if M ne N then
         return M;
      end if;
   end for;
   print "Max tries exceeded.  Trying changing B.";   
   return N;
end intrinsic;