1. No category

Polynomial Factorization with Sympy
by
Alkis Akritas
ECE/UTh/Volos
Information about the various functions used in these notes can be found in the sites
http://docs.sympy.org/dev/_modules/sympy/polys/galoistools.html
and
http://docs.sympy.org/dev/_modules/sympy/polys/factortools.html
Python plugin for TeXmacs.
Please see the documentation in Help -> plugins -> Python
Python] from sympy import *
Python] x = S(’x’)
Python]
Load the tools and look at their contents:
Python] from sympy.polys.galoistools import *
Python] from sympy.polys.factortools import *
### gf, berlekamp, zassenhaus
### dup, hensel, mignotte,
#### dup_zz_irreducible_p
Python] dir(polys.galoistools)
[’ExactQuotientFailed’, ’SYMPY_INTS’, ’__builtins__’, ’__doc__’, ’__file__’,
’__name__’, ’__package__’, ’_ceil’, ’_factor_methods’, ’_gf_pow_pnm1d2’,
’_gf_trace_map’, ’_irred_methods’, ’_raise_mod_power’, ’_sort_factors’, ’_sqrt’,
’csolve_prime’, ’division’, ’factorint’, ’gf_LC’, ’gf_Qbasis’, ’gf_Qmatrix’, ’gf_TC’,
’gf_add’, ’gf_add_ground’, ’gf_add_mul’, ’gf_berlekamp’, ’gf_cofactors’, ’gf_compose’,
’gf_compose_mod’, ’gf_crt’, ’gf_crt1’, ’gf_crt2’, ’gf_csolve’, ’gf_ddf_shoup’,
’gf_ddf_zassenhaus’, ’gf_degree’, ’gf_diff’, ’gf_div’, ’gf_edf_shoup’,
’gf_edf_zassenhaus’, ’gf_eval’, ’gf_expand’, ’gf_exquo’, ’gf_factor’, ’gf_factor_sqf’,
’gf_frobenius_map’, ’gf_frobenius_monomial_base’, ’gf_from_dict’, ’gf_from_int_poly’,
’gf_gcd’, ’gf_gcdex’, ’gf_int’, ’gf_irred_p_ben_or’, ’gf_irred_p_rabin’,
’gf_irreducible’, ’gf_irreducible_p’, ’gf_lcm’, ’gf_lshift’, ’gf_monic’, ’gf_mul’,
’gf_mul_ground’, ’gf_multi_eval’, ’gf_neg’, ’gf_normal’, ’gf_pow’, ’gf_pow_mod’,
’gf_quo’, ’gf_quo_ground’, ’gf_random’, ’gf_rem’, ’gf_rshift’, ’gf_shoup’,
’gf_sqf_list’, ’gf_sqf_p’, ’gf_sqf_part’, ’gf_sqr’, ’gf_strip’, ’gf_sub’,
’gf_sub_ground’, ’gf_sub_mul’, ’gf_to_dict’, ’gf_to_int_poly’, ’gf_trace_map’,
’gf_trunc’, ’gf_value’, ’gf_zassenhaus’, ’linear_congruence’, ’print_function’,
’prod’, ’query’, ’uniform’, ’xrange’]
Python]
We will factor the following monic polynomial f (x)∈Z[x] over the integers:
Python] f = expand( (x**3 - 7*x + 7) * (x**2 - 5) )
Python] f
x**5 - 12*x**3 + 7*x**2 + 35*x - 35
Python]
1
Caveat:: If we wanted to factor the polynomial, say p(x), where
Python] p = 10*x**6 - 9*x**5 - 7*x**4 + 40*x**2 - 36*x - 28
Python] p
10*x**6 - 9*x**5 - 7*x**4 + 40*x**2 - 36*x - 28
Python]
x
we would first make it monic by replacing x by 10 , 10 being the leading coefficient of p(x), and we would
multiply p(x) times 105. That is, the monic equivalent would be the polynomial P , where P is defined below:
Python] P = 10**5 * p.subs(x, x/10)
Python] P
x**6 - 9*x**5 - 70*x**4 + 40000*x**2 - 360000*x - 2800000
Python]
At the end, we perform the inverse transformation on the 4 factors of P over the integers, in order to find
the factors of the original polynomial p(x). That is, we replace x by 10 · x in each factor and divide out their
content.
Python] factors4P = [x + 5, x - 14, x**2 - 20*x + 200, x**2 + 20*x + 200]
Python] factors4P
[x + 5, x - 14, x**2 - 20*x + 200, x**2 + 20*x + 200]
Python] trueFactors4p = []
for i in range(len(factors4P)):
con = content(factors4P[i].subs(x, 10*x)) ### content of each factor
trueFactor = quo( factors4P[i].subs(x, 10*x), con, x)
trueFactors4p.append(trueFactor)
Python] trueFactors4p
[2*x + 1, 5*x - 7, x**2 - 2*x + 2, x**2 + 2*x + 2]
Python] prod( trueFactors4p ).expand()
10*x**6 - 9*x**5 - 7*x**4 + 40*x**2 - 36*x - 28
Python]
Therefore, to make life easy, we will only factor monic polynomials! Below are the 4 steps A - D::
A. The first thing we do in the factorization algorithm is to check whether f (x) is irreducible.
Using the Schoenemann-Eisenstein irreducibility criterion, implemented by dup_zz_irreducible_p, we see
that our polynomial is reducible. The thing to note is that except for case B below, all the other functions
used in the sequel are tools “under the hood” and the polynomials are entered with their coefficients.
Python] dup_zz_irreducible_p( Poly(f).all_coeffs(), ZZ )
None
Python]
For an irreducible polynomial the output is 1.
Python] dup_zz_irreducible_p( Poly( x**2 - 5 ).all_coeffs(), ZZ )
1
Python]
2
To be 100% sure that the polynomial is reducible we might try for a few primes — not dividing the leading
coefficient of f (x) — the functions gf_irred_p_ben_or, or gf_irred_p_rabin:
Python] gf_irred_p_ben_or( Poly(f).all_coeffs(), 2, ZZ )
0
Python] gf_irred_p_rabin( Poly(f).all_coeffs(), 3, ZZ )
0
Python]
Caveat:: If there is a prime p, for which the polynomial is irreducible in Z p[x], then the polynomial f is
irreducible over the integers, i.e. in Z[x].
B. Next, we check to see whether the polynomial f is square-free.
The function sqf_list returns the list of square-free factors; in our case we see that p(x) is square-free since
it has only itself as the single factor raised to the power 1.1
Python] sqf_list( f )
(1, [(x**5 - 12*x**3 + 7*x**2 + 35*x - 35, 1)])
By contrast, the following polynomial is not square-free and, it so happens that in this case sqf_list returns
its complete factorization 2 · (x + 1)2 · (x + 2)3.
Python] sqf_list( 2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16 )
(2, [(x + 1, 2), (x + 2, 3)])
Python]
If f (x) were not square-free, then we would decompose each square-free factor.
To find an appropriate prime we start with 2 and using the function nextprime we obtain the next one —
until a suitable prime is found. See also help(nextprime).
C. Next, using the function gf_sqf_p, we pick a suitable primeto work with; that is a prime that does not
divide the leading coefficient of f (x) and leaves it square free. For example 11, 13 and others are suitable,
while 7 is not!
Python] gf_sqf_p( Poly(f).all_coeffs(), 11, ZZ )
1
Python] gf_sqf_p( Poly(f).all_coeffs(), 7, ZZ )
0
Python]
We will proceed to factor the polynomial mod 11. We can do this by working either with the original
polynomial f (x) — in which case we take its coefficients — or with the corresponding one in Z11[x] obtained
using the function gf_from_int_poly.
Python] ff = gf_from_int_poly( Poly(f).all_coeffs(), 11)
Python] ff
### a poly in GF(p)
[1, 0, 10, 7, 2, 9]
1. The first 1 that appears after the opening parenthesis is the content of p(x), i.e. the gcd of the coefficients.
3
Python]
We now factor the original polynomial f (x) in Z11[x] using Berlekamp’s algorithm gf_berlekamp, we obtain
3 factors :
Python] gf_berlekamp( Poly(f).all_coeffs(), 11, ZZ)
[[1, 4], [1, 7], [1, 0, 4, 7]]
Python]
Likewise, we can use the corresponding polynomial in Z11[x]:
Python] gf_berlekamp( ff, 11, ZZ)
[[1, 4], [1, 7], [1, 0, 4, 7]]
Python]
Factorization modulo a prime can be also performed with the high level algorithm factor. In Sympy the
Cantor-Zassenhaus algorithm is the default for average inputs. The polynomial is entered in its normal form:
Python] factor(f, modulus=11)
(x - 4)*(x + 4)*(x**3 + 4*x - 4)
Python]
The Cantor-Zassenhaus algorithm to factor polynomials modulo a prime, p, consists of two steps:
C1:: Distinct Degree Factorization (ddf), and C2:: Equal Degree Factorization (edf), both of which are
performed only in Z p[x] with the functions gf_ddf_zassenhaus and gf_edf_zassenhaus, respectively.
Using the original polynomial f (x), and the ddf algorithm we obtain distinct factors of degrees 1 and 3.
Python] ddf = gf_ddf_zassenhaus( Poly(f).all_coeffs(), 11, ZZ)
Python] ddf
[([1, 0, 6], 1), ([1, 0, 4, 7], 3)]
Python]
Likewise, using the corresponding polynomial in Z11[x]:
Python] ddf = gf_ddf_zassenhaus( ff, 11, ZZ)
Python] ddf
[([1, 0, 6], 1), ([1, 0, 4, 7], 3)]
Python]
In other words the output of ddf is as follows:
•
Degree 1: here we have a polynomial of degree 2, which implies that we have 2 linear factors.
•
Degree 3: here we have a polynomial of degree 3, which implies that we have one polynomial of
degree three and it has already been computed.
To compute the linear factors2 we use the function gf_edf_zassenhaus on the first polynomial of degree 2
obtained above, and we have:
2. This is how we find the roots of a polynomial in Zp[x]!
4
Python] gf_edf_zassenhaus( ddf[0][0], 1, 11, ZZ)
[[1, 4], [1, 7]]
Python]
Putting the above together in one function, gf_zassenhaus, we have:
Python] fmod = gf_zassenhaus( Poly(f).all_coeffs(), 11, ZZ)
Python] fmod
[[1, 4], [1, 7], [1, 0, 4, 7]]
Python]
Likewise, using the corresponding polynomial in Z11[x]:
Python] fmod = gf_zassenhaus( ff, 11, ZZ)
Python] fmod
[[1, 4], [1, 7], [1, 0, 4, 7]]
Python]
These factors are the same as the ones obtained with Berlekamp’s algorithm.
Caveat:: Before we proceed with the final step, it is advisable to try several other small primes and pick
the one that gives you the smallest number of factors. This will make things easier in the process of finding
the true factors over the integers.
D. Having factored f (x) modulo 11, our next step is to “lift” the factors fmod over the integers. To see “how
much” we have to lift them, we use Mignotte’s theorem to compute an upper bound on the maximum value
of the coefficients of any factor of f (x).
Mignotte’s theorem is implemented by the function dup_zz_mignotte_bound and we obtain as a bound the
value 2240.
Python] dup_zz_mignotte_bound(
Poly(f).all_coeffs(), ZZ)
2240
Python]
This means that if we lift the factors modulo a power 11x, and 11x > 2240, then we are sure that any factor
of f (x) modulo 11x, will also be a factor of f (x) over the integers!
The lifting of the factors is done with the function dup_zz_hensel_lift below, where one of the inputs is
the power x . In our case, x = 3 because
Python] ceiling(
log(2240) / log(11) + 0.
)
4
Python] 11**4
14641
Python]
and lifting we obtain 3 factors over the integers:
Python] factors = dup_zz_hensel_lift(11, Poly(f).all_coeffs() , fmod, 4, ZZ)
Python] factors
5
[[1, 6582], [1, -6582], [1, 0, -7, 7]]
These factors can be easily written in their “normal” form
Python]
Python] facts = []
for poly in factors:
facts.append( Poly(poly, x).as_expr() )
Python] facts
[x + 6582, x - 6582, x**3 - 7*x + 7]
Python]
However, by trial division — using the function dup_trial_division — we see that only the third polynomial is a factor of the original f (x) over the integers.
Python] dup_trial_division(Poly(f).all_coeffs(),factors, ZZ)
[([1, -6582], 0), ([1, 6582], 0), ([1, 0, -7, 7], 1)]
Python]
Having found a true factor of f (x) we divide it out and we are left with a new polynomial f (x) = x2 − 5 of
degree 2. However, the two linear factors are still not true divisors of the updatedf (x).
Python] f = quo( f, facts[2] )
Python] f
x**2 - 5
Python] dup_trial_division(Poly(f).all_coeffs(),[[1, 6582], [1, -6582]], ZZ)
[([1, -6582], 0), ([1, 6582], 0)]
Python]
At this point we have to combine these two factors; so, using the function gf_mul, we multiply out the first
two factors of factors to obtain a second degree polynomial with coefficients [1, 0, 14636]. Notice that modulo
114 = 14641 the polynomial is equivalent to [1, 0, −5], but this is revealed by the function dup_trunc:
Python] factors
[[1, 6582], [1, -6582], [1, 0, -7, 7]]
Python] factors2 = gf_mul(factors[0], factors[1], 11**4, ZZ)
Python] factors2
[1, 0, 14636]
Python]
However, at this point we actually want the polynomials with “symmetric” coefficients; that is, we want the
negative coefficients to appear. These are obtained by running this factor through the function dup_trunc.
So, we have:
Python] factors2 = dup_trunc(factors2, 11**4, ZZ)
Python] factors2
[1, 0, -5]
6
Python]
and this time we see that we have found the true factor of the updated f (x):
Python] dup_trial_division(Poly(f).all_coeffs(),[factors2], ZZ)
[([1, 0, -5], 1)]
Python]
Hence, we have completely factored the original polynomial f (x).
How does the function dup_zz_hensel_lift work?
If you are interested, the material below explains what is “under the hood” of dup_zz_hensel_lift. More
detailed presentation in class.
The function dup_zz_hensel_lift uses the Hensel function, dup_zz_hensel_step, which can lift only 2
factors at a time, from a factorization modulo m to a factorization modulo m2. To use it we have to apply
the extended Euclidean algorithm gf_gcdex for polynomials.
Let us apply dup_zz_hensel_step, two times to lift the factorization fmod of the original polynomial f (x)
to a factorization modulo 114.
Recall the original polynomial f (x) and its three factors modulo 11:
Python] f = x**5 - 12*x**3 + 7*x**2 + 35*x - 35
Python] f
x**5 - 12*x**3 + 7*x**2 + 35*x - 35
Python] fmod
[[1, 4], [1, 7], [1, 0, 4, 7]]
Python]
To use dup_zz_hensel_step we take the product of the first two factors from fmod using the function gf_mul
Python] fact1 = gf_mul( fmod[0], fmod[1], 11, ZZ)
Python] fact1
[1, 0, 6]
Python]
and the factorization of f (x) modulo 11 becomes f (x) = fact1 · fmod[2] in Z11[x].
Python] [ fact1, fmod[2] ]
[[1, 0, 6], [1, 0, 4, 7]]
Python]
Indeed, we have:
Python] gf_mul(fact1, fmod[2], 11, ZZ)
[1, 0, 10, 7, 2, 9]
Python] gf_mul(fact1, fmod[2], 11, ZZ) == ff
7
1
Python]
We now lift this factorization from a factorization modulo 11 to one modulo 112. We first compute s, t ∈ Z11[x]
such that s · fact1+t ·fmod[2]=h, where h= gcd (fact1, fmod[2]). The underscore below means that we have
no use for the gcd.
Python] s, t, _ = gf_gcdex( fact1, fmod[2], 11, ZZ )
Python] s, t
([6, 10, 10], [5, 1])
Python] Fact1, Fact2, S, T = dup_zz_hensel_step(11, Poly(f).all_coeffs(), fact1, fmod[2],
s, t, ZZ)
Python] Fact1, Fact2, S, T
([1, 0, -5], [1, 0, -7, 7], [50, 54, 21], [-50, -54])
Python]
We have now computed Fact1, Fact2, S, T such that f (x) = Fact1 · Fact2 modulo 112 and
S · Fact1+T · Fact2=h, where h= gcd (Fact1, Fact1). Indeed, we have
Python] gf_mul(Fact1, Fact2, 11**2, ZZ)
[1, 0, 109, 7, 35, 86]
Python] gf_from_int_poly( Poly(f).all_coeffs(), 11**2 )
[1, 0, 109, 7, 35, 86]
Python]
Notice that both factors, Fact1, Fact2, are already the real factors of f (x) in Z[x]. However, as suggested
by Mignotte’s theorem, we lift once more and see that the factors do not change modulo 114.
Python] Fact1, Fact2, S, T = dup_zz_hensel_step(11**2, Poly(f).all_coeffs(), Fact1,
Fact2, S, T, ZZ)
Python] Fact1, Fact2, S, T
([1, 0, -5], [1, 0, -7, 7], [7068, -4544, 505], [-7068, 4544])
Python]
This means that we have now found the true factors of f (x).
Python] dup_trial_division(Poly(f).all_coeffs(), [Fact1, Fact2], ZZ)
[([1, 0, -5], 1), ([1, 0, -7, 7], 1)]
Python] [ Poly( Fact1, x ).as_expr() ,
Poly( Fact2, x ).as_expr() ]
[x**2 - 5, x**3 - 7*x + 7]
Python]
8