(*****************************************************************************)
(*************** Arithmetik elliptischer Kurven ueber Fp *********************)
(*
** Fuer die Punkte der Kurve y**2 = x**3 + a*x + b ueber dem
** Koerper Fp = Z/p (p prim) werden affine Koordinaten
** (x,y) mit 0 <= x < p, 0 <= y < p benuetzt. Der einzige unendlich
** ferne Punkt der Kurve ist das neutrale Element und wird durch
** (-1,-1) dargestellt. Dieser Punkt wird auch durch die
** Konstante Origin abgekuerzt.
*)
const
        Origin = (-1,-1);
end;
(*---------------------------------------------------------------------------*)
(*
** ecp_add berechnet den Punkt P+Q auf der elliptischen  Kurve
**                y**2 = x**3 + a*x + b
** ueber dem Koerper Z/p.
** Affine Koordinaten (x,y); der unendlich ferne Punkt (Nullpunkt
** bzgl. der Gruppenstruktur) wird durch (-1,-1) dargestellt.
*)
(*---------------------------------------------------------------------------*)
function ecp_add(p,a,b: integer; P,Q: array[2]) : array[2];
external
        Origin: const;
var
        m,m1,x,x1,x2,y,y1,y2: integer;
begin
        if Q = Origin then
	       return P;
        elsif P = Origin then
	       return Q; 
        end;
        x1 := P[0]; y1 := P[1];
        x2 := Q[0]; y2 := Q[1];
        if P = Q then
               m := (3 * x1*x1 + a) mod p;
               m1 := mod_inverse(2 * y1,p);
	else
               m := (y1 - y2) mod p;
               m1 := mod_inverse(x1 - x2,p);
	end;
        if m1 = 0 then
               return Origin;
	end;
        m := m*m1;
        x := (m*m - x1 - x2) mod p;
        y := (m * (x1 - x) - y1) mod p;
        return (x,y);
end.
(*---------------------------------------------------------------------------*)
(*
** Berechnet k*P auf der elliptischen Kurve y**2 = x**3+a*x+b mod p
** (k >= 0)
*)
function ecp_mult(p,a,b: integer; P: array[2]; k: integer): array[2];
external
        Origin: const;
var
        i: integer;
        Q: array[2];
begin  
        if k = 0 then
            return Origin;
        end;
        Q := P;
        for i := bit_length(k)-2 to 0 by -1 do
            Q := ecp_add(p,a,b,Q,Q);
            if bit_test(k,i) then
                Q := ecp_add(p,a,b,P,Q);
            end;
        end;
        return Q;
end.
(*--------------------------------------------------------------------------*)
(*
** Berechnet, falls x**3 + a*x + b ein Quadrat mod p,
** eine Wurzel mod p dieses Ausdrucks, i.e. die y-Koordinate 
** des Punktes (x,y) auf der elliptischen Kurve  y**2 = x**3 + a*x + b
** Ist obiger Ausdruck kein Quadrat, wird -1 zurueckgegeben
** Bemerkung: Ist (x,y) ein Punkt auf der Kurve mit y /= 0,
** so auch (x, p-y).
*)
function ecp_x2y(p,a,b,x): integer;
var
        y2: integer;
begin
        y2 := (x**3 + a*x + b) mod p
        if jacobi(y2,p) = -1 then
            return -1;
        elsif y2 = 0 then
            return 0;
        else
            return fp_sqrt(p,y2);
        end;
end;
(*----------------------------------------------------------------------*)
(*
** Sucht einen Punkt auf der elliptischen Kurve y**2 = x**3 + a*x + b
** ueber dem Koerper Fp
** Das letzte Argument x ist ein Startwert fuer die x-Koordinate
** des gesuchten Punktes. Falls dieses Argument nicht angegeben,
** oder x < 0, wird ein Zufallswert genommen
*)
function ecp_point(p,a,b: integer; x := -1): array[2];
var
        y: integer;
begin
	if x < 0 then
        	x := random(p);
	end;
        while (y := ecp_x2y(p,a,b,x)) < 0 do
            inc(x);
        end;
        if x >= p then dec(x,p); end;
        return (x,y);
end;
(*-----------------------------------------------------------------*)
(*
** berechnet die Ordnung der elliptischen Kurve
**      y*y = x*x*x + a*x + b
** ueber dem Koerper Fp (p Primzahl)
**
** ACHTUNG: Nur fuer kleine Primzahlen p < 10**4 geeignet!
*)
function ellorder(p,a,b: integer): integer;
var
    x,N: integer;
begin
    N := p + 1;
    for x := 0 to p-1 do
        N := N + jacobi(x*x*x + a*x + b,p);
    end;
    return N;
end;
(****************************************************************************)
(*
** Wurzelziehen modulo p
*)
(*------------------------------------------------------------*)
(*
** Quadratwurzel von a mod p
** Voraussetzung: p ungerade Primzahl, jacobi(a,p) = 1.
*)
function fp_sqrt(p,a: integer): integer;
begin
    if p mod 8 = 5 then
        return fp_sqrt58(p,a);
    elsif p mod 4 = 1 then
        return fp_sqrt14(p,a);
    else (* p = 3 mod 4 *)
        return a ** ((p+1) div 4) mod p;
    end;
end;
(*----------------------------------------------------------*)
(*
** Quadratwurzel von a mod p.
** Voraussetzung: p Primzahl = 5 mod 8 und jacobi(a,p) = 1.
*)
function fp_sqrt58(p,a: integer): integer;
var
    k, x0, x, j: integer;
begin
    k := p div 8;               (* p = 8*k + 5 *)
    x0 := a**k mod p;
    x := x0*a mod p;            (* x = a**(k+1) *)
    if x*x0 mod p = 1 then      (* a**(2*k+1) = 1 mod p *)
        return x;
    else
        j := 2**(2*k+1) mod p;
        return (j*x mod p);
    end;
end;
(*----------------------------------------------------------*)
(*
** Quadratwurzel von a mod p
** Voraussetzung: p Primzahl = 1 mod 4, jacobi(a,p) = 1.
*)
function fp_sqrt14(p,a: integer): integer;
var
    c: integer;
    x: array[2];
begin
    c := 1;
    while jacobi(c*c-a,p) /= -1 do inc(c); end;
    x := fp2_pow(p, c*c-a, (c,1), (p+1) div 2);
    return x[0];
end;
(*----------------------------------------------------------*)
(*
** Potenz x**ex im Koerper Fp(sqrt(D)), jacobi(D,p) = -1
*)
function fp2_pow(p,D: integer; x: array[2]; ex: integer): array[2];
var
    z0, z1, z00, x0, x1, i: integer;
begin
    if ex = 0 then
        return (1, 0);
    else
        z0 := x0 := x[0];
        z1 := x1 := x[1];
    end;
    for i := bit_length(ex)-2 to 0 by -1 do
        z00 := z0;
        z0 := (z0*z0 + D*z1*z1) mod p;
        z1 := 2*z00 * z1 mod p;
        if bit_test(ex,i) then
            z00 := z0;
            z0 := (z0*x0 + D*z1*x1) mod p;
            z1 := (z00*x1 + z1*x0) mod p;
        end;
    end;
    return (z0, z1);
end;
(*****************************************************************)