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

### Member "NZMATH-1.2.0/nzmath/poly/hensel.py" (19 Nov 2012, 17023 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 "hensel.py" see the Fossies "Dox" file reference documentation.

```    1 """
2 hensel -- Hensel Lift classes for factorization
3           of integer coefficient polynomials
4
5
6 If you need lifted factors only once, you may want to use lift_upto
7 function.  On the other hand, if you might happen to need another
8 lift, you would directly use one of the lifters and its lift method
9 for consecutive lifts.
10 """
11
12 from __future__ import division
13 import sys
14 import nzmath.arith1 as arith1
15 import nzmath.rational as rational
16 import nzmath.poly.ring as polyring
17
18
19 # module globals
20 the_ring = polyring.PolynomialRing(rational.theIntegerRing)
21 the_one = the_ring.one
22 the_zero = the_ring.zero
23
24
25 def _extgcdp(f, g, p):
26     """
27     _extgcdp(f,g,p) -> u,v,w
28
29     Find u,v,w such that f*u + g*v = w = gcd(f,g) mod p.
30     p should be a prime number.
31
32     This is a private function.
33     """
34     modp = lambda c: c % p
35     u, v, w, x, y, z = the_one, the_zero, f, the_zero, the_one, g
36     while z:
38         if zlc != 1:
39             normalizer = arith1.inverse(zlc, p)
40             x = (x * normalizer).coefficients_map(modp)
41             y = (y * normalizer).coefficients_map(modp)
42             z = (z * normalizer).coefficients_map(modp)
43         q = w.pseudo_floordiv(z)
44         u, v, w, x, y, z = x, y, z, u - q*x, v - q*y, w - q*z
45         x = x.coefficients_map(modp)
46         y = y.coefficients_map(modp)
47         z = z.coefficients_map(modp)
48     if w.degree() == 0 and w[0] != 1:
49         u = u.scalar_exact_division(w[0]) # u / w
50         v = v.scalar_exact_division(w[0]) # v / w
51         w = w.getRing().one # w / w
52     return u, v, w
53
54 def _init_cofactors(factors):
55     """
56     Return list of bi's; bi = a{i+1}*...*ar where ai's are of factors.
57
58     This is a private function.
59     """
60     bis = [the_one]
61     for ai in factors[-1:0:-1]:
62         bis.append(bis[-1] * ai)
63     bis.reverse()
64     return bis
65
66
67 class HenselLiftPair(object):
68     """
69     A class represents integer polynomial pair which will be lifted by
70     Hensel's method.
71     """
73         """
74         HenselLiftPair(f, a1, a2, u1, u2, p, q)
75
76         The parameters satisfy that:
77           a1 and a2 are monic,
78           f == a1*a2 (mod q) and
79           a1*u1 + a2*u2 == 1 (mod p),
80         with positive integers p dividing q.
81
82         If p==q, q can be omit from the argument.
83         """
84         self.f = target
85         self.a1 = factor1
86         self.a2 = factor2
89         self.p = p
90         if q is None:
91             self.q = p
92         else:
93             self.q = q
94
95     @classmethod
96     def from_factors(cls, target, factor1, factor2, p):
97         """
98         Create and return an instance of HenselLiftPair.
99
100         HenselLiftPair.from_factors(f, a1, a2, p)
101
102         The parameters satisfy that:
103           f == a1*a2 (mod p)
104         with a prime number p.
105         """
106         u1, u2 = _extgcdp(factor1, factor2, p)[:2]
107         return cls(target, factor1, factor2, u1, u2, p)
108
109     def lift_factors(self):
110         """
111         Update factors by lifted integer coefficient polynomials Ai's:
112           f == A1*A2 (mod p*q) and
113           Ai == ai (mod q) (i = 1, 2).
114         Moreover, q is updated by p*q
115
116         PRECONDITIONS (automatically satisfied):
117           f == a1*a1 (mod q)
118           a1*u1 + a2*u2 == 1 (mod p)
119           p | q
120         """
121         assert not (self.a1*self.u1 + self.a2*self.u2 - the_one).coefficients_map(lambda c: c % self.p)
122         higher = (self.f - self.a1 * self.a2).scalar_exact_division(self.q)
123         quot, rem = (self.u2 * higher).pseudo_divmod(self.a1)
124         filterq = lambda c: (c % self.p) * self.q
125         self.a1 += rem.coefficients_map(filterq)
126         self.a2 += (self.u1 * higher + self.a2 * quot).coefficients_map(filterq)
127         self.q *= self.p
128         assert not (self.a1*self.u1 + self.a2*self.u2 - the_one).coefficients_map(lambda c: c % self.p)
129
131         """
132         Update u1 and u2 with U1 and U2:
133           a1*U1 + a2*U2 == 1 (mod p**2)
134           Ui == ui (mod p) (i = 1, 2)
135         then, update p with p**2.
136
137         PRECONDITIONS (automatically satisfied):
138           a1*u1 + a2*u2 == 1 (mod p)
139         """
140         higher = (the_one - self.u1 * self.a1 - self.u2 * self.a2).scalar_exact_division(self.p)
141         quot, rem = (self.u2 * higher).pseudo_divmod(self.a1)
142         filterp = lambda c: (c % self.p) * self.p
143         self.u1 += (self.u1 * higher + self.a2 * quot).coefficients_map(filterp)
144         self.u2 += rem.coefficients_map(filterp)
145         self.p *= self.p
146
147     def lift(self):
148         """
149         The lift.
150         """
151         self.lift_factors()
153
154     def _get_factors(self):
155         """
156         getter for factors of target mod q
157         """
158         return [self.a1, self.a2]
159
160     factors = property(_get_factors)
161
162
163 class HenselLiftMulti(object):
164     """
165     A class represents integer polynomials which will be lifted by
166     Hensel's lemma.
167
168     If the number of factors is two, you should use HenselLiftPair.
169     """
170     def __init__(self, target, factors, cofactors, ladder, p, q=None):
171         """
172         HenselLiftMulti(f, factors, ladder, p, q)
173
174         The parameters satisfy that:
175          +  ai's of factors and f are all monic
176          +  bi's of cofactors are product of aj's for i < j
177          +  f == a1*...*ar (mod q)
178          +  ladder is a tuple (sis, tis) and si's of sis and ti's of tis
179             satisfy: ai*si + bi*ti == 1 (mod p),
180         with positive integers p dividing q.
181
182         If p==q, q can be omit from the argument.
183         """
184         self.f = target
185         self.a = list(factors)
186         self.r = len(self.a)
187         self.b = list(cofactors)
190         self.p = p
191         if q is None:
192             self.q = p
193         else:
194             self.q = q
195
196     @classmethod
197     def from_factors(cls, target, factors, p):
198         """
199         Create and return an instance of HenselLiftMulti.
200
201         HenselLiftMulti.from_factors(f, factors, p)
202
203         The parameters satisfy that:
204           f == a1*...*ar (mod p)
205         with a prime number p and ai's in factors.
206         """
207         bi0s = _init_cofactors(factors)
208         sis, tis = [], []
209         for ai0, bi0 in zip(factors, bi0s):
210             s, t = _extgcdp(ai0, bi0, p)[:2]
211             sis.append(s)
212             tis.append(t)
213         return cls(target, factors, bi0s, (sis, tis), p, p)
214
215     def lift_factors(self):
216         """
217         Find a lifted integer coefficient polynomials such that:
218           f = A1*A2*...*Ar (mod q*p),
219           Ai = ai (mod q) (i=1,2,...,r),
220         then, update q with q*p.
221
222         PRECONDITIONS (automatically satisfied):
223           f = a1*a2*...*ar (mod q),
224           ai*si + bi*ti = 1 (mod p) (i=1,2,...,r) and
225           p | q.
226         """
227         mini_target = self.f
228         for i in range(self.r - 1):
229             assert not (self.a[i]*self.s[i] + self.b[i]*self.t[i] - the_one).coefficients_map(lambda c: c % self.p), \
230                    "a = %s\ns = %s\nb = %s\nt = %s" %(self.a[i], self.s[i], self.b[i], self.t[i])
231             mini_pair = HenselLiftPair(mini_target,
232                                        self.a[i], self.b[i],
233                                        self.s[i], self.t[i],
234                                        self.p, self.q)
235             mini_pair.lift_factors()
236             self.a[i] = mini_pair.a1
237             mini_target = self.b[i] = mini_pair.a2
238         self.a[self.r - 1] = mini_target
239         self.q = mini_pair.q
240
242         """
243         Update ladder si's and ti's with Si's and Ti's:
244           ai*Si + bi*Ti == 1 (mod p**2) (i = 1, 2, ..., r),
245           Si == si (mod p) (i = 1, 2, ..., r) and
246           Ti == ti (mod p) (i = 1, 2, ..., r)
247         then, update p with p**2.
248
249         PRECONDITIONS (automatically satisfied):
250           ai*si + bi*ti == 1 (mod p)
251         """
252         for i in range(self.r - 1):
253             mini_pair = HenselLiftPair(self.f,
254                                        self.a[i], self.b[i],
255                                        self.s[i], self.t[i],
256                                        self.p, self.q)
258             self.s[i], self.t[i] = mini_pair.u1, mini_pair.u2
259         self.p *= self.p
260
261     def lift(self):
262         """
263         The lift.
264         """
265         self.lift_factors()
267
268     def _get_factors(self):
269         """
270         getter for factors of target mod q
271         """
272         return list(self.a)
273
274     factors = property(_get_factors)
275
276
277 class HenselLiftSimultaneously(object):
278     """
279     INVARIANTS:
280       ai's, pi's and gi's are monic
281       f == g1*g2*...*gr (mod p)
282       f == d0 + d1*p + d2*p**2 +...+ dk*p**k
283       hi == g(i+1)*...*gr
284       1 == gi*si + hi*ti (mod p) (i = 1, 2,..., r)
285       deg(si) < deg(hi), deg(ti) < deg(gi)
286       p | q
287       f == l1*l2*...*lr (mod q/p)
288       f == a1*a2*...*ar (mod q)
289       ui == ai*yi + bi*zi (mod p) (i = 1, 2,..., r)
290
291     REFERENCE:
292     G.E.Collins & M.J.Encarnaci'on.
293     Improved Techniques for factoring Univariate Polynomials,
294     J.Symb.Comp. 21, 313--327 (1996).
295     """
296
297     def __init__(self, target, factors, cofactors, bases, p):
298         """
299         Prepare attributes
300         """
301         # invariants
302         self.f = target
303         self.gis = tuple(factors)
304         self.his = tuple(cofactors)
305         self.r = len(self.gis)
306         self.sis = tuple(bases[0])
307         self.tis = tuple(bases[1])
308         self.p = p
309         self.dis = []
310
311         # q will be a power of p
312         self.q = self.p
313
314         # these are containers
315         self.ais, self.bis = {}, {} # used later
316         self.uis = self.yis = self.zis = None # used later
317
318         # functions
319         self.modp = lambda c: c % self.p
320
321     def _init_dis(self):
322         """
323         Return p-adic expansion of target polynomial.
324
325         This method is private for __init__.
326         """
327         self.dis.append(self.f.coefficients_map(self.modp))
328         rest = self.f - self.dis[-1]
329         q = self.p
330         while any(abs(x) >= q for x in self.f.itercoefficients()):
331             self.dis.append(self.f.coefficients_map(lambda c: (c // q) % self.p))
332             rest -= self.dis[-1].scalar_mul(q)
333             q *= self.p
334         if self.f != rest:
335             self.dis.append(rest.scalar_exact_division(q))
336
337     @classmethod
338     def from_factors(cls, target, factors, p, ubound=sys.maxint):
339         """
340         Create and return an instance of HenselLiftSimultaneously,
341         whose factors are lifted by HenselLiftMulti upto ubound (if it
342         is smaller than sys.maxint, or upto sys.maxint).
343
344         HenselLiftSimultaneously.from_factors(f, factors, p)
345
346         The parameters satisfy that:
347           f == a1*...*ar (mod p)
348         with a prime number p and ai's in factors.
349         """
350         lifter = HenselLiftMulti.from_factors(target, factors, p)
351         interbound = min(ubound, sys.maxint)
352         while lifter.p < interbound:
353             lifter.lift_factors()
355         new = cls(target, lifter.a, lifter.b, (lifter.s, lifter.t), lifter.p)
356         new.first_lift()
357         return new
358
359     def first_lift(self):
360         """
361         Start lifting.
362             f == l1*l2*...*lr (mod p**2)
363         Initialize di's, ui's, yi's and zi's.
364         Update ai's, bi's.
365         Then, update q with p**2.
366         """
367         assert self.p == self.q # q has not been raised yet
368         self._assertEqualModulo(self.f, arith1.product(self.gis), self.q)
369         for i in range(self.r - 1):
370             self._assertEqualModulo(self.f, arith1.product(self.gis[:i+1])*self.his[i], self.q)
371         assert self.gis[-1] == self.his[self.r - 2]
372
373         self._init_dis()
374         if len(self.dis) > 1:
375             mini_target = self.dis[0] + self.dis[1].scalar_mul(self.p)
376             self._assertEqualModulo(self.f, mini_target, self.p**2)
377         else:
378             mini_target = self.dis[0]
379         aj, bj = [], []
380         self.uis, self.yis, self.zis = [], {}, {}
381         for i in xrange(self.r - 1):
382             dividend = mini_target - self.gis[i] * self.his[i]
383             self.uis.append(dividend.scalar_exact_division(self.p))
384             self.yis[i], self.zis[i] = self._solve_yz(i)
385             aj.append(self.gis[i] + self.zis[i].scalar_mul(self.q))
386             bj.append(self.his[i] + self.yis[i].scalar_mul(self.q))
387             self._assertEqualModulo(mini_target, aj[i]*bj[i], self.q*self.p)
388             mini_target = bj[i]
389         self._assertEqualModulo(self.gis[-1], mini_target, self.q)
390         aj.append(bj[-1])
391         self.q *= self.p
392
393         for i in range(self.r - 1):
394             self._assertEqualModulo(self.f, arith1.product(aj[:i+1])*bj[i], self.q, "f != l%d * m%d" % (i, i))
395             self._assertEqualModulo(self.gis[i], aj[i], self.p)
396         self._assertEqualModulo(self.f, arith1.product(aj), self.q)
397
398         self.ais[-1], self.ais[0] = self.gis, tuple(aj)
399         self.bis[-1], self.bis[0] = self.his, tuple(bj)
400
401     def general_lift(self):
402         """
403         Continue lifting.
404             f == a1*a2*...*ar (mod p*q)
405         Update ai's, bi's, ui's, yi's and zi's.
406         Then, update q with p*q.
407         """
408         j = arith1.vp(self.q, self.p)[0]
409         if len(self.dis) > j:
410             mini_target = self.dis[j]
411         else:
412             mini_target = the_zero
413
414         self.ais[-2], self.ais[-1] = self.ais[-1], self.ais[0]
415         self.bis[-2], self.bis[-1] = self.bis[-1], self.bis[0]
416         aj, bj = [], []
417         for i in xrange(self.r - 1):
418             yi, zi = self.yis[i], self.zis[i]
419             dividend = self.ais[-2][i]*yi + self.bis[-2][i]*zi - self.uis[i]
420             v_j = dividend.scalar_exact_division(self.p)
421             if j == 2:
422                 self.uis[i] = mini_target - v_j - yi*zi
423             else:
424                 self.uis[i] = mini_target - v_j - (yi*zi).scalar_mul(self.p**(j - 2))
425             self.yis[i], self.zis[i] = self._solve_yz(i)
426             aj.append(self.ais[-1][i] + self.zis[i].scalar_mul(self.q))
427             bj.append(self.bis[-1][i] + self.yis[i].scalar_mul(self.q))
428             self._assertEqualModulo(self.f, arith1.product(aj)*bj[i], self.q*self.p, (i, j))
429             mini_target = self.yis[i]
430         aj.append(bj[-1])
431         self.q *= self.p
432
433         for i in range(self.r - 1):
434             self._assertEqualModulo(self.f, arith1.product(aj[:i+1])*bj[i], self.q, "f != a%d * b%d" % (i, i))
435             self._assertEqualModulo(self.gis[i], aj[i], self.p)
436         self._assertEqualModulo(self.f, arith1.product(aj), self.q)
437
438         self.ais[0] = tuple(aj)
439         self.bis[0] = tuple(bj)
440
441     def lift(self):
442         """
443         The lift
444         """
445         self.general_lift()
446
447     def _solve_yz(self, i):
448         """
449         Solve the equation
450             gi Y + hi Z = ui (mod p)
451         satisfying conditions:
452         (1) deg(Y) < deg(hi) and
453         (2) deg(Z) < deg(gi),
454         and return a tuple(Y, Z).
455         The method needs the following conditions:
456         (I) deg(ui) <= deg(gi)deg(hi),
457         (II)  gi si + hi ti = 1 (mod p).
458         """
459         umodp = self.uis[i].coefficients_map(self.modp)
460         quot, rem = (self.tis[i] * umodp).pseudo_divmod(self.gis[i])
461         y = (self.sis[i] * umodp + self.his[i] * quot).coefficients_map(self.modp)
462         z = rem.coefficients_map(self.modp)
463
464         assert y.degree() < self.his[i].degree()
465         self._assertEqualModulo(self.gis[i] * y + self.his[i] * z,
466                                 self.uis[i], self.p)
467
468         return y, z
469
470     def _get_factors(self):
471         """
472         getter for factors of target mod q
473         """
474         if self.ais:
475             return list(self.ais[0])
476         else:
477             return list(self.gis)
478
479     factors = property(_get_factors)
480
481     def _assertEqualModulo(self, expected, actual, modulus, message=None):
482         """
483         assert expected == actual (mod modulus)
484         """
485         if message is None:
486             for d, c in actual - expected:
487                 assert 0 == c % modulus, "degree %d\nactual = %s" % (d, str(actual))
488         else:
489             for d, c in actual - expected:
490                 assert 0 == c % modulus, "degree %d\nactual = %s\n%s" % (d, str(actual), message)
491
492
493 def lift_upto(target, factors, p, bound):
494     """
495     Hensel lift factors mod p of target upto bound and return factors
496     mod q and q.
497
498     PRECONDITIONS:
499     target == product(factors) mod p
500     POSTCONDITIONS:
501     there exist k s.t. q == p**k >= bound and
502     target == product(result) mod q
503     """
504     if len(factors) == 2:
505         lifter = HenselLiftPair.from_factors(target, factors[0], factors[1], p)
506     elif bound < sys.maxint:
507         lifter = HenselLiftMulti.from_factors(target, factors, p)
508     else:
509         lifter = HenselLiftSimultaneously.from_factors(target, factors, p)
510     while lifter.q < bound:
511         lifter.lift()
512     return lifter.factors, lifter.q
```