// A very simple demo version of the elliptic curve method
// It uses only one curve 
// A more "professional" version can be found in numlib::ecm

findcurve:=
proc(n)
// returns parameters and a point as list [a,b,x,y] 
local a,x,y;

begin
x:=random() mod (n-1) +1;
y:=random() mod (n-1) +1;
 
// this guarantees x,y \neq 0, which gives a "useful" point P 

a:=random() mod n ;
[a,(y^2-x^3-a*x) mod n,x,y]

end_proc:


  

curveadd:=
proc(n,En,P1,P2)
local x1,x2,y1,y2,d,lambda,x3;

begin
if P1=[0,0] then return (P2); end_if;
if P2=[0,0] then return (P1); end_if;
x1:=op(P1,1); y1:=op(P1,2); x2:=op(P2,1); y2:=op(P2,2);
if x1=x2 then
	if y1=-y2 then return([0,0]); end_if;
	d:=igcdex(2*y1,n);
	if op(d,1)<>1 then return(op(d,1)); end_if;
	lambda:=op(d,2)*(3*x1*x1+op(En,1));
else
	d:=igcdex(x1-x2,n);
	if op(d,1)<>1 then return(op(d,1)); end_if;
	lambda:=op(d,2)*(y1-y2);
end_if;
x3:=(lambda*lambda-x1-x2) mod n;
[x3 ,(lambda*(x1-x3)-y1) mod n] 

end_proc:

convertbase2 := 
proc(n) 
local l;
begin
   l:=[];
   while n>0 do
      l:=append(l,n mod 2);
      n:=n div 2
   end_while;
   l
end_proc:

curvemult:=
proc(n,En,k,P)
local i,l,PP;

begin
if k=0 then return([0,0]); end_if;
if k=1 then return(P); end_if;
l:=convertbase2(k);
for i from nops(l)-1 downto 1 do
	PP:=curveadd(n,En,P,P);
	if domtype(PP)=DOM_INT then return(PP); end_if;
        if l[i] mod 2 = 1 then 
		PP:=curveadd(n,En,P,PP); 
        	if domtype(PP)=DOM_INT then return(PP); end_if;
        end_if; 
end_for;
PP
end_proc:


curvemul2:=
proc(n,En,P,s,B)
local p,Q;

begin
p:=2;
result:=[0,0];
logB:=ln(float(B));
while p<=s do
        Q:=curvemult(n,En,floor(logB/ln(float(p))),P);
	if domtype(Q)=DOM_INT then return(Q); end_if;
	result:=curveadd(n,En,result,Q);
	if domtype(result)=DOM_INT then return(result); end_if; 
	p:=nextprime(p+1);
end_while;
result
end_proc:
 
factorize:=
proc(n)

begin
if isprime(n) then return(n); end_if;
c:=findcurve(n);
result:=curvemul2(n,[op(c,1..2)],[op(c,3..4)],floor(float(sqrt(n))),n);
if domtype(result)=DOM_INT then return (result,factorize(n/result)); end_if;
return(FAIL);
end_proc: