(************************************************************)
(*
** Vorlesung Algorithmische Zahlentheorie und Kryptographie 
** von Otto Forster im WS 2015/16 am Math. Inst. der LMU Muenchen
**
** Aribas Code zur Faktorisierung mit Elliptischen Kurven
**
** File author: Otto Forster,  Email:  forster AT math.lmu.de
** Date of last change: 2016-02-05
*)
(**********************************************************************)
(*
** Beispiel
**

==> f7 := 2**128 + 1.
-: 3402_82366_92093_84634_63374_60743_17682_11457

factorize(f7,16000,400).
working ................................................................
factor found after 64 curves with prime bound 10496
-: 59_64958_91274_97217

==> f8 := 2**256 + 1.
-: 115_79208_92373_16195_42357_09850_08687_90785_32699_84665_64056_40394_
57584_00791_31296_39937

==> ECfactorize(f8,16000,400).
working .......................
factor found after 23 curves with prime bound 8576
-: 1_23892_63615_52897

*)
(*****************************************************************)
(*
** Faktorisierungs-Algorithmus mit elliptischen Kurven
** N ist die zu faktorisierende Zahl, bound eine Schranke
** fuer die Primfaktoren des Multiplikators und
** anz  die Anzahl der Versuche.
** Rueckgabewert: Ein Teiler d von N  oder 0, falls erfolglos 
**
** Bemerkung: Die Faktorisierung mit elliptischen Kurven ist
** auch als eingebaute ARIBAS-Funktion  ec_factorize()  vorhanden
*)
function ECfactorize(N, bound, anz: integer): integer;
const
    anz0 = 128;
var
    k, a, x, y, d, B0, B1, s: integer;
begin
    write("working ");
    for k := 1 to anz do
        write('.');
        a := random(1000);
        x := random(N); y := 1;
        for B0 := 0 to bound-1 by anz0 do
            B1 := min(B0+anz0,bound);
            s := ppexpo(B0,B1);
            (x,y) := ecN_mult(N,a,(x,y),s);
            if y < 0 and x > 1 then
                writeln();  
                write("factor found after ",k," curves");
                writeln(" with prime bound ",B1);
                return x;
            elsif y < 0 then
                break;
            end;
        end;
    end;
    return 0;
end;
(*--------------------------------------------------------------*)
(*
** Hilfsfunktion fuer ECfactorize
**
** Multiplikation eines Punktes P mit s > 0 auf einer elliptischen
** Kurve mit affiner Gleichung  Y**2 = X**3 + a*X + b
** ueber dem Ring  Z/N
** Falls dabei ein Teiler d von N entdeckt wird, ist der
** Rueckgabewert  (d, -1); sonst s*P = (x,y) mit 0 <= x,y < N.
*)
function ecN_mult(N,a:integer; P:array[2]; s:integer): array[2];
var
    x,x0,x1,y,y0,m1,mu,Fprime,d,k: integer;
begin
    (x0,y0) := (x,y) := P;
    for k := bit_length(s)-2 to 0 by -1 do
        m1 := mod_inverse(2*y,N);
        if m1 = 0 then return (gcd(y,N),-1); end;
        Fprime := (3*x**2 + a) mod N;
        mu := (Fprime*m1) mod N;
        x1 := x;
        x := (mu**2 - 2*x) mod N;
        y := (-y - mu*(x - x1)) mod N;
        if bit_test(s,k) then
            m1 := mod_inverse(x-x0,N);
            if m1 = 0 then return (gcd(x-x0,N),-1); end;
            mu := m1*(y - y0) mod N;
            x := (mu**2 - x - x0) mod N;
            y := (-y0 - mu*(x - x0)) mod N;
        end;
    end;
    return (x,y);
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 ECfactorize
*)
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;
(*****************************************************************)