(****************************************************************)
(*
** ARIBAS code for
** several factoring routines for integers
** author: Otto Forster <forster@mathematik.uni-muenchen.de>
** date of last change: 2007-08-23
**
** This code is place under the GNU general public license (GPL)
*)
(****************************************************************)
(*
** The factoring algorithms
**     p1_factorize, pp1_factorize, ECfactorize 
** and the Aribas builtin functions
**     rho_factorize, ec_factorize, cf_factorize and qs_factorize
** should be applied only to numbers which are not prime.
** This can be tested by rab_primetest or ss_test
** Also, before applying one of these factoring algorithms,
** one should first do trial division by small primes,
** for example using the function trialdiv (below).
**
** Example calls:
**
** ==> trialdiv(10**15+1).
** ==> factorlist(2**171-1)
** ==> p1_factorize(2**67-1).
** ==> pp1_factorize(2**67-1).
** ==> ECfactorize(2**256+1, 4000, 32000).
*)
(*--------------------------------------------------------------*)
(*
** Trial division by primes p < 2**16
** Constructs an array of factors of x.
** All elements with possible exeption of the last,
** are prime factors < 2**16. If the last element
** in the array is < 2**32, it is a prime.
** The product of all elements in the array equals x.
*)
function trialdiv(x: integer): array;
var
    st: stack;
    q: integer;
begin
    q := 2;
    while q := factor16(x,q) do
        stack_push(st,q);
        x := x div q;
    end;
    stack_push(st,x);
    return stack2array(st);
end;
(*--------------------------------------------------------------*)
(*
** Constructs a list of all prime factors of x.
** If a prime p is a multiple factor of x,
** it is listed repeatedly according to its multiplicity.
** In case of failure, the empty list () is returned 
**
** Uses trial division, rho_factorize and qs_factorize.
**
** This function writes a progress report to the screen,
** which can be suppressed by setting the last argument verbose = 0.
*)
function factorlist(x: integer; verbose := 1): 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,verbose);
            count := 0;
            while (y <= 1 or y >= x) do
                if inc(count) > 2 then
                    writeln("unable to factorize ",x);
                    return ();
                end;
                if verbose then
                    writeln("trying to factorize ",x,
                            " using quadratic sieve");
                end;
                y := qs_factorize(x,verbose);
            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;
(*--------------------------------------------------------------*)
(*
** Solovay-Strassen primality test
** The argument x must be a positive integer.
** This is a probabilistic test:
** If ss_test(x) returns false, then x is certainly not a prime.
** If however the result is true, then x is only probably prime.
** To increase the probabilty, one can repeat the test.
*)
function ss_test(x: integer): boolean;
var
    b, j, u: integer;
begin
    if even(x) then return false end;
    b := 2 + random(x-2);
    j := jacobi(b,x);
    u := b ** (x div 2) mod x;
    if j = 1 and u = 1 then
        return true;
    elsif (j = -1) and (u = x-1) then
        return true;
	else
        return false;
    end;
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;
(*--------------------------------------------------------*)
(*
** Pollard's (p-1)-factoring algorithm
** In general a prime factor p of x is found, if
** p-1 is a product of prime powers q**k <= bound
*)
function p1_factorize(x: integer; bound := 16000): integer;
const
    anz0 = 128;
var
    base, d, n, n0, n1, ex: integer;
begin
    base := 2 + random(64000);
    d := gcd(base,x);
    if d > 1 then
        return d;
    end;
    writeln(); write("working ");
    for n0 := 0 to bound-1 by anz0 do
        n1 := min(n0 + anz0, bound);
        ex := ppexpo(n0,n1);
        base := base ** ex mod x;
        write('.'); flush();
        if base <= 1 then
            return 0;
        else
            d := gcd(base-1,x);
            if d > 1 then
                writeln();
                writeln("factor found with bound ",n1-1)
                return d;
            end;
        end;
    end;
    return 0;
end;
(*-----------------------------------------------------*)
(*
** (p+1)-factoring algorithm
*)
function pp1_factorize(x: integer; bound := 16000): integer;
const
    anz0 = 128;
var
    base, d, n, n0, n1, ex: integer;
begin
    base := 2 + random(64000);
    d := gcd(base,x);
    if d > 1 then
        return d;
    end;
    writeln();
    write("working ");
    for n0 := 0 to bound-1 by anz0 do
        n1 := min(n0 + anz0, bound);
        ex := ppexpo(n0,n1);
        base := mod_coshmult(base,ex,x);
        write('.'); flush();
        if base <= 1 then
            return 0;
        else
            d := gcd(base-1,x);
            if d > 1 then
                writeln();
                writeln("factor found with bound ",n1-1)
                return d;
            end;
        end;
    end;
    return 0;
end;
(*-----------------------------------------------------------------*)
(*
** Elliptic curve factorization with big prime variation.
** N is the number to be factored.
** bound and bound2 are bounds for the prime factors
** of the order of the randomly chosen elliptic curve. 
** anz is the maximal number of elliptic curves tried
** Returns a factor of N or 0 in the case of failure
*)
function ECfactorize(N: integer; bound := 1000; 
                     bound2 := 10000; anz := 200): integer;
var
    k, a, d: integer;
begin
    write("working ");
    for k := 1 to anz do
        a := 3 + random(64000);
        d := gcd(a*a-4,N);
        if d = 1 then
            write('.'); flush();
            d := ECfact0(N,a,bound);
        end;
        if d <= 0 then
            write(':'); flush();
            d := ECbigprimevar(N,a,-d,bound2);
        end;
        if d > 1 and d < N then return d; end;
    end;
    return 0;
end;
(*-----------------------------------------------------------------*)
(*
** Called by function ECfactorize, not to be called directly
**
** Faktorisierungs-Algorithmus mit der elliptischen Kurve
**      y*y = x*x*x + a*x*x + x
** bound ist Schranke fuer die Primfaktoren der Elementezahl
**      der elliptischen Kurve
*)
(*-----------------------------------------------------------------*)
function ECfact0(N,a,bound: integer): integer;
const
    anz0 = 128;
var
    x, B0, B1, s, d: integer;
    xx: array[2];
begin
    x := random(N);
    for B0 := 0 to bound-1 by anz0 do
        B1 := min(B0+anz0,bound);
        s := ppexpo(B0,B1);
        xx := mod_pemult(x,s,a,N);
        if xx[1] = 0 then
            d := xx[0];
            if d > 1 and d < N then
                writeln(); write("factor found with curve ");
                writeln("parameter ",a," and bound ",B1);
            end;
            return d;
        else
            x := xx[0];
        end;
    end;
    return -x;
end;
(*--------------------------------------------------------------*)
(*
** auxiliary function, called by ECfactorize
*)
function ECbigprimevar(N,a,x,bound: integer): integer;
const
    Maxhdiff = (22, 36, 57);
    Maxbound = (15000, 31000, 1000000);
var
    XX: array of array[2];
    maxhdiff: integer;
    c, i, q, k, d: integer;
    P,Q,R: array[2];
begin
    k := length(Maxhdiff) - 1;
    while k > 0 and bound <= Maxhdiff[k-1] do
	    dec(k);
    end;
    bound := min(bound,Maxbound[k]);
    maxhdiff := Maxhdiff[k];
    XX := alloc(array,maxhdiff+1,(0,0));

    c := ((x + a)*x + 1)*x mod N;
    P := (x,1);
    Q := ecN_dup(N,a,c,P);
    if Q[1] < 0 then return Q[0]; end;
    XX[1] := R := Q;
    for i := 2 to maxhdiff do
	    R := ecN_add(N,a,c,R,Q);
        if R[1] < 0 then return R[0]; end;
	    XX[i] := R;
    end;
    R := ecN_add(N,a,c,P,Q);	(* R = 3*P *)
    if R[1] < 0 then return R[0]; end;
    d := 0;
    q := 3; 
    while q < bound do
        k := 1; inc(q,2);
        while prime32test(q) /= 1 do
	        inc(q,2); inc(k);
	    end;
        R := ecN_add(N,a,c,R,XX[k]);
	    if R[1] < 0 then
	        d := R[0]; 
	        if d > 1 and d < N then
	            writeln();
                writeln("factor found with curve parameter ",a,
	                    ", bigprime q = ",q);
            end;
            break;
	    end;
    end;
    return d;
end;
(*--------------------------------------------------------------*)
(*
** auxiliary function, called by ECfactbpv
**
** Addition zweier Punkte P,Q auf der elliptischen Kurve
**      c*y**2 = x**3 + a*x**2 + x (modulo N)
** Falls waehrend der Rechnung durch eine nicht zu N teilerfremde
** Zahl geteilt werden muss, wird ein Paar (d,-1) zurueckgegeben,
** wobei d ein Teiler von N ist.
** Sonst Rueckgabe der Summe P+Q = (x,y) mit 0 <= x,y < N.
*)
function ecN_add(N,a,c: integer; P,Q: array[2]): array[2];
var
    x1,x2,x,y1,y2,y,m: integer;
begin
    if P = Q then
        return ecN_dup(N,a,c,P);
    end;
    x1 := P[0]; x2 := Q[0];
    m := mod_inverse(x2-x1,N);
    if m = 0 then
        return (gcd(x2-x1,N),-1);
    end;
    y1 := P[1]; y2 := Q[1];
    m := (y2 - y1)*m mod N;
    x := (c*m*m - a - x1 - x2) mod N;
    y := (- y1 - m*(x - x1)) mod N;
    return (x,y);
end;
(*-------------------------------------------------------------*)
(*
** Verdopplung eines Punktes P auf der elliptischen Kurve
**      c*y**2 = x**3 + a*x**2 + x (modulo N)
** Falls waehrend der Rechnung durch eine nicht zu N teilerfremde
** Zahl geteilt werden muss, wird ein Paar (d,-1) zurueckgegeben,
** wobei d ein Teiler von N ist.
** Sonst Rueckgabe von P+P = (x,y) mit 0 <= x,y < N.
*)
function ecN_dup(N,a,c: integer; P: array[2]): array[2];
var
    x1,x,y1,y,z,m,Pprim: integer;
begin
    x1 := P[0]; y1 := P[1];
    z := 2*c*y1;
    m := mod_inverse(z,N);
    if m = 0 then
        return (gcd(z,N),-1);
    end;
    Pprim := (((3*x1 + 2*a)*x1) + 1) mod N;
    m := Pprim*m mod N;
    x := (c*m*m - a - 2*x1) mod N;
    y := (- y1 - m*(x - x1)) mod N;
    return (x,y);
end;
(*------------------------------------------------------------------*)
(*
** Multiplication of a point P on the elliptic curve
**      c*y**2 = x**3 + a*x**2 + x (modulo N)
** by an integer s >= 1.
** If during the calculation a division by a number which is
** not coprime to N must be performed, the function returns
** immediately a pair (d,-1), where d is a divisor of N.
*)
function ecN_mult(N,a,c: integer; P: array[2]; s: integer): array[2];
var
    k: integer;
    Q: array[2];
begin
    if s = 0 then return (0,-1); end;
    Q := P;
    for k := bit_length(s)-2 to 0 by -1 do
        Q := ecN_dup(N,a,c,Q);
        if Q[1] < 0 then
            return Q;
        end;
        if bit_test(s,k) then
            Q := ecN_add(N,a,c,Q,P);
            if Q[1] < 0 then
                return Q;
            end;
        end;
    end;
    return Q;
end;
(*********************************************************************)