"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     """
   19     # keep leading coefficient
   20     lf = f.leading_coefficient()
   21 
   22     # p-adic factorization
   23     p, fp_factors = padic_factorization(f)
   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 
   44 def padic_factorization(f):
   45     """
   46     padic_factorization(f) -> p, factors
   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     """
  218     lf = f.leading_coefficient()
  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