> restart:
> protect(a,b,p,__O);
> d3:=(y2-y1)/(x2-x1):
> x3:=(d3^2-x1-x2):
> y3:=(-y1+d3*(x1-x3)):
> Esum:=unapply( ([x3,y3]), (x1,y1,x2,y2) );
# Probe:
> simplify((x2-x1)^6*(y3^2-x3^3-a*x3-b)):
> subs({seq( y1^(2*i)=(x1^3+a*x1+b)^i , i=1..4 ), seq(
> y2^(2*i)=(x2^3+a*x2+b)^i , i=1..4 ), seq(
> y1^(2*i+1)=y1*(x1^3+a*x1+b)^i , i=1..4 ), seq(
> y2^(2*i+1)=y2*(x2^3+a*x2+b)^i , i=1..4 ) }, % ):
> simplify(%);
> d3:=(3*x1^2+a)/(2*y1):
> x3:=d3^2-2*x1:
> y3:=-y1+d3*(x1-x3):
> Edub:=unapply( [x3,y3], (x1,y1) );
> unassign('d3','x3','y3');
> Esmooth:=proc(E)
>   eval( subs(E, 4*a&^3+27*b&^2 mod p))<>0;
>   evalb(%);
> end:
> `type/ellcurve` := proc(E)
>   type(E,set(name=algebraic)) and Esmooth(E);
> end:
# Prüfe, ob P auf der Kurve E liegt.
> Echeck:=proc(E::ellcurve,P::{name,list})
>   if P=__O then RETURN(true); fi;
>   evalb(eval(subs( E union { x=P[1], y=P[2] }, y^2 mod p=(x^3+a*x+b)
> mod p)))
> end:
> `type/point` := proc(P)
>   if P=__O then RETURN(true); fi;
>   if not type(P,list) then RETURN(false); fi;
>   if nops(P)<>2 then RETURN(false); fi;
>   true;
> end:
> `type/pointon` := proc(P,E)
>   if not type(P,point) then RETURN(false); fi;
>   if type(P,list(numeric)) then RETURN( Echeck(E,P) );  fi;
>   true;
> end:
# Addiere zwei Punkte auf E: P1 + P2.
> Eadd:=proc(E::ellcurve,P1::point,P2::point)
>   local x1,y1,x2,y2;
>   if not type(P1,pointon(E)) then ERROR("This is no point on E; ".P1."
> notin ".E."!" ); fi;
>   if not type(P2,pointon(E)) then ERROR("This is no point on E; ".P2."
> notin ".E."!" ); fi;
>   if P1=__O then RETURN(P2); fi;
>   if P2=__O then RETURN(P1); fi;
>   x1:=P1[1]; y1:=P1[2];
>   x2:=P2[1]; y2:=P2[2];
>   if x1<>x2 then RETURN(eval(subs(E,Esum(x1,y1,x2,y2) mod p))); fi;
>   if y1<>y2 then RETURN(__O); fi;
>   if y1=0 then RETURN(__O); fi;
>   eval(subs(E,Edub(x1,y1) mod p));
> end:  
# Skaliere einen Punkt auf E: K P.
> Escale:=proc(E::ellcurve,K::integer,P::point)
>   local Q, k;
>   if not type(P,pointon(E)) then ERROR("This is no point on E; ".P."
> notin ".E."!" ); fi;
>   k:=K;
>   if k<0 then
>     Escale(E,-k,P);
>     if %=__O then RETURN(%); fi;
>     RETURN(eval(subs(E,[%[1],-%[2]] mod p)));
>   fi;
>   if k=1 then RETURN(P); fi;
>   Escale( E, floor(k/2), P);
>   Q:=Eadd(E,%,%);
>   if k mod 2<>0 then Q:=Eadd(E,Q,P); fi;
>   Q;
> end:
> points:=proc(E::ellcurve)
>   local N, x, y2, L, i, A, B, P;
>   P:=subs(E,p);
>   A:=subs(E,a);
>   B:=subs(E,b);
>   N:=1;
>   L:={__O};
>   for x from 0 to P-1 do
>     y2:=x^3+A*x+B mod P;
>     if y2=0 then
>       N:=N+1;
>       [x,0];
>     elif numtheory[jacobi]( y2, P ) > 0 then
>       N:=N+2;
>       [msolve(y^2=y2,P)]; seq(subs(%[i],[x,y]),i=1..nops(%));
>     else __O;
>     fi;
>     #print(%);
>     L:=L union {%};
>   od;
>   L;
> end:
> count:=proc(E::ellcurve)
>   local N, x, y2, i, A, B, P;
>   P:=subs(E,p);
>   A:=subs(E,a);
>   B:=subs(E,b);
>   N:=1;
>   for x from 0 to P-1 do
>     y2:=x^3+A*x+B mod P;
>     if y2=0 then
>       N:=N+1;
>     elif numtheory[jacobi]( y2, P ) > 0 then
>       N:=N+2;
>     fi;
>   od;
>   N;
> end:
> save "elliptic-curve.m";
> 
> restart:
> protect(p,a,b,__O);
> read "elliptic-curve.m";
> 
# Beispiele
> E:={ p=7, a=1, b=2 }:
> Esmooth(E);
> E:={ p=7, a=1, b=1 }:
> Esmooth(E);
> Echeck(E,[0,1]);
> Echeck(E,[0,0]);
> points(E);
> Eadd(E,[0,1],__O);
> Eadd(E,[0,1],[0,1]);
> Escale(E,4,[0,1]);
> Eadd(E,[0,1],[0,0]);
> Escale(E,4,[0,0]);
> E:={ p=7, a=1, b=2 }:
> Esmooth(E);
> Echeck(E,[0,1]);
> points(E);
> Eadd(E,[0,1],[0,1]);
> Escale(E,4,[0,1]);
>