(****************************************************************)
(*
** ARIBAS code for
** several factoring routines for integers
** author: Otto Forster <forster@rz.mathematik.uni-muenchen.de>
** date: 99-04-30
**
**
** The factoring algorithms
**     p1_factorize, pp1_factorize, ec_factorize, ecfactor 
** and the Aribas builtin functions
**     rho_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 factors (below).
**
** Example calls:
**
** ==> factors(10**15+1).
** ==> primelist(10**9,10**9 + 1000).
** ==> p1_factorize(2**67-1).
** ==> pp1_factorize(2**67-1).
** ==> ecfactor(2**79-1).
*)
(*--------------------------------------------------------------*)
(*
** Trial division by primes p < 2**16
** Writes the found prime divisiors of x to the terminal and
** returns the last cofactor. 
** If this cofactor is less than 2**32, it is a prime
*)
function factors(x: integer)
var
        q0, q;
begin
        q0 := 2;
        while q := factor16(x,q0) do
                writeln(q);
                x := x div q;
                q0 := q;
        end;
        return(x);
end;
(*--------------------------------------------------------*)
(*
** Erstellt ein Array von Faktoren von x
** Alle Elemente, bis auf moeglicherweise das letzte, sind
** Primfaktoren < 2**16. Ist das letzte Element < 2**32,
** so ist es ebenfalls prim
*)
function factorlist(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;
(*--------------------------------------------------------------*)
(*
** Returns a list of all primes p with
**      first <= p <= last
** The argument last must satisfy last <= 2**32 and 
** last - first <= 2**16.
** If only the argument first is given,
** last := first + 1000 by default
**
** Example calls: 
**	primelist(2,100).
**	primelist(10**6).
*)
function primelist(first: integer; last := -1): array;
var
        st: stack;
        n: integer;
        vec: array;
begin
	if last < 0 then
		last := first + 1000;
	elsif last > first + 0x10000 then
		last := first + 1000;
	end;
        if last > 2**32 + 2**17 then
                writeln("upper bound <= 2**32 expected");
                return;
        end;
        if first <= 2 and last > 2 then
                stack_push(st,2);
                n := 3;
        else
                n := trialprim(first);
        end;
        while n <= last do
                stack_push(st,n);
                n := trialprim(n+2);
        end;
        return stack2array(st);
end.
(*--------------------------------------------------------------*)
(*
** Solovay-Strassen primality test
** The argument x must be > 2**16
** 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(64000);
        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 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;
(*--------------------------------------------------------*)
(*
** 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;
(*-----------------------------------------------------------------*)
(*
** Faktorisierungs-Algorithmus mit elliptischen Kurven.
** bound ist Schranke fuer Primfaktoren der Elementezahl
** der (zufaellig gewaehlten) elliptischen Kurve
** anz ist die maximale Anzahl der Versuche
*)
function ec_factorize(N: integer; bound := 500; anz := 100): integer;
var
    k, a, d: integer;
begin
    write("working ");
    for k := 1 to anz do
        a := random(64000);
        d := gcd(a*a-4,N);
        if d = 1 then
            write('.'); flush();
            d := ec_fact0(N,a,bound);
        end;
        if d > 1 and d < N then return d; end;
    end;
    return 0;
end;
(*-----------------------------------------------------------------*)
(*
** 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
** Diese Funktion wird aufgerufen von der Funktion ec_factorize
*)
(*-----------------------------------------------------------------*)
function ec_fact0(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;
(*--------------------------------------------------------------*)
(*
** Big prime variation of ec_factorize
** arguments N, bound and anz as in ec_factorize
** bound2 is a bound for the big prime
*)
function ec_factbpv(N: integer; bound := 500; 
		bound2 := 10000; anz := 100): 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 := ec_fact0(N,a,bound);
        end;
        if d <= 0 then
            write(':'); flush();
            d := ec_bigprimevar(N,a,-d,bound2);
        end;
        if d > 1 and d < N then return d; end;
    end;
    return 0;
end;
(*-------------------------------------------------------------*)
(*
** auxiliary function, called by ec_factbpv
*)
function ec_bigprimevar(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;
(*--------------------------------------------------------------*)
(*
** 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;
(*------------------------------------------------------------------*)
(*
** Tries to give a full prime factorization of x,
** using the elliptic curve method.
** If ECM fails, quadratic sieve method is used.
*)
function ecfactor(x: integer): array;
var
    st, st1: stack;
    q,y,z,count,level: integer;
    yz: array[2];
    vec: array;
    first: boolean;
begin
    writeln("Factorizing");
    writeln(x);
    if rab_primetest(x) then
	writeln("this number is prime");
	return {x};
    end;
    q := 2;
    writeln("trial division by factors < 2**16:");
    count := 0;
    while q := factor16(x,q) do
	writeln(q);
	stack_push(st,q);
        x := x div q;
	inc(count);	
    end;
    if count = 0 then
	writeln("no factors < 2**16");
    else
	writeln("remaining cofactor:");
	writeln(x);
    end;
    if rab_primetest(x) then
	writeln("this factor is prime");
	stack_push(st,x);
	return stack2array(st);	
    else
	if count > 0 then
	    writeln("this factor is composite");
	end;
	stack_push(st1,x);
    end;
    first := true;
    while not stack_empty(st1) do
        x := stack_pop(st1);
        if first then
	    level := 1;
	    first := false;
        elsif bit_length(x) < 80 then
	    level := 1;
	else
	    level := 2;
        end;
        yz := factor0(x,level);
        y := yz[0]; z := yz[1];
        writeln(y);
        if rab_primetest(y) then
	    writeln("this factor is prime");
	    stack_push(st,y);
        else
	    writeln("this factor is composite");
	    if z > 1 then
	        stack_push(st1,y);
	    else
		writeln("complete factorization failed");
	        stack_push(st,y);
	    end;
        end;
	if z <= 1 then
	    stack_push(st,z);
	else
            writeln(z);
            if rab_primetest(z) then
		writeln("cofactor is prime");
		stack_push(st,z);
    	    else
		writeln("cofactor is composite");
		stack_push(st1,z);
    	    end;
	end;
    end;
    vec := stack2array(st);
    return sort(vec);
end;
(*-----------------------------------------------------*)
(*
** auxiliary function, called by ecfactor
*)
function factor0(x: integer; level := 1): array;
var
    bound, anz, q: integer;
    first: boolean;
begin
    if rab_primetest(x) then
	return (x,-1);
    end;
    first := true;
    while level <= 2 do
        if level = 1 then
	    bound := 2000;
	    anz := 50;
        else
	    bound := 5000;
	    anz := 200;
	end;
	if first then
	    writeln("trying to factorize");
	    writeln(x);
	    writeln("by elliptic curve method, bound = ",bound);
	    first := false;
	else
	    writeln();
	    writeln("trying once more with bound = ",bound);
	end;
	q := ec_factorize(x,bound,anz);
	if q /= 0 then
	    return (q, x div q);
	elsif level >= 2 then
	    writeln();
	    writeln("now using quadratic sieve algorithm");
	    q := qs_factorize(x);
	    if q = 0 then
		writeln("factorization failed");
		return (x,1);
	    else
		return (q, x div q);
	    end;
	end;
	inc(level);
    end;
    return (x,1);
end;
(*******************************************************************)