SAGE Worksheet: GNFS Term Project

# 
# This is a Toy implementation of the General Number Field Sieve For Factoring
# 
# Author:
# 
# Math 581F Fall 2007
# Professor William Stein
# Final Project 
#
# "You haven't lived till you've sieved."
#

#  The following are the sources I used to learn about the General Number Field Sieve

#  Matthew E Briggs Master's Thesis from Virginia Polytech (This is probably the most clear and detailed)
#  The Lenstra et al book "The development of the number field sieve"
#  Pomerance's "Tale of Two Sieves"
#  Henri Cohn's "A Course in Computational Algebraic Number Theory"
#  The Math 581F Textbook / Lecture Notes (William Stein)
#  
#  Personal conversations with Peter Montgomery

#
# The General Number Field Sieve works like this
# 
# 1) Decide to use GNFS to factor N, determine algorithm parameters based on N
#    such as, the size of the rational factor base, size of the algebraic factor base
#    number of quadratic characters, degree d, and polynomial f to use
#
# 2) calculate the rational prime factor base (primes up to bound)
# 3) calculate the algebraic factor base (first degree prime ideals up to bound)
# 4) calculate the quadratic characters (first few 1st degree prime ideals beyond the algebraic factor base.)
# 5) use sieving to determine integer pairs (a,b) that have values
#    have a+b*m (m a root of f mod N) and N(a+b*theta) (theta complex root of f)
#    smooth with respect to both the algebraic and rational prime factor bases.
#    while sieving, also determine the degree mod 2 over each prime (rational and algebraic)
#    and save this information.
#    Also, record the power (mod 2) of -1 determining the legendre symbol (a+b*s / q) of the
#    first degree prime ideal (s,q) in the quadratic character base.
# 6) Once there are enough such primes more than the sum of the lengths of the bases
#    use linear algebra to determine the kernel of the transformation.
# 7) Use the result of step 6 to calculate a product of (a+b*theta)'s =  Beta^2 in the number field
# 8) Use the CRT to determine Beta from Beta^2 in the number field
# 9) replace theta in Beta^2 with m to obtain v^2 = phi(m)
# 10) calculate u^2 = product of (a+b*m)'s mod N, calculate u from (a+b*m)
# 11) u^2 - v^2 = 0 mod N, so use difference of squares trick to try to factor N.
# 12) repeat until N is factored.
#

#
# Major Weaknesses of this project:
# The major weaknesses of this project are that I relied a lot on sage
# to do math.  Although it served its purpose in teaching me about the GNFS.
#
# WARNING!!!
# On My Laptop This Notebook can not handle in the mid to high 8 digits.
#
# Other than that:
# 1) To calculate square roots, I should use the CRT, but instead
#    I let sage functions do the work for me.
#    I did this because, the CRT is a technical detail at the end of the
#    algorithm, and also, I understand that math already.
#    I wanted to focus on learning the stuff I didn't understand.
# 2) I rely on sage for all linear algebra, this wouldn't work at limits where
#    10s or 100s of thousands of elements may be in the factor base.
# 3) This program uses is_prime_power() to try and determine if an N is a prime power
#    (because the GNFS does not work well on prime powers.)
#    This function is actually the limit on sizes of N that can be factored.
# 4) If a number is a prime power, I use the ECM to factor it, in actuality
#    an algorithm would be selected heuristically.
#

#
# Global Flags: These can be set to change the behavior of the sieving algorithm
# Used for sort of anecdotal observation and experimentation of the algorithm
#

# if this value is set to false, the algorithm will not use quadratic characters
use_quadratic_characters = True;

#
# This cell implements simple heuristics for the
# parameters used in the GNFS based on given N
#

rational_factor_base_bound_ratio = 5;    # this factor seemed to work well
algebraic_factor_base_bound_ratio = 10;  # so did this
quadratic_character_base_ratio = 0.10;   # as a percentage of the algebraic factor base

degree = 3;     # hard coded in this toy implementation

a_bound_ratio = 0.004;

def determine_GNFS_bounds(N):
    B1 = max(10, ceil(rational_factor_base_bound_ratio*log(N)));
    B2 = max(20, ceil(algebraic_factor_base_bound_ratio*log(N)));

    if (use_quadratic_characters):
        Quad_character_count = max(10, ceil((1.0+quadratic_character_base_ratio)*B2))
    else:
        Quad_character_count = 0

    B3 = B2 + Quad_character_count;

    a_bound = min(max(50, ceil(a_bound_ratio * N)), (2^16 - 1) );

    print 'Parameter Selection for GNFS:'

    print 'upper bound on rational prime base (all primes under):', B1
    print 'upper bound on algebraic prime base (all first degree primes under):', B2
    print 'upper bound on quadratic character prime base (all subsequent primes under):', B3

    print 'search bounds for unit term: [', -a_bound, ',', a_bound, ']'

    return (B1, B2, B3, degree, a_bound);

#
# Note about the representation of prime base elements
# rational primes are represented as (m mod p, p) p prime in QQ
# algebraic primes are represented as (r, p) p prime in QQ, r root of f mod p.
# storing in this form helps sieving and allows common code to be used
# for both
#

# This cell defines the function that
# generates the rational factor base
# (rational primes for sieving)

def generate_rational_factor_base(m, Bound):

    rational_factor_base = [];

    for p in xrange(1,Bound):
        if(is_prime(p)):
            rational_factor_base.append((m % p, p));

    return rational_factor_base;

# This cell generate the algebraic factor base
# (first degree prime ideals of order Z[theta] in number field)

# helper functions:

# uses sage to determine roots of f mod p
def get_roots_mod_p(f, p):
    return f.change_ring(GF(p)).roots(None, False);

# generates the algebraic factor base
def generate_algebraic_factor_base(f, Bound):

    algebraic_factor_base = [];

    for p in xrange(Bound):
        if(is_prime(p)):
            roots_fmodp = get_roots_mod_p(f, p);

            for j in xrange(len(roots_fmodp)):
                algebraic_factor_base.append((roots_fmodp[j].lift(),p))
                
    return algebraic_factor_base;
# generates quadratic characters
def generate_quadratic_character_base(f, Lower_Bound, Upper_Bound):

    quadratic_character_base = [];

    for p in xrange(Lower_Bound, Upper_Bound):
        if(is_prime(p)):
            roots_fmodp = get_roots_mod_p(f, p);

            for j in xrange(len(roots_fmodp)):
                quadratic_character_base.append((roots_fmodp[j].lift(),p))
                
    return quadratic_character_base;

#
# this cell defines helper functions for polynomials
# 
#

#
# determines m-arry expansion of N
# and the coefficients of powers of m
def get_polynomial_coefficients(N, m):

    coefficient_list = [];

    temp = N;

    while(0 != temp):
        r = temp % m;
        temp = (temp - r) / m;
        coefficient_list.append(r);    

    return coefficient_list;
    

#
# function that tries to generate a
# monic, irreducible polynomial f with root m mod N
# with degree d.
# 
def generate_polynomial(polynomial_ring, N, d):
    # determine and lower and upper bounds for m

    q = 0;

    while (1 != q):

        u = floor(pow(N, 1.0/d));
        l = ceil(pow(N, (2.0*d - 3.0)/(2.0*d*(d-1))));
    
        m = floor(round(uniform(l, u)));

        (q,r) = (N).quo_rem(m^d);
    

    f = polynomial_ring(get_polynomial_coefficients(N, m));
   
    return (f, m)

#
# checks if a polynomial is monic and irreducible
#
def is_good_polynomial(f):
    return (f.is_monic() and f.is_irreducible())

#
# retries a bunch of times to generate a suitable random polynomial
#
def get_polynomial(polynomial_ring, N, d, retry_count):

    found_polynomial = False;

    for j in xrange(retry_count):
        (f, m) = generate_polynomial(polynomial_ring, N, d);
        if is_good_polynomial(f):
            found_polynomial = True;
            break;

    if (found_polynomial):
        return (f, m)
    else:
        return (0,0)

#
# This cell contains the sieving code
#
quo_rem = Integer.quo_rem
#
# Common code that will take a factor base and sieve it
#
def sieve_interval(interval_to_sieve, a_lowerbound, b, factor_base):
factor_base_count = len(factor_base); exponent_lists = [ [0 for j in
xrange(factor_base_count+1)] for k in xrange(len(interval_to_sieve))] #
# determine the sign exponents of the elements in each interval # for j
in xrange(len(interval_to_sieve)): if (interval_to_sieve[j] < 0):
exponent_lists[j][factor_base_count] = 1; interval_to_sieve[j] =
abs(interval_to_sieve[j]); # # do the factor base sieving # for j in
xrange(factor_base_count): p = factor_base[j][1]; bt =
-b*factor_base[j][0] % p; # # Much of the power of the sieving
algorithms # come from the regular access pattern and # advancing
through the array p elements at a time # # in this case we have our
factor bases stored as a pair (t,p) # such that the value in the sieve
in that position # (for the rational base this is a+b*m and for the
algebraic base it is N(a+b*theta) # is only divisible by p if a = -b*t
+ k*p, as such we # use the lower bound of the a range to determine the
first # position in the array to start sieving at advance p positions
from there for k in xrange((bt - a_lowerbound) % p,
len(interval_to_sieve), p): # # This loop goes through and calculates
the largest integer e # such that p^e divides the array entry # e = 0;
r = 0; while (0 == r): if (interval_to_sieve[k] != 0): (q,r) =
quo_rem(interval_to_sieve[k] , p); if(0 == r): r = 0; e = e + 1;
interval_to_sieve[k] = q; else: r = 1; else: r = 1;
exponent_lists[k][j] = e; return (interval_to_sieve, exponent_lists)
#
# determines the quadratic character exponent vector
#
def determine_quad_char_expons(a,b, quad_char_base): qc_expons = [0 for
j in xrange(len(quad_char_base))]; for j in
xrange(len(quad_char_base)): q = quad_char_base[j][1]; s =
quad_char_base[j][0]; qc_expons[j] = Integer((1 -
legendre_symbol(a+s*b, q))/2); return qc_expons;
#
# function: search_for_smooth_pairs
# input: 6 parameters:
# lower bound of a's
# b
# list of a+bm, sieved with respect to rational factor base
# The exponent vector with respect to the rational factor base
# list of N(a+b*theta) sieved with respect to the algebraic factor base
# The exponent vector with respect to the algebraic factor base
# output: a list of triplets of:
# (a,b) pairs of integers that are smooth with respect to
# both the algebraic and rational factor bases.
# the exponent list with respect to the rational factor base
# the exponent list with respect to the algebraic factor base
#
def search_for_smooth_pairs(a_lowerbound, b, rat_sieve, rat_exp_list,
alg_sieve, alg_exp_list, quad_char_base): smooth_list = []; for j in
xrange(len(rat_sieve)): if ( (1 == rat_sieve[j]) and (1 ==
alg_sieve[j]) ): qc_expons = determine_quad_char_expons(a_lowerbound+j,
b, quad_char_base);
smooth_list.append(((a_lowerbound+j,b),rat_exp_list[j],alg_exp_list[j],
qc_expons)); return smooth_list;
#
# Do the actual sieving
# returns a list of (a,b) pairs along with the rational, algebraic, and
quadratic exponents
#
def perform_sieving(f, m, a_bound, rat_factor_base, alg_factor_base,
quad_char_base): # determine the dimension of the space we are working
over. l1 = len(rat_factor_base); l2 = len(alg_factor_base); l3 =
len(quad_char_base); l = l1 + l2 + l3 + 2; # add two for rational and
algebraic sign exponents; d = f.degree() smooth_element_list = []; #
loop until we have more smooth elements # than we have dimension b = 1;
while(len(smooth_element_list) <= l): rat_interval = [];
alg_interval = []; for j in xrange(2*a_bound + 1): a = -a_bound + j;
rat_interval.append(a + b*m); # sieve the value a + b*m
alg_interval.append(Integer(((-b)^d)*f(-a/b))); # sieve the value
N(a+b*theta) = (-b)^d * f(-a/b) (rat_interval, rat_exp_list) =
sieve_interval(rat_interval, -a_bound, b, rat_factor_base);
(alg_interval, alg_exp_list) = sieve_interval(alg_interval, -a_bound,
b, alg_factor_base); smooth_pairs_found =
search_for_smooth_pairs(-a_bound, b, rat_interval, rat_exp_list,
alg_interval, alg_exp_list, quad_char_base);
smooth_element_list.extend(smooth_pairs_found); print 'found',
len(smooth_element_list), 'simultaneously smooth pairs.' b = b+1;
return smooth_element_list;

#
#    This cell contains a function that computes the product of the integer pairs
#    as an integer polynomial, and reduces it module the ideal generated by f
#

#
# computes product of linear polynomials mod Ideal generated by f
# (equivalent to doing math in number field
#
def compute_numberfield_product_of_pairs(polynomial_ring, f, integer_pairs, vector):

    I = polynomial_ring.ideal(f);

    product_polynomial = polynomial_ring([1]);

    for j in xrange(len(vector)):

        if (1 == vector[j]):

            linear_poly = polynomial_ring([integer_pairs[j][0], integer_pairs[j][1]]);
            product_polynomial = product_polynomial * linear_poly;
            product_polynomial.mod(I);

    return product_polynomial;

#
# computes the product of (a + b*m)'s mod N
# of the (a + b*m)'s corresponding to 1s in vector
#
def compute_integer_product_of_pairs(polynomial_ring, f, m, N, integer_pairs, vector):

    prod  = 1;

    for j in xrange(len(vector)):
        if (1 == vector[j]):
            prod  = prod*(integer_pairs[j][0] + m*integer_pairs[j][1]) % N;

    return prod;
    
#
# computes the difference of squares u^2 - v^2
# and returns (u+v) and (u-v)
#
def compute_difference_of_squares(polynomial_ring, f, m, N, integer_pairs, vector):

    found_squares = False;

    u_plus_v = 0; 
    u_minus_v = 0;

    ints_mod_N = ZZ.quo(ZZ.ideal(N));

    vsquared = ints_mod_N(compute_integer_product_of_pairs(polynomial_ring, f, m, N, integer_pairs, vector));

    if ( is_square(vsquared)):

        beta_squared = compute_numberfield_product_of_pairs(polynomial_ring, f, integer_pairs, vector);
        beta = find_square(f, beta_squared, polynomial_ring);
    
        if (False != beta):
        
            u = (beta(m)) % N;
            v = vsquared.sqrt().lift();
            u_plus_v = (u + v) % N; 
            u_minus_v = (u - v) % N;

            found_squares = True;

        else:
            print 'failed to find a square root in number field.'
    else:
        print 'integer was not square'

    return (found_squares, u_plus_v, u_minus_v);

#
#    Cheater Helper Functions
#
#    These are functions that I was too lame to implement myself
#    And Hence relied on sage's built in functionality.
#

ecm_instance = ECM();

def factor_prime_power(N):
    # The GNFS doesn't work one prime powers
    # so we use ecm to factor,
    # in reality, 
    # we would apply a heuristic on the best way to determine the prime and exponent    
    
    ecm_factors = ecm_instance.factor(N);

    prime_power_factorization = [(ecm_factors[0], len(ecm_factors))]

    return Factorization(prime_power_factorization);

#
#    This simple function does all the linear algebra (thank you sage!)
#
def find_kernel_vectors(M):
    return M.kernel().basis_matrix().rows();

#
# In reality this would be done with the CRT.
#
def find_square(f, beta_squared, polynomial_ring):
    # computes the square root of beta_squared
    # (a polynomial over ZZ of degree less than d, that represents the element in the order ZZ[theta])
    # and returns the polynomial representing the square root.
    
    K.<a> = NumberField(f);
    beta2 = beta_squared(a);
    if(not beta2.is_square()):
        print 'The supposedly found square is not square!'
        return False;
    return polynomial_ring(beta2.sqrt().list());

#
# This is the factoring algorithm
# It factors N, and returns a sage factorization object of N
#

def GNFS_Factoring_Algorithm(N):

    # determine if this is a good number to use the GNFS on
    # if it is not (i.e. it is a prime power)
    # call the specialized prime power factorization function

    if(is_prime_power(N)):

        print N, 'is a prime power, reducing to using ECM to factor'

        return factor_prime_power(N);

    # determine the algorithm bounds
    (B1, B2, B3, d, a_bound) = determine_GNFS_bounds(N);

    R.<x> = ZZ[];

    (f, m) = get_polynomial(R, N, d, 20);

    if((0 == f) and (0 == m)):
    
        # if the polynomial selection code did not find a degree d
        # polynomial in 20 tries, just use the built-in factoring
        # algorithm, because N is too small for the GNFS to work
        # at that point

        return factor(N);

    # get the rational factor base
    rational_factor_base = generate_rational_factor_base(m, B1);

    # get the algebraic factor base
    algebraic_factor_base = generate_algebraic_factor_base(f, B2);

    # get the quadratic character base
    quadratic_character_base = generate_quadratic_character_base(f, B2, B3);

    #
    # get a list of smooth elements
    #

    smooth_element_list = perform_sieving(f, m, a_bound, rational_factor_base, algebraic_factor_base, quadratic_character_base);

    #
    #    parse out the list of smooth elements into a list of integer pairs,
    #

    integer_pairs = [smooth_element_list[k][0] for k in xrange(len(smooth_element_list))];
    exponent_vectors = [];

    #
    # use the smooth element list vectors to make a matrix
    #

    for k in xrange(len(smooth_element_list)):
        exponent_vectors.extend(smooth_element_list[k][1]);
        exponent_vectors.extend(smooth_element_list[k][2]);
        exponent_vectors.extend(smooth_element_list[k][3]);

    num_expons = len(smooth_element_list[0][1]) + len(smooth_element_list[0][2]) + len(smooth_element_list[0][3]);

    M = matrix(GF(2), len(smooth_element_list), exponent_vectors )
    
    solutions = find_kernel_vectors(M);

    nontrivial_first_factor = False;
    nontrivial_second_factor = False;

    k = 0

    #
    #    Loop over the vectors that sum to zero to try and find one
    #    That leads to a non trivial factorization
    #
    while ((k < len(solutions)) and (not nontrivial_first_factor) and (not nontrivial_second_factor)):
        
        (found_squares, u_plus_v, u_minus_v) = compute_difference_of_squares(R, f, m, N, integer_pairs, solutions[k]);

        if(found_squares):
                
            first_factor_of_N = gcd(u_plus_v, N);

            second_factor_of_N = gcd(u_minus_v, Integer(N/first_factor_of_N));

            if((1 != first_factor_of_N) and (N != first_factor_of_N)):
                nontrivial_first_factor = True;
            if((1 != second_factor_of_N) and (N != second_factor_of_N)):
                nontrivial_second_factor = True;
            else:
                print 'failed to find a nontrivial factor.'
    
        else:
            print 'failed to find squares.'

        k = k + 1;


    factors_of_N = [];

    N_remaining = N;

    if(nontrivial_first_factor):
        factors_of_N.append(first_factor_of_N);
        N_remaining = Integer(N_remaining/first_factor_of_N)

    if(nontrivial_second_factor):
        remaining_second_factor = gcd(second_factor_of_N, N_remaining);
        factors_of_N.append(remaining_second_factor);
        N_remaining = Integer(N_remaining/second_factor_of_N);

    if(1 != N_remaining):
        factors_of_N.append(N_remaining);

    factorization_of_N = factor(1);

    #
    #    Recursively factor non trivial factors
    #    Or N, if we completely failed to find a factor
    #

    for j in xrange(len(factors_of_N)):
        factorization_of_N = factorization_of_N * GNFS_Factoring_Algorithm(factors_of_N[j]);

    return factorization_of_N;

#
# Some test cases follow:
#

time GNFS_Factoring_Algorithm(45113);

Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
54
upper bound on algebraic prime base (all first degree primes under):
108
upper bound on quadratic character prime base (all subsequent primes
under):                                       227
search bounds for unit term: [ -181 , 181 ]
found 25 simultaneously smooth pairs.
found 33 simultaneously smooth pairs.
found 60 simultaneously smooth pairs.
found 65 simultaneously smooth pairs.
229 is a prime power, reducing to using ECM to factor
197 is a prime power, reducing to using ECM to factor
197 * 229
CPU time: 0.60 s,  Wall time: 1.24 s

Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):                                       54
upper bound on algebraic prime base (all first degree primes under):                                       108
upper bound on quadratic character prime base (all subsequent primes under):                                       227
search bounds for unit term: [ -181 , 181 ]
found 25 simultaneously smooth pairs.
found 33 simultaneously smooth pairs.
found 60 simultaneously smooth pairs.
found 65 simultaneously smooth pairs.
229 is a prime power, reducing to using ECM to factor
197 is a prime power, reducing to using ECM to factor
197 * 229
CPU time: 0.60 s,  Wall time: 1.24 s

197 * 229

GNFS_Factoring_Algorithm(99999);

Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
58
upper bound on algebraic prime base (all first degree primes under):
116
upper bound on quadratic character prime base (all subsequent primes
under):                                       244
search bounds for unit term: [ -400 , 400 ]
found 34 simultaneously smooth pairs.
found 45 simultaneously smooth pairs.
found 82 simultaneously smooth pairs.
failed to find a nontrivial factor.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
30
upper bound on algebraic prime base (all first degree primes under):
60
upper bound on quadratic character prime base (all subsequent primes
under):                                       126
search bounds for unit term: [ -50 , 50 ]
found 18 simultaneously smooth pairs.
found 27 simultaneously smooth pairs.
found 45 simultaneously smooth pairs.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
25
upper bound on algebraic prime base (all first degree primes under):
49
upper bound on quadratic character prime base (all subsequent primes
under):                                       103
search bounds for unit term: [ -50 , 50 ]
found 17 simultaneously smooth pairs.
found 24 simultaneously smooth pairs.
found 44 simultaneously smooth pairs.
failed to find a nontrivial factor.
failed to find a nontrivial factor.
failed to find a nontrivial factor.
3 is a prime power, reducing to using ECM to factor
41 is a prime power, reducing to using ECM to factor
3 is a prime power, reducing to using ECM to factor
271 is a prime power, reducing to using ECM to factor
3^2 * 41 * 271

Parameter Selection for GNFS:
upper bound on rational prime base (all primes under): 58
upper bound on algebraic prime base (all first degree primes under): 116
upper bound on quadratic character prime base (all subsequent primes under): 244
search bounds for unit term: [ -400 , 400 ]
found 34 simultaneously smooth pairs.
found 45 simultaneously smooth pairs.
found 82 simultaneously smooth pairs.
failed to find a nontrivial factor.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under): 30
upper bound on algebraic prime base (all first degree primes under): 60
upper bound on quadratic character prime base (all subsequent primes under): 126
search bounds for unit term: [ -50 , 50 ]
found 18 simultaneously smooth pairs.
found 27 simultaneously smooth pairs.
found 45 simultaneously smooth pairs.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under): 25
upper bound on algebraic prime base (all first degree primes under): 49
upper bound on quadratic character prime base (all subsequent primes under): 103
search bounds for unit term: [ -50 , 50 ]
found 17 simultaneously smooth pairs.
found 24 simultaneously smooth pairs.
found 44 simultaneously smooth pairs.
failed to find a nontrivial factor.
failed to find a nontrivial factor.
failed to find a nontrivial factor.
3 is a prime power, reducing to using ECM to factor
41 is a prime power, reducing to using ECM to factor
3 is a prime power, reducing to using ECM to factor
271 is a prime power, reducing to using ECM to factor
3^2 * 41 * 271

GNFS_Factoring_Algorithm(784352)

Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
68
upper bound on algebraic prime base (all first degree primes under):
136
upper bound on quadratic character prime base (all subsequent primes
under):                                       286
search bounds for unit term: [ -3138 , 3138 ]
found 29 simultaneously smooth pairs.
found 68 simultaneously smooth pairs.
found 83 simultaneously smooth pairs.
found 123 simultaneously smooth pairs.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
44
upper bound on algebraic prime base (all first degree primes under):
88
upper bound on quadratic character prime base (all subsequent primes
under):                                       185
search bounds for unit term: [ -50 , 50 ]
found 20 simultaneously smooth pairs.
found 47 simultaneously smooth pairs.
found 58 simultaneously smooth pairs.
16 is a prime power, reducing to using ECM to factor
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
30
upper bound on algebraic prime base (all first degree primes under):
60
upper bound on quadratic character prime base (all subsequent primes
under):                                       126
search bounds for unit term: [ -50 , 50 ]
found 19 simultaneously smooth pairs.
found 45 simultaneously smooth pairs.
failed to find a nontrivial factor.
2 is a prime power, reducing to using ECM to factor
193 is a prime power, reducing to using ECM to factor
127 is a prime power, reducing to using ECM to factor
2^5 * 127 * 193

Parameter Selection for GNFS:
upper bound on rational prime base (all primes under): 68
upper bound on algebraic prime base (all first degree primes under): 136
upper bound on quadratic character prime base (all subsequent primes under): 286
search bounds for unit term: [ -3138 , 3138 ]
found 29 simultaneously smooth pairs.
found 68 simultaneously smooth pairs.
found 83 simultaneously smooth pairs.
found 123 simultaneously smooth pairs.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under): 44
upper bound on algebraic prime base (all first degree primes under): 88
upper bound on quadratic character prime base (all subsequent primes under): 185
search bounds for unit term: [ -50 , 50 ]
found 20 simultaneously smooth pairs.
found 47 simultaneously smooth pairs.
found 58 simultaneously smooth pairs.
16 is a prime power, reducing to using ECM to factor
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under): 30
upper bound on algebraic prime base (all first degree primes under): 60
upper bound on quadratic character prime base (all subsequent primes under): 126
search bounds for unit term: [ -50 , 50 ]
found 19 simultaneously smooth pairs.
found 45 simultaneously smooth pairs.
failed to find a nontrivial factor.
2 is a prime power, reducing to using ECM to factor
193 is a prime power, reducing to using ECM to factor
127 is a prime power, reducing to using ECM to factor
2^5 * 127 * 193

#
#    Factoring test code
#    This function will randomly generate an integer and factor it with the GNFS
#
#    It fails sometimes, due to stuff like low memory conditions, and bugs
#

def test_GNFS():

    N = Integer(ceil(uniform(1, 1000000)));

    print 'Factoring N =', N

    GNFS_Factorization = GNFS_Factoring_Algorithm(N);    
    
    if (GNFS_Factorization.value() != N):
        print 'Test Case N =', N, ' FAILED!!!'
        print 'says factorization is:', GNFS_Factorization
        return False

    print GNFS_Factorization

    print 'Test Case PASSED!'

    return True

test_GNFS()

Factoring N = 606063
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
67
upper bound on algebraic prime base (all first degree primes under):
134
upper bound on quadratic character prime base (all subsequent primes
under):                                       282
search bounds for unit term: [ -2425 , 2425 ]
found 31 simultaneously smooth pairs.
found 43 simultaneously smooth pairs.
found 85 simultaneously smooth pairs.
3 is a prime power, reducing to using ECM to factor
202021 is a prime power, reducing to using ECM to factor
3 * 202021
Test Case PASSED!
True

Factoring N = 606063
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):                                       67
upper bound on algebraic prime base (all first degree primes under):                                       134
upper bound on quadratic character prime base (all subsequent primes under):                                       282
search bounds for unit term: [ -2425 , 2425 ]
found 31 simultaneously smooth pairs.
found 43 simultaneously smooth pairs.
found 85 simultaneously smooth pairs.
3 is a prime power, reducing to using ECM to factor
202021 is a prime power, reducing to using ECM to factor
3 * 202021
Test Case PASSED!
True

#
#    At first, when reading up on the General Number Field Sieve algorithm
#    I had a hard time understanding the Quadratic character portion of
#    The algorithm, as such, I set this flag to test just what happens if
#    it is not used
#    observe that when the flag is set to false, the algorithm takes longer
#    and has many more errors about failing to find a square in the number field
#
use_quadratic_characters = False

GNFS_Factoring_Algorithm(784352)

Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
68
upper bound on algebraic prime base (all first degree primes under):
136
upper bound on quadratic character prime base (all subsequent primes
under):                                       136
search bounds for unit term: [ -3138 , 3138 ]
found 48 simultaneously smooth pairs.
found 101 simultaneously smooth pairs.
The supposedly found square is not square!
failed to find a square root in number field.
failed to find squares.
failed to find a nontrivial factor.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
42
upper bound on algebraic prime base (all first degree primes under):
84
upper bound on quadratic character prime base (all subsequent primes
under):                                       84
search bounds for unit term: [ -50 , 50 ]
found 40 simultaneously smooth pairs.
The supposedly found square is not square!
failed to find a square root in number field.
failed to find squares.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
35
upper bound on algebraic prime base (all first degree primes under):
70
upper bound on quadratic character prime base (all subsequent primes
under):                                       70
search bounds for unit term: [ -50 , 50 ]
found 24 simultaneously smooth pairs.
found 53 simultaneously smooth pairs.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under):              
32
upper bound on algebraic prime base (all first degree primes under):
63
upper bound on quadratic character prime base (all subsequent primes
under):                                       63
search bounds for unit term: [ -50 , 50 ]
found 23 simultaneously smooth pairs.
found 50 simultaneously smooth pairs.
The supposedly found square is not square!
failed to find a square root in number field.
failed to find squares.
4 is a prime power, reducing to using ECM to factor
127 is a prime power, reducing to using ECM to factor
2 is a prime power, reducing to using ECM to factor
4 is a prime power, reducing to using ECM to factor
193 is a prime power, reducing to using ECM to factor
2^5 * 127 * 193

Parameter Selection for GNFS:
upper bound on rational prime base (all primes under): 68
upper bound on algebraic prime base (all first degree primes under): 136
upper bound on quadratic character prime base (all subsequent primes under): 136
search bounds for unit term: [ -3138 , 3138 ]
found 48 simultaneously smooth pairs.
found 101 simultaneously smooth pairs.
The supposedly found square is not square!
failed to find a square root in number field.
failed to find squares.
failed to find a nontrivial factor.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under): 42
upper bound on algebraic prime base (all first degree primes under): 84
upper bound on quadratic character prime base (all subsequent primes under): 84
search bounds for unit term: [ -50 , 50 ]
found 40 simultaneously smooth pairs.
The supposedly found square is not square!
failed to find a square root in number field.
failed to find squares.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under): 35
upper bound on algebraic prime base (all first degree primes under): 70
upper bound on quadratic character prime base (all subsequent primes under): 70
search bounds for unit term: [ -50 , 50 ]
found 24 simultaneously smooth pairs.
found 53 simultaneously smooth pairs.
Parameter Selection for GNFS:
upper bound on rational prime base (all primes under): 32
upper bound on algebraic prime base (all first degree primes under): 63
upper bound on quadratic character prime base (all subsequent primes under): 63
search bounds for unit term: [ -50 , 50 ]
found 23 simultaneously smooth pairs.
found 50 simultaneously smooth pairs.
The supposedly found square is not square!
failed to find a square root in number field.
failed to find squares.
4 is a prime power, reducing to using ECM to factor
127 is a prime power, reducing to using ECM to factor
2 is a prime power, reducing to using ECM to factor
4 is a prime power, reducing to using ECM to factor
193 is a prime power, reducing to using ECM to factor
2^5 * 127 * 193

use_quadratic_characters = True