/*
 *  An implementation of a search for a non-skewed NFS polynomial, based
 *  on Brian Murphy's thesis.
 *
 *  By Contini - May, 2000.
 *
 *  Please report bugs, improvements, and suggestions to
 *  contini@maths.usyd.edu.au
 */



/*
 *  A polynomial for  n  can be reconstructed from its degree  d  and its
 *  value of  m1  (assuming  m2 = 1).    This function does that.
 */
function m1ToPoly( n, m1, d )

    rem := n;
    coeffs := [];
    for i := d to 0 by -1 do
        temp := m1^i;
        coeff := rem div temp;
        Append (~coeffs, coeff);
        rem -:= 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;
 
    P<X, Y> := PolynomialRing(Integers(), 2);
    F := P!0;
    for i in [0..d] do
        F +:= coeffs[i+1]*X^(d-i)*Y^i;
    end for;

    return F;
end function;



/*
 *  Compute  cont_p  (defined in section 3.2.1 of Murphy's thesis) of
 *  homogeneous polynomial  F  corresponding to univariate polynomial  f
 *  for the prime  p .  The following input values are required:
 *      - disc  should be the discriminant of  f .
 *      - bound is the size of the sample space to look through in
 *      order to estimate  cont_p  for the case of  p  dividing  disc .
 *      For example,  bound=10000  may be a good choice.
 */
function ComputeCont_p( p, disc, bound, F, f );
    if disc mod p eq 0 then
        /* estimate using some small sample space */
        cont_p := 0;
        upper_limit := bound^2;
        for i in [1..bound] do
            repeat
                a := Random( upper_limit );
                b := Random( upper_limit );
            until Gcd( a, b ) eq 1;
            value := Evaluate( F, [a, b] );
            while value mod p eq 0 do
                value div:= p;
                cont_p +:= 1;
            end while;
        end for;
        cont_p := cont_p / (0.0 + bound);
    else
        K := GF(p);
        roots_mod_p := Roots( f, K );
        q_p := #roots_mod_p;
        /* add in the contribution of the projective roots
         */
        d := Degree( f );
	if Coefficient(f, d) mod p eq 0 then
	    q_p +:= 1;
	end if;
            
        cont_p := q_p * p / (p * p - 1.0);
    end if;

    return cont_p;
end function;



/*
 *  Compute the  alpha  value corresponding to polynomial  F .  For
 *  primes  p  that divide the discriminant, do a random sampling of
 *  bound  elements to estimate  cont_p  (the expected exponent of
 *  prime  p  on the set of  F  values).  Include an early abort strategy:
 *  if by the time we reach  abort_p , the estimated  alpha  value is >
 *  abort_alpha  then stop the computation and return (the estimated)  alpha .
 *  If the early abort will not be used, then just set  abort_p  and
 *  abort_alpha  to 0.
 *  Function uses primes up to  max_p .
 */
function ComputeAlpha( F, bound, abort_p, abort_alpha, max_p );
 
    Q<X> := PolynomialRing(Integers());
    f := Q!Evaluate( F, [X, 1] );
 
    disc := Discriminant( f );
    a_d := Coefficient(f, Degree(f));
    primes := [p : p in [1..max_p] | IsPrime(p)];
 
    alpha := 0.0;
 
    for p in primes do
        cont_p := ComputeCont_p( p, disc, bound, F, f );
 
        alpha +:= (1.0/(p-1) - cont_p)*Log(p);
 
        /* early-abort strategy */
        if p eq abort_p and alpha gt abort_alpha then
            break;
        end if;
    end for;
 
    return alpha;
 
end function;


/*
 *  Compute  rho( x )  where  rho  is the Dickman rho function.
 *  This can be optimized by precomputing the coefficients of
 *  the polynomials and storing those coefficients for future calls.
 */
function DickmanRho( x );

    if x le 1.0 then
	return 1.0;
    end if;
    k := Ceiling( x );
    zeta := k - x;
 
 
    c := [];
    c[1] :=  1.0 - Log( 2.0 );
    for i in [1..54] do
	c[i+1] := 1.0 /(i * 2^i);
    end for;

    prev_c := [];
    l := 2;
    while l lt k do
	for i in [0..54] do
	    prev_c[i+1] := c[i+1];
	end for;

	for i in [1..54] do
	    sum := 0.0;

	    for j in [0..i-1] do
		sum := sum + prev_c[j+1]/(i * (l + 1.0)^(i-j) );
	    end for;

	    c[i+1] := sum;
	end for;

	sum := 0.0;
	for j in [1..54] do
	    sum +:= c[j+1]/( j + 1.0 );
	end for;
	c[1] := sum/l;

	l +:= 1;
    end while;
 
    rho_k := 0.0;
    for i in [0..54] do
	rho_k +:= c[i+1] * zeta^i;
    end for;

    return rho_k;

end function;


/*
 *  Quickly compute an approximation to the Dickman rho function.
 */
function ApproxDickmanRho( u );

    if u le 1.0 then
	return 1.0;
    end if;
    return u^-u;

end function;


/*
 *  Rate the polynomial pair  F1  and  F2 , where  F1  has degree  d
 *  and corresponding  alpha  value and  F2  is linear.  B1  and  B2
 *  represent upperbounds for the primes in the algebraic and rational
 *  factor bases, respectively.  High ratings tend to mean better
 *  polynomials.
 *
 *  Actually, we modify his rating method in 2 ways.
 *  First, we divide through by  K  (the number of intervals) so that
 *  the rating is not proportional to the number of intervals we choose.
 *  Second, we add in a  log B1^d  to the numerator of  u1  (which is
 *  equivalent to adding  d  to the whole quantity of  u1 ) so that the
 *  rho values are more typical of what will happen in practice.  A
 *  similar modification is done to the  u2  computation.
 */
function RatePolynomials( F1, F2, alpha, B1, B2, d )

    pi := Pi(RealField());
    K := 1000;

    rating := 0.0;
    logB1 := Log( B1 );
    logB2 := Log( B2 );
    for i in [1..K] do
	theta := pi/K*(i - 0.5);
	tmp := Abs( Evaluate( F1, [Cos( theta ), Sin( theta ) ] ) );
	u1 := ( Log( tmp ) + alpha ) / logB1 + d;
	tmp := Abs( Evaluate( F2, [Cos( theta ), Sin( theta ) ] ) );
	u2 := Log( tmp ) / logB2 + 1.0;
	rating +:= ApproxDickmanRho( u1 ) * ApproxDickmanRho( u2 );
    end for;

    return rating/K;

end function;



/*
 *  Let  X  be a sequence, where  X[i]  represents the ranking
 *  of the  i'th  best polynomial under some ranking scheme.  This
 *  function computes the correlation coefficient of the ranking scheme
 *  (i.e. it gives us an idea of how good the ranking scheme is).
 *
 *  for example, to compute the correlation coefficient of the E-rank
 *  in table 6.4 of Murphy's thesis, you would send in  X :=
 *  [ 1, 3, 6, 7, 2, 8, 4, 10, 5, 9, 15, 11, 18, 13, 12, 14, 17, 16 ] .
 *
 *  Remark: there seem to be small differences between Murphy's and my
 *  calculations.
 */
function CorrelationCoeff( X )
    /*  compute  eXY  =  expected value of  X * Y  where Y is identity
        perumtation */
    eXY := &+[X[i]*i: i in [1..#X]]/(#X + 0.0);
    eX := (#X+1.0)/2.0;
    eX2 := &+[X[i]*X[i]: i in [1..#X]]/(#X + 0.0);
    sigma := Sqrt( eX2 - eX*eX );
    return (eXY - eX*eX)/sigma^2;

end function;


/*
 *  Choose a leading coefficient for our polynomial.  We want this leading
 *  coefficient to have lots of small prime factors, and to be close to
 *  high_a_d .  This procedure uses randomness, so that we should get a
 *  new leading coefficient every time it is called (assuming it is not
 *  called too many times).
 */
function ChooseLeadingCoeff( high_a_d );

    /*  Start with some value for the leading coefficient  c  that has
     *  many small prime factors.  However, don't let this start value
     *  get too close to  high_a_d , since this can severely limit the
     *  number of unique leading coefficients we can get.
     */
    c := 32 * 27 * 25 * 49 ;
    c_limit := high_a_d div 300;
    if 10*c gt c_limit then
	c := 2 * 3 * 5 * 7;
	if 10*c gt c_limit then
	    c := 2 * 3;
	    if 10*c gt c_limit then
		c := 2;
	    end if;
	end if;
    end if;

    if high_a_d lt 1000 then
	primes := [2, 3, 5, 7, 11];
    else
	primes := [2, 3, 5, 7, 11, 13, 17, 19, 23, 29];
    end if;
    num_primes := #primes;

    min := high_a_d div 2;
    while c lt min do
	p := primes[ Random( 1, num_primes ) ];
	if c*p le high_a_d then
	    c := c * p;
	end if;
    end while;

    return c;

end function;


/*
 *  Let  n  be the number we want to factor.  This algorithm works
 *  something like Murphy's procedure 5.1.4 to find a good NFS polynomial.
 *  (though I've modified it to fit my own understanding, so there are
 *  definitely some differences from procedure 5.1.4).  Remark: this
 *  function alone randomly chooses one leading coefficient  a_d  for
 *  the polynomials being searched.  The user should try calling this
 *  several times so several leading coefficients are tried.
 *
 *  Returns a sequence of tuples which contain candidate polynomials and
 *  information about those candidates.  More precisely, the tuple has
 *  the following items (in this order):
 *	- the polynomial  F  ( F  is homogeneous).
 *	- a rating for the polynomial  F  (see 5.7 in Murphy's thesis).
 *      - the alpha value (see section 3.2.3 of Murphy's thesis) of candidate
 *	polynomial  F .
 *	- the  m1  value corresponding to polynomial  F (i.e  F(m1, m2) = n ).
 *	- the  m2  value corresponding to polynomial  F .  This function
 *	always has  m2=1 .
 *
 *  Description of input variables:
 *	n		the number we want to factor.
 *
 *	chi		only consider polynomials that have all their
 *			leading coefficients  <=  m1 * chi
 *
 *	max_k		for each choice of the leading coefficient, only
 *			consider polynomials whose  m1  lies within this
 *			range: [best_m1-max_k..best_m1+max_k] , where
 *			"best_m1" indicates the  m1  which minimizes the
 *			coefficient  a_{d-1}  of polynomial  F  (where  d
 *			is the degree of  F ).  This is used to limit the
 *			amount of time spent on a single leading coefficient.
 *			Also the number of  m1  considered may be even
 *			smaller than the range indicated above since some of
 *			these  m1  may correspond to a different leading
 *			coefficient.
 *
 *	max_alpha	The maximum value of  alpha  to allow for a polynomial
 *			to be a candidate.  Generally you want this to be
 *			a negative number, say -1.0 or smaller.  Values less
 *			than -3.0 are considered very good.
 *
 *	max_cands	maximum number of candidate polynomials to allow.
 *			(this is used to limit the amount of time spent
 *			on a single choice for the leading coefficient).
 *
 */
function PolySearch1( n, chi, max_k, max_alpha, max_cands );

    /* The set of return polynomials */
    candidates := [];

    /* select a degree according to the asymptotic formula */
    d := Round( Root( 3 * Log(n)/Log(Log(n)) , 3 ) );
print "degree ", d;
 
    /* Compute  B  - an upper bound for the factor base primes.  This
       will be used in Murphy's rating scheme for the candidate
       polynomials */
    dlogd := d*Log(d);
    temp := 1.0/d * Log( n );
    e := dlogd + Sqrt( dlogd^2 + 4*temp*Log(temp) );
    B := Round( Exp( 0.5*e ) );

    P<X, Y> := PolynomialRing(Integers(), 2);
 
    m2 := 1;
    /*  Need to choose  m1  between  n^[1/(d+1)]  and  n^[1/d] .
	Murphy's procedure 5.1.4 actually requests that  m1  is
	significantly larger than  n^[1/(d+1)].  (This is implied
	by the request that  a_d  is significantly smaller than its
	corresponding  m1 ).  We allow the lowest  m1  to be the
	geometric mean of these two values.  The use of the geomretic
	mean just seems to work well - it can be changed.
     */
    low_m1 := Iroot( Iroot( n, d+1 ) * Iroot( n, d ) , 2 );
    high_m1 := Iroot( n, d );
    low_a_d := n div high_m1^d;
    high_a_d := n div low_m1^d;
 
    c := ChooseLeadingCoeff( high_a_d );
print "using leading coeff ", c;
    a_d := c;
    if a_d lt low_a_d then
	a_d := ((low_a_d div c) + 1) * c;
    end if;

    /* The value of  m1  that causes  a_{d-1} to be minimized */
    best_m1 := Iroot( n div a_d, d );
    /* Determine  a_{d-1}  for this value of  best_m1 */
    a_dminus1 := (n - a_d*best_m1^d) div best_m1^(d-1);
    coeff_bound := Round( best_m1 * chi );
    /* we try values between values [best_m1-k...best_m1+k] */
    k := (coeff_bound - a_dminus1) div (a_d * d);

    if k gt max_k then
	k := max_k;
    end if;
print "k=",k;

    for m1 in [best_m1-k..best_m1+k] do

	rem := n;
	coeffs := [];
	for i := d to 0 by -1 do
	    temp := m1^i;
	    coeff := rem div temp;
	    Append (~coeffs, coeff);
	    rem -:= 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;
     
	/* check that the coefficients are sufficiently small */
	small_coeffs := true;
	for i in [1..#coeffs] do
	    if Abs( coeffs[i] ) gt coeff_bound then
		small_coeffs := false;
		break;
	    end if;
	end for;

	if small_coeffs eq false then
	    continue;
	end if;

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

	/* obtain a quick estimate for  alpha  */
	alpha := ComputeAlpha( F, 2500, 11, max_alpha+0.5, 100 );
print m1, alpha;

	if alpha le max_alpha then
	    /*  get a more accurate value for  alpha  */
	    alpha := ComputeAlpha( F, 10000, 0, 0, 2000 );
print "recomputed alpha: ", alpha;
	    if alpha gt max_alpha then
		continue;
	    end if;
	    /* create the linear polynomial */
	    F2 := X*m2 - Y*m1;
	    rating := RatePolynomials( F, F2, alpha, B, B, d );
	    Append( ~candidates, <F, rating, alpha, m1, m2> );
	    if #candidates eq max_cands then
		break;
	    end if;
print F;
print "";
	end if;

    end for;
    return candidates;
end function;





/*
 *  Let  n  be the number we want to factor.  This algorithm calls a
 *  function named PolySearch1() which is similar to Murphy's procedure 5.1.4
 *  to find a good NFS polynomial.  Each call to PolySearch1() uses a
 *  (hopefully) new leading coefficient, that is somewhat randomly chosen.
 *  Up to  max_leading_coeffs  will be tried.
 *
 *  Returns a sequence of tuples which contain candidate polynomials and
 *  information about those candidates.  More precisely, the tuple has
 *  the following items (in this order):
 *	- the polynomial  F  ( F  is homogeneous).
 *	- a rating for the polynomial  F  (see 5.7 in Murphy's thesis).
 *      - the alpha value (see section 3.2.3 of Murphy's thesis) of candidate
 *	polynomial  F .
 *	- the  m1  value corresponding to polynomial  F (i.e  F(m1, m2) = n ).
 *	- the  m2  value corresponding to polynomial  F .  This function
 *	always has  m2=1 .
 *
 *  Description of input variables:
 *	n		the number we want to factor.
 *
 *	chi		only consider polynomials that have all their
 *			leading coefficients  <=  m1 * chi
 *
 *	max_k		for each choice of the leading coefficient, only
 *			consider polynomials whose  m1  lies within this
 *			range: [best_m1-max_k..best_m1+max_k] , where
 *			"best_m1" indicates the  m1  which minimizes the
 *			coefficient  a_{d-1}  of polynomial  F  (where  d
 *			is the degree of  F ).  This is used to limit the
 *			amount of time spent on a single leading coefficient.
 *			Also the number of  m1  considered may be even
 *			smaller than the range indicated above since some of
 *			these  m1  may correspond to a different leading
 *			coefficient.
 *
 *	max_alpha	The maximum value of  alpha  to allow for a polynomial
 *			to be a candidate.  Generally you want this to be
 *			a negative number, say -1.0 or smaller.  Values less
 *			than -3.0 are considered very good.
 *
 *	max_cands	maximum number of candidate polynomials to allow.
 *			(this is used to limit the amount of time spent
 *			on a single choice for the leading coefficient).
 *
 */
function MurphyPolynomial1( n, chi, max_k, max_alpha, max_cands,
  max_leading_coeffs );

    candidates := [];
    for i in [1..max_leading_coeffs] do
	new_cands := PolySearch1( n, chi, max_k, max_alpha, max_cands );
	for j in [1..#new_cands] do
	    Append( ~candidates, new_cands[j] );
	end for;
	max_cands -:= #new_cands;
	if max_cands eq 0 then
	    break;
	end if;
    end for;

    return candidates;

end function;