ECPoints[b_,c_,n_]:=Module[{PtLst, i, j, lhs, rhs},
PtLst={};
For[i=0,i<n,i++,
rhs=Mod[i^3+b*i+c,n];
For[j=0,j<n,j++,
lhs=Mod[j^2,n];
If[rhs==lhs,PtLst=Append[PtLst,List[i,j]]];
]
];
PtLst
]

ECAllPoints[b_,c_,n_]:=Module[{PtLst, i, j, lhs, rhs},
PtLst={};
For[i=0,i<n,i++,
rhs=Mod[i^3+b*i+c,n];
For[j=0,j<n,j++,
lhs=Mod[j^2,n];
If[rhs==lhs,PtLst=Append[PtLst,List[i,j]]];
]
];
PtLst=Append[PtLst,List[Infinity,Infinity]];
PtLst
]

ECOrder[b_,c_,n_]:=Module[{t, i, j, lhs, rhs},
t=1;
For[i=0,i<n,i++,
rhs=Mod[i^3+b*i+c,n];
For[j=0,j<n,j++,
lhs=Mod[j^2,n];
If[rhs==lhs,t++];
]
];
t
]

ECPlot[b_,c_,n_]:=Module[{},
ListPlot[ECPoints[b,c,n]]
]

ECPointOnCurve[b_,c_,n_,pt_]:=Module[{t, i, j, lhs, rhs,res},
i=pt[[1]];
j=pt[[2]];
rhs=Mod[i^3+b*i+c,n];
lhs=Mod[j^2,n];
If[rhs==lhs,res=True,res=False];
res
]

ECPointWithXBF[b_,c_,n_,x_]:=Module[{i, lhs, rhs, res},
res={};
rhs=Mod[x^3+b*x+c,n];
i = 0;
While[i<n,
lhs=Mod[i^2,n];
If[rhs==lhs,res={x,i};i = n];
i++;
];
res
]

ECPointWithYBF[b_,c_,n_,y_]:=Module[{i, j, lhs, rhs, res},
res={};
rhs=Mod[y^2,n];
i = 0;
While[i<n,
lhs=Mod[i^3+b*i+c,n];
If[rhs==lhs,res={i,y};i = n];
i++;
];
res
]

ECGenerateCurveConstant[b_,n_,pt_]:=Module[{},
Mod[pt[[2]]^2-(pt[[1]]^3+b*pt[[1]]),n]
]

ECPointWithX[b_,c_,n_,x_]:=Module[{y,i,rhs, res,sols},
res={};
rhs=Mod[x^3+b*x+c,n];
sols=Values[Solve[y^2==rhs,Modulus->n]];
For[i=1,i≤Length[sols],i++,
y=sols[[i,1]];
res=Append[res,List[x,y]]
];
res
]

ECPointWithY[b_,c_,n_,y_]:=Module[{x,i,lhs, res,sols},
res={};
lhs=Mod[y^2,n];
sols=Values[Solve[x^3+b*x+c==lhs,Modulus->n]];
For[i=1,i≤Length[sols],i++,
x=sols[[i,1]];
res=Append[res,List[x,y]]
];
res
]

ECPointAdd[b_,c_,n_,pt1_,pt2_]:=Module[{x,y,invy,m},
If[pt1==-1,Return[-1]];
If[pt2==-1,Return[-1]];
If[pt1[[1]]==Infinity,Return[pt2]];
If[pt2[[1]]==Infinity,Return[pt1]];
If[pt1==pt2,
If[Mod[2*pt1[[2]],n]==0,Return[List[Infinity,Infinity]]];
If[GCD[2*pt1[[2]],n]>1,Return[-1]];
invy=ModularInverse[2*pt1[[2]],n];
m=Mod[(3*pt1[[1]]^2+b)*invy,n];,
If[Mod[pt1[[1]],n]==Mod[pt2[[1]],n],Return[List[Infinity,Infinity]]];
If[GCD[pt1[[1]]-pt2[[1]],n]>1,Return[-1]];
invy=ModularInverse[pt1[[1]]-pt2[[1]],n];
m=Mod[(pt1[[2]]-pt2[[2]])*invy,n];
];
x=Mod[m^2-pt1[[1]]-pt2[[1]],n];
y=Mod[m*(pt1[[1]]-x)-pt1[[2]],n];
List[x,y]
]

ECPointScalarMult[b_,c_,n_,m_,pt1_]:=Module[{retpt,newy, t, pt},
If[pt1==-1,Return[-1]];
retpt=List[Infinity,Infinity];
t = m;
pt=pt1;
If[t <0,t=-t;newy=Mod[-pt[[2]],n];pt=List[pt[[1]],newy]];
While[t >0,
If[Mod[t,2]==1,retpt=ECPointAdd[b,c,n,retpt,pt]];
pt=ECPointAdd[b,c,n,pt,pt];
t=Floor[t/2];
];
retpt
]

ECPointFactorialScalarMult[b_,c_,n_,m_,pt1_]:=Module[{pt,i},
pt=pt1;
For[i=2,i≤m,i++,pt=ECPointScalarMult[b,c,n,i,pt]];
pt
]

ECPointOrder[b_,c_,n_,pt1_]:=Module[{t, pt,i,r},
If[pt1==-1,Return[-1]];
retpt=List[Infinity,Infinity];
t = True;
i=1;
While[t ,
pt=ECPointScalarMult[b,c,n,i,pt1];
If[pt==-1,Return[-1]];
If[pt[[1]]==Infinity,r=i;t=False];
i++;
];
r
]