<<FiniteFields`;
(* Initialize Elliptic Curve Params y^2=x^3+a*x+b *)
p=39402006196394479212279040100143613805079739270465446667948293404245721771496870329047266088258938001861606973112319;
a=GF[p][{-3}];
b=GF[p][{16^^b3312fa7e23ee7e4988e056be3f82d19181d9c6efe8141120314088f5013875ac656398d8a2ed19d2a85c8edd3ec2aef}];
(* Large prime order *)
g=39402006196394479212279040100143613805079739270465446667946905279627659399113263569398956308152294913554433653942643;
(* Lucas sequence *)
Lucas[P_,Q_,k_]:=( 
delta=P*P-4*Q;
U=1;
V=P;
For[i==BitLength[k]-1,i>0,i--,
{U=Mod[U*V,p],
V=Mod[Floor[(V*V+delta*U*U)/2],p],
If[IntegerString[k,2][[i]]==1,
{U=Mod[Floor[(P*U+V)/2],p],
V=Mod[Floor[(P*V+delta*U)/2],p]}]
}
];
Return[{U,V}];
);
(* Square equation solution *)
SquareSolve[b_,p_]:=( 
z=0;
Which[
Mod[p,4]==3,{
k=Floor[p/4],
z=PowerMod[b,k+1,p]
},
Mod[p,8]==5,{
k=Floor[p/8],
z1=PowerMod[b,2*k+1,p],
If[z1==1,z=PowerMod[b,k+1,p]],
z=PowerMod[4*b,k+1,p],
If[Floor[z/2]==1,
z=Floor[(z+p)/2],
z=Floor[z/2]]
},
Mod[p,4]==1,{
k=Floor[p/4],
Q=b,P=0,
While[
z==0,
{
P=RandomInteger[{p-1}],
U=Lucas[P,Q,2*k+1][[1]],
V=Lucas[P,Q,2*k+1][[2]],
z=Floor[U/2],
If[PowerMod[z,2,p]==b,Break],
If[V==P-1,z="No square roots exist",Break]
}
]
}
];
Return[z];
);
(* generating random Point *)
randomPoint[p_]:=( 
x=RandomInteger[{p-1}];
alpha=x^3+ToElementCode[a]*x+ToElementCode[b] ;
If[alpha==0,Return[{x,0}]];
beta=SquareSolve[alpha,p ];
mu=RandomInteger[{0,1}];
y=(-1)^mu*beta;
Return[{x,y}];
);
(* adding two Points *)
addPoints[P1_,P2_]:=( 
x=Mod[PowerMod[P1[[2]]-P2[[2]],2,p]*PowerMod[P1[[1]]-P2[[1]],-2,p]-P1[[1]]-P2[[1]],p];
y=Mod[(2*P1[[1]]+P2[[1]])*(P2[[2]]-P1[[2]])*PowerMod[P2[[1]]-P1[[1]],-1,p]-PowerMod[P1[[2]]-P2[[2]],3,p]*PowerMod[P1[[1]]-P2[[1]],-3,p]-P1[[2]],p];
Return[{x,y}];
);
(* Doubling Point *)
dblPoint[G_]:=( 
lambda=Mod[(3*PowerMod[G[[1]],2,p]+ToElementCode[a])*PowerMod[2*G[[2]],-1,p],p];
x1=Mod[PowerMod[lambda,2,p]-2*G[[1]],p];
y1=Mod[lambda*(G[[1]]-x1)-G[[2]],p];
Return[{x1,y1}];
);
(* Inverting Point *)
invertPoint[P_]:=(
x=P[[1]];
y=-P[[2]];
Return[{x,y}];
);
(* Scalar multiplying Point *)
sclmulPoint[P_,k_]:=( 
If[k==0,Return[{0}]];
S=P;
For[i==BitLength[k]-1,i>1,i--,
{
S=dblPoint[S,p],
If[IntegerString[k,2][[i]]==1,S=addPoints[S,P,p]],
}
];
Return[S];
);
(* Checking if Point belongs to curve y^2 = x^3 + a*x + b in the field GF[p] *)
checkPoint[P_]:=( 
str=" ";
If[Mod[P[[2]]^2,p]==Mod[P[[1]]^3+ToElementCode[a]*P[[1]]+ToElementCode[b],p],str="Point belongs to curve",str="Point doesn't belong to curve"];
Return[str]
);

(* Generating 5 random Points *)
P1=randomPoint[p];
P2=randomPoint[p];
P3=randomPoint[p];
P4=randomPoint[p];
P5=randomPoint[p];
Print["P1:(",BaseForm[P1[[1]],16]," ",BaseForm[P1[[2]],16],")"];
checkPoint[P1]
Print["P2:(",BaseForm[P2[[1]],16]," ",BaseForm[P2[[2]],16],")"];
checkPoint[P2]
Print["P3:(",BaseForm[P3[[1]],16]," ",BaseForm[P3[[2]],16],")"];
checkPoint[P3]
Print["P4:(",BaseForm[P4[[1]],16]," ",BaseForm[P4[[2]],16],")"];
checkPoint[P4]
Print["P5:(",BaseForm[P5[[1]],16]," ",BaseForm[P5[[2]],16],")"];
checkPoint[P5]

(* Adding five pairs of generated points *)
P1P2=addPoints[P1,P2];
P2P4=addPoints[P2,P4];
P3P1=addPoints[P3,P1];
P4P5=addPoints[P4,P5];
P5P3=addPoints[P5,P3];
Print["(P1+P2):(",BaseForm[P1P2[[1]],16]," ",BaseForm[P1P2[[2]],16],")"];
checkPoint[P1P2]
Print["(P2+P4):(",BaseForm[P2P4[[1]],16]," ",BaseForm[P2P4[[2]],16],")"];
checkPoint[P2P4]
Print["(P3+P1):(",BaseForm[P3P1[[1]],16]," ",BaseForm[P3P1[[2]],16],")"];
checkPoint[P3P1]
Print["(P4+P5):(",BaseForm[P4P5[[1]],16]," ",BaseForm[P4P5[[2]],16],")"];
checkPoint[P4P5]
Print["(P5+P3):(",BaseForm[P5P3[[1]],16]," ",BaseForm[P5P3[[2]],16],")"];
checkPoint[P5P3]

(* Doubling five generated points *)
P1dbl=dblPoint[P1];
P2dbl=dblPoint[P2];
P3dbl=dblPoint[P3];
P4dbl=dblPoint[P4];
P5dbl=dblPoint[P5];
Print["2*P1:(",BaseForm[P1dbl[[1]],16]," ",BaseForm[P1dbl[[2]],16],")"];
checkPoint[P1dbl]
Print["2*P2:(",BaseForm[P2dbl[[1]],16]," ",BaseForm[P2dbl[[2]],16],")"];
checkPoint[P2dbl]
Print["2*P3:(",BaseForm[P3dbl[[1]],16]," ",BaseForm[P3dbl[[2]],16],")"];
checkPoint[P3dbl]
Print["2*P4:(",BaseForm[P4dbl[[1]],16]," ",BaseForm[P4dbl[[2]],16],")"];
checkPoint[P4dbl]
Print["2*P5:(",BaseForm[P5dbl[[1]],16]," ",BaseForm[P5dbl[[2]],16],")"];
checkPoint[P5dbl]

(* Inverting generated points *)
P1inv=invertPoint[P1];
P2inv=invertPoint[P2];
P3inv=invertPoint[P3];
P4inv=invertPoint[P4];
P5inv=invertPoint[P5];
Print["-P1:(",BaseForm[P1inv[[1]],16]," ",BaseForm[P1inv[[2]],16],")"];
checkPoint[P1inv]
Print["-P2:(",BaseForm[P2inv[[1]],16]," ",BaseForm[P2inv[[2]],16],")"];
checkPoint[P2inv]
Print["-P3:(",BaseForm[P3inv[[1]],16]," ",BaseForm[P3inv[[2]],16],")"];
checkPoint[P3inv]
Print["-P4:(",BaseForm[P4inv[[1]],16]," ",BaseForm[P4inv[[2]],16],")"];
checkPoint[P4inv]
Print["-P5:(",BaseForm[P5inv[[1]],16]," ",BaseForm[P5inv[[2]],16],")"];
checkPoint[P5inv]

(* Scalar multilying of five generated points with random numbers from field GF[p] *)
n1=RandomInteger[(p-1)];
P1mul=sclmulPoint[P1,n1];
n2=RandomInteger[(p-1)];
P2mul=sclmulPoint[P2,n2];
n3=RandomInteger[(p-1)];
P3mul=sclmulPoint[P3,n3];
n4=RandomInteger[(p-1)];
P4mul=sclmulPoint[P4,n4];
n5=RandomInteger[(p-1)];
P5mul=sclmulPoint[P5,n5];
Print["n1 = ",BaseForm[n1,16],"nP1:(",BaseForm[P1mul[[1]],16]," ",BaseForm[P1mul[[2]],16],")"];
checkPoint[P1mul]
Print["n2 = ",BaseForm[n2,16],"nP2:(",BaseForm[P2mul[[1]],16]," ",BaseForm[P2mul[[2]],16],")"];
checkPoint[P2mul]
Print["n3 = ",BaseForm[n3,16],"nP3:(",BaseForm[P3mul[[1]],16]," ",BaseForm[P3mul[[2]],16],")"];
checkPoint[P3mul]
Print["n4 = ",BaseForm[n4,16],"nP4:(",BaseForm[P4mul[[1]],16]," ",BaseForm[P4mul[[2]],16],")"];
checkPoint[P4mul]
Print["n5 = ",BaseForm[n5,16],"nP5:(",BaseForm[P5mul[[1]],16]," ",BaseForm[P5mul[[2]],16],")"];
checkPoint[P5mul]

(* Generating 5 random Points of a large prime order *)
P1lpo=randomPoint[g];
P2lpo=randomPoint[g];
P3lpo=randomPoint[g];
P4lpo=randomPoint[g];
P5lpo=randomPoint[g];
Print["P1':(",BaseForm[P1lpo[[1]],16]," ",BaseForm[P1lpo[[2]],16],")"];
checkPoint[P1]
Print["P2':(",BaseForm[P2lpo[[1]],16]," ",BaseForm[P2lpo[[2]],16],")"];
checkPoint[P1]
Print["P3':(",BaseForm[P3lpo[[1]],16]," ",BaseForm[P3lpo[[2]],16],")"];
checkPoint[P1]
Print["P4':(",BaseForm[P4lpo[[1]],16]," ",BaseForm[P4lpo[[2]],16],")"];
checkPoint[P1]
Print["P5':(",BaseForm[P5lpo[[1]],16]," ",BaseForm[P5lpo[[2]],16],")"];
checkPoint[P1]

(* Doubling point P1 *)
P12=dblPoint[P1]
(* Multiplying point P1 by 5 *)
P15=sclmulPoint[P1,5];
(*Print["2*P1:(",BaseForm[P12[[1]],16],"; ",BaseForm[P12[[2]],16],")"];*)
Print["5*P1:(",BaseForm[P15[[1]],16],"; ",BaseForm[P15[[2]],16],")"];