(* ::Package:: *)

(* :Title: Elliptic Curves *)

(* :Author: Gerald Teschl 
			http://www.mat.univie.ac.at/~gerald/
			Gerald.Teschl@univie.ac.at  *)

(* :Summary:
	This package provides some functions for working with elliptic curves over Z_p. *)

(* :Context: EllipticCurves` *)

(* :Package Version: 1.2.1 *)

(* :Copyright: Copyright 2023, G. Teschl  *)

(* :History: 
	Version 1.0 by G. Teschl
	Version 1.1 by G. Teschl
	Version 1.2 by G. Teschl *)

(* :Keywords: Elliptic Curves *)

(* :Source: Some algoritms are from R. Crandall and C. Pomerance, Prime Numbers: a Computational Perspective, Springer 2001. *)

(* :Warning: None. *)

(* :Mathematica Version: 6.0 *)

(* :Limitation: Some algortihms are not sutiable for handling large numbers. *)

(* :Discussion: TODO.
    *)


BeginPackage["EllipticCurves`"];

ECset::usage = 
	"ECset[{a,b,c,d},p] defines the elliptic curve d y^2 = x^3 + c x^2 + a x +b over Z_p.";
ECget::usage = 
	"ECget[] retuns the curve parameters {a,b,c,d,p}.";
ECadd::usage = 
	"ECadd[P,Q] add two points.";
ECdbl::usage = 
	"ECdbl[P] double a point.";
ECsub::usage = 
	"ECsub[P,Q] subtract two points.";
ECneg::usage = 
	"ECneg[P] negative point.";
ECmlt::usage = 
	"ECmlt[n,P] compute the multiple of a point.";
ECaffine::usage = 
	"ECaffine[P] converte a point from modified projective to affine coordinates.";
ECmpc::usage = 
	"ECmpc[P] converte a point from affine to modified projective coordinates.";
ECPointQ::usage = 
	"ECPointQ[P] test if a point is on the curve.";
ECRandomPoint::usage = 
	"ECRandomPoint[] generate a random point on the curve.";
ECAllPoints::usage = 
	"ECAllPoints[] return all points on the curve.";
ECPointShort::usage = 
	"ECPointShort[P] convert a point to short form.";
ECPointLong::usage = 
	"ECPointLong[P] convert a point to standard form.";
ECHexToPoint::usage = 
	"ECHexToPoint[string] convert a point in hexadecimal format to standard form.";
ECPointToHex::usage = 
	"ECPointToHex[P] convert a point to hexadecimal format.";
ECOrder::usage = 
	"ECOrder[P] compute the order of P.";
ECGroupOrder::usage = 
	"ECGroupOrder[] compute the order of the group.";
ECPsi::usage =
	"ECPis[n_,x_,y_] division polynomial.";
ECLog::usage =
	"ECLog[a,b,p] computes the discrete logarithm of b with base a modulo p.";
ECLogRho::usage =
	"ECLogRho[a,b,p] computes the discrete logarithm of b with base a modulo p using Pollard's rho method.";
ECLogBSGS::usage =
	"ECLogBSGS[a,b,p] computes the discrete logarithm of b with base a modulo p using Shank's baby-step/giant-step method.";



Begin["`Private`"];


(* Set the curve *)

ECset::singular="Error! Singular curve.";
ECset::noprime="Error! p must be an odd prime.";
ECset::pequal3="Warning! No reduced Weierstrass form for p=3. Division polynomials are unsupported.";
ECset::cnonzero="Warning! Modified projective coordinates will not work if c is nonzero.";

ECa=0;ECb=1;ECc=0;ECd=1;ECaa=ECa;ECbb=ECb;ECp=5;ECn=6;ECnf=0;

ECset[{a_Integer,b_Integer,c_Integer:0,d_Integer:1},p_Integer,n_Integer:0]:=Block[{\[CapitalDelta]=Mod[4a^3+27 b^2+c(18 a b - a^2 c + 4 b c^2),p]},
	If[\[CapitalDelta]==0,Message[ECset::singular];Return[$Failed,Block]];
	If[p<3||!PrimeQ[p],Message[ECset::noprime];Return[$Failed,Block]];
	If[Mod[d,p]==0,Message[ECset::dzero];Return[$Failed,Block]];
	ECp=p;ECa=Mod[a,p];ECb=Mod[b,p];ECc=Mod[c,p];ECd=Mod[d,p];ECdi=PowerMod[ECd,-1,p];
	If[ECp!=3,
		Block[{co3=Mod[PowerMod[3,-1,p] ECd ECc,ECp]},ECaa=Mod[ECd^2 ECa- 3 co3^2,ECp];ECbb=Mod[ECd^3ECb+co3(2 co3^2-ECd^2 ECa),ECp];],
		ECaa=ECa;ECbb=ECb;
		Print["Warning: Modified Projective Coordinates not supported!"];
	];
	ECn=n;
	ResetECpsi[];
	If[n>0,ECnf=FactorInteger[ECn]];
];
ECset[]:=Block[{},
	Print["p=",ECp];
	Print["Curve: "<>If[ECd>1,ToString[ECd],""]<>"\!\(\*SuperscriptBox[\(y\), \(2\)]\) = \!\(\*SuperscriptBox[\(x\), \(3\)]\)"<>If[ECc>0," + "<>If[ECc>1,ToString[ECc],""]<>"\!\(\*SuperscriptBox[\(x\), \(2\)]\)",""]<>If[ECa>0," + "<>If[ECa>1,ToString[ECa],""]<>"x",""]<>If[ECb>0," + "<>ToString[ECb],""]];
	If[ECn>0,Print["Group order: ",ECn]];
	{ECa,ECb,ECc,ECd,ECp}
];
ECget[]:={ECa,ECb,ECc,ECd,ECp};
ECget[1]:={ECaa,ECbb,ECp};



ECset["secp256k1"]:=ECset[{0,7},2^256-2^32-2^9-2^8-2^7-2^6-2^4-1,115792089237316195423570985008687907852837564279074904382605163141518161494337];
ECset["secp192r1"]:=ECset[{-3,2455155546008943817740293915197451784769108058161191238065},6277101735386680763835789423207666416083908700390324961279,6277101735386680763835789423176059013767194773182842284081];
ECset["secp256r1"]:=ECset[{-3,41058363725152142129326129780047268409114441015993725554835256314039467401291},115792089210356248762697446949407573530086143415290314195533631308867097853951,115792089210356248762697446949407573529996955224135760342422259061068512044369];
ECset["Curve25519"]:=ECset[{1,0,486662},2^255-19,8(2^252+27742317777372353535851937790883648493)];


(* Basic arithmetic / affine coordinates *)

ECneg[{\[Infinity],\[Infinity]}]:={\[Infinity],\[Infinity]};
ECneg[{x_Integer,y_Integer}]:={x,Mod[-y,ECp]};
ECdbl[{x1_,y1_}]:=Which[{x1,y1}=={\[Infinity],\[Infinity]},{\[Infinity],\[Infinity]},
	y1==0,{\[Infinity],\[Infinity]},
	True,With[{s=Mod[ECdi(3x1^2+2 ECc x1+ECa)PowerMod[2y1,-1,ECp],ECp]},With[{x3=Mod[ECd s^2-2x1-ECc,ECp]},{x3,Mod[s(x1-x3)-y1,ECp]}]]];
ECadd[{x1_,y1_},{x2_,y2_}]:=Which[{x1,y1}=={\[Infinity],\[Infinity]},{x2,y2},
	{x2,y2}=={\[Infinity],\[Infinity]},{x1,y1},
	Mod[x1-x2,ECp]==0&&Mod[y1+y2,ECp]==0,{\[Infinity],\[Infinity]},
	Mod[x1-x2,ECp]==0&&Mod[y1-y2,ECp]==0,ECdbl[{x1,y1}],
	True,Block[{s=Mod[(y2-y1)PowerMod[x2-x1,-1,ECp],ECp]},
	With[{x3=Mod[ECd s^2-x1-x2-ECc,ECp]},{x3,Mod[s(x1-x3)-y1,ECp]}]]];
ECadd[P1_List,P2_List,P3__]:=ECadd[ECadd[P1,P2],P3];
ECsub[P1_List,P2_List]:=ECadd[P1,ECneg[P2]];
ECmlt[0,{x_,y_}]:={\[Infinity],\[Infinity]};
ECmlt[n_Integer,P_List]:=ECmlt[-n,ECneg[P]]/; n<0;
ECmlt[n_Integer,P_List]:=Module[{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];
		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]]]
	,{b,hl-1,2,-1}];
	P2];

(* Basic arithmetic / modified projective coordinates *)

ECneg[{0,1,0}]:={0,1,0};ECneg[{X_Integer,Y_Integer,Z_Integer}]:=Mod[{X,-Y,Z},ECp];
ECdbl[{X_Integer,Y_Integer,Z_Integer}]:=Block[{M,S,X2},
	If[Mod[X Y,ECp]==0,Return[{0,1,0},Block]];
	M=Mod[3PowerMod[X,2,ECp]+ECa PowerMod[Z,4,ECp],ECp];S=Mod[4 X PowerMod[Y,2,ECp],ECp];
	X2=Mod[PowerMod[M,2,ECp]-2S,ECp];
	Mod[{X2,M(S-X2)-8 PowerMod[Y,4,ECp],2Y Z},ECp]
	];
ECadd[{X1_Integer,Y1_Integer,Z1_Integer},{X2_Integer,Y2_Integer,Z2_Integer}]:=Block[{U1,U2,S1,S2,W,R,T,M,X3},
	If[Z1==0,Return[{X2,Y2,Z2},Block]];
	If[Z2==0,Return[{X1,Y1,Z1},Block]];
	U1=Mod[X2 PowerMod[Z1,2,ECp],ECp];U2=Mod[X1 PowerMod[Z2,2,ECp],ECp];
	S1=Mod[Y2 PowerMod[Z1,3,ECp],ECp];S2=Mod[Y1 PowerMod[Z2,3,ECp],ECp];
	W=Mod[U1-U2,ECp];R=Mod[S1-S2,ECp];
	If[W==0,
		If[R==0,Return[ECdbl[{X1,Y1,Z1}],Block],Return[{0,1,0},Block]]
	];
	T=Mod[U1+U2,ECp];M=Mod[S1+S2,ECp];
	X3=Mod[R^2-T W^2,ECp];
	Mod[{X3,PowerMod[2,-1,ECp]((T W^2-2 X3)R-M W^3),Z1 Z2 W},ECp]
	];
ECmlt[0,{X_,Y_,Z_}]:={0,1,0};

(* Conversion to affine coordinates and back *)

ECaffine[{0,Y_,0}]:={\[Infinity],\[Infinity]};
ECaffine[{X_,Y_,Z_}]:=Mod[{X PowerMod[Z,-2,ECp],Y PowerMod[Z,-3,ECp]},ECp] /; Z!=0;

ECmpc[{\[Infinity],\[Infinity]}]:={0,1,0};
ECmpc[{x_,y_}]:=Mod[{x,y,1},ECp];

(* Functions to test/generate points *)

ECPointQ[{x_,y_}]:=If[{x,y}=={\[Infinity],\[Infinity]},True,Mod[ECd y^2-x(x^2+ECc x+ECa) -ECb,ECp]==0];
ECPointQ[{X_,Y_,Z_}]:=If[{X,Z}=={0,0},True,ECPointQ[ECaffine[{X,Y,Z}]]];

ECRandomPoint[n_:1]:=Block[{x,y2,y},
	x=RandomInteger[ECp-1];
	y2=ECdi(x(x^2+ECc x+ECa) +ECb);
	While[JacobiSymbol[y2,ECp]==-1,
		x=RandomInteger[ECp-1];
		y2=ECdi(x(x^2+ECc x+ECa) +ECb);
		];
	y=PowerMod[y2,1/2,ECp];
	If[RandomInteger[1]==1,{x,y},ECneg[{x,y}]]
];
ECRandomPoint[n_Integer]:=Table[ECRandomPoint[],{n}];

ECAllPoints[]:=Block[{tmp={{\[Infinity],\[Infinity]}},y2,y},
Do[
	y2=Mod[ECdi(x(x^2+ECa ) +ECb),ECp];
	If[y2==0,tmp=Append[tmp,{x,0}],
		If[JacobiSymbol[y2,ECp]==1,y=PowerMod[y2,1/2,ECp];tmp=Append[tmp,{x,y}];tmp=Append[tmp,ECneg[{x,y}]]]
	]
	,{x,0,ECp-1}];
	tmp];

ECPointShort[{x_Integer,y_Integer}]:={x,Mod[y,2]};
ECPointShort[{X_Integer,Y_Integer,Z_Integer}]:=ECPointShort[ECaffine[{X,Y,Z}]];
ECPointLong[{x_Integer,s_Integer}]:=Block[{y2=Mod[x(x^2+ECc x+ECa)+ECb,ECp],y},
	y=PowerMod[y2,1/2,ECp];
	If[Mod[y-s,2]==0,{x,y},ECneg[{x,y}]]
];
ECHexToPoint[str_String]:=Block[{str1=StringTake[str,2],str2=StringDelete[StringDrop[str,2]," "]},
	If[str1=="04",
		FromDigits[StringPartition[str2,StringLength[str2]/2],16],
		ECPointLong[{FromDigits[str2,16],If[str1=="02",0,1]}]
	]
];
ECPointToHex[P_List]:=Block[{PX,s},
	{PX,s}=ECPointShort[P];
	PX=IntegerString[PX,16];
	PX=ToUpperCase[StringInsert[PX," ",Range[9,StringLength[PX],8]]];
	If[s==0,"02 ","03 "]<>PX
];

(* Compute the order of a point *)

ECOrder[P_List]:=Block[{p,k,o,Pinf},
	If[ECn==0,ECGroupOrder[]];
	If[ECnf==0,ECnf=FactorInteger[ECn]];
	o=ECn;
	If[Length[P]==2,Pinf={\[Infinity],\[Infinity]},Pinf={0,1,0}];
	Do[
		{p,k}=ECnf[[j]];
		i=1;While[i<=k&&ECmlt[ECn/p^i,P]==Pinf,i++];
		o=o/p^(i-1),{j,Length[ECnf]}];
	o];

Options[ECGroupOrder] = {Method -> "Automatic", Verbose -> 0};
Options[ECSchoof] = {Verbose -> 0};

ECGroupOrder[n_Integer:0,opts:OptionsPattern[ECGroupOrder]]:=Block[{},
	If[n>0,
		ECn=n;
		ECnf=0;
	];
	If[ECn==0,
		If[ECp==3||(ECp< 10^4&&OptionValue[Method]=="Automatic")||OptionValue[Method]=="Jacobi",
			If[OptionValue[Verbose]>0,Print["Summing Jacobi symbol."]];

			ECn=ECp+1+\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(x = 0\), \(ECp - 1\)]\(JacobiSymbol[PowerMod[x, 3, ECp] + ECc\ PowerMod[x, 2, ECp] + ECa\ x\  + ECb, ECp]\)\),
			If[OptionValue[Verbose]>0,Print["Using Schoof algorithm."]];
			ECn=ECSchoof[FilterRules[{opts}, Options[ECSchoof]]]
		];
		ECnf=0;
	];
	ECn
]

(* Division polynomials, we store the values, but need to reset them if a new curve is selected *)

ResetECpsi[]:=Block[{},
	Clear[ECpsi];
	ECpsi[-1,x_] = -1;
	ECpsi[0,x_] = 0;
	ECpsi[1,x_] = 1;
	ECpsi[2,x_] = 2;
	ECpsi[3,x_] = 3 x^4+ 6 ECaa x^2 + 12 ECbb x - ECaa^2;
	ECpsi[4,x_] = Expand[4(x^6 + 5 ECaa x^4 + 20 ECbb x^3 - 
		5 ECaa^2 x^2 - 4 ECaa ECbb x - 8 ECbb^2 - ECaa^3)];
	ECpsi[m_Integer?(#>4 &),x_]:=ECpsi[m,x]=PolynomialMod[If[EvenQ[m], 
		ECpsi[m/2,x]/2 *(ECpsi[m/2+2,x] ECpsi[m/2-1,x]^2 - ECpsi[m/2-2,x] ECpsi[m/2+1,x]^2),
		If[EvenQ[(m-1)/2],(x^3+ECaa x +ECbb)^2,1]ECpsi[(m-1)/2 + 2,x] ECpsi[(m-1)/2,x]^3 -If[OddQ[(m-1)/2],(x^3+ECaa x +ECbb)^2,1]
		ECpsi[(m-1)/2 + 1,x]^3 ECpsi[(m-1)/2 -1,x]
		],ECp];
];
ECPsi[m_Integer,{x_Integer,y_Integer}]:=Mod[If[EvenQ[m],y,1]ECpsi[m,x],ECp];
ECPsi[m_Integer,{x_,y_}]:=If[EvenQ[m],y,1]If[ECc==0,ECpsi[m,x],PolynomialMod[ECpsi[m,x+PowerMod[3,-1,ECp]ECc],ECp]];



(* Schoof algortihm to compute the group order *)

(* Internal functions *)

reduce[pol_]:=PolynomialMod[pol,ss,Modulus->ECp];
polneg[{X_,Y_,Z_}]:=PolynomialMod[{X,-Y,Z},ECp];
poldbl[{X_,Y_,Z_}]:=Block[{M,S,X2,Y2},
	If[X==0||Y==0,Return[{0,1,0},Block]];
	M=reduce[Z^2];M=reduce[M^2];Y2=reduce[Y^2];
	M=reduce[3X^2+ECaa M];S=reduce[4 X (x^3+ECaa x +ECbb)];S=reduce[S Y2];
	X2=reduce[M^2-2S];
	reduce[{X2 (x^3+ECaa x +ECbb)Y2,(x^3+ECaa x +ECbb)reduce[Y Y2](reduce[M(S-X2)]-8(x^3+ECaa x +ECbb)^2 reduce[Y2^2]),2(x^3+ECaa x +ECbb)Y2 Z}]
];
poladd[{X1_,Y1_,Z1_},{X2_,Y2_,Z2_}]:=Block[{U1,U2,S1,S2,W,R,T,M,W2,X3,Y3,Z3},
	If[Z1==0,Return[{X2,Y2,Z2},Block]];
	If[Z2==0,Return[{X1,Y1,Z1},Block]];
	U1=reduce[ Z1^2];U2=reduce[Z2^2];
	S1=reduce[U1 Z1];S2=reduce[U2 Z2];
	S1=reduce[Y2 S1];S2=reduce[Y1 S2];
	U1=reduce[X2 U1];U2=reduce[X1 U2];
	W=reduce[U1-U2];R=reduce[S1-S2];
	If[W===0,
		If[R===0,Return[poldbl[{X1,Y1,Z1}],Block],Return[{0,1,0},Block]]
	];
	T=reduce[U1+U2];M=reduce[S1+S2];
	W2=reduce[W^2];
	Y3=reduce[T W2];
	X3=reduce[(x^3+ECaa x +ECbb)R R-Y3];
	Z3=reduce[Z1 Z2];
	W2=reduce[W2 W];
	{X3,reduce[PowerMod[2,-1,ECp]((Y3-2 X3)R-M W2)],reduce[Z3 W]}
];
(* polsub[{X1_,Y1_,Z1_},{X2_,Y2_,Z2_}]:=poladd[{X1,Y1,Z1},polneg[{X2,Y2,Z2}]];
polaffine[{X_,Y_,Z_}]:=PolynomialMod[{X/reduce[Z^2],Y/reduce[Z^3]},ECp]; *)
polpow[pol_, exp_] := 
	Module[{b = IntegerDigits[exp,2], acc = pol},
		Do[
			 acc=reduce[acc*acc];
			If[b[[i]] == 1,
				acc=reduce[acc * pol]
			],
			{i,2,Length[b]}
		];
		acc
	];

polmult[n_,{X_,Y_}]:=Block[{c,Y2,XX,YY,ZZ,zz,GG,psi,psip,psipp,psim,psimm},
	psi[x_] = reduce[ECpsi[n, x]];
	psip[x_] = reduce[ECpsi[n + 1, x]];
	psipp[x_] = reduce[ECpsi[n + 2, x]];
	psim[x_] = reduce[ECpsi[n - 1, x]];
	psimm[x_] = reduce[ECpsi[n - 2, x]];
	c=PowerMod[4,-1,ECp];
	Y2=reduce[(x^3+ECaa x +ECbb)Y^2];
	ZZ=reduce[If[EvenQ[n],Y2,1]psi[X]];
	XX=reduce[X ZZ^2-Y2 psip[X]psim[X]];
	YY=reduce[c Y If[EvenQ[n],Y2,1](psipp[X]psim[X]^2-psimm[X]psip[X]^2)];
	{XX,YY,ZZ}
]

(* Main function *)

ECSchoof[opts:OptionsPattern[ECSchoof]]:=Block[{primes={2},prod=2,pr=2,res,gcd,schoof,t},
	While[prod < 4Sqrt[ECp],
		pr=NextPrime[pr];
		primes=Append[primes,pr];
		prod *=pr];
	Do[If[prod/primes[[j]]>4Sqrt[ECp],
		prod=prod/primes[[j]];primes=Delete[primes,j]],{j,Length[primes]-1,2,-1}];
	If[OptionValue[Verbose]>0,Print["Schoof: primes=",primes]];
	gcd=PolynomialGCD[x^ECp-x, x^3 + ECaa x + ECbb,Modulus-> ECp];
	If[ Exponent[gcd, x] >= 1, schoof={0}, schoof={1}];
	If[OptionValue[Verbose]>0,Print["Schoof: prime=2, l=",schoof[[1]]]];
	Do[
		l = primes[[j]];
		ss =PolynomialMod[PowerMod[l,-1,ECp]ECpsi[l,x],ECp] ;
		t=Block[{Q,QZ2,QZ3,P0,P1,P2,P3,PZ},
			P0 = {polpow[x, ECp],  polpow[ x^3 + ECaa x + ECbb,( ECp-1)/2],1};
			P1 = {polpow[P0[[1]],ECp],  polpow[ P0[[2]],ECp+1],1};
			P2=polmult[Mod[ECp, l],{x,1}];
			Q=poladd[P1,P2];
			If[Q[[3]]===0, Return[0,Block]];
			P3=P0;
			QZ2=reduce[Q[[3]]^2];QZ3=reduce[Q[[3]] QZ2];
			Do[
				PZ=reduce[P3[[3]]^2];
				If[reduce[P3[[1]]QZ2-Q[[1]]PZ]===0,PZ=reduce[P3[[3]] PZ];If[reduce[P3[[2]]QZ3-Q[[2]]PZ]===0,Return[k,Block],Return[l-k,Block]]];
				P3=poladd[P3,P0]
			,{k,(l-1)/2}]
		];
		If[OptionValue[Verbose]>0,Print["Schoof: prime=",l,", t=",t]];
		If[IntegerQ[t],schoof=Append[schoof,t],Break[]]
	,{j,2,Length[primes]}];
	If[Length[schoof]<Length[primes],Return[$Failed,Block]];
	t=ChineseRemainder[schoof,primes];
	ECp + 1 - If[t^2 > 4ECp, t - prod, t]];




(* Compute Logarithms *)

ECLog::nosol="Warning! `2` is not a multiple of `1`";
ECLog::noinv="Warning! `1` is not invertible modulo `2`";
ECLog::method="Error! Unkonw method.";
ECLog::failedrho="Failed! Could not find a working collission.";

Options[ECLog] = {Method -> "Automatic", Verbose -> 0};

(* Choose the method *)

ECLogSwitch[P_List,Q_List,m_Integer,opts:OptionsPattern[ECLog]]:=Switch[OptionValue[Method],
	"Automatic",If[m<10^2,ECLogNaive[P,Q,m,opts],ECLogBSGS[P,Q,m,opts]],
	"Naive",ECLogNaive[P,Q,m,opts],
	"BSGS",ECLogBSGS[P,Q,m,opts],
	"Rho",If[m<5,ECLogNaive[P,Q,m,opts],ECLogRho[P,Q,m,opts]], (* Pollard rho tends to fail for small values *)
	_,Message[ECLog::method];Null];

(* Naive method for small numbers *)
       
ECLogNaive[P_List,Q_List,m_Integer,opts:OptionsPattern[ECLog]]:=Block[{j=0,PP},
	If[Length[P]==2,PP={\[Infinity],\[Infinity]},PP={0,1,0}];
	If[OptionValue[Verbose]>1,Print["ECLogNaive: P=",P," Q=",Q,", m=",m]];
	While[PP != Q && j < m, j++; PP = ECadd[PP, P]];
	If[j<m,j,Message[ECLog::nosol,a,b,p];Null]
];
ECLogNaive[P_List,Q_List,opts:OptionsPattern[ECLog]]:=ECLogNaive[P,Q,ECOrder[P],opts];

(* Pollard rho algorithm *)

ECLogRho[P_List,Q_List,m_Integer,opts:OptionsPattern[ECLog]]:=Block[{trys=10,j0,j=1,maxj=Floor[10Sqrt[m]],f,x0,y0,x,\[Alpha]x,\[Beta]x,y,\[Alpha]y,\[Beta]y,d,mm},
	If[OptionValue[Verbose]>1,Print["ECLogRho: P=",P," Q=",Q," m=",m]];
	f[{x_,\[Alpha]_,\[Beta]_}]:=Block[{r},
		r=If[x[[1]]!=\[Infinity],Mod[x[[1]],3],1];
		Switch[r,
			0,{ECadd[x, P],Mod[\[Alpha]+1,m],\[Beta]},
			1,{ECmlt[2,x],Mod[2\[Alpha],m],Mod[2\[Beta],m]},
			2,{ECadd[x, Q],\[Alpha],Mod[\[Beta]+1,m]}]];
	While[trys>0,
		If[OptionValue[Verbose]>0,Print["ECLogRho: Searching for a collission"]];
		j0=RandomInteger[{1,m-1}];
		x0=ECadd[Q,ECmlt[j0,P]];y0=x0;
		{x,\[Alpha]x,\[Beta]x}=f[{x0,j0,1}];{y,\[Alpha]y,\[Beta]y}=f[{x,\[Alpha]x,\[Beta]x}];
		While[j<maxj&&x!= y,
			{x,\[Alpha]x,\[Beta]x}=f[{x,\[Alpha]x,\[Beta]x}];
			{y,\[Alpha]y,\[Beta]y}=f[f[{y,\[Alpha]y,\[Beta]y}]];
			j++
		];
		If[OptionValue[Verbose]>0,If[x==y,Print["ECLogRho: Collission after ",j," steps."],Print["ECLogRho: No collission after ",j," steps."]]]
		If[x==y&&Mod[\[Beta]y-\[Beta]x,m]!=0,
			d=GCD[\[Beta]y-\[Beta]x,m];mm=m/d;
			If[Mod[\[Alpha]x-\[Alpha]y,d]!=0,Message[ECLog::nosol,P,Q];Return[Null,Block]];
			l=Mod[(\[Alpha]x-\[Alpha]y)/d PowerMod[(\[Beta]y-\[Beta]x)/d,-1,mm],mm];
			While[ECmlt[l,P]!=Q,l=l+mm];
			Return[l,Block]
			];
		trys--;
	];
	Message[ECLog::failedrho];
	Null
];
ECLogRho[P_List,Q_List,opts:OptionsPattern[ECLog]]:=ECLogRho[P,Q,ECOrder[P],opts];

(* Shanks baby-step/giant-step algorithm *)

ECLogBSGS[P_List,Q_List,q_Integer,opts:OptionsPattern[ECLog]]:=Block[{m=Ceiling[Sqrt[q]],bsteps=CreateDataStructure["HashTable"],Pinf,bs=P,c,gs=Q},
	If[OptionValue[Verbose]>1,Print["ECLogBSGS: P=",P," Q=",Q," m=",m]];
	If[OptionValue[Verbose]>0,Print["ECLogBSGS: Computing baby steps."]];
	If[Length[P]==2,Pinf={\[Infinity],\[Infinity]},Pinf={0,1,0}];
	bsteps["Insert",Pinf->0];
	bsteps["Insert",P->1];
	Do[bs=ECadd[P,bs];bsteps["Insert",bs->j],{j,2,m-1}];
	If[OptionValue[Verbose]>0,Print["ECLogBSGS: Computing giant steps."]];
	c=ECmlt[-m,P];
	Do[
		If[bsteps["KeyExistsQ",gs],Return[bsteps["Lookup",gs]+m j,Block]];
		gs=ECadd[gs,c]
	,{j,0,m-1}]
	Message[ECLog::nosol,a,b,p];
];
ECLogBSGS[P_List,Q_List,opts:OptionsPattern[ECLog]]:=ECLogBSGS[P,Q,ECOrder[P],opts];

(* Pohlig-Hellman reduction for prime powers q^k *)

ECLogPH[P_List,Q_List,q_Integer,k_Integer,opts:OptionsPattern[ECLog]]:=Block[{P1,Pj,Qj,xj,yj},
	If[OptionValue[Verbose]>0,Print["ECLogPH: P=",P," Q=",Q," q=",q," k=",k]];
	If[k==1||q^k<10^4,Return[ECLogSwitch[P,Q,q^k,opts]]];
	Pj=ECmlt[q^(k-1),P];P1=Pj;
	Qj=ECmlt[q^(k-1),Q];
	xj=ECLogSwitch[P1,Qj,q,opts];
	Do[
		Pj=ECmlt[q^(k-j),P];
		Qj=ECmlt[q^(k-j),Q];
		yj=ECLogSwitch[P1,ECadd[Qj,ECmlt[-xj,Pj]],q,opts];
		If[yj==Null,Return[Null,Block]];
		xj=Mod[xj+q^(j-1) yj,q^j];
	,{j,2,k}];
	xj
];

(* Pohlig-Hellman reduction for coprime factos using CRT *)

ECLog[P_List,Q_List,opts:OptionsPattern[ECLog]]:=Block[{m,mf,Pinf,Pj,Qj,pj,kj,l,ll={},mm={}},
	m=ECOrder[P];
	If[OptionValue[Verbose]>0,Print["ECLog: ord(P)=",m]];
	If[Length[P]==2,Pinf={\[Infinity],\[Infinity]},Pinf={0,1,0}];
	If[ECmlt[m,Q]!=Pinf,Return[Message[ECLog::nosol,P,Q];Null,Block]];
	mf=FactorInteger[m];
	If[Length[mf]==1,Return[ECLogPH[P,Q,mf[[1,1]],mf[[1,2]],opts]]];
	Do[
		{pj,kj}=mf[[j]];
		Pj=ECmlt[m/pj^kj,P];
		Qj=ECmlt[m/pj^kj,Q];
		l=ECLogPH[Pj,Qj,pj,kj,opts];
		If[l==Null,Return[Null,Block]];
		mm=Append[mm,pj^kj];
		ll=Append[ll,l];
	,{j,Length[mf]}];
	ChineseRemainder[ll,mm]
];

End[]


Protect[{ECset,ECget,ECneg,ECadd,ECdbl,ECsub,ECmlt,ECaffine,ECmpc,ECOrder,ECGroupOrder,ECPointQ,ECPointShort,ECPointLong,ECPointToHex,ECHexToPoint,ECRandomPoint,ECPsi,ECLog,ECLogBSGS}];

EndPackage[]
