Help:=proc(): print(`EC(a,b,p), AddE(a,b,p,A,B) `):end:

#EC(a,b,p): all the points on the elliptic curve y^2=x^3+a*x+b
EC:=proc(a,b,p) local x,y,S:
 if not isprime(p) then
   RETURN(FAIL):
 fi:

S:={}:

for x from 0 to p-1 do
 for y from 0 to p-1 do
  if y^2-x^3-a*x-b mod p=0 then
    S:=S union {[x,y]}:
  fi:
 od:
od:

S:

end:


#AddE(a,b,p,A,B) adding the two points A and B on the elliptic curve y^2=x^3+a*x+b
AddE:=proc(a,b,p,A,B) local s,C,xC,yC:

if A[1]<>B[1] then
 s:=(A[2]-B[2])*(A[1]-B[1])&^(-1) mod p:

  xC:=s^2-A[1]-B[1] mod p:
  yC:=(A[2]-s*(A[1]-B[1])) mod p:
    RETURN([xC,yC]):
else
  s:=((3*A[1]^2+a)+a)*(2*A[2])&^(-1) mod p:
  xC:=s^2 -2*A[2] mod p:
  yC:=(A[2]-s*(A[1]-xC)) mod p:
  RETURN([xC,yC]):
fi:

end: