/* Functions for elliptic curves over the finite field Z/pZ; could be
extended to more general finite fields one day.
Author: John Cremona (john.cremona@nottingham.ac.uk)
Last edited: 13/10/04
*/
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ A random point on an elliptic curve defined over Z/pZ
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
ellzprandompoint(ep)=
local(p,xp,yp);
p=ep.a1.mod; yp=[];
while(#yp==0,yp=ellordinate(ep,xp=Mod(random(p),p)));
[xp,yp[1]];
}
/* Example:
e=ellinit(Mod(1,101)*[0,0,1,-7,6]);
lift(vector(5,j,ellzprandompoint(e)))
returns something like
[[53, 6], [44, 6], [63, 54], [21, 5], [99, 3]]
*/
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Reduction of rational points mod p
\\ [0] -> [0]
\\ [x,y] -> [0] if val(x),val(y)<0 else [Mod(x,p),Mod(y,p)]
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
ellreducepoint(P,p)=if(P==[0],[0],if(valuation(P[1],p)<0,[0],Mod(1,p)*P));
}
/* Examples:
ellreducepoint([0],5)
\\ returns [0]
ellreducepoint([1/25,1/125],5)
\\ returns [0]
ellreducepoint([1/25,1/125],3)
\\ returns [Mod(1, 3), Mod(2, 3)]
*/
\\ Extension to built-in data type: the return vector from ellinit()
\\ on input of type t_INTMOD has length 19 but entries 14-19 are unused
\\ and set to 0. We use these as follows (provisional):
\\
\\ E[14] : the group order n = #E(Fp)
\\ E[15] : factorization of the group order n
\\ E[16] : the group structure: [n] if cyclic, or
\\ [n1,n2] with n=n1*n2 and n2|n1
\\ E[17] : the group generators: [P] is cyclic, or
\\ [P1,P2] with order(Pi)=ni
\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Input can be either a vector of INT_MODs (in which case the
\\ parameter p is ignored so can be omitted) or a vector of INTs, in
\\ which case we need to also give the prime p. No check is made that p
\\ is prime. Singular data will cause a run-time error
\\
\\ Output is an extended ellinit structure as explained above
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
ellzpinit(ai,p)=
local(ep,gp);
ai=vecextract(ai,[1,2,3,4,5]);
if(type(ai[1])=="t_INTMOD",
ep=ellinit(ai);
p=ai[1].mod;
ap=ellap(ellinit(lift(ep)),p),
if(p==0,
print("Error in ellzpinit: p=0");
return(0),
ep=ellinit(Mod(1,p)*ai);
ap=ellap(ellinit(ai,1),p)));
ep[15]=factor(ep[14]=1+p-ap);
gp=ellzpgroup(ep);
if(#gp<2,print("Problem finding group structure"),
ep[16]=gp[1];
ep[17]=gp[2]);
ep;
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Extracting the group order / factored group order:
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
ellzp.grouporder=ellzp[14];
ellzp.factoredgrouporder=ellzp[15];
/* Example:
e=ellzpinit([0,0,1,-7,6],101);
e.grouporder
e.factoredgrouporder
returns
98
[2 1]
[7 2]
*/
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Extracting the group structure:
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
ellzp.groupstr=[ellzp[16],ellzp[17]];
ellzp.iscyclic=(#ellzp[16])==1;
ellzp.isotype=ellzp[16];
ellzp.generators=ellzp[17];
/* Example 1 :
e=ellzpinit([0,0,1,-7,6],101);
e.groupstr
e.iscyclic
e.isotype
e.generators
returns (possibly with a different generator)
[[98], [[Mod(38, 101), Mod(8, 101)]]]
1
[98]
[[Mod(38, 101), Mod(8, 101)]]
/* Example 2 :
e=ellzpinit([0,0,1,-7,6],103);
e.groupstr
e.iscyclic
e.isotype
e.generators
returns (possibly with a different first generator)
[[52, 2], [[Mod(88, 103), Mod(25, 103)], [Mod(83, 103), Mod(51, 103)]]]
0
[52, 2]
[[Mod(88, 103), Mod(25, 103)], [Mod(83, 103), Mod(51, 103)]]
*/
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Order of a point (using the factored group order):
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
ellzppointorder(ellzp,t)=
local(ans,fm,p,m,s);
\\print("Finding order of ",t);
\\print("Group order = ",ellzp.grouporder);
ans=1; fm=(ellzp.factoredgrouporder)~;
for(j=1,#fm,p=fm[1,j];
\\print("p=",p);
m=factorback(vecextract(fm,concat("^",Str(j)))~);
\\print("m=",m);
s=ellpow(ellzp,t,m);
\\print(s);
while(s!=[0],s=ellpow(ellzp,s,p);ans*=p));
ans;
}
/* Example:
ellzp=ellzpinit([0,0,1,-7,6],101);
vector(10,j,ellzppointorder(ellzp,ellzprandompoint(ellzp)))
returns something like
[7, 49, 98, 49, 98, 98, 14, 98, 98, 49]
which shows (assuming that 98 does appear in the list, as it does with
probability more than 99.6%) that the group is cyclic. We can then
get a generator:
P=[0];while(ellzppointorder(ellzp,P)<98,P=ellzprandompoint(ellzp));P
returns (perhaps)
Mod(3, 101), Mod(3, 101)]
*/
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Group structure:
\\ returns either
\\ [[n],[P]] (if cyclic of order n) where P is a generator
\\ or
\\ [[n1,n2],[P1,P2]] is (Z/n1)x(Z/n2) with n2|n1
\\ and P1, P2 generators of order n1,n2
\\
\\ Internal function used in ellzpinit, need/should not be called by user
\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
debug_group=0;
\\ Utility function: given positive integers m,n returns [l,m',n']
\\ where gcd(m',n')=1, lcm(m',n')=m'n'=lcm(m,n), m'|m, n'|n:
\\
\\ (Nothing to do with elliptic curves, may be useful elsewhere)
{
tidy_lcm(m,n)=
local(g,l,m0,n0);m0=m;n0=n;
g=gcd(m,n);
l=m*n/g; \\ = lcm(m,n)
g=gcd(m,n/g); \\ divisible by primes dividing n to a higher power than m
while(g!=1, m/=g; g=gcd(m,g));
n=l/m;
if((l==m*n)&&
(gcd(m,n)==1)&&
(m0%m==0)&&
(n0%n==0),,
print("error in tidy_lcm(",m0,",",n0,")"));
[l,m,n];
}
\\ Utility function: given a positive integer n returns the
\\ largest d such that d^2|n.
\\
\\ (Nothing to do with elliptic curves, may be useful elsewhere)
{
maxsquarefreediv(n)=
local(fn);
fn=factor(n);
fn[,2]\=2;
factorback(fn)
}
\\ Utility function: given P1 of order n1 and P2 of order n2,
\\ construct a point Q of order n3=lcm(n1,n2); returns [Q,n3]
\\
\\ (algorithm may be useful elsewhere)
{
ellzpmerge1(ellzp,P1,n1,P2)=
local(n2,abc,m,m1,m2);
if(ellpow(ellzp,P2,n1)==[0],return([P1,n1]));
n2=ellzppointorder(ellzp,P2);
if(n2%n1==0,return([P2,n2]));
abc=tidy_lcm(n1,n2);
m=abc[1];
m1=n1/abc[2];
m2=n2/abc[3];
[elladd(ellzp,ellpow(ellzp,P1,m1),ellpow(ellzp,P2,m2)), m];
}
\\ Utility function: given independent points P1,P2 of order n1,n2,
\\ and a further point Q, replace P2 by a point still independent of
\\ P1 whose order is the lcm of n2 and the order of Q mod .
\\ n2target is n/n1 (n=group order) so n2|n2target and we are aiming
\\ to achieve n2=n2target
\\
\\ All of P1,P2,n1,n2 may be changed; as we have no writeable
\\ variable parameters in gp we return [P1,n1,P2,n2]
\\
\\ (algorithm may be useful elsewhere)
{
ellzpmerge2(ellzp,P1,n1,P2,n2,n2target,Q)=
local(Q1,Q2,P1n1,a,w,m,l,ordQ,lmn);
if(debug_group>2,print("merge2 called with ",[P1,n1,P2,n2,Q,n2target]));
if(n2==0,print("Error: n2=0 at place A");1/0);
Q1 = ellpow(ellzp,Q,n2);
if(Q1==[0], if(debug_group>1,print("Order(Q) divides n2"));
return([P1,n1,P2,n2])
);
Q2 = ellpow(ellzp,Q1,n1/n2);
if(Q2!=[0],
\\ P1 needs updating and we discard P2
if(debug_group>1, print("Order(Q) does not divide n1, updating P1"));
P1n1=ellzpmerge1(ellzp,P1,n1,Q);
P1=P1n1[1]; n1=P1n1[2];
if(debug_group>1, print("New P1 has order ",n1));
if(n2>1, P2 = [0]; n2=1);
return([P1,n1,P2,n2])
);
\\ Now we find a multiple a*P1 such that Q-a*P1 is killed by
\\ n2target so we can apply the Weil Pairing of order n2target
Q1 = ellpow(ellzp,Q,n2target);
Q2 = ellpow(ellzp,P1,n2target); \\ has exact order n1/n2target
a = bg_algorithm(ellzp,Q2,Q1,0,n1/n2target);
if(debug_group,
print("Dlog of ",Q1," w.r.t. ",Q2," (order ",n1/n2target,") is ",a);
print("Check: a*Q2-Q1 = ",ellsub(ellzp,ellpow(ellzp,Q2,a),Q1))
\\a*Q2-Q1
);
Q = elladd(ellzp,Q,ellpow(ellzp,P1,-a)); \\ = Q-a*P1
if(Q==[0],
if(debug_group, print("Q-a*P1 = ",Q,", no use"));
return([P1,n1,P2,n2])
);
if(debug_group,
print("Replacing Q by Q-a*P1 = ",Q," where a = ", a);
if(ellpow(ellzp,Q,n2target)==[0],
print1("order divides n2target, OK; computing Weil pairing of ");
print("order " ,n2target),
print("order does not divide n2target -- WRONG!")
));
\\ compute Weil pairing of P1 and Q, which is an n2target'th root of unity:
Q1 = ellpow(ellzp,P1,(n1/n2target));
if(debug_group,
print("order((n1/n2target)*P1) = ",Q1," is ",ellzppointorder(ellzp,Q1));
print("order(Q) = ",Q," is ",ellzppointorder(ellzp,Q))
);
if(debug_group,print("Computing Weil pairing of those two points..."));
w = ellweilpairing(ellzp,Q1,Q,n2target);
m = znorder(w);
if(debug_group, print("w = ", w," of order ", m ));
\\ Compare this with n2 to see if we have gained:
l = lcm(n2,m);
if(l==n2, return([P1,n1,P2,n2]) ); \\ no gain
ordQ = ellzppointorder(ellzp,Q);
if(debug_group, print("ordQ = ", ordQ, "; ordQ/m = ",ordQ/m));
Q1 = ellpow(ellzp,Q,ordQ/m); \\ of order m
if(l==m, \\ replace P2, n2 by Q1, m
P2 = Q1;
n2 = m;
return([P1,n1,P2,n2])
);
\\ Now P2,Q1 have orders n2,m & both are independent of P1,
\\ so we combine them to get a point of order l=lcm(n2,m)
\\ still independent of P1
lmn = tidy_lcm(n2,m);
n2=lmn[2]; m=lmn[3]; if(n2==0,print("Error: n2=0 at place B");1/0);
P2 = elladd(ellzp,ellpow(ellzp,P2,l/n2),ellpow(ellzp,Q,l/m));
n2 = l; if(n2==0,print("Error: n2=0 at place C");1/0);
if(debug_group,
print("Changed P2 = ",P2,":\t order(P2) = ",n2);
if(ellzppointorder(ellzp,P2)!=n2, print("that's wrong!"))
);
return([P1,n1,P2,n2]);
}
{
ellzpgroup(ellzp)=
local(N,p,P1,n1,n2,n,Q,Pm,n2target);
p=ellzp.a1.mod;
if(debug_group,
print("Curve is ", lift(vecextract(ellzp,[1,2,3,4,5])), " mod ",p);
print("Factorization of p-1 is ",factor(p-1)));
fN=ellzp.factoredgrouporder; \\ factorization of the order
N=ellzp.grouporder; \\ the order
if(debug_group,print("N = ", N, " is the group order"));
if(N==1,return([[],[]]));
fN[,2]\=2;
n2max = gcd(factorback(fN),p-1);
iscyclic=(n2max==1); \\ 0 means we don't yet know
if(debug_group,
print("A priori maximal possible n2 = ", n2max);
if(iscyclic,print("Group must be cyclic"))
);
\\ Start finding random points...
P1 = ellzprandompoint(ellzp);
n1 = ellzppointorder(ellzp,P1);
if(debug_group,print("P1 = ", P1, ", order ",n1));
if(n1==N, return([[N],[P1]])); \\ we hit the jackpot, quick exit
n2max = gcd(n2max,N/n1);
if(debug_group, print("maximal possible n2 = ", n2max));
iscyclic=(n2max==1); \\ 0 means we don't yet know
if(debug_group, if(iscyclic,print("Group must be cyclic")));
\\ use more random points to increase the order; after a few we are
\\ highly likely to be holding a point of maximal order (hence a
\\ generator when E is cyclic, which happens often)
n=0;
while((n<10)||((n1*gcd(n1,p-1))%N!=0),n+=1;
Q = ellzprandompoint(ellzp);
Pm=ellzpmerge1(ellzp,P1,n1,Q);
P1=Pm[1]; n1=Pm[2];
if(debug_group,print("P1 = ", P1, ", order ",n1));
if(n1==N, return([[N],[P1]])); \\ cyclic, generator P1, exit
n2max = gcd(n2max,N/n1);
if(debug_group, print("maximal possible n2 = ", n2max));
iscyclic=(n2max==1); \\ 0 means we don't yet know
if(debug_group, if(iscyclic,print("Group must be cyclic")));
);
if(n1==N, return([[N],[P1]])); \\ cyclic, generator P1, exit
\\ See whether [n1,N/n1] is possible group structure...
if(N%n1,print("problem 1");return([])); \\contradicts Lagrange...
n2=N/n1;
if(n1%n2,print("problem 2");return([])); \\ should not happen
\\ since once we get here, (n1*gcd(n1,p-1))%N==0
if((p-1)%n2,print("problem 3");return([])); \\ should not happen
if(debug_group,
print1("ellzpgroup found probable group structure, now finding second ");
print("generator of order ",n2)
);
P2=[0];
n2target=n2;
n2=1; \\ holds order(P2) always
while(n2