(*---------------------------------------------------------------*)
(*
** Stellt fest, ob x eine Primitivwurzel mod p ist.
** Es wird vorausgesetzt, dass p prim ist,
** qvec ein Array, dessen Komponenten die Primfaktoren von p-1 sind,
** und x eine nicht durch p teilbare ganze Zahl.
*)
function check_primroot(x,p: integer; qvec: array): boolean;
var
    i,m: integer;
begin
    for i := 0 to length(qvec)-1 do
        m := (p-1) div qvec[i];
        if x**m mod p = 1 then
            return false;
        end;
    end;
    return true;
end;
(*---------------------------------------------------------------*)
(*
** Liefert fuer eine ganze Zahl x einen Vektor,
** dessen Komponenten die Primfaktoren von x sind.
** Mehrfache Faktoren werden nur einmal aufgezaehlt.
** Funktioniert nur, wenn x hoechstens einen Primfaktor > 2**16
** hat, der aber < 2**32 sein muss.
*)
function factor_vec(x: integer): array;
var
    q, x0: integer;
    st: stack;
begin
    q := 2; x0 := x;
    while q := factor16(x,q) do
        stack_push(st,q);
        x := x div q;
        while x mod q = 0 do x := x div q; end;
    end;
    if x > 2**32 then
        writeln("WARNING: factors of ",x0," too big");
    end;
    if x > 1 then stack_push(st,x); end;
    return stack2array(st);
end;
(*---------------------------------------------------------------*)
(*
** Liefert fuer eine ungerade Primzahl p
** die kleinste Primitivwurzel mod p, die >= g ist.
** Setzt voraus, dass p-1 von der Funktion factor_vec
** vollstaendig faktorisiert werden kann.
*)
function primroot(p: integer; g := 2): integer;
var
    k, x: integer;
    qvec: array;
begin
    qvec := factor_vec(p-1);
    for k := 0 to p-2 do
        x := (g+k) mod p;
        if x /= isqrt(x)**2 and check_primroot(x,p,qvec) then
            return x;
        end;
    end;
end;
(*---------------------------------------------------------------*)