"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:
   37         zlc = z.leading_coefficient()
   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     """
   72     def __init__(self, target, factor1, factor2, ladder1, ladder2, p, q=None):
   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
   87         self.u1 = ladder1
   88         self.u2 = ladder2
   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 
  130     def lift_ladder(self):
  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()
  152         self.lift_ladder()
  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)
  188         self.s = list(ladder[0])
  189         self.t = list(ladder[1])
  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 
  241     def lift_ladder(self):
  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)
  257             mini_pair.lift_ladder()
  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()
  266         self.lift_ladder()
  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()
  354             lifter.lift_ladder()
  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