/*
 *  An implementation of algorithm 12.1 in the paper
 *  "Factoring Integers with the Number Field Sieve" by J.P. Buhler,
 *  H.W. Lenstra, and Carl Pomerance.  This is the variant of the
 *  Number Field Sieve that uses a homogenous polynomial.
 *
 *  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.  It is
 *  recommended to try this on numbers of no bigger than 30-digits.
 *  There are various problems that can happen for larger numbers,
 *  such as running out of memory due to not taking advantage of
 *  the sparseness of the matrix!
 *
 *  Three methods for selecting the polynomial are included here.
 *  None of them are optimized or complete.  However, the polynomials
 *  should work most of the time.  This implementation of NFS takes  n
 *  along with the polynomial  f, m1, m2,  and  d  for input.
 *
 *  "print" statements that can be ommitted are left justified.
 *
 *  by Contini - March 23, 2000.
 *  updated: November 12, 2003.
 */



/*
 *  This method of selecting the polynomial takes  m2 = 1  but allows
 *  c_d != 1 .
 */
function SelectPolynomial1(n);
    /* select a degree according to the asymptotic formula */
    d := Round( Root( 3 * Log(n)/Log(Log(n)) , 3 ) );
print "degree ", d;

    P<X,Y> := PolynomialRing( IntegerRing(), 2 );

    m2 := 1;
    m1 := Iroot( n, d+1 ) ;

    coeffs := [];
    for i := d to 0 by -1 do
        temp := m1^i;
        coeff := n div temp;
        Append (~coeffs, coeff);
        n -:= coeff*temp;
    end for;
 
    m1div2 := m1 div 2;
    /* adjust the coefficients to be between  -m1div2  and  m1div2  */
    for i := d+1 to 2 by -1 do
	if coeffs[i] gt m1div2 then
	    coeffs[i] -:= m1;
	    coeffs[i-1] +:= 1;
	end if;
    end for;

    f := P!0;
    for i in [0..d] do
	f +:= coeffs[i+1]*X^(d-i)*Y^i;
    end for;

print f, m1, m2;

    return f, m1, m2, d;
end function;


/*
 *  This method of selecting the polynomial takes  c_d = 1  but allows
 *  m2 != 1 .  See (12.3) in the paper of Buhler, Lenstra, and Pomerance.
 */
function SelectPolynomial2(n);
    /* select a degree according to the asymptotic formula */
    d := Round( Root( 3 * Log(n)/Log(Log(n)) , 3 ) );
print "degree ", d;

    P<X,Y> := PolynomialRing( IntegerRing(), 2 );

    /* A not very optimized algorithm for choosing  m1  and  m2  */
    m1 := Iroot( n, d+1 );
    closeness_ratio := 0.90;
    min_m2 := Round( m1 * closeness_ratio );
    max_m2 := Round( m1 * 1/closeness_ratio );

    while true do
	m1 +:= 1;
	tmp := n - m1^d;
	m2 := 0;
	for j in Divisors( tmp ) do
	    if j ge min_m2 and j le max_m2 then
		if Gcd( m1, j ) ne 1 then
		    continue;
		end if;
		if Abs( m1 - j ) lt Abs( m1 - m2 ) then
		    m2 := j;
		end if;
	    end if;
	end for;
	if m2 ne 0 then
	    break;
	end if;
    end while;

    f := P!0;
    quot := (n-m1^d) div m2;
    rem := quot;
    if false then
	m1div2 := m1 div 2;
	for i in [0..d-2] do
	    tmp := m2^(d-1-i);
	    c_i := (rem * Modinv( tmp, m1 ) ) mod m1;
	    if c_i gt m1div2 then
		c_i -:= m1;
	    end if;
	    rem := (rem - c_i*tmp) div m1;
	    f +:= c_i*X^i*Y^(d-i);
	end for;
	f +:= rem*X^(d-1)*Y;
    else
	for i in [0..d-2] do
	    tmp := m1^(d-1-i);
	    c_i := (rem * Modinv( tmp, m2 ) ) mod m2;
	    rem := (rem - c_i*tmp) div m2;
	    f +:= c_i*X^(d-1-i)*Y^(i+1);
	end for;
	f +:= rem*Y^d;
    end if;
    f := f + X^d;

print f, m1, m2;

    return f, m1, m2, d;
end function;



/*
 *  This method of selecting the polynomial allows  c_d != 1  and
 *  m2 != 1 .  See (12.4) in the paper of Buhler, Lenstra, and Pomerance.
 */
function SelectPolynomial3(n);
    /* select a degree according to the asymptotic formula */
    d := Round( Root( 3 * Log(n)/Log(Log(n)) , 3 ) );
print "degree ", d;

    P<X,Y> := PolynomialRing( IntegerRing(), 2 );
    m1 := Iroot( n, d+1 );
    m2 := m1 + 1;

    L := RMatrixSpace( IntegerRing(), d+2, d+3 )!0;
    L[1][1] := n^2;
    L[1][2] := 1;
    for i in [0..d] do
	L[i+2][1] := n*m1^i * m2^(d-i);
	L[i+2][i+3] := 1;
    end for;
    RL := LLL(L);

    j := 0;
    for i := 1 to d+2 do
	if RL[i][1] eq 0 and RL[i][2] ne 0 then
	    j := i;
	    break;
	end if;
    end for;
    if j eq 0 then
	print "ERROR- LLL failed";
	return 0;
    end if;

    f := 0;
    for i in [0..d] do
	f +:= RL[1][i+3]*X^i*Y^(d-i);
    end for;
 
    /*  need to put in checks to make sure we get a good polynomial */

print f, m1, m2;

    return f, m1, m2, d;

end function;


/*
 *  Find the values of  a*m1 - b*m2  that are smooth
 */
procedure RationalSieve( ~sieve_array, y, FB, log_primes, ~prev_roots,
  m1_times_m2_inv, sieve_len, b );
    /*  
     *  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 := m1_times_m2_inv[i];
	if loc eq -1 then
	    /*  m2  is divisible by p */
	    if b mod p ne 0 then
		/*  p  does not divide any of the values in the sieve array */
		continue;
	    end if;
	    /*  if we get here,  p  divides every location in the sieve
		array */
	    for loc in [1..sieve_len] do
		sieve_array[loc] +:= logp;
	    end for;
	    continue;
	end if;

	loc +:= prev_roots[i];
	if loc ge p 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;


/*
 *  Determine values of  a  such that  f(a, b)  is divisible by the
 *  primes in the factor base.  Sieve on these values to find smooth
 *  elements.
 */
procedure AlgebraicSieve( ~nf_sieve_array, y, ~nf_prev_roots, log_primes,
  R_p_set, sieve_len, b );
    /*  
     *  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];

	if r eq -1 then
	    /*  -1 represents the infinity element */
	    if b mod p ne 0 then
		/*  p  does not divide anywhere in the sieve array */
		continue;
	    end if;
	    /*  if we get here,  p  divides every location in the sieve
		array */
	    for loc in [1..sieve_len] do
		nf_sieve_array[loc] +:= logp;
	    end for;
	    continue;
	end if;

	loc := nf_prev_roots[i] + r;
	if loc ge p 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, f, m2, m1, b, 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, cols, c_d );

    bm1 := b*m1;
    threshold := Round( Log( Abs( m2 - bm1 ) ) ) - error_term;
    zero := [0: i in [1..d]];
    nf_threshold := Round( Log( Abs( Evaluate(f, [1, b]) ) ) ) - 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 rational side really is smooth */
		relation := [];
		rational_side := a*m2 - bm1;
		for i in [1..#FB] do
		    p := FB[i];
		    pr := prev_roots[i];
		    if a mod p eq pr or (pr eq -1 and b mod p eq 0) 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 */
		f_ab := IntegerRing()!Evaluate(f, [a, b]);

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

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

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

		/*  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];
		    ls := LegendreSymbol( (c_d * (a - b * r) ) mod p, p );
		    if ls eq -1 then
			M[num_rels, i + char_offset] := 1;
		    end if;
		end for;

		/*  store the algebraic element */
		x := zero;
		x[1] := c_d * a; x[2] := -b;
		alg_elt := J!x;
		Append( ~smooth_alg_elts, alg_elt );

		/*  store the sign */
		if rational_side eq -1 then
		    M[num_rels, cols-1] := 1;
		end if;
		/*  this column is to assure that an even  number of
		    elements are used in the final relation */
		M[num_rels, cols] := 1;

		if num_rels eq enough_rels then
		    break;
		end if;

		/*  update the thresholds to minimize false reports */
		f_ab := IntegerRing()!Evaluate(f, [a, b]);
		nf_threshold := Round( Log( Abs( f_ab ) ) );
		nf_threshold -:= error_term;
		threshold := Round( Log( Abs( a*m2 - bm1 ) ) ) - error_term;
	    end if;
	end if;
    end for;
end procedure;



/*
 *  This version requires you to send in the polynomial  f
 *  along with  m1, m2,  and  d .
 */
function NFS( n, f, m1, m2, d )

    f_X := Derivative( f, 1 );

    /*  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];
    m1_times_m2_inv := [];
    i := 0;
    for p in FB do
	i +:= 1;
	if m2 mod p eq 0 then
	    /*  flag no solution */
	    Append( ~m1_times_m2_inv, -1 );
	    prev_roots[i] := -1;
	else
	    Append( ~m1_times_m2_inv, (Modinv(m2, p) * m1) mod p );
	end if;
    end for;
print "primes in FB: ",#FB;

    /*  prepare data for the algebraic sieve */
    Q<X> := PolynomialRing( Integers() );
    f1 := Q!Evaluate( f, [X,1] );
    c_d := Coefficient( f1, d );
    f_c_d := Q!Evaluate( f, [X, c_d] );
    J<omega> := NumberField( f_c_d div c_d );

    cols := #FB;
    R_p_set := [];
    nf_prev_roots := [];
    for p in FB do
	K := GF(p);
	roots_mod_p := Roots( f1, 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;
	if c_d mod p eq 0 then
	    /*  add the infinity element, which will be represented by -1 */
	    Append( ~R_p_set, [p, -1, Round(Log(p))] );
	    Append( ~nf_prev_roots, -1 );
	    cols +:= 1;
	end if;
    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
	while c_d mod p eq 0 do
	    p := NextPrime( p );
	end while;
	K := GF(p);
	roots_mod_p := Roots( f1, 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;

    /*  Add one column for the sign of  a*m2 - b*m1  and one column
    which will always be 1 to make sure that the number of relations
    in any dependency is even */
    cols +:= 2;

    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;
    sieve_len := y;

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

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

	ScanSieveArray( J, f, m2, m1, b, 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, cols, c_d );

        b +:= 1;

print "b = ",b,",         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;

	S_size := 0;
	list := [];
	for i in [1..num_rels] do
	    if soln[i] eq 1 then
		Append( ~list, smooth_alg_elts[i] );
		S_size +:= 1;
	    end if;
	end for;
	Append( ~list, (Evaluate( f_X, [omega, c_d] )/c_d)^2 );

	/*  Compute the number field element of the relation.
	    In a real implementation, we wouldn't actually multiply this
	    element out.
	 */
	prod := &*list;

	/*  Compute the square root of the algebraic element */
	beta := Root(prod, 2);

	/* map  beta  to an integer modulo n */
	m2inv := Modinv( m2, n );
	seq := Eltseq(beta);
	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( (c_d * m1 * m2inv) mod n, k, n ) ) mod n;
	end for;
	v := (elt * Modexp( m2, d-1, n) ) mod n;

	/*  compute the integer element of the relation.
	 */
	w := 1;
	for i in [1..#FB] do
	    e := exps[i];
	    if e gt 0 then
		w := (w * ( FB[i]^(e div 2) ) ) mod n;
	    end if;
	end for;

	h := Modexp( c_d, d-2+(S_size div 2), n ) * Evaluate( f_X, [m1, m2] )
	  mod n;
	l := Modexp( m2, S_size div 2, n );

	tmp := (h*w-l*v) mod n;
	gcd := Gcd(tmp, n);
print "Gcd of ",tmp," and ",n," 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;