/*
* 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 := 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 := 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 := 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, );
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;