Elliptic Curves

The purpose of this notebook is to show how to use the Mathematica packages EllipticCurves', which is available via
	 http://www.mat.univie.ac.at/pub/teschl/notebooks/EllipticCurves.m
and needs to be  stored in the AddOns/Applications subfolder of your Mathematica folder. On multi-user systems, you can install an add-on either in the central Mathematica directory (provided you have access to it) or else in your individual user's Mathematica directory (usually ~/.Mathematica/). (You can also put them  somewhere else, as long as you change the $PATH variable in Mathematica such that the files can be found by Mathematica or evaluate the packages before using this  notebook.) 

The latest version of this notebook is available from
	 http://www.mat.univie.ac.at/pub/teschl/notebooks/EllipticCurves.nb

©2023 by Gerald Teschl  <mailto:Gerald.Teschl@univie.ac.at>, <http://www.mat.univie.ac.at/~gerald/>

SetDirectory[NotebookDirectory[]];

Needs["EllipticCurves`"]

The code is mainly for educational purpose and you are invited to look at the implementation. The computation of the group order uses a basic version of Schoof’s algorithm and only works for small groups (around maximal 12 decimal digits, depending on your hardware and your patience). The computation of the order of a group element is done via the group order and inherits this limitation (unless the group order is known and can be factorized). Computation of discrete logs (DLP) is done using a Pohlig-Hellman reduction and Shank’s baby-step/giant-step algorithm (or optionally Pollard’s rho algorithm). Naturally this only works for small orders.

Basic operations

Use the curve associated with (y^3)= (x^3)+3x+1 modulo p=11.

ECset[{3,1},11]

The order of the group is

ECGroupOrder[]

The computation only works for small modulus p. For large groups it can be specified as a third (optional) argument in case it is known.

In case you forgot what group you are working with (note that it will remember the group order)

ECset[]

Choose two random points on the curve

{P,Q}=ECRandomPoint[2]

Add them

ECadd[P,Q]

or compute a multiple

ECmlt[3,P]

A point plus its negative is the identity

ECadd[P,ECneg[P]]

The order of an group element can be computed using

ECOrder[P]

We can check if it is on our curve via

ECPointQ[P]

The differential polynomial of order n will annihilate all points of order n

ECPsi[ECOrder[P],{x,y}]

ECPsi[ECOrder[P],P]

ECPsi[8,P]

The differential polynomials also satisfy the following relation

{a,b,c,d,p}=ECget[];
n=5;
Mod[{P[[1]]-ECPsi[n-1,P]ECPsi[n+1,P]PowerMod[((ECPsi[n,P])^2),-1,p],(ECPsi[n+2,P]((ECPsi[n-1,P])^2)-ECPsi[n-2,P]((ECPsi[n+1,P])^2))PowerMod[4 P[[2]] ((ECPsi[n,P])^3),-1,p]},p]==ECmlt[n,P]

Solve the discrete log problem (DLP) on the curve

ECLog[P,ECmlt[5,P]]

Here is a list of all points on the curve

ec=ECAllPoints[]

Create a Cayley table

Table[ECadd[ec[[i]],ec[[j]]],{i,Length[ec]},{j,Length[ec]}]//MatrixForm

Check associativity via the Cayley table (Light’s associativity test)

tmp=True;Do[If[i<=j<=k,tmp=And[tmp,ECadd[ec[[i]],ec[[j]],ec[[k]]]==ECadd[ec[[i]],ec[[k]],ec[[j]]]]],{i,Length[ec]},{j,Length[ec]},{k,Length[ec]}];tmp

Not all functions are listed here. Use the following to get the full list:

Names["EllipticCurves`*"]

Lenstra’s elliptic curve factorization method

The idea of H. Lenstra’s method is that, if p is not prime but composite, the group operations will fail. Moreover, the reason for this failure will give us a factor of p.

Here are the definitions of the group operations from the package. However, modified to catch the case, when the multiplicative inverse modulo p does not exists:

ecneg[{x1_,y1_}]:={x1,-y1};
ecdbl[{x1_,y1_}]:=Which[
	y1==0,{$Failed,0},
	True,Block[{s},
	s=GCD[p,2y1];
	If[s>1,Return[{"Success!",s},Block]];
	s=Mod[(3x1^2+a)PowerMod[2y1,-1,p],p];
	Mod[{s^2-2x1,s(x1-s^2+2x1)-y1},p]]];
ecadd[{x1_,y1_},{x2_,y2_}]:=Which[
	Mod[x1-x2,p]==0&&Mod[y1+y2,p]==0,{∞,∞},
	Mod[x1-x2,p]==0&&Mod[y1-y2,p]==0,ecdbl[{x1,y1}],
	True,Block[{s},
	s=GCD[p,x2-x1];
	If[s>1,Return[{"Success!",s},Block]];
	s=Mod[(y2-y1)PowerMod[x2-x1,-1,p],p];
	Mod[{s^2-x1-x2,s(x1-s^2+x1+x2)-y1},p]]];
ecmlt[n_Integer,P_List]:=Block[{Pn=ecneg[P],P2=P,h,hb,hl,k,kl,kb},
	h=Reverse[IntegerDigits[3n,2]];hl=Length[h];
	k=Reverse[IntegerDigits[n,2]];kl=Length[k];
	Do[
		P2=ecdbl[P2];
	If[P2[[1]]=="Success!",Return[P2,Block]];
		hb=h[[b]];
		If[b<=kl,kb=k[[b]],kb=0];
		If[{hb,kb}=={1,0},P2=ecadd[P2,P],If[{hb,kb}=={0,1},P2=ecadd[P2,Pn]]];
	If[P2=="Success!",Return[$Failed,Block]];
	,{b,hl-1,2,-1}];
	P2];

As a test, let us choose a Mersenne number

Mersenne[k_]:=(2^k)-1;
p=Mersenne[59];

We need to guess (a multiple of) the order of a randomly chosen point P0. To this end we take the product of all prime powers up to a given bound

PrimeProduct[S_]:=Block[{p=2,prod=1},While[p<= S,prod=prod*(p^(Floor[Log[p,S]])); p=NextPrime[p]];prod];

Computing the corresponding multiple of P0, the group operations will fail in case we have correctly guessed a multiple of the order of P0 with respect to the elliptic curve associates with a prime factor of p.

As a first curve we try

a=3;P0={1,1};b=Mod[((P0[[2]])^2)-(((P0[[1]])^3)+a P0[[1]] ),p];s=PrimeProduct[80];
ecmlt[s,P0]

This one did not work. Indeed the smallest prime factor of p is

p1=FactorInteger[p][[1,1]]

and the corresponding group order of our elliptic curve is

ECset[{a,b},p1];FactorInteger[ECGroupOrder[]]

Unfortunately our starting point has the largest factor contained in its order

ECOrder[P0]//FactorInteger

But this large prime factor was missed by our bound 80! So one option is to change P0 and hope that its order does not contain the largest prime factor. Here the chances are slim and it is difficult to find such a point! Here is one (note that finding the point requires knowledge of the factorization of p, so without this knowledge your best option is to choose a random point):

Until[Mod[ECOrder[P0],383]!=0,P0=ECRandomPoint[]];P0

With this point we succeed:

ecmlt[s,P0]

Alternatively, we can keep our original starting point and increase the bound:

a=3;P0={1,1};b=Mod[((P0[[2]])^2)-(((P0[[1]])^3)+a P0[[1]] ),p];s=PrimeProduct[400];
ecmlt[s,P0]

However, the real power of Lenstra’s method is that we can change the curve:

a=4;P0={1,1};b=Mod[((P0[[2]])^2)-(((P0[[1]])^3)+a P0[[1]] ),p];s=PrimeProduct[80];
ecmlt[s,P0]

This time our original bound was sufficient to cover all prime factors:

ECset[{a,b},p1];FactorInteger[ECGroupOrder[]]

ECOrder[P0]//FactorInteger

The bitcoin curve and the ECDSA

There are a few built-in cryptographic curves

ECset["secp256k1"]

Note: In case of large primes, the order must be specified, as the package can only handle small sized groups.

Here are the parameters:

ECset[];

PrimeQ[ECGroupOrder[]]

Here is the base point used by bitcoin (https://en.bitcoin.it/wiki/Secp256k1)

G=ECHexToPoint["02 79BE667E F9DCBBAC 55A06295 CE870B07 029BFCDB 2DCE28D9 59F2815B 16F81798"]

Check that it is on the curve

ECPointQ[G]

The order of the base point is of course the group order (since the latter is prime) and its bit length is

Ceiling[Log[2,ECGroupOrder[]]]

Create a random private key

SeedRandom[Method->"OpenSSL"];
a=RandomInteger[{2,ECOrder[G]-1}]

The corresponding public key is (https://en.bitcoin.it/wiki/Secp256k1)

A=ECmlt[a,G]

or in hex format

ECPointToHex[A]

Let us sign a text using the ECDSA (https://doi.org/10.6028/NIST.FIPS.186-5)

h=Hash["The Magic Words are Squeamish Ossifrage","SHA256"]

Since the bit length of the hash equals the bit length of the order of the base point, there is no need to cut the hash (which otherwise would be required by the standard).

BitAnd[(2^(Ceiling[Log[2,ECOrder[G]]]))-1,h]==h

Now we can compute the signature

m=ECOrder[G];
s=0;
While[s==0,
Block[{r},
ρ=0;
While[ρ==0,
r=RandomInteger[m-1];
ρ=Mod[ECmlt[r,G][[1]],m]
];
s=Mod[PowerMod[r,-1,m]Mod[h+ρ a,m],m];
]
]
{ρ,s}

Check the signature:

v1=Mod[h PowerMod[s,-1,m],m];
v2=Mod[ρ PowerMod[s,-1,m],m];
Mod[ECadd[ECmlt[v1,G],ECmlt[v2,A]][[1]],m]==ρ

Note that (ρ,-s) also works (since we are only checking the x-coordinate during verification):

v1=Mod[h PowerMod[-s,-1,m],m];
v2=Mod[ρ PowerMod[-s,-1,m],m];
Mod[ECadd[ECmlt[v1,G],ECmlt[v2,A]][[1]],m]==ρ

We could also create a bitcoin address as follows (https://en.bitcoin.it/wiki/Technical_background_of_version_1_Bitcoin_addresses):

We start with the public key in hex:

AA=ToLowerCase[StringDelete[ECPointToHex[A]," "]]

Make sure the x coordinate is 32 bytes long:

AA=StringTake[AA,2]<>StringRepeat["0",66-StringLength[AA]]<>StringTake[AA,{3,-1}]

To compute the hash (first SHA256, then RIPEMD160) we need to convert the hex value to bytes:

HexToBytes[h_]:=ExportByteArray[IntegerDigits[FromDigits[h,16],256,Ceiling[StringLength[h]/2]],"Byte"];

h=Hash[HexToBytes[AA],"SHA256","HexString"]

h="00"<>Hash[HexToBytes[h],"RIPEMD160","HexString"]

Finally we need the checksum (the first eight hex digits from twice SHA265):

check=Hash[HexToBytes[Hash[HexToBytes[h],"SHA256","HexString"]],"SHA256","HexString"]

h=h<>StringTake[check,8]

This result needs to be base58 encoded (https://en.bitcoin.it/wiki/Base58Check_encoding)

Base58CheckEncode[s_String]:=Module[{o="",ss=s},
While[StringTake[ss,2]=="00",
o=o<>"1";ss=StringDrop[ss,2]];
o<>StringJoin[IntegerDigits[FromDigits[ss,16],58]/.{0->"1",1->"2",2->"3",3->"4",4->"5",5->"6",6->"7",7->"8",8->"9",9->"A",10->"B",11->"C",12->"D",13->"E",14->"F",15->"G",16->"H",17->"J",18->"K",19->"L",20->"M",21->"N",22->"P",23->"Q",24->"R",25->"S",26->"T",27->"U",28->"V",29->"W",30->"X",31->"Y",32->"Z",33->"a",34->"b",35->"c",36->"d",37->"e",38->"f",39->"g",40->"h",41->"i",42->"j",43->"k",44->"m",45->"n",46->"o",47->"p",48->"q",49->"r",50->"s",51->"t",52->"u",53->"v",54->"w",55->"x",56->"y",57->"z"}]
]

Base58CheckEncode[h]

Mathematica can also create a key pair for you:

key=GenerateAsymmetricKeyPair["EllipticCurve"]

You can access the public/private key as follows:

ECmlt[key["PrivateKey"]["PrivateMultiplier"],G]==key["PrivateKey"]["PublicCurvePoint"]

In addition, the key contains byte arrays of both keys:

ExportByteArray[IntegerDigits[key["PrivateKey"]["PrivateMultiplier"],256],"Byte"]==key["PrivateKey"]["PrivateByteArray"]

tmp="04"<>IntegerString[key["PrivateKey"]["PublicCurvePoint"][[1]],16]<>IntegerString[key["PrivateKey"]["PublicCurvePoint"][[2]],16];
HexToBytes[tmp]==key["PrivateKey"]["PublicByteArray"]

Mathematica can also compute the Bitcoin address (note that it uses the long form of the public key for the computation)

BlockchainKeyEncode[key["PublicKey"],"Address",BlockchainBase->"Bitcoin"]

BarcodeImage[%,"QR"]

The Montgomery curve Curve25519

You can also play with the Montgomery curve (https://en.wikipedia.org/wiki/Curve25519)

ECset["Curve25519"]

The recommended generator is

G=ECPointLong[{9,0}]

On a Montgomery curve multiplication can be done without the y coordinate if only the x coordinate is of interest (like in elliptic curves signatures).

You can check that its order is given by

ECOrder[G]==2^252+27742317777372353535851937790883648493

PrimeQ[%]

ECPointQ[{0,0}]

Note that the point (0,0) is of order 2

ECmlt[2,{0,0}]

ECPointLong[{0,1}]

and the point (1,*) is of order 4

ECPointLong[{1,0}]

ECOrder[ECPointLong[{1,0}]]

The first point, which is on the curve, and has prime order is

j=1;While[!PrimeQ[ECOrder[ECPointLong[{j,0}]]],j++];j