An Elliptic Curve Package for Mathematica

John McGee, Wolfram Research, Technical Services

Join the Conversation #WolframTechConf

Initialization

Needs["EllipticCurves`"]

What is an Elliptic Curve?

An elliptic curve E over the field 𝔽 is the zero set of a non-singular cubic polynomial in two variables.

	E(𝔽)={(x,y)∈(OverBar[𝔽])|p(x,y)=0;p(x,y)∈𝔽[x,y],deg(p)=3}
	
So E is an algebraic set (variety), having an equation the general form

	(Subscript[a, 9]) (x^3)+(Subscript[a, 8]) (x^2) y+(Subscript[a, 7]) (x^2)+(Subscript[a, 6]) x (y^2)+(Subscript[a, 5]) x y+(Subscript[a, 4]) x+(Subscript[a, 3]) (y^3)+(Subscript[a, 2]) (y^2)+(Subscript[a, 1]) y+(Subscript[a, 0])==0
	
Using appropriate changes of variables for fields of characteristic other than 2 or 3, we can transform the general curve into short Weierstrass form as:

	(y^2)=(x^3)+A x + B   where 4(A^3)+27 (B^2)≠0

Some general cubic curves

Why Elliptic Curves?

* Mathematically Interesting
	number theory, group theory, commutative ring theory, complex variables ...
	
* Algebraic Number Theory
	the congruent number problem, Fermat’s Last Theorem
	
* Computation Number Theory
	Elliptic curve integer factorization
	
* Cryptography
	Massey-Omura, El Gamal

Chord and Tangent Point Arithmetic

We can create a geometric arithmetic of points by defining:

	P⊕Q is the reflection about the x-axis, of the third point of intersection 
		of the line through P-Q with the curve.  We use the tangent line if P=Q

Algebraic Point Arithmetic

Chord and Tangent Arithmetic takes on a simple algebraic form for the curve in short Weierstrass form.  We find for P=((Subscript[x, 1]),(Subscript[y, 1])),Q=((Subscript[x, 2]),(Subscript[y, 2])) on (y^2)=(x^3)+A x+B

	λ=({{{, ({{(((Subscript[y, 2])-(Subscript[y, 1]))/((Subscript[x, 2])-(Subscript[x, 1]))), (Subscript[x, 2])-(Subscript[x, 1])≠0}, {((3 ((Subscript[x, 1])^2)+A)/(2 (Subscript[y, 1]))), P=Q}})}}) 	     ({{(Subscript[x, 3])=(λ^2)-(Subscript[x, 1])-(Subscript[x, 2])       }, {(Subscript[y, 3])=λ ((Subscript[x, 1])-(Subscript[x, 3]))-(Subscript[y, 1])}})
	
This gives an arithmetic of points, but we do not have a group yet,

	we need an identity and inverses

Derivation of solution for intersection point (Subscript[x, 3])

The intersection of the line y=λ x+r with the curve (y^2)=(x^3)+A x +B gives a cubic in x only

PolynomialForm[Collect[(((λ x+r))^2)-((x^3)+a x+b),x],TraditionalOrder->True]

But we know that (Subscript[x, 1]),(Subscript[x, 2]) are solutions so we have

PolynomialForm[Collect[(((λ x+r))^2)-((x^3)+a x+b)==Expand[(x-(Subscript[x, 1]))(x-(Subscript[x, 2]))(x-(Subscript[x, 3]))],x],TraditionalOrder->True]

Equating the coefficients of (x^2) and solving for (Subscript[x, 3]) gives

Solve[Coefficient[(((λ x+r))^2)-((x^3)+a x+b)-Expand[(x-(Subscript[x, 1]))(x-(Subscript[x, 2]))(x-(Subscript[x, 3]))],(x^2)]==0,(Subscript[x, 3])]

Projective Space and the Point at Infinity

We can think of the projective space (ℙ^2) as the set of all directions in (ℝ^3), i.e. the set of equivalence classes of nonzero points in (ℝ^3), where P∼Q⇔Q=k P for some non-zero constant k.

	(ℙ^2)={[x,y,z]∈(ℝ^3)|(x,y,z)∼(a,b,c)⇔(a,b,c)=(k x,k y,k z)}
	
We have an injective copy of (ℝ^2)->(ℙ^2)  where (x,y)|->[x,y,1], but (ℙ^2) has extra elements, points at infinity of the form [x,y,0].

We can transform our Weierstrass equation into an equation of the complete elliptic curve as a projective variety with the map

	f(x,y)|->(z^3)f(x/z,y/z)  in our case giving F(x,y,z)=(z^3)(((((y/z)))^2)-((((x/z)))^3)-A (x/z)-B) = (y^2)z-(x^3)-A x (z^2)-B (z^3)
	
	Since this polynomial is homogeneous, if F(x,y,z)=0 then F(k x,k y,k z)=(k^3)F(x,y,z)=0
	
	When z=1 we obtain the original curve F(x,y,1)=f(x,y)=(y^2)-(x^3)-A x-B=0
	
	When z=0, we obtain F(x,y,z)=-(x^3)=0, which gives us a new solution in projective space [0,1,0].
	
	It turns out that this is the identity we need.  We define 𝒪=[0,1,0].  It corresponds to the vertical direction.

A representation of projective space

pnts=Flatten[Table[{x,y,1},{x,-5,5},{y,-5,5}],1];

Show[ContourPlot3D[{z==1,z==0},{x,-6,6},{y,-6,6},{z,-1,2},Mesh->None,ContourStyle->{{Opacity[0.3],Green}},Boxed->False,Axes->False,ImageSize->1000],Graphics3D[{PointSize[0.01],Point[{0,0,0}],Red,Point/@pnts,Magenta,Line[{{0,0,0},2 #}]&/@pnts,Line[{{0,0,0},{25 Cos[#],25 Sin[#],0}}]&/@Range[0,2 π,π/16]}]]

The Elliptic Curve Group Law

With the algebraic definition of point addition above, along with the point at infinity 𝒪, we have the group law

	(1)	P+Q=R			using the addition laws above
	
	(2)	P+𝒪=P			𝒪 is the identity
	
	(3)	P+(-P)=\[ScriptZero]			inverses P=((Subscript[x, 0]),(Subscript[y, 0])),-P=((Subscript[x, 0]),-(Subscript[y, 0]))
	
	(4)	P+Q=Q+P			commutativity
	
	(5)	P+(Q+R)=(P+Q)+R	associativity
	
The proof of associativity is non-trivial.

Since proof of the group law depends only of the field properties, these same arithmetic gives an elliptic curve for any field, ℚ,ℂ and the finite fields (Subscript[𝔽, (p^n)])

Examples

Names["EllipticCurves`*"]

?ecJinvariant

?ecHasseBounds

randomECPoint[]

p=Prime[12345]

{a,b}=RandomInteger[p,2]

isECQ[a,b,p]

ecJinvariant[]

setEC[a,b,p]

pntP=randomECPoint[]

pntQ=randomECPoint[]

pntP⊕pntQ

14⊗((3⊗pntP)⊕(4⊗pntQ)⊕pntP)

NestList[(#⊕pntP)&,pntP,10-1]

Table[k⊗pntP,{k,1,10}]

25⊗pntP

The Order of a Point

A point P has finite order n if n P=𝒪 and n is the smallest positive integer that makes this true.

We can compute n P efficiently using the binary expansion of n = (Subscript[b, 0])+(Subscript[b, 1])2+(Subscript[b, 2])(2^2)+...+(Subscript[b, m])(2^m).
    then n P=(Subscript[b, 0])P+(Subscript[b, 1])2 P+(Subscript[b, 2])2(2P)+(Subscript[b, 3])2(4P)+...+(Subscript[b, m]) 2 ((2^(m-1))P).

Torsion Points: E[n]={P∈E((OverBar[K]))|n P=𝒪}.     E[n]∼(Subscript[ℤ, n])×(Subscript[ℤ, n])     |E[n]|=(n^2)

If |P| is small we can find it by direct search, else Shank’s Baby Step - Giant Step algorithm is effective.

	Given P∈E((Subscript[𝔽, q])), compute Q=(q+1)P then find a match  j P=Q+k (2 m P} for some j,k
	
		this match usually happens quickly due to the birthday paradox.
		
	Then (q+1-2 k m-j)P=M P=𝒪 so that |P| is a divisor of M.
	
Shank’s Baby Step - Giant Step algorithm is 𝒪((q^(1/4))), so not feasible for large q.

Difficult to parallelize

$tbljP=NestList[ecAdd[#,pntP]&,stpnt,jblksize-1]; (* {P,2P,...,nP} *)

$tbljP=AssociationThread[$tbljP[[All,1]],Range[Length[$tbljP]]-1];   (* <|(Subscript[kP, x])->k|>*)

q[{a_,b_,c_}]:={ecAdd[a,pnt2mP],ecAdd[b,pntn2mP],c+1};

{kpnt2mP,kpntn2mP,k}=NestWhile[q,{kpnt2mP,kpntn2mP,1},Not[KeyExistsQ[$tbljP,#[[1,1]]]||KeyExistsQ[$tbljP,#[[2,1]]]]&,1,m];

A Theorem of Gauss

In general it is not easy to determine |E((Subscript[𝔽, q]))|. For the special case of the Fermat Equation (x^3)+(y^3)=1 over (Subscript[𝔽, p])

A Theorem of Gauss
	Let (Subscript[M, p]) be the number of projective solutions to the homogeneous equation (x^3)+(y^3)+(z^3)=0 over (Subscript[𝔽, p])
	
	(a) If p≢1 (mod 3) then (Subscript[M, p])=p+1  with 1 point at infinity [1,-1,0]
	
	(b) if p≡1(mod 3) then (Subscript[M, p])=p+1+A 
		where A,B are integer solutions to 4p=(A^2)+27(B^2)
		with 3 points at infinity [1,-1,0],[1,-1 u,0],[1,-1 (u^2),0] where ord(u)=3
		Since |(Subscript[𝔽, p]^×)|=p-1=3 k so (Subscript[𝔽, p]^×) has a element of order 3 by Cauchy’s Theorem
		
For the more general case we need to do more work ...

Example

(* Return the number of projective solutions to (x^3)+(y^3)+(z^3) = 0 *)
(* return {N,I,Z) where N is the number of projective points, I is the number of points at infinity, Z = p(mod 3) *)
GaussFermatCurvePoints[p_]:=Module[{a,b,r,s},
If[Mod[p,3]≠1,Return[{p+1,1,Mod[p,3]}]];
{a,b}={r,s}/.Flatten@Solve[4p==(r^2)+27(s^2)∧r>0∧s>0,{r,s},Integers];
If[Mod[a,3]==2,a*=-1];
{p+1+a,3,Mod[p,3]}
]

GaussFermatCurvePoints[31]

({{X, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, {Y, 1, 1, 1, 0, 0, 0, 1, 1, 1, 5, 5, 5, 6, 11, 11, 11, 17, 17, 17, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 27, 27, 27, 30}, {Z, 6, 26, 30, 6, 26, 30, 11, 24, 27, 11, 24, 27, 0, 1, 5, 25, 17, 22, 23, 17, 22, 23, 17, 22, 23, 1, 5, 25, 11, 24, 27, 0, 1, 5, 25, 0}})

GaussFermatCurvePoints[29]

({{X, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, {Y, 1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28}, {Z, 28, 28, 3, 24, 1, 13, 8, 19, 9, 5, 7, 10, 26, 15, 4, 21, 12, 27, 23, 20, 6, 18, 14, 25, 17, 2, 22, 11, 16, 0}})

A Simple Approach to Group Order

A theorem of Lagrange tells us that ord(P)=|<P>| divides |E((Subscript[𝔽, q]))|

A theorem due to Hasse bounds the group order:  q+1-2(Sqrt[q])≤|E((Subscript[𝔽, q]))|≤q+1+2(Sqrt[q])

If we know the order of a set of points {(Subscript[O, 1]),(Subscript[O, 2]),...(Subscript[O, k])} where M=LCM[(Subscript[O, 1]),(Subscript[O, 2]),...,(Subscript[O, K])] is the unique

	integer q+1-2(Sqrt[q])≤ M ≤q+1+2(Sqrt[q]) then M=|E((Subscript[𝔽, q]))|
	
Lemma:  If P∈E((Subscript[𝔽, p])) with ord(P)>4 (Sqrt[p]) then k=ord(P) has a unique multiple M within the Hasse bounds.

Lemma:  If E((Subscript[𝔽, p])) contains one or more points P with ord(P)>4 (Sqrt[p]) then a randomly chosen point is likely
	   to have large order.  (|E[n]|=(n^2) )

Experiments of Point Orders

For the curve given by (y^3)=(x^3)+A x+B over (Subscript[𝔽, p]) with
      A=75875415081265986964805, B=51748887602164205111353,P= 81155059153329437714899
      
 We have the following point order distribution for 1000 randomly chosen points.  Each of these orders >4 (Sqrt[p]).

Division Polynomials

For the elliptic curve f(x,y)=(y^2)-(x^3)-A x-B=0, the (n^th) division polynomial (Subscript[ψ, n])(x,y)=0∀P∈E[n]

These are defined recursively as:

	(Subscript[ψ, 0])=0, (Subscript[ψ, 1])=1, (Subscript[ψ, 2])=2 y                   	Note: E[0]=E((OverBar[K])), E[1]={𝒪}, E[2]={((Subscript[e, 1]),0},{(Subscript[e, 2]),0},{(Subscript[e, 3]),0}}
	
	(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[ψ, 2 m+1])=(Subscript[ψ, m+2]) (Subscript[ψ, m]^3)-(Subscript[ψ, m-1])(Subscript[ψ, m+1]^3)		(Subscript[ψ, 2m])=(((2y))^(-1))(Subscript[ψ, m])((Subscript[ψ, m+1])(Subscript[ψ, m+1]^2)-(Subscript[ψ, m-2])(Subscript[ψ, m+1]^2))
	
Notes: the degree of these polynomials grows quadratically with n  (E[n]∼(Subscript[ℤ, n])×(Subscript[ℤ, n]))

	we can reduce the exponents of y to zero or one with (y^(2n))->((((x^3)+A x+B))^n)  i.e. (Subscript[ψ, n])(x,y)∈K[x,y]/(f)

Direct ‘Computation’ of n P

Using the division polynomials (Subscript[ψ, n])(x,y) we define two related polynomials

	(Subscript[ϕ, m])=x (Subscript[ψ, m]^2)-(Subscript[ψ, m+1])(Subscript[ψ, m-1])		(Subscript[ω, m])=(((4 y))^(-1))((Subscript[ψ, m+2])(Subscript[ψ, m-1]^2)-(Subscript[ψ, m-2])(Subscript[ψ, m+1]^2))
	
Then it can be shown (by induction on n) that

	n P=((((Subscript[ϕ, n])(x))/((Subscript[ψ, n]^2)(x))),(((Subscript[ω, n])(x,y))/(((Subscript[ψ, m])(x,y))^3)))
	
Notes:  the rational function for the x-coordinate is a polynomial in x only.

	this is not a good way to compute n P for a specific n,P

Schoof’s Algorithm for |E((Subscript[𝔽, q]))|

Theorem of Hasse:  for E((Subscript[𝔽, q])) if N=|E((Subscript[𝔽, q]))| then

	q+1-2(Sqrt[q])≤N ≤q+1+2(Sqrt[q])     so |E((Subscript[𝔽, q]))|=q+1-a  with |a|≤2 (Sqrt[q])
	
(1) Choose a set of primes S={2,3,5,...,L} such that (Subscript[Π, 𝓅∈S])𝓅>4(Sqrt[q])

(2) For each prime 𝓅∈S determine τ (mod 𝓅) using the division polynomials (Subscript[ψ, 𝓅])
	(can be run in parallel)

	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, 19]),(Subscript[f, l]))≠1 then ((Subscript[ϕ, p])^2)+k≡±τ (Subscript[ϕ, p]) (mod l)
	13.		    for some point in E[l] so we test
	14.			If gcd((Subscript[p, 19]),(Subscript[f, l]))≠1 then  t≡τ(mod l) else t≡-τ (mod l) 
	15.	    Next τ

(3) Determine τ using the Chinese remainder theorem

Problem: the size of the division polynomials.  Solution: Elkies, Atkins (SEA)

	choose special primes 𝓅, such that the division polynomials have factors ∼ (𝓅-1)/2, and compute mod these factors
	
	Gives approximatley 𝒪((log^6)(n)) performance.

Example

r1=Timing[ecGroupOrderSchoof[ec8]]

some time passes ....

FactorInteger[Last[r1]]

The Discrete Logarithm Problem

Given two points P,Q=k P, determine k.  This problem is hard for a general abelian group G.
If |G|=𝓅 then G∼(Subscript[𝔽, 𝓅]) and the index calculus can be used to solve for k.

The (SuperStar[MOV]) attack is more general, it uses the Weil pairing to convert the problem in E((Subscript[𝔽, q])) to one in (Subscript[𝔽, (q^m)]^×)

	where E[n]⊆E((Subscript[𝔽, (q^m)]^×)), then (Subscript[e, n]):E[n]×E[n]->(Subscript[μ, n]) where (P,Q)↦ζ; (ζ^n)=1.
	
(SuperStar[ ])Menezes, Okamoto, and Vanstone

In any case, knowing the size of the group E((Subscript[𝔽, q])) and ord(P) are critical for security of EC cryptography

Future Work  (as time permits)

Elliptic Curves over ℂ

	The Weil Paring      (Subscript[ℯ, n]):E[n]×E[n]->(Subscript[μ, n])
	
	Divisors
	
	Finding Elkies Primes

Thank you! - Questions?

References

R. Crandall and Carl Pomerance. Prime Numbers a Computational Perspective (2001)
	PrimeKit (Mathematica)
	
R. Schoof. Elliptic curves over finite fields and the computation of square roots mod p. (Math. Comp. 1985)

J. H. Silverman and J. Tate. Rational Points on Elliptic Curves (1992)

L. Washington. Elliptic Curves Number Theory and Cryptography (2008)