(*-------------------------------------------------------*)
(*
** Diskreter Logarithmus
** siehe Funktion dlog(p,a,x) weiter unten
*)
(*-------------------------------------------------------*)
(*
** Sei X = (x1, x2, ..., xn) und d = gcd(x1,...,xn)
** Die Funktion bezout berechnet ein Array (v1,...vn),
** so dass v1*x1 + v2*x2 + ... + vn*xn = d
*)
function bezout(X: array): array;
var
    lambda: array;
    d, i, j, u, v, len: integer;
begin
    len := length(X);
    if len = 0 then
        return ();
    end;
    lambda := alloc(array,len);
    lambda[0] := 1;
    d := X[0];
    for i := 1 to len-1 do
        d := gcdx(d,X[i],u,v);
        for j := 0 to i-1 do
            lambda[j] := lambda[j]*u;
        end;
        lambda[i] := v;
    end;
    return lambda;
end;
(*-----------------------------------------------------*)
(*
** The argument zz must be an integer >= 2
** trialprime(zz) returns the least integer x >= zz,
** which has no prime divisor p < min(2**16,x)
*)
function trialprim(zz)
begin
        if even(zz) then
                inc(zz);
        end;
        while factor16(zz) do
                inc(zz,2);
        end;
        return zz;
end.
(*----------------------------------------------------*)
(*
** Returns the smallest (probable) prime >= zz
*)
function nextprime(zz)
begin
        zz := trialprim(zz);
        if zz <= 2**32 + 2**17 then
                return zz;
        end;
        while not rab_primetest(zz) do
                zz := trialprim(zz+2);
                write('.'); flush();
        end;
        writeln(" probable prime:");
        return zz;
end.
(*--------------------------------------------------------------*)
(*
** Erstellt ein Array mit den Primpotenz-Faktoren von x.
** Falls letztes Element des Arrays > 2**32, ist dies
** nicht notwendig eine Primzahlpotenz
*)
function primepowlist(x: integer): array;
var
    st: stack;
    q, y, t: integer;
begin
    y := abs(x);
    q := 2;
    while q := factor16(y,q) do
        t := q;
        y := y div q;
        while y mod q = 0 do
            y := y div q;
            t := t*q;
        end;
        stack_push(st,t);
    end;
    if y > 1 then
        stack_push(st,y);
    end;
    return stack2array(st);
end;
(*---------------------------------------------------------------*)
(*
** Erstellt ein Array aus den Primfaktoren von x, wobei jeder
** Primfaktor nur einmal aufgezaehlt wird. 
** Falls das letzte Element des Arrays > 2**32 ist,
** ist es nicht notwendig prim.
*)
function primefactors(x: integer): array;
var
    st: stack;
    q, y: integer;
begin
    y := abs(x);
    q := 2;
    while q := factor16(y,q) do
        stack_push(st,q);
        y := y div q;
        while y mod q = 0 do
            y := y div q;
        end;
    end;
    if y > 1 then
        stack_push(st,y);
    end;
    return stack2array(st);
end;
(*----------------------------------------------------*)
(*
** Die Funktion berechnet eine Primitivwurzel mod p
** durch Suche von a ab.
*)
function primroot(p: integer; a := 2): integer;
var
    i, k, b: integer;
    found: boolean;
    qvec: array;
begin
    qvec := primefactors(p-1);
    k := -1;
    while inc(k) < p do
        b := (a+k) mod (p-1);
        found := true;
        for i := 0 to length(qvec)-1 do
            if b**((p-1) div qvec[i]) mod p = 1 then
                found := false;
                break;
            end;
        end;
        if found then
            return b;
        end;
    end;
    return 0;   (* this case should not happen *)
end;
(*-----------------------------------------------------------------*)
function logtab(p,g: integer): array;
var
    vec: array[p];
    x,k: integer;
begin
    vec[0] := -1;
    x := 1;
    for k := 0 to p-2 do
        vec[x] := k;
        x := x*g mod p;
    end;
    return vec;
end;
(*-----------------------------------------------------------------*)
function balken(n)
var
    str: string;
begin
    return(alloc(string,n,'*'));
end;
(*-----------------------------------------------------------------*)
function logplot(p)
var
    g,k,y: integer;
    vec: array;
begin
    g := primroot(p);
    vec := logtab(p,g);
    writeln("Diskreter Logarithmus mod ",p," zur Basis ",g);
    for k := 1 to p-1 do
        y := vec[k];
        writeln(k:4,y:4,"  ",balken(y));
    end;
end;
(*-----------------------------------------------------------------*)
(*
** Hilfsfunktion fuer dlog_h0
*)
function preptab(var Gtab: array; p,b,anz,collinc: integer): integer;
var
    y, i, k, hlen, hidx, hidx1: integer;
    collcount := 0;
begin    
    hlen := length(Gtab);
    Gtab[1] := (0,1);
    y := 1
    for k := 1 to anz-1 do
        y := y*b mod p;
        hidx := y mod hlen;
        if Gtab[hidx] = (0,0) then
            Gtab[hidx] := (k,y);
        else
            for i := 1 to hlen do
                inc(collcount);
                hidx1 := (hidx + i*collinc) mod hlen;
                if Gtab[hidx1] = (0,0) then
                    Gtab[hidx1] := (k,y);
                    break;
                end;
            end;
        end;
    end;
    return collcount;
end;
(*----------------------------------------------------*)
(*
** p ist eine ungerade Primzahl, N ein Teiler von p-1.
** a ist ein Element von (Z/p)* von der Ordnung N
** und x ein Element in der von a erzeugten Untergruppe.
** Die Funktion gibt den Logarithmus von x bzgl. a zurueck.
** Der Algorithmus benutzt die Giant-Step-Baby-Step-Methode
** und Hashing
*)
function dlog_h0(p,a,x,N: integer): integer;
external
    Gtab: array of array[2];
var
    gstep, hlen, collinc, b, y, i, j, k, hidx, hidx1, coll: integer;
begin
    gstep := isqrt(N) + 1;
    hlen := nextprime(2*gstep);
    collinc := nextprime((hlen div 10) + 1);
    Gtab := alloc(array,hlen,(0,0));
    b := a**gstep mod p;
    coll := preptab(Gtab,p,b,gstep,collinc);
    b := mod_inverse(a,p);
    y := (x mod p);
    for i := 0 to gstep-1 do
        hidx := y mod hlen;
        for j := 0 to hlen-1 do
            hidx1 := (hidx + j*collinc) mod hlen;
            if Gtab[hidx1] = (0,0) then
                break;
            elsif Gtab[hidx1][1] = y then
                k := Gtab[hidx1][0];
                return (k*gstep + i);
            end;
        end;
        y := y*b mod p;
    end;
    writeln("dlog_h0: this case should not happen");
    halt(-1);
end;
(*----------------------------------------------------------*)
(*
** Zufallsgenerator
**      x --> (i,j) in (Z/r) x (Z/r)
*)
function randij(x, r: integer; var i,j: integer);
var
    y,z: integer;
begin
    y := x mod 8093;
    z := x mod 12347;
    i := (y*y + 2) mod r;
    j := z mod r;
end;
(*----------------------------------------------------------------*)
(*
** Monte-Carlo-Algorithmus zur Berechnung des diskreten Logarithmus.
** q muss ein Teiler von p-1 sein, moeglichst prim,
** a ein Element von (Z/p)* der Ordnung q
** und x ein Element der von a erzeugten Untergruppe
*)
function dlog_rq(p,a,x,q: integer): integer;
const
    rlen = 13;
var
    Aarr, Xarr: array[rlen];
    aexp, xexp: array[rlen];
    i,j,k,I,J,y,z,nu: integer;
    found: boolean;
begin

    for i := 0 to rlen-1 do
        k := random(q);
        aexp[i] := k;
        Aarr[i] := a**k mod p;
        k := random(q);
        xexp[i] := k;
        Xarr[i] := x**k mod p;
    end;

    found := false;
    I := J := 0;
    z := y := 1 + random(q-1);
    for nu := 1 to 10*isqrt(q) do

        randij(y,rlen,i,j);
        dec(I,xexp[i]); inc(J,aexp[j]);
        y := y * Xarr[i] * Aarr[j] mod p;

        randij(z,rlen,i,j);
        inc(I,xexp[i]); dec(J,aexp[j]);
        z := z * Xarr[i] * Aarr[j] mod p;
        randij(z,rlen,i,j);
        inc(I,xexp[i]); dec(J,aexp[j]);
        z := z * Xarr[i] * Aarr[j] mod p;

        if nu mod 512 = 1 then
            write('.'); flush();
        end;
        if z = y then
            writeln(" (",nu," steps)");
            found := true;
            break;
        end;
    end;
    if found then
        I := I mod q; J := J mod q;
        if gcd(I,q) /= 1 then
            return -1;
        else
            return (mod_inverse(I,q)*J) mod q;
        end;
    else
        return -2;
    end;
end;
(*----------------------------------------------------*)
(*
** Berechnet den diskreten Logarithmus mod p
** von x zur Basis a.
** x muss in der von a erzeugten Untergruppe von (Z/p)* liegen
*)
function dlog(p,a,x): integer;
var
    pvec, Nvec, lvec, lambda: array;
    i, len, k, A, X, q, n: integer;
    verb: boolean;
begin
    if not rab_primetest(p) then
        writeln("prime number expected: ",p);
        halt(-1);
    end;
    pvec := primepowlist(p-1);
    len := length(pvec);
    if (q := pvec[len-1]) > 2**32 then
        writeln("warning: calculation will take a long time");
    end;
    Nvec := alloc(array,len);
    for i := 0 to len-1 do
        Nvec[i] := (p-1) div pvec[i];
    end;
    lambda := bezout(Nvec);
    lvec := alloc(array,len);
    write("working ");
    for i := len-1 to 0 by -1 do
        q := pvec[i];
        write(q," "); flush();
        A := a**Nvec[i] mod p;
        X := x**Nvec[i] mod p;
        if q > 2**16 then
            n := dlog_rq(p,A,X,q);
            if n < 0 then
                writeln("random walk failed");
                halt(-1);
            end;
            lvec[i] := n;
        else
            lvec[i] := dlog_h0(p,A,X,q);
        end;
    end;
    k := 0;        
    for i := 0 to len-1 do
        k := k + lvec[i]*Nvec[i]*lambda[i];
    end;
    writeln();
    return k mod (p-1);
end;
(*************************************************************)