Package of all the EC group procedures
by Marion Scheepers
This worksheet covers the procedures for computing order of certain EC groups, and of points in such groups. It also covers the Hellman-Pohlig-Silver discrete logarithm algorithm for elliptic curve groups with an example.
> with(numtheory): with(GaussInt):
Warning, the protected name order has been redefined and unprotected
Warning, the name GIgcd has been redefined
Cubic integers library
CIConjugate
> CIConjugate:=proc(X) local a, b, answer; a:=X[1]; b:=X[2]; answer:=[a-b,-b]; end:
CINorm
> CINorm:=proc(X) local a, b, answer; a:=X[1]; b:=X[2]; answer:=(a^2-a*b+b^2); end:
CISum
> CISum:=proc(X,Y) local answer; answer:=[X[1]+Y[1],X[2]+Y[2]]; end:
CIProduct
> CIProduct:=proc(X,Y) local answer; answer:=[X[1]*Y[1] - X[2]*Y[2],X[1]*Y[2] + Y[1]*X[2] - X[2]*Y[2]]; end:
CIExponent
> CIExponent:=proc(Point,exponent) local spi, spj, spm, spbitscalar, spbitlength, spW, spPt, spQ; if (scalar = 0) then spW:=[1,0]; else spbitscalar:=convert(exponent,binary,decimal); spbitlength:=length(spbitscalar); for spj from 1 to spbitlength do spm[spbitlength-spj]:= (spbitscalar mod (10^(spbitlength-spj+1))- (spbitscalar mod (10^(spbitlength-spj))))/(10^(spbitlength-spj)); od: if (spm[0] = 0) then spW:=[1,0]; else spW:= Point; fi: spPt[0]:=Point; for spj from 1 to (spbitlength-1) do spPt[spj]:= CIProduct(spPt[spj-1],spPt[spj-1]); if (spm[spj] = 1) then spW:=CIProduct(spW,spPt[spj]) fi: od: RETURN(spW); fi: end:
CIisprime
> CIisprime:=proc(X) local a, b, c, d, result; a:=X[1]; b:=X[2]; c:= CINorm(X); if (b=0 and isprime(a) = true and a mod 3 = 2) then result:=TRUE; elif (a = 1 and b = -1) then result:=TRUE; elif (isprime(c)) then result:=TRUE; else result:=FALSE; fi: RETURN(result); end:
CIMod
> CIMod:=proc(a,b) local r, q, x, y, u, v, m, n, P, w; P:=CIProduct(a,CIConjugate(b)); x:=op(1,P)/CINorm(b); y:=op(2,P)/CINorm(b); u:=floor(x); m:=ceil(x); if (x-u > m-x) then q[1]:=-m; else q[1]:=-u; fi: w:=floor(y); n:=ceil(y); if (y-w > n-y) then q[2]:=-n; else q[2]:=-w; fi: r:=CISum(a,CIProduct(q,b)); RETURN(r); end:
CIModSum
> CIModSum:=proc(A,B,P) local R, S; R:=CISum(A,B); S:=CIMod(R,P); RETURN(S); end:
CIModProduct
> CIModProduct:=proc(A,B,P) local R, S; S:=CIProduct(A,B); R:=CIMod(S,P); RETURN(R); end:
CIModExponent
> CIModExponent:=proc(Point,exponent,Modulus) local spi, spj, spm, spbitscalar, spbitlength, spW, spPt; if (scalar = 0) then spW:=[1,0]; else spbitscalar:=convert(exponent,binary,decimal); spbitlength:=length(spbitscalar); for spj from 1 to spbitlength do spm[spbitlength-spj]:= (spbitscalar mod (10^(spbitlength-spj+1))- (spbitscalar mod (10^(spbitlength-spj))))/(10^(spbitlength-spj)); od: if (spm[0] = 0) then spW:=[1,0]; else spW:= Point; fi: spPt[0]:=Point; for spj from 1 to (spbitlength-1) do spPt[spj]:= CIModProduct(spPt[spj-1],spPt[spj-1],Modulus); if (spm[spj] = 1) then spW:= CIModProduct(spW,spPt[spj],Modulus); fi: od: fi: RETURN(spW); end:
CIRepresent
> CIRepresent:=proc(p) local a, j, q, r0, r1, r2, u, v, x, y, Form, repr, result; if (isprime(p) = false) then printf("Oops! This is not a prime number\n"); elif (p mod 3 = 2) then printf("Oops! This prime number is not of the correct form\n"); else y:= TonSh(-3,p); if ( y < p-y) then x:=y; else x:=p-y; fi: ############################# Initialize values of r-sequence r0:=p; ## r1:=p; ## r2:=x; ## ############################# for j from 1 to p do ## r0:=r1; ## a:= r2; ## r2:= r1 mod r2; ## r1:=a; ## ############################# Form:=[r0 mod 3, r1 mod 3, r2 mod 3]; if (r1^2 < p) then if (Form[1]*Form[3] > 0) then if (Form = [1,2,1] or Form = [2,1,2]) then u:= r1; v:= (r1+r2)/3; else u:= r1; v:=(r0-r1)/3; fi: else if (Form[3] > 0) then u:=r1; v:=r0/3; elif (Form[1] > 0) then u:=r1; v:=r2/3; else if (r0^2 < 3*p) then u:=r1; v:=r0/3; else u:=r1; v:=r2/3; fi: fi: fi: break; fi: od: fi: result:=[u+v,u-v]; RETURN(result); end:
Procedures for EC arithmetic over GF(P), P>3 prime
Representation of curves and points
For all these procedures we are working on elliptic curve groups over fields of the form GF(p) where p is a prime number larger than 3.
Curve
The curve will be represented by a triple [A,B,P] where:
1) P is a prime number 2) A is the coefficient of x and
3) B is the constant term in the equation
y^2 = x^3 + A*x + B mod P.
Point
A point on the curve will be represented as a triplet [e,x,y] where: 1) e = 0 means the point is the identity of the EC-group 2) e = 1 means the point is a solution to the cubic, mod P (the prime) 3) x is the x-coordinate of a solution when e = 1, arbitrary when e = 0 4) y is the y-coordinate of a solution when e = 1, arbitrary when e = 0
Tonelli-Shanks for square roots mod a prime
> ModSqrt:=proc(TSa,Prime) local TSRoot, TSb, TSi, TSj, TSm, TSmm, TSn, TSq, TSr, TSt, TSx, TSy, TSz, TSe, TSK, TSTwoPow, TSOdd; if(Prime mod 4 = 3) then TSx := TSa&^((Prime+1)/4) mod Prime; elif(legendre(TSa,Prime)=-1) then printf("%d does not have a square root mod %d \n", TSa, Prime); else TSe:=0; TSq:=Prime-1; for TSj from 1 to (Prime-1) do if(TSq mod 2 = 0) then TSe:= TSe+1; TSq:=TSq/2; else break; fi: od: for TSi from 1 to 100 do TSn:=TSi; if(legendre(TSn,Prime) = -1) then break; fi: od: TSz:=TSn&^TSq mod Prime; TSy:=TSz; TSr:=TSe; TSx:=TSa&^((TSq-1)/2) mod Prime; TSb:=TSa*(TSx^2) mod Prime; TSx:=TSa*TSx mod Prime; for TSi from 1 to TSe do if (TSb = 1) then break; else for TSm from 1 to infinity do TSt:=TSb&^(2^TSm) mod Prime; if (TSt = 1) then TSmm:= TSm; break; fi: od: TSt:=TSy&^(2^(TSr-TSmm-1)) mod Prime; TSy:=TSt&^2 mod Prime; TSr:=TSmm; TSx:=TSx*TSt mod Prime; ### WARNING: `Prime` might conflict with Maple's meaning of that name TSb:=TSb*TSy mod Prime; fi: od: fi: TSRoot:=TSx; end:
Search
> Search:=proc(lowerbound,upperbound,Curve) local sj, sx, sy, sResult, sA, sB, sP; sA:=Curve[1]; sB:=Curve[2]; sP:=Curve[3]; if (4*sA^3+27*sB^2 mod sP = 0) then printf("This is not an elliptic curve. \n"); else for sj from lowerbound to upperbound do if ( legendre(sj^3+sA*sj+sB,sP)=1 ) then sx:=(sj^3+sA*sj+sB) mod sP; sy:= ModSqrt(sx,sP); sResult:=[1,sj,sy]; break; fi: od: fi: RETURN(sResult); end:
Weil's Theorems on order of EC group
Modular Exponentiation for Gaussian Integers
> GIModExp:=proc(Input,Exp,Mod) local lpj, lgm, powpt, expbitscalar, expbitlength, lcarry, lcarry1, lbuild, lbuild1; if(Exp=0) then powpt:=1; else ### Convert exponent to binary ### expbitscalar:=convert(Exp,binary,decimal); expbitlength:=length(expbitscalar); ### Then reverse the bits from least significant to most ### for lpj from 1 to expbitlength do lgm[expbitlength-(lpj-1)]:= (expbitscalar mod (10^(expbitlength-lpj+1))- (expbitscalar mod (10^(expbitlength-lpj))))/(10^(expbitlength-lpj)); od: ### Now set up the iterated computation. ### lcarry:=Input; if(lgm[1]=0) then lbuild:=1; else lbuild:=lcarry; fi: for lpj from 2 to expbitlength do lcarry1:=GIrem(lcarry^2, Mod); lcarry:=lcarry1; if(lgm[lpj]=0) then lbuild1:= lbuild; else lbuild1:=GIrem( lcarry*lbuild, Mod); fi: lbuild:=lbuild1; od: powpt:=lbuild; fi: end:
Quartic Legendre symbol for Gaussian primes p with Norm(p) mod 4 = 1 and prime GIQLegendre(n,p);...
> GIQLegendre:=proc(n,p) local Norm, e, L: #### First, get the norm of p #### ### WARNING: `Norm` might conflict with Maple's meaning of that name Norm:=GInorm(p); #### Check if it is prime, and that p mod 4 = 1 #### if ( Norm mod 4 = 1 and isprime(Norm)) then ### WARNING: `Norm` might conflict with Maple's meaning of that name e:= (Norm - 1)/4; L:=GIModExp(n,e,p); else printf("Oops, the modulus is not a Gaussian prime!"); fi: end:
Fermat's great theorem
> FermatGreat:=proc(N) local A, B, a, k, n, X, Z; if (isprime(N)) then if (N mod 4 = 1) then Z:=sum2sqr(N); X:=[op(1,op(1,Z)),op(2,op(1,Z))]; RETURN(X); break; else printf("%d is prime, but %d mod 4 = 3\n", N, N); fi: else printf("%d is not prime\n",N); fi: end:
Weil's Theorem for =1 mod 4 primes and Y^2 = X^3 - D*X Weil4mod1(D,n);...
It is know that if n mod 4 = 3 and n is prime, then for these curves the group order is n+1. These are examples then of super singular curves over GF(n).
> Weil4mod1:=proc(D,n) local Y, a, j, p, L, X: ## First check that n is the right kind of prime ## if (n mod 4 = 1) then Y:=FermatGreat(n); a[1]:=Y[1]+I*Y[2]; a[2]:=Y[1]-I*Y[2]; a[3]:=Y[2]+I*Y[1]; a[4]:=Y[2]-I*Y[1]; a[5]:=-a[1]; a[6]:=-a[2]; a[7]:=-a[3]; a[8]:=-a[4]; for j from 1 to 8 do if (GIrem(a[j],2+2*I) = 1) then p:=a[j]; break; fi: od: L:=GIQLegendre(D,p); X:= n+1 - conjugate(L)*p - L*conjugate(p); RETURN(X); else printf("%d is not the right type of prime\n",n); fi: end:
Weil's Theorem for =1 mod 3 primes and Y^2 = X^3 + D Weil3mod1(D,n);...
It is known that for n mod 3 = 2 and prime these curves are super singular over GF(n), and the corresponding group has order n+1.
> Weil3mod1:=proc(d,n) local D, i, j, X1, x1, X2, x2, X3, x3, X4, x4, X5, x5, X6, x6, P, Num, LS, Three, N, ANSW0, ANSW1, ANSW2, ANSW3, ANSW4; D:=[d,0]; N:=[n+1,0]; if (n mod 3 = 2) then printf("This prime is not =1 mod 3\n"); else Three :=[3,0]; X1:=CIRepresent(n); x1:=CIMod(X1,Three); X2:=[-X1[1],-X1[2]]; x2:=CIMod(X2,Three); X3:=CIProduct([0,1],X1); x3:=CIMod(X3,Three); X4:=CIProduct([0,1],X2); x4:=CIMod(X4,Three); X5:=CIProduct([0,1],X3); x5:=CIMod(X5,Three); X6:=CIProduct([0,1],X4); x6:=CIMod(X6,Three); if ((x1[1] = -1) and (x1[2] = 0)) then P:=X1; Num:=CIProduct([4,0],D); LS :=CIModExponent(Num,((CINorm(P)-1)/6),P); ANSW0:=CIConjugate(LS); ANSW1:=CIProduct(ANSW0,P); ANSW2:=CIConjugate(ANSW1); ANSW3:=CISum(ANSW1,ANSW2); ANSW4:=CISum(N,ANSW3); elif ((x2[1] = -1) and (x2[2] = 0)) then P:=X2; Num:=CIProduct([4,0],D); LS :=CIModExponent(Num,((CINorm(P)-1)/6),P); ANSW0:=CIConjugate(LS); ANSW1:=CIProduct(ANSW0,P); ANSW2:=CIConjugate(ANSW1); ANSW3:=CISum(ANSW1,ANSW2); ANSW4:=CISum(N,ANSW3); elif ((x3[1] = -1) and (x3[2] = 0)) then P:=X3; Num:=CIProduct([4,0],D); LS :=CIModExponent(Num,((CINorm(P)-1)/6),P); ANSW0:=CIConjugate(LS); ANSW1:=CIProduct(ANSW0,P); ANSW2:=CIConjugate(ANSW1); ANSW3:=CISum(ANSW1,ANSW2); ANSW4:=CISum(N,ANSW3); elif ((x4[1] = -1) and (x4[2] = 0)) then P:=X4; Num:=CIProduct([4,0],D); LS :=CIModExponent(Num,((CINorm(P)-1)/6),P); ANSW0:=CIConjugate(LS); ANSW1:=CIProduct(ANSW0,P); ANSW2:=CIConjugate(ANSW1); ANSW3:=CISum(ANSW1,ANSW2); ANSW4:=CISum(N,ANSW3); elif ((x5[1] = -1) and (x5[2] = 0)) then P:=X5; Num:=CIProduct([4,0],D); LS :=CIModExponent(Num,((CINorm(P)-1)/6),P); ANSW0:=CIConjugate(LS); ANSW1:=CIProduct(ANSW0,P); ANSW2:=CIConjugate(ANSW1); ANSW3:=CISum(ANSW1,ANSW2); ANSW4:=CISum(N,ANSW3); else # ((x6[1] = 0) and (x6[2] = -1)) then P:=X6; Num:=CIProduct([4,0],D); LS :=CIModExponent(Num,((CINorm(P)-1)/6),P); ANSW0:=CIConjugate(LS); ANSW1:=CIProduct(ANSW0,P); ANSW2:=CIProduct(LS,CIConjugate(P)); #CIConjugate(ANSW1); ANSW3:=CISum(ANSW1,ANSW2); ANSW4:=CISum(N,ANSW3); fi: RETURN(ANSW4); fi: end:
Order of a group GPOrd(Curve);...
> GPOrd:=proc(Curve) local A, B, P, Answer; A:=Curve[1]; B:=Curve[2]; P:=Curve[3]; if (A = 0) then if ( P mod 3 = 2) then Answer:=P+1; elif ( P mod 3 = 1) then Answer:=Weil3mod1(B,P)[1]; else printf("The curve is not of the right sort (1)\n"); fi: elif (B = 0) then if ( P mod 4 = 3) then Answer:=P+1; elif ( P mod 4 = 1) then Answer:=Weil4mod1(-A,P); else printf("The curve is not of the right sort (2)\n"); fi: else printf("This curve is not of the right sort (3)\n"); fi: end:
Elliptic Curve Addition
> ECPlus:=proc(Point1,Point2,Curve) local pe1, px1, py1, pe2, px2, py2, ps, px, py, pR, pS, pA, pB, pP, pe; pe1:=Point1[1]; pe2:=Point2[1]; pA:=Curve[1]; px1:=Point1[2]; px2:=Point2[2]; pB:=Curve[2]; py1:=Point1[3]; py2:=Point2[3]; pP:=Curve[3]; if (4*pA^3 + 27*pB^2 mod pP = 0) then printf("This is not an elliptic curve. \n"); elif ( pe1 = 1 and (py1^2 mod pP <> px1^3 +pA*px1 + pB mod pP)) then printf("Point 1 is not on the elliptic curve. \n"); elif ( pe2 = 1 and (py2^2 mod pP <> px2^3 +pA*px2 + pB mod pP)) then printf("Point 2 is not on the elliptic curve. \n"); else if( pe1 = 0) then pR:= Point2; elif(pe2= 0) then pR:= Point1; else if (px1 = px2 and 0 = (py1+py2) mod pP) then pe:=0; px:=0; py:=0; elif(px1= px2 and py1 = py2) then pe:=1; ps:=(3*px1^2+pA)/(2*py1) mod pP; px:=(ps^2 - (px1+px2)) mod pP; py:=ps*(px1-px) - py1 mod pP; else pe:=1; ps:=(py1-py2)/(px1-px2) mod pP; px:=(ps^2 - (px1+px2)) mod pP; py:=ps*(px1-px) -py1 mod pP; fi: pR:=[pe,px,py]; fi: fi: RETURN(pR); end:
Scalar Multiple of a point on a curve
> ScalarProduct:=proc(scalar,Point,Curve) local spi, spj, spm, spbitscalar, spbitlength, spA, spB, spP, spW, spPt, spQ, spResult; if (scalar = 0) then spW:=[0,0,0]; else spbitscalar:=convert(scalar,binary,decimal); spbitlength:=length(spbitscalar); for spj from 1 to spbitlength do spm[spbitlength-spj]:= (spbitscalar mod (10^(spbitlength-spj+1))- (spbitscalar mod (10^(spbitlength-spj))))/(10^(spbitlength-spj)); od: if (spm[0] = 0) then spW:=[0,0,0]; else spW:= Point; fi: spPt[0]:=Point; for spj from 1 to (spbitlength-1) do spPt[spj]:= ECPlus(spPt[spj-1],spPt[spj-1],Curve); od: for spi from 1 to spbitlength-1 do if (spm[spi] = 1) then spQ[spi]:=spPt[spi]; else spQ[spi]:=[0,0,0]; fi: od: for spi from 1 to spbitlength-1 do spW:= ECPlus(spW,spQ[spi],Curve); od: if (spW[1] = 0) then spW:=[0,0,0]; fi: fi: spResult:=spW; end:
Order of a point on a curve
> OrderOfPointOnCurve:=proc(Point,Curve) local ordi, ordj, ordz, ordL, ordN, ordP, GpOrd, ordPt, ordprimes, ordpowers; ############################################### ordP:=Curve[3]; if (Curve[2] = 0) then if (ordP mod 4 = 3) then GpOrd:=ordP+1; else GpOrd:=Weil4mod1(-Curve[1],ordP); fi: else if (Curve[1] = 0 and ordP mod 3 = 2) then GpOrd:=ordP +1; elif(Curve[1] = 0 and ordP mod 3 = 1) then GpOrd:=Weil3mod1(Curve[2],ordP)[1]; else printf("This group does not fit the requirements of the procedure\n"); fi: fi: ################################################ ordL:=ifactor(GpOrd); for ordi from 1 to nops(ordL) do ordprimes[ordi]:=op(1,op(1,op(ordi,ordL))); if( nops(op(ordi,ordL)) = 1) then ordpowers[ordi]:=1; else ordpowers[ordi]:=op(2,op(ordi,ordL)); fi: od: ################################################# ordz:=GpOrd; ordPt:=ScalarProduct(ordz,Point,Curve); for ordi from 1 to nops(ordL) do for ordj from 1 to ordpowers[ordi]+1 do if (ordPt[1] = 0) then if(ordj = ordpowers[ordi]+1) then break; else ordz:= ordz/ordprimes[ordi]; ordPt:=ScalarProduct(ordz,Point,Curve); fi: elif (ordj = 1) then break; else ordz:= ordz*ordprimes[ordi]; ordPt:=ScalarProduct(ordz,Point,Curve); break; fi: od: od: ordN:=ordz; end:
Hellman-Pohlig-Silver discrete log attack on EC groups
> HPSonEC:=proc(Generator,Point,Curve) local dla, dld, dli, dlj, dlk, dln, dls, dlt, dlx, dly, dlz, dlK, dlN, dlN2, dlNN, dlP, dlPt, dlX, dlY, dlZZ, dlchi, dlprimes, dlpowers, dlUnknown, dlGpOrder; ###################### First compute the order of the group ###### if (Curve[2] = 0) then if (Curve[3] mod 4 = 3) then dlGpOrder:=Curve[3]+1; else dlGpOrder:=Weil4mod1(Curve[1],Curve[3]); fi: else if (Curve[3] mod 3 = 2 and Curve[1] = 0) then dlGpOrder:=Curve[3]+1; fi: fi: ####################### Then compute the order of the generator ### dlN:= OrderOfPointOnCurve(Generator,Curve); dlN2:= OrderOfPointOnCurve(Point,Curve); if (gcd(dlN,dlN2)<dlN2) then printf("The point is not a multiple of the generator\n"); RETURN; else dlK:=ifactor(dlN); for dli from 1 to nops(dlK) do dlprimes[dli]:=op(1,op(1,op(dli,dlK))); if (nops(op(dli,dlK))=1) then dlpowers[dli]:=1; else dlpowers[dli]:=op(2,op(dli,dlK)); fi: od: for dli from 1 to nops(dlK) do dlz:= dlN/dlprimes[dli]; dlchi:=ScalarProduct(dlz,Generator,Curve); for dlj from 0 to (dlpowers[dli] -1) do if (dlj = 0) then dla[dlj]:=Point; else dlZZ:=(dlN-1)*dld[dli,dlj-1]*(dlprimes[dli]^(dlj-1)) mod dlN; dlPt:=ScalarProduct(dlZZ,Generator,Curve); dla[dlj]:=ECPlus(dla[dlj-1],dlPt,Curve); fi: dly:=dlN/(dlprimes[dli]^(dlj+1)); dlt[dlj]:=ScalarProduct(dly,dla[dlj],Curve); dls:=[0,0,0]; for dlk from 0 to (dlprimes[dli]-1) do if (dls = dlt[dlj]) then dld[dli,dlj]:= dlk; break; fi: dls:=ECPlus(dls,dlchi,Curve); od: od: #################################################################### dln[dli]:=0; for dlj from 0 to dlpowers[dli]-1 do dln[dli]:= dln[dli] + dld[dli,dlj]*(dlprimes[dli]^dlj) mod dlN; od: od: #################################################################### ############### Applying the Chinese Remainder Theorem ############# #################################################################### dlUnknown:=0; for dli from 1 to nops(dlK) do dlNN[dli]:=dlN/(dlprimes[dli]^dlpowers[dli]); dlX[dli] :=1/dlNN[dli] mod (dlprimes[dli]^dlpowers[dli]); dlUnknown:= dlUnknown+ dlNN[dli]*dln[dli]*dlX[dli] mod dlN; od: fi: end:
An example
> P:=17628652901592224077691838623064378402878020122363779179718011376257597439: A:=394857398475: B:=0: Curve:=[A,B,P]:
A point from the curve
> Q1:=Search(623457623,28374682736482,Curve);
A scalar product of a point from the curve
> ### WARNING: `Scalar` might conflict with Maple's meaning of that name Scalar:=1188229207: Q:=ScalarProduct(Scalar,Q1,Curve);
A sum of points on the curve
> QQ:=ECPlus(Q,Q1,Curve);
Order of a point on the curve
> xx:=OrderOfPointOnCurve(Q1,Curve);
Discrete Log of a point on the curve
> yy:=HPSonEC(Q1,Q,Curve);
>