## "Fossies" - the Fresh Open Source Software Archive

### Member "NZMATH-1.2.0/nzmath/poly/factor.py" (19 Nov 2012, 9241 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "factor.py" see the Fossies "Dox" file reference documentation.

```    1 from __future__ import division
2 import nzmath.arith1 as arith1
3 import nzmath.prime as prime
4 import nzmath.rational as rational
5 import nzmath.combinatorial as combinatorial
6 import nzmath.intresidue as intresidue
7 import nzmath.finitefield as finitefield
8 import nzmath.poly.uniutil as uniutil
9 import nzmath.poly.ring as poly_ring
10 import nzmath.poly.hensel as hensel
11
12 def zassenhaus(f):
13     """
14     zassenhaus(f) -> list of factors of f.
15
16     Factor a squarefree monic integer coefficient polynomial f with
17     Berlekamp-Zassenhaus method.
18     """
21
24     if len(fp_factors) == 1:
25         return [f]
26
27     # purge leading coefficient from factors
28     for i,g in enumerate(fp_factors):
29         if g.degree() == 0:
30            del fp_factors[i]
31            break
32
33     # lift to Mignotte bound
34     blm = upper_bound_of_coefficient(f)
35     bound = p**(arith1.log(2*blm,p)+1)
36
37     # Hensel lifting
38     lf_inv_modq = intresidue.IntegerResidueClass(lf, bound).inverse()
39     fq = f.coefficients_map(lambda c: (lf_inv_modq*c).minimumAbsolute()) # fq is monic
40     fq_factors, q = hensel.lift_upto(fq, fp_factors, p, bound)
41
42     return brute_force_search(f, fq_factors, bound)
43
45     """
47
48     Return a prime p and a p-adic factorization of given integer
49     coefficient squarefree polynomial f. The result factors have
50     integer coefficients, injected from F_p to its minimum absolute
51     representation. The prime is chosen to be 1) f is still squarefree
52     mod p, 2) the number of factors is not greater than with the
53     successive prime.
54     """
55     num_factors = f.degree()
56     stock = None
57     for p in prime.generator():
58         fmodp = uniutil.polynomial(
59             f.terms(),
60             finitefield.FinitePrimeField.getInstance(p))
61         if f.degree() > fmodp.degree():
62             continue
63         g = fmodp.getRing().gcd(fmodp,
64                                 fmodp.differentiate())
65         if g.degree() == 0:
66             fp_factors = fmodp.factor()
67             if not stock or num_factors > len(fp_factors):
68                 stock = (p, fp_factors)
69                 if len(fp_factors) == 1:
70                     return stock
71                 num_factors = len(fp_factors)
72             else:
73                 break
74     p = stock[0]
75     fp_factors = []
76     for (fp_factor, m) in stock[1]:
77         assert m == 1 # since squarefree
78         fp_factors.append(minimum_absolute_injection(fp_factor))
79     return (p, fp_factors)
80
81 def brute_force_search(f, fp_factors, q):
82     """
83     brute_force_search(f, fp_factors, q) -> [factors]
84
85     Find the factorization of f by searching a factor which is a
86     product of some combination in fp_factors.  The combination is
87     searched by brute force.
88     """
89     factors = []
90     d, r = 1, len(fp_factors)
91     while 2*d <= r:
92         found, combination = find_combination(f, d, fp_factors, q)
93         if found:
94             factors.append(found)
95             for picked in combination:
96                 fp_factors.remove(picked)
97             f = f.pseudo_floordiv(found).primitive_part()
98             r -= d
99         else:
100             d += 1
101     factors.append(f.primitive_part())
102     return factors
103
104 def padic_lift_list(f, factors, p, q):
105     """
106     padicLift(f, factors, p, q) -> lifted_factors
107
108     Find a lifted integer coefficient polynomials such that:
109       f = G1*G2*...*Gm (mod q*p),
110       Gi = gi (mod q),
111     from f and gi's of integer coefficient polynomials such that:
112       f = g1*g2*...*gm (mod q),
113       gi's are pairwise coprime
114     with positive integers p dividing q.
115     """
116     ZpZx = poly_ring.PolynomialRing(
117         intresidue.IntegerResidueClassRing.getInstance(p))
118     gg = arith1.product(factors)
119     h = ZpZx.createElement([(d, c // q) for (d, c) in (f - gg).iterterms()])
120     lifted = []
121     for g in factors:
122         gg = gg.pseudo_floordiv(g)
123         g_mod = ZpZx.createElement(g)
124         if gg.degree() == 0:
125             break
126         u, v, w = extgcdp(g, gg, p)
127         if w.degree() > 0:
128             raise ValueError("factors must be pairwise coprime.")
129         v_mod = ZpZx.createElement(v)
130         t = v_mod * h // g_mod
131         lifted.append(g + minimum_absolute_injection(v_mod * h - g_mod * t)*q)
132         u_mod = ZpZx.createElement(u)
133         gg_mod = ZpZx.createElement(gg)
134         h = u_mod * h + gg_mod * t
135         lifted.append(g + minimum_absolute_injection(h)*q)
136     return lifted
137
138 def extgcdp(f, g, p):
139     """
140     extgcdp(f,g,p) -> u,v,w
141
142     Find u,v,w such that f*u + g*v = w = gcd(f,g) mod p.
143     """
144     zpz = intresidue.IntegerResidueClassRing.getInstance(p)
145     f_zpz = uniutil.polynomial(f, zpz)
146     g_zpz = uniutil.polynomial(g, zpz)
147     zero, one = f_zpz.getRing().zero, f_zpz.getRing().one
148     u, v, w, x, y, z = one, zero, f_zpz, zero, one, g_zpz
149     while z:
150         q = w // z
151         u, v, w, x, y, z = x, y, z, u - q*x, v - q*y, w - q*z
152     if w.degree() == 0 and w[0] != zpz.one:
153         u = u.scalar_exact_division(w[0]) # u / w
154         v = v.scalar_exact_division(w[0]) # v / w
155         w = w.getRing().one # w / w
156     u = minimum_absolute_injection(u)
157     v = minimum_absolute_injection(v)
158     w = minimum_absolute_injection(w)
159     return u, v, w
160
161 def minimum_absolute_injection(f):
162     """
163     minimum_absolute_injection(f) -> F
164
165     Return an integer coefficient polynomial F by injection of a Z/pZ
166     coefficient polynomial f with sending each coefficient to minimum
167     absolute representatives.
168     """
169     coefficientRing = f.getCoefficientRing()
170     if isinstance(coefficientRing, intresidue.IntegerResidueClassRing):
171         p = coefficientRing.m
172     elif isinstance(coefficientRing, finitefield.FinitePrimeField):
173         p = coefficientRing.getCharacteristic()
174     else:
175         raise TypeError("unknown ring (%s)" % repr(coefficientRing))
176     half = p >> 1
177     g = {}
178     for i, c in f.iterterms():
179         if c.n > half:
180             g[i] = c.n - p
181         else:
182             g[i] = c.n
183     return uniutil.polynomial(g, rational.theIntegerRing)
184
185
186 def upper_bound_of_coefficient(f):
187     """
188     upper_bound_of_coefficient(polynomial) -> int
189
190     Compute Landau-Mignotte bound of coefficients of factors, whose
191     degree is no greater than half of the given polynomial.  The given
192     polynomial must have integer coefficients.
193     """
194     weight = 0
195     for c in f.itercoefficients():
196         weight += abs(c)**2
197     weight = arith1.floorsqrt(weight) + 1
198     degree = f.degree()
199     lc = f[degree]
200     m = (degree >> 1) + 1
201     bound = 1
202     for i in range(1, m):
203         b = combinatorial.binomial(m - 1, i) * weight + \
204             combinatorial.binomial(m - 1, i - 1) * lc
205         if bound < b:
206             bound = b
207     return bound
208
209 def find_combination(f, d, factors, q):
210     """
211     find_combination(f, d, factors, q) -> g, list
212
213     Find a combination of d factors which divides f (or its
214     complement).  The returned values are: the product g of the
215     combination and a list consisting of the combination itself.
216     If there is no combination, return (0,[]).
217     """
219     ZqZX = poly_ring.PolynomialRing(
220         intresidue.IntegerResidueClassRing.getInstance(q))
221
222     if d == 1:
223         for g in factors:
224             product = minimum_absolute_injection(ZqZX.createElement(lf*g))
225             if divisibility_test(lf*f, product):
226                 return (product.primitive_part(), [g])
227     else:
228         for idx in combinatorial.combinationIndexGenerator(len(factors), d):
229             picked = [factors[i] for i in idx]
230             product = lf * arith1.product(picked)
231             product = minimum_absolute_injection(ZqZX.createElement(product))
232             if divisibility_test(lf*f, product):
233                 return (product.primitive_part(), picked)
234     return 0, [] # nothing found
235
236 def divisibility_test(f, g):
237     """
238     Return boolean value whether f is divisible by g or not, for polynomials.
239     """
240     if g[0] and f[0] % g[0]:
241         return False
242     if isinstance(f, uniutil.FieldPolynomial) and f % g:
243         return False
244     elif isinstance(f, uniutil.UniqueFactorizationDomainPolynomial) and f.pseudo_mod(g):
245         return False
246     return True
247
248 def integerpolynomialfactorization(f):
249     """
250     integerpolynomialfactorization -> list of (factors,index) of f.
251
252     Factor a integer coefficient polynomial f with
253     Berlekamp-Zassenhaus method.
254     """
255     cont = f.content()
256     prim = f.primitive_part()
257     F = [prim]
258     G = prim
259     c = 0
260     one = G.getRing().one
261     while (G.differentiate() and F[c] != one):
262         deriv = G.differentiate()
263         F.append(F[c].subresultant_gcd(deriv))
264         c = c + 1
265         G = F[c]
266     sqfree_part = F[0].pseudo_floordiv(F[0].subresultant_gcd(F[1])).primitive_part()
267     N = zassenhaus(sqfree_part)
268
269     if cont != 1:
270         result = [(one.scalar_mul(cont) ,1)]
271     else:
272         result = []
273
274     F.reverse()
275     e = len(F)
276     for factor in N:
277         for deg, deriv in enumerate(F):
278             if not (deriv.pseudo_mod(factor)):
279                 result.append((factor, (e-deg)))
280                 break
281     return result
```