(**********************************************************************)
(*
    file ecpordschoof.ari

    ARIBAS Code 
    zur Berechnung der Punktezahl elliptischer Kurven
    ueber dem Koerper Z/p nach Schoof 
    
    File author: Otto Forster <forster@math.lmu.de>
    Date of last change: 2012-07-11

------------------
Beispielsitzung:

==> p := next_prime(random(2**96)).

working .,,,,,,,,,, probable prime:
-: 4973_67861_26822_48580_45464_43217

==> a := random(100).
-: 31

==> b := random(p).
-: 1689_13746_46130_81806_36284_86952

==> ecp_discr(p,a,b).
-: 4676_37285_58602_61300_73702_16221

==> N := ecp_ordschoof(p,a,b).
Schoof method: using primes L up to 17

    ....
(dazwischen weitere Meldungen)
    ....

trace modulo 510510 is 322131
now entering Pollard's lambda algorithm
two kangaroos: up to 177348 jumps
.,.,.....,,,,.,,,,.,.,.,...,,..,.,.......,,.,.,.,,,.,,, kangaroos escaped
two kangaroos: up to 177348 jumps
.,..,,,,,,...,,,.,.,,,.,,.., trapped after 122931 jumps
-: 4973_67861_26822_38149_75967_64477

==> P := ecp_point(p,a,b).
-: (1174_38176_86365_65233_20816_66449, 4731_38925_57873_90733_21145_48632)

==> ecp_mult(p,a,P,N).
-: (0, -1)

*)
(*********** 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
** (0,-1) dargestellt. Dieser Punkt wird durch Origin abgekuerzt.
** Die elliptische Kurve wird den Funktionen entweder durch die
** Parameter p,a,b oder, wenn b implizit gegeben ist, nur durch
** die Parameter p,a uebergeben.
*)
const
        Origin = (0,-1);
end;
(*--------------------------------------------------------------------*)
(**** Funktionen:

    ecp_add(p,a: integer; P,Q: array[2]) : array[2];
        Addition P+Q
    ecp_dup(p,a: integer; P: array[2]) : array[2];
        Verdopplung 2*P
    ecp_mult(p,a: integer; P: array[2]; k: integer): array[2];
        Multiplikation k*P
    ecp_point(p,a,b: integer; x := -1): array[2];
        Sucht Punkt auf elliptischer Kurve
    ecp_ord0(p,a,b: integer): integer;
        bestimmt Ordnung der elliptischen Kurve fuer kleine p
    ecp_discr(p,a,b: integer): integer;
        berechnet Diskriminante
    factorlist(N: integer): array;
        berechnet Liste der Primfaktoren von N mit Vielfachheiten
*)
(*--------------------------------------------------------------------*)
(*
** ecp_add berechnet den Punkt P+Q auf der elliptischen  Kurve
**                y**2 = x**3 + a*x + b
** ueber dem Koerper Z/p. Der Parameter b kommt nicht als Argument
** der Funktion ecp_add vor, er ist durch die Punkte auf
** der Kurve implizit gegeben.
** Affine Koordinaten (x,y); der unendlich ferne Punkt (Nullpunkt
** bzgl. der Gruppenstruktur) wird durch (0,-1) dargestellt.
*)
(*--------------------------------------------------------------------*)
function ecp_add(p,a: integer; P,Q: array[2]) : array[2];
external
    Origin: const;
var
    m,m1,d,x,x1,x2,y,y1,y2: integer;
begin
    if Q = Origin then
        return P;
    elsif P = Origin then
	   return Q; 
    elsif P = Q then
        return ecp_dup(p,a,P);
    end;
    x1 := P[0]; y1 := P[1];
    x2 := Q[0]; y2 := Q[1];
    d := (x1 - x2) mod p;
    if d = 0 then
        return Origin;
    end;
    m1 := mod_inverse(d,p);
    m := (y1 - y2)*m1 mod p;
    x := (m*m - x1 - x2) mod p;
    y := (m * (x1 - x) - y1) mod p;
    return (x,y);
end.
(*--------------------------------------------------------------------*)
(*
** Berechnet 2*P
*)
function ecp_dup(p,a: integer; P: array[2]) : array[2];
external
    Origin: const;
var
    m,m1,d,x,x1,y,y1: integer;
begin
    if P = Origin then
	    return P;
    end;
    x1 := P[0]; y1 := P[1];
    m := (3 * x1*x1 + a) mod p;
    d := 2*y1 mod p;
    if d = 0 then
        return Origin;
    end;
    m1 := mod_inverse(d,p);
    m := m*m1 mod p;
    x := (m*m - 2*x1) 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
*)
function ecp_mult(p,a: integer; P: array[2]; k: integer): array[2];
external
    Origin: const;
var
    i: integer;
    Q: array[2];
begin  
    if (k = 0) or (P = Origin) then
        return Origin;
    elsif k < 0 then
        P[1] := (-P[1]) mod p;
        k := abs(k);
    end;
    Q := P;
    for i := bit_length(k)-2 to 0 by -1 do
        Q := ecp_dup(p,a,Q);
        if bit_test(k,i) then
            Q := ecp_add(p,a,Q,P);
        end;
    end;
    return Q;
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, y2: integer;
begin
    if x < 0 then
        x := random(p);
    end;
    while true do
        y2 := (x**3 + a*x + b) mod p;            
        if jacobi(y2,p) = -1 then
            inc(x);
        elsif y2 = 0 then
            y := 0;
            break;
        else
            y := gfp_sqrt(p,y2);
            break;
        end;
    end;
    if x >= p then x := x mod 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)
**
** Nur fuer kleine Primzahlen p < 10**5 geeignet!
*)
function ecp_ord0(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;
(*-----------------------------------------------------------------*)
function ecp_discr(p,a,b: integer): integer;
begin
    return (4*a**3 + 27*b**2) mod p;
end;
(*****************************************************************)
(*
** (used by ecp_ordPlam)
*)
function ranidx(x,r1,r0: integer): integer;
begin
    return (x**2 mod r1) mod r0;
end;
(*----------------------------------------------------------*)
(*
** auxiliary function for the lambda method
** used in ecp_ordPlam
*)
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,0);
    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;
(*--------------------------------------------------------*)
(*
** Calculates the order of the elliptic curve
**      y*y = x*x*x + a*x + b  over Z/p
** using the lambda method
** The last argument tr is a pair (c,n); 
** the hypothesis is that the trace modulo n equals c
** Used in connection with Schoof's method
*)
function ecp_ordlamTr(p,a,b: integer; tr: array[2]): integer;
const
    maxtries = 8;
var
    i, k, d, N, N0, r: integer;
    P: array[2];
begin
    N0 := 1;
    d := isqrt(4*p);
    for i := 1 to maxtries do
        P := ecp_point(p,a,b);
        N := ecp_ordPlamTr(p,a,P,tr);
        if N < 0 then
            continue;
        end;
        r := N0 div gcd(N0,N);
        if r > 1 then
            write(';');
            N := N*r;
        end;
        if N > 2*d then
            k := (p+1+d) div N;
            return k*N;
        else
            write(';');
            N0 := N;
        end;
    end;
    return -1;
end;
(*--------------------------------------------------------*)
type
    kanglist = pointer to kangdata;
    kangdata = record
       P: array[2];
       mu: integer;
       next: kanglist;
    end;
end;
(*--------------------------------------------------------------------*)
function kangstore(P: array[2]; mu: integer; var KL: kanglist): kanglist;
var
    KL1: kanglist;
begin
    new(KL1);
    KL1.P := P;
    KL1.mu := mu;
    KL1.next := KL;
    return KL := KL1;
end;
(*--------------------------------------------------------------------*)
function kangsearch(P: array[2]; var KL: kanglist): integer;
var
    ptr: kanglist;
begin
    ptr := KL;
    while ptr /= nil do
        if P = ptr.P then
            return ptr.mu;
        end;
        ptr := ptr.next;
    end;
    return -1;
end;
(*--------------------------------------------------------------------*)
(*
** Calculates a multiple of the order of a point P on the elliptic curve
**    y**2 = x**3 + a*x + b  over Z/p. 
** using the lambda method.
** The third argument is a pair (c1,n1) and it is supposed
** that the order N of the elliptic curve satisfies
**     c mod n1 = c1, where c := (p+1) - N is the trace
** Return value: order of the point or -1 in case of failure
**
** This implementation uses two simultaneously jumping kangaroos
*)
function ecp_ordPlamTr(p,a: integer; P: array[2]; tr: array[2]): integer;
external
    Origin: array;
var
    kang1, kang2: kanglist;
    KK, SS: array;
    n, i, k, r0, r1, s0, nu, mu: integer;
    range, range0, tsteps, c1, n1, mask: integer;
    P1, T, W: array[2];
    found: boolean;
begin
    (c1,n1) := tr;
    range0 := isqrt(4*p);
    if range0 < n1 then
        return (p + 1 - c1)
    end;
    range := 1 + isqrt(4*p) div n1;
    r0 := ini_steps(KK,range);
    SS := alloc(array,r0,(0,0));
    P1 := ecp_mult(p,a,P,n1);
    for i := 0 to r0-1 do
        SS[i] := ecp_mult(p,a,P1,KK[i]);
        KK[i] := n1*KK[i];
    end;
    r1 := next_prime(range + random(range),0);
    tsteps := max(128,3*isqrt(range));
    k := max(bit_length(tsteps)-4,1);
    k := min(k,15);
    mask := 2**k - 1;
    s0 := random(n1);
    nu := p + 1 + s0 - c1;
    T := ecp_mult(p,a,P,nu);
    mu := s0;
    W := ecp_mult(p,a,P,mu);
    writeln("two kangaroos: up to ",2*tsteps," jumps");
    found := false;
    for i := 1 to 2*tsteps do	
        k := ranidx(T[0],r1,r0);
        T := ecp_add(p,a,T,SS[k]);
        nu := nu + KK[k];
        (* always T = nu*P *)
        if bit_and(T[0],mask) = 1 then
            write('.'); flush();
            kangstore(T,nu,kang1);
            n := kangsearch(T,kang2);
            if n >= 0 then
                n := nu - n;
                found := true;
                break;
            end;
        end;
        k := ranidx(W[0],r1,r0);
        W := ecp_add(p,a,W,SS[k]);
        mu := mu + KK[k];
        (* always W = mu*P *)
        if bit_and(W[0],mask) = 1 then
            write(','); flush();
            kangstore(W,mu,kang2);
            n := kangsearch(W,kang1);
            if n >= 0 then
                n := n - mu;
                found := true;
                break;
            end;
            if abs(nu - mu - p - 1) > range0 then
                write('!'); flush();
                break;
            end;
        end;
    end;
    if found then
        writeln(" trapped after ",i," jumps");
        return n;
    else
        writeln(" kangaroos escaped");
        return -1;
    end;
end;
(*--------------------------------------------------------*)
(*
** N must be a multiple of the order of the point P
** on the elliptic curve  y**2 = x**3 + a*x + b  over Z/p.
** The function returns the exact order of P
*)
function ecp_exactord(p,a: integer; P: array[2]; N: integer): integer;
external
    Origin: array;
var
    ff: array;
    i, q, q0: integer;
begin
    ff := factorlist(N);
    q0 := 0;
    for i := 0 to length(ff)-1 do
        q := ff[i];
        if q = q0 then
            continue;
        elsif ecp_mult(p,a,P,N div q) = Origin then
            N := N div q;
            q0 := 0;
        else
            q0 := q;
        end;
    end;
    return N;
end;
(*****************************************************************)
(*****************************************************************)
(*
** Erstellt ein Array der Primfaktoren von x,
** wobei jeder Primfaktor so oft aufgezaehlt wird,
** wie es seiner Vielfachheit entspricht
** Verwendet rho_factorize und qs_factorize;
** gut geeignet fuer Zahlen mit bis zu 30-40 Dezimalstellen
**
** Wird das optionale Argument verbose mit einem
** Wert /= 0 angegeben, erfolgen Meldungen
*)
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 := bit_length(x);
            bound := 4*bound**2;
            if verbose then
                writeln("trying to factorize ",x," using Pollard rho")
            end;
            y := rho_factorize(x,bound,0);
            count := 0;
            while (y <= 1 or y >= x) do
                if inc(count) > 2 then
                    writeln("unable to factorize ",x);
                    halt(-1);
                end;
                if verbose then
                    writeln("trying to factorize ",x,
                            " using quadratic sieve");
                end;
                y := qs_factorize(x,0);
            end;
            if verbose then
                writeln("found factor ",y);
            end;
            stack_push(st1,x div y);
            stack_push(st1,y);
        end;
    end;
    vec := stack2array(st);
    return sort(vec);
end;
(**************************************************************)
(*
** Polynomial arithmetic in Fp[X]
** Fp = Z/pZ, p an odd prime number
** A polynomial F(X) = a0 + a1*X + ... + an*X**n
** is represented by an array (a0, a1, ..., an), 0 <= ak < p
** The coefficient an is supposed to be /= 0 modulo p
** The zero polynomial is represented by the empty array ().
*)
(**************************************************************)
(*
** reduces all coefficients of the polynomial F(X) modulo p
** and deletes leading zeroes.
*)
function fpX_red(p: integer; F: array): array;
var
    d: integer;
begin
    F := F mod p;
    d := length(F)-1;
    while d >= 0 and F[d] = 0 do
        dec(d);
    end;
    return F[0..d];
end;
(*------------------------------------------------------------*)
(*
** Addition of two polynomials in Fp[X]
*)
function fpX_add(p: integer; F,G: array): array
begin
    return fpX_red(p,F+G);
end;
(*------------------------------------------------------------*)
(*
** Subtraction of two polynomials in Fp[X]
*)
function fpX_sub(p: integer; F,G: array): array
begin
    return fpX_red(p,F-G);
end;
(*------------------------------------------------------------*)
(*
** multiplies F(X) by a scalar lambda
*)
function fpX_scal(p: integer; F: array; lambda: integer): array;
begin
    if lambda mod p = 0 then
        return ();
    end;
    return ((lambda*F) mod p);
end;
(*------------------------------------------------------------*)
(*
** Product of two polynomials in Fp[X]
*)
function fpX_mult(p:integer; F,G: array): array;
var
    n,m,i,i0,i1,k: integer;
    x: integer;
    R: array;
begin
    n := length(F)-1;
    m := length(G)-1;
    if n<0 or m<0 then
        return ();
    end;
    R := alloc(array,n+m+1);
    for k := 0 to n+m do
        x := 0;
        i0 := max(0,k-m);
        i1 := min(n,k);
        for i := i0 to i1 do
            x := x + F[i]*G[k-i];
        end;
        R[k] := x mod p;
    end;
    return R;
end;
(*------------------------------------------------------------*)
(*
** Square of a polynomial in Fp[X]
*)
function fpX_square(p:integer; F: array): array;
var
    n,m,i,i0,k: integer;
    x: integer;
    R: array;
begin
    n := length(F)-1;
    if n<0 then
        return ();
    end;
    R := alloc(array,2*n+1);
    for k := 0 to 2*n do
        i0 := max(0,k-n);
        m := (k+1) div 2;
        x := 0;
        for i := i0 to m-1 do
            x := x + F[i]*F[k-i];
        end;
        x := 2*x;
        if even(k) then
            x := x + F[m]**2;
        end;
        R[k] := x mod p;
    end;
    return R;
end;
(*------------------------------------------------------------*)
(*
** Divides F(X) by G(X) in Fp[X]
** The leading coefficient of G is supposed to be invertible
** Returns a pair (Q(X),R(X)) such that
**     F(X) = Q(X)*G(X) + R(X)
** and deg R(X) < deg G(X)
*)
function fpX_divide(p: integer; F,G: array): array[2] of array;
var
    Q: array;
    i,j,k,m,a,c,c1: integer;
begin
    k := length(F)-1;
    m := length(G)-1;
    if k < m then 
        return ((),F); 
    end;
    c := G[m];
    if c /= 1 then
        c1 := mod_inverse(c,p);
        if c1 = 0 then
            writeln("error in fpX_divide: leading coeff not invertible");
            halt(-1);
        end;
        G := fpX_scal(p,G,c1); 
    end;
    Q := alloc(array,k-m+1);
    for j := k to m by -1 do
        a := F[j] mod p;
        Q[j-m] := (a*G[m]) mod p;
        if a /= 0 then
            for i := 1 to m do
                F[j-i] := F[j-i] - a*G[m-i];
            end;
        end;
    end;
    if c /= 1 then
        Q := fpX_scal(p,Q,c1);
    end;
    k := m-1;
    for j := 0 to k do
        F[j] := F[j] mod p;
    end;
    while k >= 0 and F[k] = 0 do
        dec(k);
    end;
    return (Q,F[0..k]);
end;
(*------------------------------------------------------------*)
(*
** Divides F(X) by G(X) in Fp[X]
** The leading coefficient of G is supposed to be invertible
** The result Q(X) satisfies
**     deg(F(X) - Q(X)*G(X)) < deg G(X)
** Stripped down version of fpX_divide
*)
function fpX_div(p: integer; F,G: array): array;
var
    Q: array;
    i,j,k,m,a,c,c1: integer;
begin
    k := length(F)-1;
    m := length(G)-1;
    if k < m then 
        return (); 
    end;
    c := G[m];
    if c /= 1 then
        c1 := mod_inverse(c,p);
        if c1 = 0 then
            writeln("error in fpX_div: leading coeff not invertible");
            halt(-1);
        end;
        G := fpX_scal(p,G,c1); 
    end;
    Q := alloc(array,k-m+1);
    for j := k to m by -1 do
        a := F[j] mod p;
        Q[j-m] := (a*G[m]) mod p;
        if a /= 0 then
            for i := 1 to m do
                F[j-i] := F[j-i] - a*G[m-i];
            end;
        end;
    end;
    if c /= 1 then
        Q := fpX_scal(p,Q,c1);
    end;
    return Q;
end;
(*------------------------------------------------------------*)
(*
** Calculates F(X) mod G(X) in Fp[X]
** The leading coefficient of G is supposed to be invertible
** Stripped down version of fpX_divide
*)
function fpX_mod(p: integer; F,G: array): array;
var
    i,j,k,m,a,c: integer;
begin
    k := length(F)-1;
    m := length(G)-1;
    if k < m or m < 0 then 
        return F; 
    end;
    c := G[m];
    if c /= 1 then
        c := mod_inverse(c,p);
        if c = 0 then
            writeln("error in fpX_mod: leading coeff not invertible");
            halt(-1);
        end;
        G := fpX_scal(p,G,c);
    end;
    for j := k to m by -1 do
        a := F[j] mod p;
        if a /= 0 then
            for i := 1 to m do
                dec(F[j-i],a*G[m-i]);
            end;
        end;
    end;
    return fpX_red(p,F[0..m-1]);
end;
(*------------------------------------------------------------*)
function fpX_gcd(p: integer; P,Q: array): array;
var
    F: array;
    m, k, c: integer;
begin
    m := length(Q);
    while m > 0 and (Q[m-1] mod p = 0) do
        dec(m);
    end;
    while length(Q) > 0 do
        F := Q;
        Q := fpX_mod(p,P,Q);
        P := F;
    end;
    k := length(P);
    if k > 0 then
        c := mod_inverse(P[k-1],p);
        P := fpX_scal(p,P,c);
    end;
    return P;
end;
(*------------------------------------------------------------*)
(*
** modular power F(X)**ex mod Q(X)
*)
function fpX_modpow(p: integer; F: array; ex: integer; Q: array): array;
var
    i: integer;
    R: array;
begin
    if ex < 0 then
        writeln("error");
        writeln("power: exponent must be non-negative: ",ex);
        return ();
    elsif ex = 0 then
        return {1};
    else
        R := F;
    end;
    for i := bit_length(ex)-2 to 0 by -1 do
        R := fpX_mod(p,fpX_square(p,R),Q);
        if bit_test(ex,i) then
            R := fpX_mod(p,fpX_mult(p,R,F),Q);
        end;
    end;
    return R;
end;
(*------------------------------------------------------------*)
(*
** Evaluates F(x) mod p.
*)
function fpX_val(p: integer; F: array; x: integer): integer;
var
    i,z: integer;
begin
    z := 0;
    for i := length(F)-1 to 0 by -1 do
        z := (z*x + F[i]) mod p;
    end;
    return z;
end;
(*-------------------------------------------------------------*)
(*
** Calculates the gcd of F(X) and F'(X)
*)
function fpX_discr(p: integer; F: array): array;
var
    F1: array;
    k, deg: integer;
begin
    deg := length(F)-1;
    if deg <= 1 then
        return {1};
    end;
    F1 := alloc(array,deg);
    for k := 1 to deg do
        F1[k-1] := k*F[k];
    end;
    return fpX_gcd(p,F,F1);
end;
(**************************************************************)
(*
** Arithmetic in the fields Fpn = GF(p**n), p odd prime
*)
(**************************************************************)
(*
** A field Fpn can be defined as Fp[X]/(Q(X))
** where Q(X) is an irreducible polynomial of degree n.
** The elements of the field are represented by
** polynomials of degree <= n-1
** A polynomial F(X) = a0 + a1*X + a2*X**2 + ... + ak*X**k
** is represented by the vector (a0,a1,...,ak)
** The highest coefficient ak is supposed to be /= 0
** The zero polynomial is represented by the empty array ().
*)
(*------------------------------------------------------------*)
(*
** Q is the polynomial defining the field.
** If Q is not irreducible, only a ring is obtained.
** If in this latter case, during calculating an inverse,
** a factor of Q is found, it is stored in Q1
*)
type
    _Fieldpn = record
        p,deg: integer;
        Q: array;
        Q1: array;
    end;
end;
(*------------------------------------------------------------*)
(*
** Initialize the field Fp[X]/(Q(X))
** The highest coefficient of Q must be invertible mod p
** If Q(X) is not irreducible, the result
** is not a field, but only a ring
** Returns the degree
*)
function fpn_init(var field: _Fieldpn; p: integer; Q: array): integer;
var
    deg: integer;
    lambda: integer;
begin
    deg := length(Q)-1;
    field.p := p;
    field.deg := deg; 
    lambda := Q[deg];
    if lambda /= 1 then
        lambda := mod_inverse(lambda,p);
        Q := fpX_scal(p,Q,lambda);
    end;
    field.Q := Q;
    field.Q1 := {1};
    return deg;
end;
(*------------------------------------------------------------*)
function fpn_update(var field: _Fieldpn): integer;
var
    Q1: array;
    p,c: integer;
    deg1: integer;
begin
    p := field.p;
    Q1 := field.Q1;
    deg1 := length(Q1) - 1;
    field.Q1 := {1};
    if (deg1 <= 0) or (deg1 >= field.deg) or
        fpX_mod(p,field.Q,Q1) /= () then
        return field.deg;
    end;
    if deg1 > (field.deg div 2) then
        Q1 := fpX_div(p,field.Q,Q1);
        deg1 := length(Q1)-1;
    end;
    c := Q1[deg1];
    Q1 := fpX_scal(p,Q1,mod_inverse(c,p));
    field.Q := Q1;
    field.deg := deg1;
    return deg1;
end;
(*------------------------------------------------------------*)
(*
** Retrieve the defining polynomial of the field
*)
function fpn_fieldpol(var field: _Fieldpn): array;
begin
    return field.Q; 
end;
(*------------------------------------------------------------*)
function fpn_mult(var field: _Fieldpn; F,G: array): array;
var
    p: integer;
    R: array;
begin
    p := field.p;
    R := fpX_mult(p,F,G);
    return fpX_mod(p,R,field.Q);
end;
(*------------------------------------------------------------*)
function fpn_square(var field: _Fieldpn; F: array): array;
var
    p: integer;
    R: array;
begin
    p := field.p;
    R := fpX_square(p,F);
    return fpX_mod(p,R,field.Q);
end;
(*------------------------------------------------------------*)
(*
** Calculates the inverse of xi in the field Fpn,
** using the extended Euclidean algorithm
*)
function fpn_inverse(var field: _Fieldpn; F: array): array;
var
    G, Q1, Q2, Qtemp, P, R: array;
    p, c: integer;
    deg1: integer;
begin
    G := field.Q;
    p := field.p;
    Q1 := (); Q2 := {1}; 
    while F /= () do
        (P,R) := fpX_divide(p,G,F);
        Qtemp := Q2;
        Q2 := fpX_sub(p,Q1,fpX_mult(p,P,Q2));
        Q1 := Qtemp;
        G := F;
        F := R;
    end;
    deg1 := length(G) - 1;
    if deg1 /= 0 then
        if deg1 >= 1 and deg1 < field.deg then
            writeln("error in fpn_inverse: field polynomial is reducible");
            field.Q1 := G;
        end;
        return ();
    else
        c := mod_inverse(G[0],p);
        return fpX_scal(p,Q1,c);
    end;
end;
(*------------------------------------------------------------*)
function fpn_power(var field: _Fieldpn; F: array; ex: integer): array;
begin
    return fpX_modpow(field.p,F,ex,field.Q);
end;
(*------------------------------------------------------------*)
(*
** Evaluates F(xi) in Fpn, where F in Fp[X], xi in Fpn
*)
function fpn_polval(var field: _Fieldpn; F,xi: array): array;
var
    p, ak: integer;
    k, deg: integer;
    Z: array;
begin
    p := field.p;
    deg := length(F)-1;
    if deg < 0 then
        return ();
    elsif length(xi) <= 1 then
        if xi = () then
            return {F[0]};
        else
            return fpX_val(p,F,xi[0]);
        end;
    end;
    Z := {F[deg]};
    for k := deg-1 to 0 by -1 do
        Z := fpn_mult(field,Z,xi);
        ak := F[k];
        if Z = () then
            Z := {ak};
        else 
            Z[0] := Z[0] + ak;
        end;
    end;
    return fpX_red(p,Z);
end;
(***************************************************************)
(*
** Irreducibility test of polynomials over Fp
*)
(*-------------------------------------------------------------*)
(*
** Tests whether the polynomial F is irreducible over Fp = GF(p)
** Supposes that the highest coefficient of F is invertible mod p.
**
** This function uses the following theorem:
** A polynomial F(X) in Fp[X] of degree n is irreducible
** if and only if
**    gcd(X**(p**m)-X,F) = 1 for all m with 1 <= m <= n/2
*)
function fpX_testirred(p: integer; F: array): boolean;
var
    n,m,lambda: integer;
    X,Xpm,G: array;
begin
    n := length(F)-1;
    lambda := F[n];
    if lambda /= 1 then
        lambda := mod_inverse(lambda,p);
        F := fpX_scal(p,F,lambda);
    end;
    X := (0,1);
    Xpm := X;
    for m := 1 to (n div 2) do
        Xpm := fpX_modpow(p,Xpm,p,F);
        G := fpX_sub(p,Xpm,X);
        G := fpX_gcd(p,G,F);
        if G /= {1} then
            return false;
        end;
        write('.'); flush();
    end;
    return true;
end;
(**************************************************************)
type
    pnpoint = array[2] of array;
end;
const
    Error = (0,0,-2);
end;
(*------------------------------------------------------------*)
(*
** Addition of two points on the elliptic curve
**     y**2 = x**3 + a*x + b
** The coefficients are in the ground field Fp,
** whereas the points have coordinates in Fpn
**
** The coefficient b is implicitely given by the points P,Q
*)
function ecpn_add(var field: _Fieldpn; a: integer; P,Q: pnpoint): pnpoint;
external
    Origin, Error: const;
var
    p: integer;
    m,m1,m2,d,x,x1,x2,y,y1,y2: array;
begin
    if Q = Error or P = Error then
        return Error;
    end;
    if Q = Origin then
        return P;
    elsif P = Origin then
        return Q; 
    elsif P = Q then
        return ecpn_dup(field,a,P);
    end;
    x1 := P[0]; y1 := P[1];
    x2 := Q[0]; y2 := Q[1];
    p := field.p;
    d := fpX_sub(p,x1,x2);
    if d = () then
        return Origin;
    end;
    m1 := fpn_inverse(field,d);
    if m1 = () then
        return Error;
    end;
    m := fpn_mult(field,m1,fpX_sub(p,y1,y2));
    m2 := fpn_square(field,m);
    x := fpX_sub(p,m2,fpX_add(p,x1,x2));
    y := fpn_mult(field,m,fpX_sub(p,x1,x));
    y := fpX_sub(p,y,y1);
    return (x,y);
end;
(*--------------------------------------------------------------------*)
(*
** Calculates 2*P
*)
function ecpn_dup(var field: _Fieldpn; a: integer; P: pnpoint) : pnpoint;
external
    Origin, Error: const;
var
    p: integer;
    m,m1,m2,d,x,x1,y,y1: array;
begin
    if P = Origin or P = Error then
	    return P;
    end;
    x1 := P[0]; y1 := P[1];
    if y1 = () then
        return Origin;
    end;
    p := field.p;
    m := fpn_square(field,x1);
    m := fpX_scal(p,m,3);
    m := fpX_add(p,m,{a});	(* m = 3*x1**2 + a *)
    d := fpX_scal(p,y1,2);
    m1 := fpn_inverse(field,d);
    if m1 = () then
        return Error;
    end;
    m := fpn_mult(field,m,m1);
    m2 := fpn_square(field,m);
    x := fpX_sub(p,m2,fpX_scal(p,x1,2));
    y := fpn_mult(field,m,fpX_sub(p,x1,x));
    y := fpX_sub(p,y,y1);
    return (x,y);
end;
(*--------------------------------------------------------------------*)
(*
** Calculates k*P, where k is an integer
*)
function ecpn_mult(var field: _Fieldpn; a: integer; P: pnpoint; k: integer): 
				pnpoint;
external
    Origin, Error: const;
var
    i: integer;
    Q: pnpoint;
begin  
    if P = Error then
        return Error;
    end;
    if (k = 0) or (P = Origin) then
        return Origin;
    elsif k < 0 then
        k := abs(k);
        P[1] := fpX_scal(field.p,P[1],-1);
    end;
    Q := P;
    for i := bit_length(k)-2 to 0 by -1 do
        Q := ecpn_dup(field,a,Q);
        if bit_test(k,i) then
            Q := ecpn_add(field,a,Q,P);
        end;
    end;
    return Q;
end;
(*--------------------------------------------------------------------*)
function ecpn_frob(var field: _Fieldpn; P: pnpoint): pnpoint;
var
    p: integer;
    x, y: array;
begin
    p := field.p;
    x := fpn_power(field,P[0],p);
    y := fpn_power(field,P[1],p);
    return (x,y);
end;
(*--------------------------------------------------------------------*)
(*
** Calculates -P
*)
function ecpn_neg(var field: _Fieldpn; P: pnpoint): pnpoint;
external
    Origin, Error: const;
begin
    if P = Origin or P = Error then
        ;
    else
        P[1] := fpX_scal(field.p,P[1],-1);
    end;
    return P;
end;
(**************************************************************)
(*
** Twisted versions of ecpn_add, ecpn_dup, ecpn_mult, ecpn_frob
** These functions apply to the elliptic curve
**      t*y**2 = x**3 + a*x + b,
** where a,b in Fp and t in Fpn
** If t doesnt have a square root in Fpn, this is equivalent
** to consider the curve y**2 = x**3 + a*x + b over the
** field Fpn[sqrt(t)] and points (x,y*sqrt(t))
*)
(*------------------------------------------------------------*)
(*
** Addition of two points on the elliptic curve
**     t*y**2 = x**3 + a*x + b
** The coefficients are in the ground field Fp,
** whereas the points have coordinates in Fpn
**
** The coefficient b is implicitly given by the points P,Q
*)
function ecpn_tadd(var field: _Fieldpn; var t: array; a: integer; 
                   P,Q: pnpoint): pnpoint;
external
    Origin, Error: const;
var
    p: integer;
    m,m1,m2,d,x,x1,x2,y,y1,y2: array;
begin
    if P = Error or Q = Error then
        return Error;
    end;
    if Q = Origin then
        return P;
    elsif P = Origin then
        return Q; 
    end;
    if P = Q then
        return ecpn_tdup(field,t,a,P);
    end;
    x1 := P[0]; y1 := P[1];
    x2 := Q[0]; y2 := Q[1];
    p := field.p;
    d := fpX_sub(p,x1,x2);
    if d = () then
        return Origin;
    end;
    m1 := fpn_inverse(field,d);
    if m1 = () then
        return Error;
    end;
    m := fpn_mult(field,m1,fpX_sub(p,y1,y2));
    m2 := fpn_square(field,m);
    m2 := fpn_mult(field,m2,t);     
    (* this is the only place where the twist t appears *)
    x := fpX_sub(p,m2,fpX_add(p,x1,x2));
    y := fpn_mult(field,m,fpX_sub(p,x1,x));
    y := fpX_sub(p,y,y1);
    return (x,y);
end;
(*--------------------------------------------------------------------*)
(*
** Calculates 2*P
*)
function ecpn_tdup(var field: _Fieldpn; var t: array; a: integer; 
                   P: pnpoint) : pnpoint;
external
    Origin, Error: const;
var
    p: integer;
    m,m1,m2,d,x,x1,y,y1: array;
begin
    if P = Origin or P = Error then
	    return P;
    end;
    x1 := P[0]; y1 := P[1];
    if y1 = () then
        return Origin;
    end;
    p := field.p;
    m := fpX_scal(p,fpn_square(field,x1),3);
    m := fpX_add(p,m,{a});          (* m = 3*x1**2 + a *)
    d := fpX_scal(p,y1,2);
    d := fpn_mult(field,d,t);       (* twist *)
    m1 := fpn_inverse(field,d);
    if m1 = () then                 (* inverse does not exist *)
        return Error;
    end;
    m := fpn_mult(field,m,m1);      (* m = (3*x1**2 + a)/(2*t*y1) *)
    m2 := fpn_square(field,m);
    m2 := fpn_mult(field,m2,t);     (* twist *)
    x := fpX_sub(p,m2,fpX_scal(p,x1,2));
    y := fpn_mult(field,m,fpX_sub(p,x1,x));
    y := fpX_sub(p,y,y1);
    return (x,y);
end;
(*--------------------------------------------------------------------*)
(*
** Calculates k*P, where k is an integer
*)
function ecpn_tmult(var field: _Fieldpn; var t: array; a: integer; 
                    P: pnpoint; k: integer): pnpoint;
external
    Origin, Error: const;
var
    i: integer;
    Q: pnpoint;
begin  
    if P = Error then
        return Error;
    end;
    if (k = 0) or (P = Origin) then
        return Origin;
    elsif k < 0 then
        k := abs(k);
        P[1] := fpX_scal(field.p,P[1],-1);
    end;
    Q := P;
    for i := bit_length(k)-2 to 0 by -1 do
        Q := ecpn_tdup(field,t,a,Q);
        if bit_test(k,i) then
            Q := ecpn_tadd(field,t,a,Q,P);
        end;
    end;
    return Q;
end;
(*--------------------------------------------------------------------*)
function ecpn_tfrob(var field: _Fieldpn; var t: array;
                    P: pnpoint): pnpoint;
var
    p: integer;
    x, y, y2, t1: array;
begin
    p := field.p;
    x := fpn_power(field,P[0],p);
    y2 := fpn_square(field,P[1]);
    t1 := fpn_mult(field,y2,t);
    t1 := fpn_power(field,t1,(p-1) div 2);
    y := fpn_mult(field,t1,P[1]);
    return (x,y);
end;
(**************************************************************)
(*
** Division polynomials for elliptic curves
*)
(**************************************************************)
type
    _PsiData = record
        pab: array[3];
        psivec: array of array;
    end;
end;
var
    PsiMax := 97;
    PsiData: _PsiData;
end;
(*-------------------------------------------------------------*)
function Psi_prepare(p,a,b,N: integer): integer;
external
    PsiData: _PsiData;
    PsiMax: integer;
begin
    if PsiData.pab /= (p,a,b) then
        PsiData.pab := (p,a,b);
        PsiData.psivec := alloc(array,PsiMax+1,());
        PsiData.psivec[1] := {1};
        PsiData.psivec[2] := {1};
    end;
    if N > 0 and N <= PsiMax and PsiData.psivec[N] /= () then
        return N;
    else
        return 0;
    end;
end;
(*-------------------------------------------------------------*)
function Psi_retrieve(N: integer): array;
external
    PsiData: _PsiData;
    PsiMax: integer;
begin
    if N < 0 or N > PsiMax then
        return ();
    else
        return PsiData.psivec[N];
    end;
end;
(*-------------------------------------------------------------*)
function Psi_store(N: integer; psi: array): integer;
external
    PsiData: _PsiData;
    PsiMax: integer;
begin
    if N < 0 or N > PsiMax then
        return -1;
    else
        PsiData.psivec[N] := psi;
        return N;
    end;
end;
(*-------------------------------------------------------------*)
(*
** Polynom fuer 3-Teilungspunkte auf elliptischer Kurve
**	       y**2 = x**3 + a*x + b
** ueber dem Koerper Fp
*)
function Psi3(p,a,b: integer): array;
var
    psi: array[5];
begin
    if Psi_prepare(p,a,b,3) = 3 then
        return Psi_retrieve(3);
    end;
    psi[0] := (-a*a) mod p;
    psi[1] := 12*b mod p;
    psi[2] := 6*a mod p;
    psi[3] := 0;
    psi[4] := 3;
    Psi_store(3,psi);
    return psi;
end;
(*------------------------------------------------------------*)
(*
** Das Polynom fuer die Vierteilungspunkte ist
**	Psi4(X,Y) := 2*Y*Psi4bar(X)
*)
function Psi4bar(p,a,b: integer): array;
var
    psi: array[7];
    a2,b2: integer;
begin
    if Psi_prepare(p,a,b,4) = 4 then
        return Psi_retrieve(4);
    end;
    a2 := a*a mod p;
    b2 := b*b mod p;
    psi[0] := (-16*b2-2*a2*a) mod p;
    psi[1] := (-8*a*b) mod p;
    psi[2] := (-10*a2) mod p;
    psi[3] := 40*b mod p;
    psi[4] := 10*a mod p;
    psi[5] := 0;
    psi[6] := 2;
    Psi_store(4,psi);
    return psi;
end;
(*------------------------------------------------------------*)
function Psi5(p,a,b: integer): array;
var
    psi, psi3, psi4bar, Y2, Y4, F, G: array;
begin
    if Psi_prepare(p,a,b,5) = 5 then
        return Psi_retrieve(5);
    end;
    psi3 := Psi3(p,a,b);
    G := fpX_square(p,psi3);
    G := fpX_mult(p,G,psi3);
    psi4bar := Psi4bar(p,a,b);
    Y2 := (b,a,0,1);
    Y2 := fpX_scal(p,Y2,4);
    Y4 := fpX_square(p,Y2);
    F := fpX_mult(p,Y4,psi4bar);
    psi := fpX_sub(p,F,G);
    Psi_store(5,psi);
    return psi;
end;
(*------------------------------------------------------------*)
(*
** Das Polynom fuer die 6-Teilungspunkte ist
**	Psi6(X,Y) := 2*Y*Psi6bar(X)
*)
function Psi6bar(p,a,b: integer): array;
var
    psi, psi3,psi5,psi4bar,G: array;
begin
    if Psi_prepare(p,a,b,6) = 6 then
        return Psi_retrieve(6);
    end;
    psi5 := Psi5(p,a,b);
    psi4bar := Psi4bar(p,a,b);
    G := fpX_square(p,psi4bar);
    G := fpX_sub(p,psi5,G);
    psi3 := Psi3(p,a,b);
    psi := fpX_mult(p,psi3,G);
    Psi_store(6,psi);
    return psi;
end;
(*------------------------------------------------------------*)
(*
** Returns the N-th division polynomial if N is odd
** If N is even, and the returned polynomial is denoted by psi(X),
** then the N-th division polynomial is 2*Y*psi(X).
*)
function PsiN(p,a,b,N: integer): array;
var
    psi: array;
begin
    if Psi_prepare(p,a,b,N) = N then
        return Psi_retrieve(N);
    end;
    if N <= 6 then
        if N = 3 then
            return Psi3(p,a,b);
        elsif N = 4 then
            return Psi4bar(p,a,b);
        elsif N = 5 then
            return Psi5(p,a,b);
        elsif N = 6 then
            return Psi6bar(p,a,b);
        else
            return ();	(* this case should not happen *)
        end;
    end;
    if even(N) then
        psi := PsiNev(p,a,b,N);
    else
        psi := PsiNodd(p,a,b,N);
    end;
    Psi_store(N,psi);
    return psi;
end;
(*------------------------------------------------------------*)
function PsiNev(p,a,b,N: integer): array;
var
    psim,psim1,psim2,psim_1,psim_2,F,G: array;
    m: integer;
begin
    if N < 8 then
        writeln("PsiNev: N must be >= 8");
        halt(-1);
    end;
    m := N div 2;
    psim := PsiN(p,a,b,m);
    psim1 := PsiN(p,a,b,m+1);
    psim2 := PsiN(p,a,b,m+2);
    psim_1 := PsiN(p,a,b,m-1);
    psim_2 := PsiN(p,a,b,m-2);
    F := fpX_square(p,psim_1);
    F := fpX_mult(p,F,psim2);
    G := fpX_square(p,psim1);
    G := fpX_mult(p,G,psim_2);
    F := fpX_sub(p,F,G);
    return fpX_mult(p,F,psim);
end;
(*------------------------------------------------------------*)
function PsiNodd(p,a,b,N: integer): array;
var
    psim,psim1,psim2,psim_1,F,G,Y2,Y4: array;
    m: integer;
begin
    if N < 7 then
        writeln("PsiNev: N must be >= 7");
        halt(-1);
    end;
    m := N div 2;
    psim := PsiN(p,a,b,m);
    psim1 := PsiN(p,a,b,m+1);
    psim2 := PsiN(p,a,b,m+2);
    psim_1 := PsiN(p,a,b,m-1);
    F := fpX_square(p,psim);
    F := fpX_mult(p,F,psim);
    F := fpX_mult(p,F,psim2);
    G := fpX_square(p,psim1);
    G := fpX_mult(p,G,psim1);
    G := fpX_mult(p,G,psim_1);
    Y2 := (b,a,0,1);
    Y2 := fpX_scal(p,Y2,4);
    Y4 := fpX_square(p,Y2);
    if even(m) then
        F := fpX_mult(p,F,Y4);
    else
        G := fpX_mult(p,G,Y4);
    end;
    return fpX_sub(p,F,G);
end;
(**************************************************************)
(*
** The order of an elliptic curve y**2 = x**3 + a*x + b over
** the field Fp = Z/p can be written as
**               N = (p + 1) - c,
** where c is the trace of the Frobenius endomorphism of the curve
** The function calculates (c mod L) for a prime L < p
*)
function ecp_traceL(p,a,b,L: integer): integer;
external
    Origin, Error: const;
    PsiMax: integer;
var
    Field: _Fieldpn;
    update: boolean;
    deg, deg0, deg1, k, y2, x, pL, pL1, c1, c2: integer;
    psi, F3, Y2, F: array;
    P0, P1, P2, qP, S, R, x0p, y0p, xP, x1P, zP: array;
begin
    if L < 2 then
        writeln("L must be a prime >= 2");
        return -1;
    elsif L > PsiMax then
        writeln("L too big: ",L);
        return -1;
    elsif prime32test(L) /= 1 then
        writeln("L not prime: ",L);
        return -1;
    end;

    if L = 2 then
        return ecp_ordpar(p,a,b);
    end;
    F3 := (b,a,0,1);
    psi := PsiN(p,a,b,L);
    deg := fpn_init(Field,p,psi);
    writeln(L,"-division polynomial has degree ",deg);
    P0 := ((0,1),{1});
    write("calculating x1 = x-coordinate of P1 = frob(P0) "); flush();
    x0p := fpn_power(Field,(0,1),p);
    writeln("... done");
    Y2 := fpn_polval(Field,F3,(0,1))
    pL := p mod L;
    if pL <= L div 2 then
        pL1 := pL;
    else
        pL1 := L - pL;
    end;

    write("comparing x1 with x-coordinate of k*P0 ");
    update := false;
    R := Origin;
    for k := 1 to L div 2 do
        R := ecpn_tadd(Field,Y2,a,R,P0);
        xP := R[0];
        write('.'); flush();
        zP := fpX_sub(p,x0p,xP);
        F := fpX_gcd(p,zP,psi);
        if F /= {1} then
            writeln(" done");
            Field.Q1 := F;
            deg := fpn_update(Field);
            writeln("found factor of division polynomial of degree ",deg);
        else
            if k = pL1 then
                if pL = pL1 then
                    qP := R;
                else
                    qP := ecpn_neg(Field,R);
                end;
            end;
            continue;
        end;
        if deg = 1 then
            F := fpn_fieldpol(Field);
            x := p - F[0];
            y2 := fpX_val(p,F3,x);
            if jacobi(y2,p) = 1 then
                return (p+1) mod L;
            else
                return (-p-1) mod L;
            end;
        end;

        P1 := ecpn_tfrob(Field,Y2,P0);
        R := ecpn_tmult(Field,Y2,a,P0,k);
        if P1 = R then
            c1 := k;
        elsif P1 = ecpn_neg(Field,R) then
            c1 := L - k;
        else
            update := true;
            break;
        end;
        c2 := mod_inverse(c1,L);
        c2 := c2*pL mod L;
        return (c1 + c2) mod L;
    end;
    if update = false then
        writeln(" done");
        write("calculating y-coordinate of P1 = frob(P0) "); flush();
        y0p := fpn_power(Field,Y2,(p-1) div 2);
        writeln("... done");
        P1 := (x0p, y0p);
    else
        qP := ecpn_tmult(Field,Y2,a,P0,pL);
    end;
    write("calculating P2 = frob(P1) "); flush();
    P2 := ecpn_tfrob(Field,Y2,P1);
    writeln("... done");
    pL := p mod L;
    writeln("p mod ",L," = ",pL);
    write("calculating P2 + ",pL,"*P0 "); flush();
    S := ecpn_tadd(Field,Y2,a,P2,qP);
    if S = Error then
        writeln("this case should not happen");
        return -1;            
    end;
    writeln("... done");
    writeln("comparing P2 + ",pL,"*P0 and k*P1 for k = 0,1,2,...");
    R := Origin;
    if S = R then
        return 0;
    end;
    for k := 1 to (L div 2) do
        R := ecpn_tadd(Field,Y2,a,R,P1);
        write('.'); flush();
        if R = Error then
            writeln("this case should not happen");
            return -1;            
        end;
        if R[0] = S[0] then
            writeln();
            if S[1] = R[1] then
                return k;
            else
                return L - k;
            end;
        end;
    end;
    return -1;
end;
(*------------------------------------------------------------*)
(*
** returns 0, if the order of the elliptic curve is even
** and 1, if the order is odd.
** The order is even iff the polynomial 
**     F(X) = X**3 + a*X + b
** has a zero in Fp. 
** This is tested by calculating the gcd of F(X) and X**p - X.
*)
function ecp_ordpar(p,a,b: integer): integer;
var
    F,G,X,Xp: array;
begin
    F := (b,a,0,1);
    X := (0,1);
    Xp := fpX_modpow(p,X,p,F);
    G := fpX_sub(p,Xp,X);
    G := fpX_gcd(p,G,F);
    if G = {1} then
        return 1;
    else
        return 0;
    end;
end;
(*------------------------------------------------------------*)
(***************************************************************)
(*
** The argument cpvec must be a vector of pairs
**    ((c1,p1),(c2,p2),...,(ck,pk)),
** where p1, p2, ... , pk are relatively prime.
** The function returns a pair (c,N), 
** where N := p1*p2*...*pk and c mod pi = ci
** If the pi are not relatively prime, () is returned.
*)
function chin_rem(cpvec: array of array[2]): array;
var
    k,len,u,v,a,c,d,p,N: integer;
begin
    len := length(cpvec);
    if len = 0 then
        return (0,1);
    elsif len = 1 then
        return cpvec[0];
    end;
    (a,N) := cpvec[0];
    for k := 1 to len-1 do
        (c,p) := cpvec[k];
        d := gcdx(N,p,u,v);
        if d /= 1 then
            return ();
        end;
        a := (a*v*p + c*u*N);
        N := p*N;
        a := a mod N;
    end;
    return (a,N);
end.
(*--------------------------------------------------------*)
(*
** Calculates the order of an elliptic curve by a combination 
** of Schoof's method and Pollard's probabilistic version 
** of Shanks' giant step baby step method
*)
function ecp_ordschoof(p,a,b: integer; Lmax := 0): integer;
var
    Lvec := ( 2, 3, 5, 7,11,13, 17, 19, 23, 29, 31, 37, 41, 43);
    bvec := (32,40,50,64,72,84,100,112,128,140,150,163,172,180);
    L, cL, c, k, N, n, n1, len, hasse: integer;
    tr: array;
    st: stack;
begin
    if p < 1000 then
        return ecp_ord0(p,a,b);
    end;
    len := length(Lvec);
    if Lmax <= 0 then
        n := bit_length(p);
        k := len-1;
        while k > 0 and n <= bvec[k-1] do
            dec(k);
        end;
    else
        k := 0;
        while k < len-1 and Lvec[k] < Lmax do
            inc(k);
        end;
    end;
    Lvec := Lvec[0..k];
    Lmax := Lvec[k];
    writeln("Schoof method: using primes L up to ",Lmax);
    for k := 0 to length(Lvec)-1 do
        L := Lvec[k];
        writeln("(************** trace mod ",L," **************)");
        cL := ecp_traceL(p,a,b,L);
        if cL >= 0 then
            stack_push(st,(cL,L));
            writeln("trace mod ",L," is ",cL);
        end;
    end;
    (c,n1) := tr := chin_rem(stack2array(st));
    writeln("(*************************************************)");
    writeln("trace modulo ",n1," is ",c);
    hasse := 2*isqrt(4*p);
    if n1 > hasse then
        if 2*c > hasse then
            c := c - n1;
        end;
        N := (p+1) - c;
    else
        writeln("now entering Pollard's lambda algorithm");
        N := ecp_ordlamTr(p,a,b,tr);
    end;
    return N;
end;
(***************************************************************)