Elliptic Curve Algorithms in Mathematica

René Schoof's Algorithm for #E((Subscript[𝔽, p]))

John J. McGee

Mathematics and Statistics Department
Radford University, Virginia

A bad day for Alice and Bob [8]

The Discrete Logarithm Problem

Solve for k∈ℤ

	(2^k)=2048					k=ln(2048)/ln(2)

	(17^k)≡16958 (mod  17499)			Index Calculus O((e^(Sqrt[2 ln p ln ln p])))  [9]
	
	k x = y  (in an Abelian Group)
	
Basis for many cryptographic applications  (Massey-Omura, ElGamal) [9]

	Structure of the group determines the difficulty of finding k

Elliptic Curves

An elliptic curve is the solution set of a nonsingular cubic polynomial equation in two unknowns

E(F)={(x,y)∈F×F|f(x,y)=0}    F a field,  f(x,y)∈F[x,y]

The general equation of a cubic in two variables is given by

a (x^3)+b (y^3)+c (x^2)y+d x (y^2)+ e (x^2)+f (y^2)+g x y+h x+i y+j = 0

When char(F)∉{2,3} by a change of variables we can convert to Weierstrass form

(y^2)=(x^3)+A x+B

Elliptic Curves in (ℝ^2) - we can draw the graph of E(ℝ)

Code

PlotEc[a_,b_]:=Module[{p1,p2,p3,t1},
p1=Plot[(x^3)+a x+b,{x,-4,5},PlotStyle->{ColorData[1,1],Thickness[0.003]}];
p2=Plot[Sqrt[(x^3)+a x+b],{x,-5,9},PlotStyle->{ColorData[1,2],Thickness[0.005]}];
p3=Plot[-Sqrt[(x^3)+a x+b],{x,-5,9},PlotStyle->{ColorData[1,2],Thickness[0.005]}];
Return[{p1,p2,p3}];
];
ManupulateEc[]:=Manipulate[
Show[PlotEc[a,b],ImageSize->Large,PlotRange->{{-5,9},{-30,30}},AxesLabel->{"x","y"},PlotLabel->Style[TraditionalForm[(y^2) == (x^(3 ))+  a x +b],20]],
Row[Map[Style[#,FontSize->18]&,{"Weierstrass Normal Form Elliptic Curve: ",TraditionalForm[ (y^2)== (x^3)+A x +B]}]],
{{a,0,"A"},-25,25},{{b,0,"B"},-25,25},Deployed->True];

Interactive Example

ManupulateEc[]

Elliptic Curve "Chord and Tangent" Point Addition

Every line between two points on the curve intersects a third point (Bezout)

If we define A⊕B as the reflection of the third point of intersection about the axis of symmetry
then the points form an algebraic group under ⊕

Code

This function performs the elliptic curve group addition (Subscript[P, 1])+(Subscript[P, 2]) over the finite field ℝ.  The parameters are {a,b}, (Subscript[P, 1])={(Subscript[x, 1]),(Subscript[y, 1])},(Subscript[P, 2])={(Subscript[x, 2]),(Subscript[y, 2])},  The function accepts and returns {∞,∞} for the group identity 𝒪.  It returns {} if either point does not lie on the elliptic curve.

EcAdd[{a_,b_},p1_,p2_]:=Module[{m,x1,y1,x2,y2,x3,y3,w},
{x1,y1}=p1; {x2,y2}=p2;
(* Handle identity cases *)
If[x1==∞,Return[p2]];
(* (Subscript[p, 1]) + (-(Subscript[p, 1])) = ∞ *)
If[x1==x2&& y1+y2==0,Return[{∞,∞}]];
(* If we are doubling a point *)
If[p1==p2,
(* Check for vertical tangent *)
If[y1==0,Return[{∞,∞}]];
(* Compute the slope of the tangent *)
m=(3 (x1^2)+a)/(2 y1);
,
(* else compute the slope of the chord *)
m=(y2-y1)/(x2-x1)
];
x3=(m^2)-x1-x2;
y3=m(x1-x3)-y1;
Return[{x3,y3}];
];

EcY[{a_,b_},x_]:=Sqrt[(x^3)+a x+b];

PlotEcAdd[pt1_,pt2_]:=Module[{pt3,pt4,g},
pt4=EcAdd[{a,b},pt1,pt2];
pt3={pt4[[1]],-pt4[[2]]};
g={ContourPlot[(y^2)==(x^3)+a x+b,{x,-5,5},{y,-4,4},Axes->True,PerformanceGoal->"Quality",PlotPoints->{15,20}],Graphics[{PointSize[0.02],Point[{pt1,pt2,pt3}],ColorData[1,2],Point[pt4],Line[{pt1,pt3}],ColorData[1,3],Line[{pt3,pt4}]}]};
Return[g];
];

a=-3;b=-1;
x1=-1.2;

pt1={x1,EcY[{a,b},x1]}

ManupulateEcAdd[]:=Manipulate[Show[PlotEcAdd[pt1,{pt[[1]],EcY[{a,b},pt[[1]]]}],ImageSize->440],Style[Row[{"  Elliptic Curve Point Addition  ",TraditionalForm[(y^2)==(x^3)+a x +b]}],FontSize->18],{pt,{-1.2,0},{-0.38,0},Locator,Appearance->Style["|",20,Red]}]

Example

ManupulateEcAdd[]

Elliptic Curve Group Law

First we compute the slope of the chord (or tangent)

λ={({{ ((Subscript[y, 2])-(Subscript[y, 1]))((((Subscript[x, 2])-(Subscript[x, 1])))^(-1)), (Subscript[x, 1])≠(Subscript[x, 2]) }, {(3((Subscript[x, 1])^2)+A)(((2 (Subscript[y, 1])))^(-1)),     otherwise }})

We then intersect this line with the elliptic curve to find ((Subscript[x, 3]),-(Subscript[y, 3]))

(Subscript[x, 3])=(λ^2)-(Subscript[x, 1])-(Subscript[x, 2]),       (Subscript[y, 3])=λ((Subscript[x, 1])-(Subscript[x, 3]))-(Subscript[y, 1])

If the line is vertical ((Subscript[x, 1])=(Subscript[x, 2]),(Subscript[y, 1])=-(Subscript[y, 2])) we define the intersection point to be 𝒪
	the point at infinity, derived from examining the curve in projective space.
	
𝒪 acts as the identity in the group.

Proof of group properties depends only on the algebraic form of the equations, valid for any field.

Order of the group E((Subscript[𝔽, p]))

When F is a finite field, size of the solution set is finite 

By Lagrange's theorem each point has order dividing #E((Subscript[𝔽, p]))=N
	
	 ({{i.e   k P =}, { }}) ({{P⊕P⊕...⊕P}, {k times}})({{=𝒪   for some integer k|N}, { }}) 

Counting the number of points in E((Subscript[𝔽, p])) or E((Subscript[𝔽, (p^k)])) is important in many E.C. applications.

In 1985 René Schoof published an efficient algorithm ( O((log^8)(p)) [5]  ) to accomplish this.

"A four-year-old child could understand that.  Run out and find me a four-year-old child, 
I can't make head or tail out of it." - Groucho Marx (Duck Soup-1933)

The Elliptic Curve Discrete Log Problem - ECDLP

We compute an integer multiple of a point P∈E(K) as

	 ({{ k P =}, { }}) ({{P⊕P⊕...⊕P}, {k times}})  
Given two points P,Q, where Q=k P, determine k.

How hard is the ECDLP?  It depends on the group structure.

For an elliptic curve over (Subscript[𝔽, p]) with p=564538252084441556247016902735257 ∼ (10^33)

Mississauga, Ontario – (November 6, 2002) – Certicom Corporation (TSX: CIC), a leading provider of wireless security solutions, today announced that Chris Monico and his team of mathematicians from Notre Dame have successfully solved the Certicom ECCp-109 Challenge. The challenge was solved utilizing a massive amount of computing power including 10,000 computers (mostly PCs) running 24 hours a day for 549 days.

“Without a doubt, this is one of the world’s largest single math computations ever completed,” said Chris Monico, a post doctorate researcher in mathematics at Notre Dame. “Cracking this code was a massive undertaking in terms of time, resources and CPU cycles. It has been the most difficult elliptic curve discrete algorithm problem ever computed.”  [10]

Hasse's Theorem

This theorem bounds the number of points in E((Subscript[𝔽, p])) very close to p

#E((Subscript[𝔽, p]))=p+1-t, with |t|≤2(Sqrt[p]).

Example:  for p=29996224275907

	29996213322146≤#E((Subscript[𝔽, p]))≤29996235229670

The Frobenius Endomorphism

Maps points with coordinates in the algebraic closure of (Subscript[𝔽, p])

(Subscript[ϕ, p]):E((Subscript[(OverBar[𝔽]), p]))->E((Subscript[(OverBar[𝔽]), p]))  with (Subscript[ϕ, p])((x,y))=((x^p),(y^p))

It has the following property crucial to Schoof's algorithm

(Subscript[ϕ, p]^2)-t (Subscript[ϕ, p])+p=0	∀ P∈E((Subscript[(OverBar[𝔽]), p]))          characteristic equation of the Frobenius

In terms of a point (x,y) on the curve this means

((x^(p^2)),(y^(p^2)))⊕ t ((x^p),(y^p))⊕ p( x,y) = 𝒪

The Division Polynomials

The division polynomial (Subscript[ψ, k])(x,y) goes to zero precisely for all points
	(x,y)∈E((Subscript[𝔽, p])) with k P=𝒪,  i.e. points with order dividing k
	
They are defined recursively as follows in terms of the E.C. parameters a,b

(Subscript[ψ, 0])=0,(Subscript[ψ, 1])=1,(Subscript[ψ, 2])=2 y,
		(Subscript[ψ, 3])=3 (x^4)+6 a (x^2)+12 b x-(a^2),
		(Subscript[ψ, 4])=4 y((x^6)+5 a (x^4)+20 b (x^3)-5 (a^2) (x^2)-4 a b x-8 (b^2)-(a^3))
		(Subscript[ψ, 2n+1])=(Subscript[ψ, n])((Subscript[ψ, n+2])(Subscript[ψ, n-1]^2)-(Subscript[ψ, n-2])(Subscript[ψ, n+1]^2))   n∈ℤ,n>2
		(Subscript[ψ, 2n+1])=(Subscript[ψ, n+2])(Subscript[ψ, n]^3)-(Subscript[ψ, n+1]^3)(Subscript[ψ, n-1])                n∈ℤ,n>1

Example   (y^2)=(x^3)+5098808 x +14586785 over (Subscript[𝔽, 15518471])

(Subscript[ψ, 10])=(8 (x^30)+13675911 (x^28)+10352566 (x^27)+11965774 (x^26)+6804971 (x^25)+9232487 (x^24)+6127877 (x^23)+
8693774 (x^22)+3909418 (x^21)+2723979 (x^20)+8608469 (x^19)+8304216 (x^18)+1548915 (x^17)+1989531 (x^16)+2294803 (x^15)+8149372 (x^14)+8007390 (x^13)+12194582 (x^12)+599787 (x^11)+13749587 (x^10)+11196271 (x^9)+9692489 (x^8)+819246 (x^7)+3960226 (x^6)+10786594 (x^5)+3956883 (x^4)+8684493 (x^3)+7657384 (x^2)+14848542 x+3814154) y

Point Multiplication Formula

The division polynomials can be used to derive a formula for n (x,y) as:

n P=(x-(((Subscript[ψ, n-1])(Subscript[ψ, n+1]))/(Subscript[ψ, n]^2)),(((Subscript[ψ, n+2])(Subscript[ψ, n-1]^2)-(Subscript[ψ, n-2])(Subscript[ψ, n+1]^2))/(4 y (Subscript[ψ, n]^3))))

This allows us to write the Frobenius formula explicitly as a rational equation

((x^(p^2)),(y^(p^2)))⊕ t ((x^p),(y^p))⊕ p( x,y) = 𝒪

The only unknown parameter in this equation is t.

Schoof's Algorithm

The essence of Schoof's algorithm is to solve the Frobenius equation for a set of small primes ℓ, whose product is at least 4(Sqrt[p]).  Focusing on points only in E[ℓ], we arrive at

((x^(p^2)),(y^(p^2)))⊕ k( x,y) ≡ - t ((x^p),(y^p)) (mod (Subscript[ψ, ℓ]), p);  k≡p(mod ℓ)

The order of the polynomials is kept small by performing the polynomial algebra
modulo the (ℓ^th) division polynomial, with coefficients computed modulo p

In this way we can determine t (mod ℓ) for a set of small primes and use the Chinese Remainder
Theorem to recover t (mod M) where M=(∏_(i=1)^m)(Subscript[ℓ, i]).   Since M>4(Sqrt[p]) we determine
the unique t that falls in Hasse's bound, thus we obtain #E((Subscript[𝔽, p]))

#E((Subscript[𝔽, p]))=p+1-t, with |t|≤2(Sqrt[p]).

Equation 19y

There are several special cases and we handle the x and y components of the Frobenius equation separately.   Further, we can work with polynomials in x only, by handling even and odd powers distinctly and replacing (y^(2k)) with ((((x^3)+a x+b))^k) .  One particularly complex case is the y coordinate test of Schoof (19).

(Subscript[p, 19y])(x,y)=4(Subscript[f, τ]^(3 p))(y^p) (((2 (x^(p^2))+x) α (β^2)-(β^3) (y^(p^2))-(α^3)) (Subscript[f, k]^2)-α (β^2)(Subscript[f, k-1]) (Subscript[f, k+1]) )
		 -(β^3) (Subscript[f, k]^2) ((((Subscript[f, τ-1]^2) (Subscript[f, τ+2])-(Subscript[f, τ-2]) (Subscript[f, τ+1]^2)))^p)
where
	α= (Subscript[ψ, k+2])(Subscript[ψ, k-1]^2)-(Subscript[ψ, k-2])(Subscript[ψ, k+1]^2)-4(y^((p^2)+1))(Subscript[ψ, k]^3)

	β=4y (Subscript[ψ, k]) ((x-(x^(p^2)))(Subscript[ψ, k]^2)-(Subscript[ψ, k-1])(Subscript[ψ, k+1]))

Example for(y^2)=(x^3)+5098808 x +14586785 over (Subscript[𝔽, 15518471])

(Subscript[f, 19Y])=15100019 (x^23)+8577465 (x^22)+12352612 (x^21)+4217845 (x^20)+1801075 (x^19)+13143806 (x^18)+6382365 (x^17)+
	11576028 (x^16)+2549073 (x^15)+5504927 (x^14)+2759619 (x^13)+10279388 (x^12)+4113150 (x^11)+11229691 (x^10) +
	1265164 (x^9)+15017167 (x^8)+6700100 (x^7)+9033889 (x^6)+1079185 (x^5)+3338292 (x^4)+10079001 (x^3)+
	3606043 (x^2)+1507954 x+7919516

Schoof's Algorithm Summary

1. If gcd((x^3)+a x+b, (x^p)-x)=1 then t≡0 (mod 2), else t≡1 (mod 2)
2. Create a set of small primes S={(Subscript[l, i])} such that (∏_(i=1)^L)(Subscript[l, i])>4 (Sqrt[p]).
3. Compute the first (Subscript[l, L])+2 division polynomials (Subscript[ψ, k]).
4. For each l∈S, compute k≡p ( mod l)
5. 	If gcd((Subscript[p, 16]),(Subscript[f, l]))≠1 then there exists P∈E[l] such that  (Subscript[ϕ, l]^2)(P)=±k P.
6.		If k is not a quadratic residue mod l, then t≡0(mod l) else
7.		Compute w such that (w^2)≡k (mod l)
8.		If gcd((Subscript[p, 17]),(Subscript[f, l]))=1 then  t≡0(mod l), else
9.		If gcd((Subscript[p, 18]),(Subscript[f, l]))≠1 then  t≡2w(mod l), else  t≡-2w(mod l). 
10.	else we are in case two
11.	    For each τ≤(l+1)/2
12.		If gcd((Subscript[p, 19X]),(Subscript[f, l]))≠1 then ((Subscript[ϕ, p])^2)(P)+k P≡±τ (Subscript[ϕ, p])(P) (mod l)
13.		    for some point in E[l] so we test
14.			If gcd((Subscript[p, 19 Y]),(Subscript[f, l]))≠1 then  t≡τ(mod l) else t≡-τ (mod l) 
15.	    Next τ
16. Next l
17. At this point we have computed t (mod (Subscript[l, i])) for all (Subscript[l, i])∈S, 
18.	so we use the Chinese Remainder Theorem to compute
19.	     T≡t (mod N) where N=(∏_(i=1)^L)(Subscript[l, i]).
20. If T is within Hasse's bounds then t=T, else t≡-T(mod N) and
21.	#E((Subscript[𝔽, p]))=p+1-t.

Results

Table 1 summarizes the results of running Schoof's algorithms on several curves.  The first column describes the curve (y^2)=(x^3)+A x+B over (Subscript[𝔽, p]) by giving A,B,p.  The second column show the results of computing (Subscript[τ, i])≡t (mod (Subscript[l, i])) as pairs {(Subscript[τ, i]),(Subscript[l, i])}.  The third column gives the unique t within Hasses' bounds that satisfies the congruences of column two.   #E((Subscript[𝔽, p])) is then computed as p+1-t. 
 
                                                Table 1 - Results of Schoof's Algorithm
   ({{A,B,p, t(mod (Subscript[l, i])), t, #E((Subscript[𝔽, p])), Correct}, {13,215,229, {0,2},{0,5},{4,7}, -10, 240, yes}, {106,166,197, {0,2},{2,3},{0,5},{1,11}, -10, 208, yes}, {31,16,137, {1,2},{0,3},{4,5},{2,7}, 9, 129, yes}, {503,367,523, {1,2},{0,5},{6,7},{7,11}, -15, 539, yes}, {1333,1129,3571, {1,2},{2,11},{0,13}, -13, 3559, yes}})

A large example

Largest example run to date  #E((Subscript[𝔽, p])): (y^2)=(x^3)+A x+B  with  p∼(10^30)

p=2184972043155134786254507044209   A=812066892826496630200758496181, B=393895671555005480312164618491 

"Hasse's bounds "2184972043155131829924320373490" ≤ #E(\!\(\*SubscriptBox[\(𝔽\), \(p\)]\)) ≤ "2184972043155137742584693714930

τ (mod (Subscript[ℓ, i]))={{2,0},{3,0},{5,2},{7,1},{11,10},{13,12},{17,10},{23,8},{29,17},{31,18},{37,4},{41,2},{43,18},{47,17}}

τ=-2235306593865918  via Chinese Remainder Theorem

#E((Subscript[𝔽, p]))=2184972043155137021561100910128	 12.48 hours runtime

FactorInteger[2184972043155137021561100910128]

Representative Mathematica Functions for Elliptic Curves over (Subscript[𝔽, p])

SquareRootModPShanksTonelli[ a,p ]	 	PolynomialPowerMod[f,n,{g,k}]

ComputeDivisionPolynomials[ {a,b,p},k ]		EcAddMod[ {a,b,p},(Subscript[P, 1]),(Subscript[P, 2]) ]

EcPowerMod[ {a,b,p},Q,k ]				EcPointOrderBabyGiant[ {a,b,p},Q ]

PollardRhoEcMod[ec,PX1,QX1,{a0,b0},p]		LenstraECMFactor[ n,t ]

EcGroupOrder[ {a,b,p} ]				EcGroupOrderSchoof[{a, b, p}]

Future Work

● Convert existing code into a Mathematica 6.0 package

● Extend the code to the Schoof - Elkies - Atkin (SEA) algorithm

● Implement code for elliptic curves over ℂ making use of Weierstrass ℘-function

Thank you!  - Questions?

References

[1] Brown, Ezra - Square Roots from 1; 24, 51, 10 to Dan Shanks, College Math Journal (1999).

[2] Dummit, D., Foote, R. Abstract Algebra, John Wiley and Sons (2004)

[3] Lercier, R. and Morain, F., Counting the number of points on elliptic curves over finite fields: strategies and performances, Advances in Cryptology, Proc. Eurocrypt'95, LNCS 921, L.C. Guillou and J.J. Quisquater, Eds., Springer-Verlag, 1995, pp. 79--94.

[4] Rosen, Kenneth H. -  Elementary Number Theory and its Applications, Addison and Wesley (2000)

[5] Schoof, René - Elliptic curves over finite fields and the computation of square roots mod p, Mathematics of Computation, Vol. 44, No 170 (1985), 482-494

[6] Schoof, René - Counting points on elliptic curves over finite fields, Journal de Téorie des Nombres de Bordeaux 7 (1995), 219-254

[7] Silverman, Joseph - Rational Points on Elliptic Curves (1992)

[8] Singh, Simon - The Code Book - The Science of Secrecy from Ancient Egypt to Quantum Cryptography (1999)

[9] Washington, Lawrence - Elliptic Curves - Number Theory and Cryptography, Chapman & Hall (2003)

[10] Certicom (http://www.certicom.com)