NZMATH  1.2.0 About: NZMATH is a Python based number theory oriented calculation system.   Fossies Dox: NZMATH-1.2.0.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)
hensel.py
Go to the documentation of this file.
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
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
nzmath.poly.hensel.HenselLiftSimultaneously.lift
def lift(self)
Definition: hensel.py:441
nzmath.poly.hensel.HenselLiftMulti.lift
def lift(self)
Definition: hensel.py:261
nzmath.poly.hensel.HenselLiftPair.__init__
Definition: hensel.py:72
nzmath.bigrange.range
def range(start, stop=None, step=None)
Definition: bigrange.py:19
nzmath.poly.hensel.HenselLiftSimultaneously._get_factors
def _get_factors(self)
Definition: hensel.py:470
nzmath.poly.hensel.HenselLiftMulti.lift_factors
def lift_factors(self)
Definition: hensel.py:215
nzmath.poly.hensel.HenselLiftSimultaneously.uis
uis
Definition: hensel.py:316
nzmath.poly.hensel.HenselLiftPair.p
p
Definition: hensel.py:89
nzmath.poly.hensel.HenselLiftPair.a1
a1
Definition: hensel.py:85
nzmath.poly.hensel._init_cofactors
def _init_cofactors(factors)
Definition: hensel.py:54
nzmath.poly.hensel.HenselLiftMulti.a
a
Definition: hensel.py:185
nzmath.poly.hensel.HenselLiftMulti.b
b
Definition: hensel.py:187
Definition: hensel.py:130
nzmath.poly.hensel.HenselLiftPair.a2
a2
Definition: hensel.py:86
nzmath.poly.hensel.HenselLiftMulti.r
r
Definition: hensel.py:186
nzmath.poly.ring
Definition: ring.py:1
nzmath.poly.hensel.HenselLiftSimultaneously.yis
yis
Definition: hensel.py:316
nzmath.poly.hensel.HenselLiftSimultaneously.bis
bis
Definition: hensel.py:315
nzmath.rational
Definition: rational.py:1
nzmath.poly.hensel.HenselLiftMulti.s
s
Definition: hensel.py:188
nzmath.poly.hensel.HenselLiftMulti.p
p
Definition: hensel.py:190
nzmath.poly.hensel.HenselLiftSimultaneously._solve_yz
def _solve_yz(self, i)
Definition: hensel.py:447
nzmath.poly.hensel.HenselLiftPair.from_factors
def from_factors(cls, target, factor1, factor2, p)
Definition: hensel.py:96
nzmath.poly.hensel.HenselLiftSimultaneously.zis
zis
Definition: hensel.py:316
nzmath.poly.hensel.HenselLiftSimultaneously.__init__
def __init__(self, target, factors, cofactors, bases, p)
Definition: hensel.py:297
nzmath.poly.hensel.HenselLiftMulti.t
t
Definition: hensel.py:189
nzmath.poly.hensel.HenselLiftPair.u2
u2
Definition: hensel.py:88
nzmath.poly.hensel.HenselLiftPair.lift
def lift(self)
Definition: hensel.py:147
nzmath.poly.hensel.HenselLiftSimultaneously.his
his
Definition: hensel.py:304
nzmath.poly.hensel.HenselLiftPair.u1
u1
Definition: hensel.py:87
nzmath.poly.hensel.HenselLiftMulti.__init__
def __init__(self, target, factors, cofactors, ladder, p, q=None)
Definition: hensel.py:170
nzmath.poly.hensel.lift_upto
def lift_upto(target, factors, p, bound)
Definition: hensel.py:493
nzmath.poly.hensel.HenselLiftMulti._get_factors
def _get_factors(self)
Definition: hensel.py:268
nzmath.arith1
Definition: arith1.py:1
nzmath.poly.hensel.HenselLiftSimultaneously._assertEqualModulo
def _assertEqualModulo(self, expected, actual, modulus, message=None)
Definition: hensel.py:481
nzmath.poly.hensel.HenselLiftSimultaneously.modp
modp
Definition: hensel.py:319
nzmath.poly.hensel.HenselLiftPair
Definition: hensel.py:67
nzmath.poly.hensel.HenselLiftMulti.q
q
Definition: hensel.py:192
nzmath.poly.hensel.HenselLiftPair.q
q
Definition: hensel.py:91
nzmath.poly.hensel.HenselLiftSimultaneously._init_dis
def _init_dis(self)
Definition: hensel.py:321
nzmath.poly.hensel.HenselLiftPair.f
f
Definition: hensel.py:84
nzmath.poly.hensel.HenselLiftMulti.from_factors
def from_factors(cls, target, factors, p)
Definition: hensel.py:197
nzmath.poly.hensel.HenselLiftMulti
Definition: hensel.py:163
nzmath.poly.hensel._extgcdp
def _extgcdp(f, g, p)
Definition: hensel.py:25
nzmath.poly.hensel.HenselLiftSimultaneously.general_lift
def general_lift(self)
Definition: hensel.py:401
nzmath.poly.hensel.HenselLiftSimultaneously.sis
sis
Definition: hensel.py:306
nzmath.poly.hensel.HenselLiftSimultaneously.first_lift
def first_lift(self)
Definition: hensel.py:359
nzmath.poly.hensel.HenselLiftPair._get_factors
def _get_factors(self)
Definition: hensel.py:154
nzmath.poly.hensel.HenselLiftSimultaneously.tis
tis
Definition: hensel.py:307
Definition: hensel.py:241
nzmath.poly.hensel.HenselLiftMulti.f
f
Definition: hensel.py:184
nzmath.poly.hensel.HenselLiftSimultaneously.gis
gis
Definition: hensel.py:303
nzmath.poly.hensel.HenselLiftSimultaneously.from_factors
def from_factors(cls, target, factors, p, ubound=sys.maxint)
Definition: hensel.py:338
nzmath.poly.hensel.HenselLiftPair.lift_factors
def lift_factors(self)
Definition: hensel.py:109
nzmath.poly.hensel.HenselLiftSimultaneously.dis
dis
Definition: hensel.py:309