/*
 *  An implementation of the Number Field Sieve as described in
 *  "Factoring Integers with the Number Field Sieve" by J.P. Buhler,
 *  H.W. Lenstra, and Carl Pomerance.  (Algorithm 11.1)
 *
 *  This code is not at all optimized - it is demonstration of the
 *  Magma computer algebra package, and to for people who want to
 *  learn the Number Field Sieve algorithm or how to code it.  Note
 *  that the algebraic square root code is available only for the
 *  most recent release of Magma.
 *
 *  This code is best demonstrated for numbers around 20-digits.
 *  Larger numbers will cause the program to take up too much memory,
 *  and numbers that are too small do not work well with this algorithm.
 *
 *  "print" statements that can be ommitted are left justified.
 *
 *  by Contini - Fri March 3rd, 2000.
 */

function SelectPolynomial(n);
    /* select a degree according to the asymptotic formula */
    d := Round( Root( 3 * Log(n)/Log(Log(n)) , 3 ) );
print "degree ", d;

    m := Iroot( n, d );
print "m := ", m, ";";

    coeffs := [];
    for i := d to 0 by -1 do
	temp2 := m^i;
	coeff := n div temp2;
	/* Append (~coeffs, coeff); */
	coeffs[i+1] := coeff;
	n -:= coeff*temp2;
    end for;

    P<X> := PolynomialRing( IntegerRing() );
    f := P!coeffs;
    return f, m, d;

end function;


procedure RationalSieve( ~sieve_array, y, FB, log_primes, ~prev_roots,
  m_mod_p, sieve_len);
    /*  
     *  RATIONAL SIEVE
     */
    sieve_array := [0 : i in [1..sieve_len]];

    for i in [1..#FB] do
	p := FB[i];
	logp := log_primes[i];
	/*
	 *  Compute  loc  - the first location in the sieve array
	 *  where  p  divides the corresponding residue
	 */
	loc := prev_roots[i] - m_mod_p[i];
	if loc lt 0 then
	    loc +:= p;
	end if;
	prev_roots[i] := loc;
	/*  sieve begins at location 1, so skip location 0 */
	if loc eq 0 then
	    loc := p;
	end if;

	/*  MAIN SIEVE LOOP  */
	while loc le sieve_len do
	    sieve_array[loc] +:= logp;
	    loc +:= p;
	end while;

    end for;

end procedure;


procedure AlgebraicSieve( ~nf_sieve_array, y, ~nf_prev_roots, log_primes,
  R_p_set, sieve_len );
    /*  
     *  NUMBER FIELD SIEVE
     */
    nf_sieve_array := [0 : i in [1..sieve_len]];
    for i in [1..#R_p_set] do
	R_p := R_p_set[i];
	p := R_p[1];
	r := R_p[2];
	logp := R_p[3];

	loc := nf_prev_roots[i] - r;
	if loc lt 0 then
	    loc +:= p;
	end if;
	nf_prev_roots[i] := loc;
	if loc eq 0 then
	    loc := p;
	end if;

	/*  MAIN SIEVE LOOP  */
	while loc le sieve_len do
	    nf_sieve_array[loc] +:= logp;
	    loc +:= p;
	end while;

    end for;
end procedure;


procedure ScanSieveArray( J, b, bm, y, d, sieve_array, nf_sieve_array,
  error_term, ~num_rels, sieve_len, ~M, FB, prev_roots, nf_prev_roots,
  R_p_set, enough_rels, ~smooth_alg_elts, char_cols, char_offset );
	/*  scan sieve array
	 */
    threshold := Round( Log( bm ) ) - error_term;
    zero := [0: i in [1..d]];
    x := zero;
    x[1] := 1; x[2] := b;
    nf_threshold := Round( Log( Abs( Norm( J!x ) ) ) ) - error_term;
    for a in [1..sieve_len] do
	if sieve_array[a] ge threshold then
	    /*  rational side appears to be smooth */
	    if nf_sieve_array[a] ge nf_threshold then
		/*  number field side appears to be smooth */
		if Gcd(b, a) ne 1 then
		    /*  relation is redundant */
		    continue;
		end if;

		/*  see if it really is smooth */
		relation := [];
		rational_side := a + bm;
		for i in [1..#FB] do
		    p := FB[i];
		    if a mod p eq prev_roots[i] then
			e := 0;
			quot, rem := Quotrem( rational_side, p );
			repeat
			    rational_side := quot;
			    e +:= 1;
			    quot, rem := Quotrem( rational_side, p );
			until rem ne 0;
			Append( ~relation, [i, e] );
		    end if;
		end for;
		if Abs(rational_side) ne 1 then
		    /*  rational side is not smooth :-( */
		    continue;
		end if;

		/*  try algebraic side */
		x := zero;
		x[1] := a; x[2] := b;
		alg_elt := J!x;
		norm := IntegerRing()!Norm(alg_elt);

		for i in [1..#R_p_set] do
		    R_p := R_p_set[i];
		    p := R_p[1];
		    if a mod p eq nf_prev_roots[i] then
			e := 0;
			quot, rem := Quotrem( norm, p );
			repeat
			    norm := quot;
			    e +:= 1;
			    quot, rem := Quotrem( norm, p );
			until rem ne 0;
			Append( ~relation, [i+#FB, e] );
		    end if;
		end for;
		if Abs(norm) ne 1 then
		    /*  algebraic side is not smooth :-( */
		    continue;
		end if;

		/*  The relation is a keeper! */
print "a := ", a, "; b:= ", b,";";
		num_rels +:= 1;

		/*  Compute the character columns corresponding to
		    alg_elt  */
		for i in [1..#char_cols] do
		    p := char_cols[i][1];
		    r := char_cols[i][2];
		    if LegendreSymbol( (a + b*r) mod p, p ) eq -1 then
			M[num_rels, i + char_offset] := 1;
		    end if;
		end for;

		/*  store the algebraic element */
		Append( ~smooth_alg_elts, alg_elt );

		/*  store relation in matrix */
		for j in [1..#relation] do
		    M[num_rels, relation[j, 1]] := relation[j, 2];
		end for;

		if num_rels eq enough_rels then
		    break;
		end if;

		/*  update the  nf_threshold  - since the norms
		    grow somewhat quickly, we waste a lot of time
		    on false reports if this update is not made */
		x := zero;
		x[1] := a; x[2] := b;
		nf_threshold := Round( Log( Abs( Norm( J!x ) ) ) );
		nf_threshold -:= error_term;
	    end if;
	end if;
    end for;
end procedure;


function NFS( n )

    f, m, d := SelectPolynomial(n);
    fprime := Derivative( f );

    J<alpha> := NumberField( f );

    /*  choose  y  - the bound for smoothness  */
    dlogd := d*Log(d);
    temp := 1.0/d * Log( n );
    e := dlogd + Sqrt( dlogd^2 + 4*temp*Log(temp) );
    y := Round( Exp( 0.5*e ) );
    print "smoothness bound is ", y;

    /*  prepare data for rational sieve */
    FB := [p : p in [2..y] | IsPrime( p ) ];
    log_primes := [Round(Log(p)) : p in FB];
    prev_roots := [0: p in FB];
    m_mod_p := [m mod p: p in FB];

    /*  prepare data for the algebraic sieve */
    cols := #FB;
    R_p_set := [];
    nf_prev_roots := [];
    for p in FB do
	K := GF(p);
	roots_mod_p := Roots( f, K );
	for i in [1..#roots_mod_p] do
	    r := IntegerRing()!roots_mod_p[i][1];
	    /*  also include the rounded logarithm of  p  in the storage */
	    Append( ~R_p_set, [p, r, Round(Log(p))] );
	    Append( ~nf_prev_roots, 0 );
	    cols +:= 1;
	end for;
    end for;
    smooth_alg_elts := [];

    /*  create the character columns */
    char_offset := cols;
    num_char_cols := 3 * Ilog2( n );
    char_cols := [];
    p := NextPrime( y );
    i := 0;
    repeat
	K := GF(p);
	roots_mod_p := Roots( f, K );
	for j in [1..#roots_mod_p] do
	    r := IntegerRing()!roots_mod_p[j][1];
	    Append( ~char_cols, [p, r] );
	    i +:= 1;
	    if i eq num_char_cols then
		break;
	    end if;
	end for;

	p := NextPrime( p );
    until i eq num_char_cols;
    cols +:= num_char_cols;

    enough_rels := cols+50;
print "will quit after", enough_rels," relations";

    /*  allow for an error in the sieve threshold */
    error_term := 3;

    /*  set up matrix */
    M := RMatrixSpace (IntegerRing(), enough_rels, cols)!0;

    num_rels := 0;
    b := 1;
    /*  bm = b*m */
    bm := m;
    sieve_len := y * 25;

    while true do
	RationalSieve( ~sieve_array, y, FB, log_primes, ~prev_roots,
	  m_mod_p, sieve_len );

	AlgebraicSieve( ~nf_sieve_array, y, ~nf_prev_roots, log_primes,
	  R_p_set, sieve_len );

	ScanSieveArray( J, b, bm, y, d, sieve_array, nf_sieve_array,
	  error_term, ~num_rels, sieve_len, ~M, FB, prev_roots,
	  nf_prev_roots, R_p_set, enough_rels, ~smooth_alg_elts,
	  char_cols, char_offset );

	bm +:= m;
	b +:= 1;
print "Number of relations: ", num_rels;

	if num_rels ge enough_rels then
	    break;
	end if;
    end while;


print "matrix size: ",num_rels," by ", cols;
print "finding null space";

    gf2 := GF(2);
    Mprime := RMatrixSpace (gf2, num_rels, cols)!M;
    time null_space := NullSpace( Mprime );

    /*  "free" Mprime */
    Mprime := 0;

    for j in  [1..Dimension(null_space)] do
print "trying null space vector #",j;

	/*  get the first element of the solution vector */
	soln := ChangeUniverse (Eltseq (null_space.j), IntegerRing());

	/*  compute the exponents of a relation that involves an
	    algebraic square and an integer square
	 */

	exps := RSpace(IntegerRing(), cols)!0;
	for i in [1..num_rels] do
	    if soln[i] eq 1 then
		exps +:= M[i];
	    end if;
	end for;

	/*  compute the square root  x  of the integer element of the relation
	 */
	x := 1;
	for i in [1..#FB] do
	    e := exps[i];
	    if e gt 0 then
		x := (x * ( FB[i]^(e div 2) ) ) mod n;
	    end if;
	end for;
	x := ( x * Evaluate( fprime, m ) ) mod n;
print x,"^2 mod n = ", x^2 mod n;

	list := [];
	for i in [1..num_rels] do
	    if soln[i] eq 1 then
		Append( ~list, smooth_alg_elts[i] );
	    end if;
	end for;
	Append( ~list, Evaluate( fprime, alpha )^2 );

	/*  Compute the square root  beta  of the algebraic element */
	ok, l := Root(2, list, []);
        if ok eq false then
            print "Element is not a square... continuing!";
            continue;
        end if;

	/*  map the sqrt to an integer  v  modulo  n  */
	v := 1;
	for i in [1..#l] do
	    list := l[i];
	    prod_list := 1;
	    for j in [1..#list] do
		seq := Eltseq( list[j] );
		elt := 0;
		for k in [0..d-1] do
		    num := Numerator( seq[k+1] ) mod n;
		    den := Denominator( seq[k+1] ) mod n;
		    elt +:= ( ( num * Modinv( den, n ) mod n) *
		      Modexp( m, k, n ) ) mod n;
		end for;
		prod_list := (prod_list * elt) mod n;
	    end for;
	    v := v * Modexp( Modinv( prod_list, n ), 2^(i-1), n ) mod n;
	end for;
	
print v,"^2 mod n = ", v^2 mod n;

	gcd := Gcd(x-v, n);
print "Gcd of ",x-v," and ",v," is ", gcd;
	if gcd ne 1 and gcd ne n then
	    print "FACTOR: ", gcd;
	    return gcd;
	end if;
    end for;

print "no factor found :-( ";
    return 0;

end function;