Elliptic Curves<Text-field layout="Heading 1" style="Heading 1">Cubic curves</Text-field>plot([sqrt(x^3-2*x+4),-sqrt(x^3-2*x+4)],x=-2..3);<Text-field layout="Section" style="Section">The group law and addition formulas</Text-field>ecadd:=proc(A,B,b1,b2)
if A[1]=B[1]
then
if A[2]=-B[2]
then
"Identity";
else
[((3*A[1]^2+2*b2*A[1]+b1)/(2*A[2]))^2-b2-2*A[1],
-(((3*A[1]^2+2*b2*A[1]+b1)/(2*A[2]))*(((3*A[1]^2+2*b2*A[1]+b1)/(2*A[2]))^2-b2-3*A[1])+A[2])];
end if;
else
[((B[2]-A[2])/(B[1]-A[1]))^2-b2-A[1]-B[1],
-(((B[2]-A[2])/(B[1]-A[1]))*(((B[2]-A[2])/(B[1]-A[1]))^2-b2-2*A[1]-B[1])+A[2])];
end if;
end proc:ecadd([-1,4],[2,5],0,0);ecadd([2,-3],[2,-3],1,-2);ecadd([2,-3],[2,3],1,-2);<Text-field layout="Heading 1" style="Heading 1"><Font executable="false">Elliptic curves mod p</Font></Text-field>for x from 0 to 10 do
[x,modp(x^3+1,11)]
end do;countpointsmodp:=proc(b0,b1,b2,p)
local c, t;
c:=1;
for t from 0 to p-1 do
c:=c+1+numtheory[legendre](t^3+b2*t^2+b1*t+b0,p);
end do;
c;
end proc:countpointsmodp(1,0,0,11);countpointsmodp(1,1,0,541);<Text-field layout="Section" style="Section">Application: Encryption via elliptic curves</Text-field>countpointsmodp(1,2,0,1377359);converttopoint:=proc(m,b0,b1,b2,p)
local x0, f;
x0:=1000*m;
f:=x->x^3+b2*x^2+b1*x+b0;
while not numtheory[legendre](f(x0),p)= 1 do
x0:=x0+1;
end do;
[x0,modp(f(x0)&^((p+1)/4),p)]
end proc:converttopoint(1234,1,2,0,1377359);gcd(11111,1375269);gcd(22222,1375269);modp(1/11111,1375269);modp(1/22222,1375269);<Text-field layout="Section" style="Section">Application: Elliptic curve method of factorization</Text-field>gcd(modp(2&^(100!),137703491)-1,137703491);pollardfactor:=proc(m)
local a, n;
a:=2;
n:=1;
while gcd(a-1,m)=1 do
a:=modp(a&^n,m);
n:=n+1;
end do;
gcd(a-1,m);
end proc: pollardfactor(137703491);modp(2/5,7);<Text-field layout="Section" style="Section">Fermat's last theorem</Text-field>sort(expand(q*product((1-q^(4*n))^2*(1-q^(8*n))^2,n = 1..2)),q,ascending);series(q*product((1-q^(4*n))^2*(1-q^(8*n))^2,n = 1..2),q = 0,10);series(q*product((1-q^(4*n))^2*(1-q^(8*n))^2,n = 1 .. 7),q = 0,30);for t from 2 to 10 do
[ithprime(t),countpointsmodp(0,-1,0,ithprime(t))]
end do;coeff(expand(q*product((1-q^(4*n))^2*(1-q^(8*n))^2,n = 1 .. 16)),q^65);<Text-field layout="Section" style="Section">Notes: Elliptic Curve Calculations</Text-field>ecaddprojective:=proc(A,B,b1,b2)
if A[3]=0 then return(B) else
if B[3]=0 then return(A) else
if A[1]=B[1] then
if A[2]=-B[2] then
return([0,1,0]) else
return([((3*A[1]^2+2*b2*A[1]+b1)/(2*A[2]))^2-b2-2*A[1],
-(((3*A[1]^2+2*b2*A[1]+b1)/(2*A[2]))*(((3*A[1]^2+2*b2*A[1]+b1)/(2*A[2]))^2-b2-3*A[1])+A[2]),1])
end if;
else return([((B[2]-A[2])/(B[1]-A[1]))^2-b2-A[1]-B[1],
-(((B[2]-A[2])/(B[1]-A[1]))*(((B[2]-A[2])/(B[1]-A[1]))^2-b2-2*A[1]-B[1])+A[2]),1])
end if;
end if;
end if;
end proc: ecpointmultiple:=proc(A,b1,b2,k)
local q, r, l, i;
l:=convert(k,base,2);
q:=[0,1,0];
r:=[A[1],A[2],A[3]];
for i from 1 to nops(l) do
if l[i]= 1 then
q:=ecaddprojective(q,r,b1,b2)
end if;
r:=ecaddprojective(r,r,b1,b2);
end do;
return(q);
end proc:ecpointmultiple([-1,4,1],0,0,3);ecaddprojectivemodp:=proc(A,B,b1,b2,p)
if A[3]=0 then return(B) else
if B[3]=0 then return(A) else
if modp(A[1],p)=modp(B[1],p) then
if modp(A[2],p)=modp(-B[2],p) then
return([0,1,0]) else
return([modp(((3*A[1]^2+2*b2*A[1]+b1)/(2*A[2]))^2-b2-2*A[1],p),
modp(-(((3*A[1]^2+2*b2*A[1]+b1)/(2*A[2]))*(((3*A[1]^2+2*b2*A[1]+b1)/(2*A[2]))^2-b2-3*A[1])+A[2]),p),1])
end if;
else return([modp(((B[2]-A[2])/(B[1]-A[1]))^2-b2-A[1]-B[1],p),
modp(-(((B[2]-A[2])/(B[1]-A[1]))*(((B[2]-A[2])/(B[1]-A[1]))^2-b2-2*A[1]-B[1])+A[2]),p),1])
end if;
end if;
end if;
end proc: ecpointmultiplemodp:=proc(A,b1,b2,k,p)
local q, r, l, i;
l:=convert(k,base,2);
q:=[0,1,0];
r:=[A[1],A[2],A[3]];
for i from 1 to nops(l) do
if l[i]= 1 then
q:=ecaddprojectivemodp(q,r,b1,b2,p)
end if;
r:=ecaddprojectivemodp(r,r,b1,b2,p);
end do;
return(q);
end proc:ecpointmultiplemodp([1234005, 349433, 1],2,0,11111,1377359);ecpointmultiplemodp([1114312, 498654, 1],2,0,22222,1377359);ecpointmultiplemodp([710108, 1324551, 1],2,0,283322,1377359);ecpointmultiplemodp([1075576, 1307157, 1],2,0,141661,1377359);ecfactor:=proc(m)
local k, kbinary, b0, b1, b2, n, a, q, r;
for a from 1 to 40 do
k:=(a+10)!;
b2:=0;
b1:=a;
b0:=-a;
if 1<gcd(-4*a^3-27*a^2,m) then
return(gcd(-4*a^3-27*a^2,m)); end if;
q:=[1,1,1];
r:=[0,1,0];
kbinary:=convert(k,base,2);
for n from 1 to nops(kbinary) do
if kbinary[n]=1 then
r:=ecaddprojective(r,q,b1,b2);
if gcd(denom(r[1]),m)>1 then
return(gcd(denom(r[1]),m)) else
r[1]:=modp(r[1],m) end if;
if gcd(denom(r[2]),m)>1 then
return(gcd(denom(r[2]),m)) else
r[2]:=modp(r[2],m) end if;
end if;
q:=ecaddprojective(q,q,b1,b2);
if gcd(denom(q[1]),m)>1 then
return(gcd(denom(q[1]),m)) else
q[1]:=modp(q[1],m) end if;
if gcd(denom(q[2]),m)>1 then
return(gcd(denom(q[2]),m)) else
q[2]:=modp(q[2],m) end if;
end do;
end do;
return([]);
end proc: ecfactor(137703491);