(***************************************************************)
(*
** ARIBAS code for
** Arithmetic on elliptic curves over GF(2**n)
**
** author: (C) 2004 Otto Forster
**   Email: forster [AT] mathematik [DOT] uni-muenchen [DOT] de
**   WWW:   http://www.mathematik.uni-muenchen.de/~forster
** date of last change:
**   2004-08-16
**
** This code is placed under the GNU general public licence
*)
(*-------------------------------------------------------------*)
(*
Example code:
We work in the field GF(2**53), choose two elements
a,b in GF(2**53) to define the elliptic curve
    y**2 + x*y = x**3 + a*x**2 + b
and construct two random points P,Q an this elliptic
curve and add them.
Then the order of the elliptic curve is calculated.

==> gf2n_init(53).
-: 9_00719_92547_41063

==> a := 2x1.
-: 2x1

==> b := gf2nint(random(2**53)).
-: 2x3_5037_6F62_DDF5

==> P := ec2n_findpoint(a,b).
-: (2x19_128C_D804_AE6F, 2x18_4E93_3052_7869)

==> Q := ec2n_findpoint(a,b).
-: (2x1B_AE2F_2C3E_C8A2, 2x1B_1327_1161_7C4E)

==> R := ec2n_add(a,P,Q).
-: (2x1D_3382_5ECA_E672, 2x14_2DD8_B914_C49E)

==> N := ec2n_order(a,b).
,,
-: 9_00719_91888_99670

==> ec2n_pow(a,P,N).
-: ()

==> ec2n_pow(a,P,N div 2).
-: (2x0, 2x1D_CF4F_B6B4_6B29)

==> ec2n_pow(a,P,N div 5).
-: (2x1E_E1DA_FF64_60ED, 2x18_09EC_8FD0_2E0A)

*)
(***************************************************************)
(*
** A curve over GF(2**n) is given by an equation
**    y**2 + x*y = x**3 + a*x**2 + b
** plus a point at infinity, which represents the
** neutral element of the group. This point is
** denoted by ECOrig
*)
type
    ECPoint = array of gf2nint;
const
    ECOrig = ();
end;
(***************************************************************)
(*
** Addition of two points P,Q on the elliptic curve
**    y**2 + x*y = x**3 + a*x**2 + b
** The coefficient b is implicitly determined by P.
** Q must lie on the same elliptic curve as P.
*)
function ec2n_add(a: integer; P,Q: ECPoint) : ECPoint;
external
    ECOrig: const;
var
    lambda,x,x1,x2,y,y1,y2: gf2nint;
begin
    if Q = ECOrig then
        return P;
    elsif P = ECOrig then
        return Q;
    elsif P = Q then
        return ec2n_dup(a,P);
    end;
    x1 := P[0]; y1 := P[1];
    x2 := Q[0]; y2 := Q[1];
    if x1 = x2 then
        return ECOrig;
    end;
    lambda := (y1 - y2)/(x1 - x2);
    x := lambda**2 + lambda + x1 + x2 + a;
    y := y2 + (x + x2)*lambda + x;
    return (x,y);
end.
(*---------------------------------------------------------*)
(*
** Calculates P + P on the elliptic curve
**    y**2 + x*y = x**3 + a*x**2 + b
*)
function ec2n_dup(a: integer; P: ECPoint) : ECPoint;
external
    ECOrig: const;
var
    lambda,x,x1,y,y1: integer;
begin
    if P = ECOrig then
       return P;
    end;
    x1 := P[0]; y1 := P[1];
    if x1 = 2x0 then
        return ECOrig;
    end;
    lambda := x1 + y1/x1;
    x := lambda**2 + lambda + a;
    y := y1 + (x + x1)*lambda + x;
    return (x,y);
end.
(*---------------------------------------------------------*)
(*
** returns the negative of the point P
*)
function ec2n_neg(P: ECPoint): ECPoint;
external
    ECOrig: const;
var
    x,y: integer;
begin
    if P /= ECOrig then
        P[1] := P[1] + P[0];
    end;
    return P;
end;
(*---------------------------------------------------------*)
(*
** Calculates k*P on the elliptic curve
**    y**2 + x*y = x**3 + a*x**2 + b
** Hypothesis: k >= 0.
*)
function ec2n_pow(a: integer; P: ECPoint; k: integer): ECPoint;
external
    ECOrig: const;
var
    i: integer;
    R: ECPoint;
begin
    if k = 0 then
        return ECOrig;
    elsif k < 0 then
        P := ec2n_neg(P);
        k := -k;
    end;
    R := P;
    for i := bit_length(k)-2 to 0 by -1 do
        R := ec2n_dup(a,R);
        if bit_test(k,i) then
            R := ec2n_add(a,R,P);
        end;
    end;
    return R;
end;
(*-------------------------------------------------------------*)
(*
** returns a solution x in GF(2**n) of the quadratic equation
**    x*x + x = alfa
** Hypothesis: gf2n_trace(alfa) = 0
** Algorithm:
** a) find  beta in GF(2**n)  such that  trace(beta) = 1
** b) calculate the trace polynomial T(beta*X) mod (X**2 + X + alfa)
**    where T(Z) = Z + Z**2 + Z**4 + ... + Z**(2**(n-1))
** The polynomials T(beta*X) and (X**2 + X + alfa) have exactly
** one common zero.
*)
function gf2n_quadsolve(alfa: integer): integer;
var
    beta, a, b, b2: gf2nint;
    i, deg, trc: integer;
begin
	deg := gf2n_degree();
    trc := 0;
    while trc = 0 do
        beta := gf2nint(random(2**deg));
        trc := gf2n_trace(beta);
    end;
    a := 2x0; b := beta;
    for i := 1 to deg-1 do
        b2 := b**2;
        a := a**2 + b2*alfa;
        b := b2 + beta;
    end;
    return a;
end;
(*-------------------------------------------------------------*)
(*
** Given a,b,x, tries to find a coordinate y such that (x,y) lies on
** the elliptic curve
**    y**2 + x*y = x**3 + a*x**2 + b
** In case of success, returns true and the variable parameter y is set
** accordingly;
** if y does not exist, the return value is false
**
** Method: Solve the equation
**    (y/x)**2 + (y/x) = x + a + (b/x**2)
*)
function ec2n_x2y(a,b,x: gf2nint; var y: gf2nint): boolean;
var
    deg: integer;
    x2, alfa, yx: gf2nint;
begin
    if x = 2x0 then
        deg := gf2n_degree();
        y := b ** (2**(deg-1));
        return true;
    end;
    (* here x /= 2x0 *)
    alfa := x + a + b/x**2;
    if gf2n_trace(alfa) /= 0 then
        return false;
    end;
    yx := gf2n_quadsolve(alfa);
    y := yx * x;
    return true;
end;
(*---------------------------------------------------------*)
(*
** Searches a point on the elliptic curve
**    y**2 + x*y = x**3 + a*x**2 + b
** The third optional argument xi determines the x-coordinate
** where the search begins. When xi is not given,
** a random value is used.
** In case of success, a point P is returned;
** in case of failure the empty array () is returned.
*)
function ec2n_findpoint(a,b: gf2nint; xi := -1): ECPoint;
const
    maxtrials = 128;
var
	found: boolean;
    x, y: gf2nint;
    n, deg: integer;
begin
    if xi < 0 then
        deg := gf2n_degree();
        xi := random(2**deg);
    end;
    for n := 1 to maxtrials do
        x := gf2nint(xi);
        found := ec2n_x2y(a,b,x,y);
        if found then
            return (x,y);
        end;
        inc(xi);
    end;
    return ();
end;
(*-------------------------------------------------------------*)
(*
** Counts the number of points on the elliptic curve
**    y**2 + x*y = x**3 + a*x**2 + b
** over GF(2**n)
** Works only for small n
*)
function ec2n_count(a,b: gf2nint): integer;
var
    deg, count: integer;
    x, xi, alfa: integer;
begin
    deg := gf2n_degree();
    count := 2;     (* (0,y) and ECOrig *)
    for xi := 1 to 2**deg-1 do
        x := gf2nint(xi);
        alfa := x + a + b/x**2;
        if gf2n_trace(alfa) = 0 then
            inc(count,2);
        end;
    end;
    return count;
end;
(*----------------------------------------------------------*)
(*
** auxiliary function for the lambda method
** used in ec2n_ordaux
*)
function ini_steps(var KK: array; range: integer): integer;
var
    anz, anz0, i, averg, x: integer;
begin
    averg := isqrt(range);
    anz0 := 2; x := 4;
    while x <= anz0*averg do
        inc(anz0);
        x := x*2;
    end;
    anz := next_prime(2*anz0 - 1);
    KK := alloc(array,anz);
    KK[0] := 1;
    for i := 1 to anz0 do
        KK[i] := KK[i-1]*2;
    end;
    for i := anz0+1 to anz-1 do
        KK[i] := 1 + 2*random(averg);
    end;
    return anz;
end;
(*--------------------------------------------------------*)
(*
** used in ec2n_ordaux
*)
function xhash(x: gf2nint; r: integer): integer;
const
    Modul0 = 7777801;    (* prime *)
var
    xi: integer;
begin
    xi := integer(x) mod Modul0;
    xi := (xi*xi + 2) mod Modul0;
    return (xi mod r);
end;
(*---------------------------------------------------------------------*)
(*
** Calculates a multiple of the order of a point P on the elliptic curve
**    E: y**2 + x*y = x**3 + a*x**2 + b
** over GF(2**n).
** The last optional argument mm is an integer, (default = 2),
** and must satisfy the hypothesis that ord(E) = 0 modulo mm
** Uses lambda (cangaroo) method; complexity grows as 2**(n/4)
*)
function ec2n_ordaux(a: integer; P: ECPoint; mm := 2): integer;
external
    ECOrig: const;
var
    KK: array;
    SS: array of ECPoint;
    deg, q, q1, i, k, anz, nu, mu, range, tsteps, wsteps: integer;
    R0, P2, T, W: ECPoint;
begin
    deg := gf2n_degree();
    q := 2**deg;
    range := isqrt(4*q) div mm;
    tsteps := 2*isqrt(range);
    wsteps := 2*tsteps;
    anz := ini_steps(KK, range);
    P2 := ec2n_pow(a,P,mm);
    if P2 = ECOrig then
        return mm;
    end;
    SS := alloc(array,anz,());
    for i := 0 to anz-1 do
        SS[i] := ec2n_pow(a,P2,KK[i]);
        KK[i] := KK[i]*mm;
    end;
    q1 := (q+1) div mm;
    R0 := ec2n_pow(a,P2,q1);
    T := R0;
    nu := q1 * mm;
    if T = ECOrig then
        return nu;
    end;
    for i := 1 to tsteps do
        k := xhash(T[0],anz);
        nu := nu + KK[k];
        T := ec2n_add(a,T,SS[k]);
        if T = ECOrig then
            return nu;
        end;
    end;
    W := P2;
    mu := mm;
    for i := 1 to wsteps do
        k := xhash(W[0],anz);
        mu := mu + KK[k];
        W := ec2n_add(a,W,SS[k]);
        if W = ECOrig then
            return mu;
        end;
        if W[0] = T[0] then
            if W[1] = T[1] then
                return nu - mu;
            else
                return nu + mu;
            end;
        end;
    end;
    return -1;
end;
(*-----------------------------------------------------*)
(*
** Calculates the order of a point P on an elliptic
** curve over GF(2**n)
*)
function ec2n_pointorder(a: integer; P: ECPoint): integer;
external
    ECOrig: array;
var
    q,i,N: integer;
    vec: array;
begin
    N := ec2n_ordaux(a,P);
    if N < 0 then
        return N;
    end;
    vec := primefactors(N);
    for i := 0 to length(vec)-1 do
        q := vec[i];
        if (q < N) and (ec2n_pow(a,P,N div q) = ECOrig) then
            N := N div q;
            while (q < N) and (N mod q = 0) do
                if ec2n_pow(a,P,N div q) = ECOrig then
                    N := N div q;
                else
                    break;
                end;
            end;
        end;
    end;
    return N;
end;
(*----------------------------------------------------------*)
(*
** Calculates the order of an elliptic curve
**    E: y**2 + x*y = x**3 + a*x**2 + b
** over GF(2**n).
** Runs reasonably fast for n <= 64
*)
function ec2n_order(a,b: gf2nint): integer;
const
    maxtries = 6;
var
    n,q,range,N,Nmax,m,m0,i,d: integer;
    P: ECPoint;
begin
    n := gf2n_degree();
    if n <= 10 then
        return ec2n_count(a,b);
    end;
    q := 2**n;
    range := isqrt(4*q);
    Nmax := q + 1 + range;
    for i := 1 to maxtries do
        P := ec2n_findpoint(a,b);
        m := ec2n_pointorder(a,P);
        if m > 2*range+1 then
            N := (Nmax div m) * m;
            return N;
        elsif m > 0 then
            write('!');
        else
            write(',');
        end;
    end;
    return -1;
end;
(*****************************************************************)
(*
** Returns a list of all prime factors of x,
** where each prime factor is listed one or several
** times according to its multiplicity.
** The function uses the Pollard rho method and
** the quadratic sieve.
** If the function fails, the first element of
** the returned list equals 1
**
** If the optional argument verbose is given a value /=0,
** then the function writes progress reports to the screen
**
** This function works reasonably fast for integers up to
** 40-50 decimal places
*)
function factorlist(x: integer; verbose := 0): array;
var
    st, st1: stack;
    q, y, bound: integer;
    vec: array;
    count: integer;
begin
    x := abs(x);
    if x < 2 then
        return ();
    end;
    q := 2;
    while q := factor16(x,q) do
        stack_push(st,q);
        x := x div q;
        if verbose then
            writeln(q);
        end;
    end;
    if x < 2**32 then
        stack_push(st,x);
        if verbose then
            writeln(x);
        end;
    else
        stack_push(st1,x);
    end;
    while not stack_empty(st1) do
        x := stack_pop(st1);
        if rab_primetest(x) then
            stack_push(st,x);
            if verbose then
                writeln(x);
            end;
        else
            bound := 4*bit_length(x)**2;
            if verbose then
                writeln("trying to factorize ",x," using Pollard rho")
            end;
            y := rho_factorize(x,bound,verbose);
            if y = 0 then
                if verbose then
                    writeln("trying to factorize ",x,
                            " using quadratic sieve");
                end;
                y := qs_factorize(x,verbose);
            end;
            if y <= 1 or y >= x then
                if verbose then
                    writeln("unable to factorize ",x);
                end;
                stack_push(st,1);
                stack_push(st,x);
            else
                if verbose then
                    writeln("found factor ",y);
                end;
                stack_push(st1,x div y);
                stack_push(st1,y);
            end;
        end;
    end;
    vec := stack2array(st);
    return sort(vec);
end;
(*---------------------------------------------------------*)
(*
** Returns a list of all prime factors of x,
** where each prime factor is listed only once.
** If the function fails, the first element of
** the returned list equals 1.
** If the optional argument verbose is given a value /=0,
** then the function writes progress reports to the screen
*)
function primefactors(x: integer; verbose := 0): array;
var
    i,p,p0: integer;
    vec: array;
    st: stack;
begin
    vec := factorlist(x,verbose);
    p0 := 0;
    for i := 0 to length(vec)-1 do
        p := vec[i];
        if p /= p0 then
            stack_push(st,p);
            p0 := p;
        end;
    end;
    return stack2array(st);
end;
(*****************************************************************)