(******************************************************************)
(*
** Vorlesung Algorithmische Zahlentheorie und Kryptographie 
** von Otto Forster im WS 2015/16 am Math. Inst. der LMU Muenchen
**
** ARIBAS Beispiel-Code zu Paragraph 10
** Die (p-1)- und die Rho-Faktorisierungs-Methode
** 
** File author: Otto Forster,  Email:  forster AT math.lmu.de
** Date of last change: 2015-12-17
*)
(*******************************************************************)
(* Beispiel

==> f6 := 2**64 + 1.
-: 18446_74407_37095_51617

==> p1_factorize(f6).
working
factor found with bound 128
-: 274177

==> poll_rho(f6).
working .
factor found after
512 iterations
-: 274177

*)
(*******************************************************************)
(*
** Stellt alle Primfaktoren von x, die kleiner als 2**16 sind,
** fest und gibt sie aus.
** Funktionswert: letzter Primfaktor oder letzter Cofaktor
** Ist dieser < 2**32, so ist er prim
*)
function trialdiv(x: integer): integer;
var
    q: integer;
begin
    q := 2;
    while q := factor16(x,q) do
        writeln(q);
        x := x div q;
    end;
    return x;
end;
(*----------------------------------------------------------*)
(*
** Produkt aller Primzahlen B0 < p <= B1
** und aller ganzen Zahlen isqrt(B0) < n <= isqrt(B1)
** Diese Funktion wird gebraucht von den Funktionen
** p1_factorize, pp1_factorize und ec_factorize
*)
function ppexpo(B0,B1: integer): integer;
var
    x, m0, m1, i: integer;
begin
    x := 1;
    m0 := max(2,isqrt(B0)+1); m1 := isqrt(B1);
    for i := m0 to m1 do
        x := x*i;
    end;
    if odd(B0) then inc(B0) end;
    for i := B0+1 to B1 by 2 do
        if prime32test(i) > 0 then x := x*i end;
    end;
    return x;
end;
(*-----------------------------------------------------*)
(*
** (p-1)-Faktorisierungs-Algorithmus von Pollard
** Findet i.a. Primfaktoren p von x, wenn p-1 Produkt
** kleiner Primzahlen q <= bound ist.
** Defaultwert fuer bound = 32000
*)
function p1_factorize(N: integer; bound := 32000): integer;
const
    anz0 = 128;
var
    base, d, B0, B1, ex, count: integer;
begin
    base := 2 + random(N-2);
    if (d := gcd(base,N)) > 1 then return d end;
    write("working ");
    count := 0;
    for B0 := 0 to bound-1 by anz0 do
        B1 := min(B0+anz0, bound);
        ex := ppexpo(B0,B1);
        base := base ** ex mod N;
        if base = 1 then return 0 end;
        d := gcd(base-1,N);
        if d > 1 then
            writeln();
            writeln("factor found with bound ",B1);
            return d;
        end;
        inc(count);
        if bit_and(count,7) = 0 then write('.'); end;
    end;
    return 0;
end;
(*--------------------------------------------------------*)
(*
** Versucht, die Zahl N mittels der Pollardschen Rho-Methode
** zu faktorisieren. Der Parameter anz ist eine Schranke fuer 
** die Anzahl der Iterationen.
** Ein Faktor p von N wird im allgemeinen gefunden, falls
** anz groesser als die Quadratwurzel von p ist.
**
** Bemerkung: In ARIBAS als eingebaute Funktion mit dem Namen
** rho_factorize vorhanden
*)
function poll_rho(N: integer; anz := 32000): integer;
const
    anz0 = 256;
var
    x, y, i, d, P: integer;
begin
    y := x := random(N);
    write("working ");
    for i := 0 to (anz-1) div anz0 do
        P := accumdiff(x,y,N,anz0);
        d := gcd(P,N);
        if d > 1 and d < N then
            writeln(); writeln("factor found after ");
            writeln((i+1)*anz0," iterations");
            return(d);
        end;
        if bit_and(i,7) = 0 then write('.'); end;
    end;
    return 0;
end;
(*----------------------------------------------------------------*)
function accumdiff(var x,y: integer; N, anz: integer): integer;
var
    i, P: integer;
begin
    P := 1;
    for i := 1 to anz do
        x := (x*x + 2) mod N;
        y := (y*y + 2) mod N;
        y := (y*y + 2) mod N;
        P := P * (y-x) mod N;
    end;
    return P;
end;
(********************************************************************)