# Ticket #2179: mpoly-factor.patch

File mpoly-factor.patch, 11.5 KB (added by jbmohler, 15 years ago)
• ## sage/rings/polynomial/multi_polynomial_element.py

```# HG changeset patch
# User Joel B. Mohler <jbm5@lehigh.edu>
# Date 1203193950 18000
# Node ID e260778cd2581075aede8855c50dd8ee2b925504
# Parent  fe7c1be520fa34a1ceb6d88133e7f9c5c5d2be64
kroneckers trick factoring algorithm for multi-variate polynomials over ZZ

diff -r fe7c1be520fa -r e260778cd258 sage/rings/polynomial/multi_polynomial_element.py```
 a import polydict import polydict from sage.structure.factorization import Factorization from multi_polynomial_factor import factor_ZZx from sage.rings.polynomial.polynomial_singular_interface import Polynomial_singular_repr class MPolynomial_polydict(Polynomial_si sage: k.factor() (s^2 + 2/3) * (x + s*y)^2 * (x + (-s)*y)^5 * (x^2 + s*x*y + s^2*y^2)^5 """ if self.parent().base() is ZZ: return factor_ZZx(self) # I do not think this applied anymore.  Or at least it's # more relevant to optimizing the NTL build. #\note{Singular multi-variate polynomial factorization is very
• ## new file sage/rings/polynomial/multi_polynomial_factor.py

`diff -r fe7c1be520fa -r e260778cd258 sage/rings/polynomial/multi_polynomial_factor.py`
 - """ Factorization algorithm for multivariate polynomials over ZZ. The algorithm uses Kronecker's Trick -- that is specialize a polynomial in n variables by substituting a large (prime) integer to get a polynomial in n-1 variables.  Iterate until we have a bi-variate or uni-variate and factor directly. AUTHORS: -- Joel B. Mohler: Initial implementation TESTS: sage: R. = ZZ[] sage: f = (1-x)*(z-y) sage: f.factor() (-1) * (-1*y + z) * (x - 1) sage: R.=ZZ[] sage: t=-p10^170*X1^10*X2^10+p10^130*X1^10*X2^5+p10^130*X1^5*X2^10-p10^90*X1^5*X2^5+p10^80*X1^5*X2^5-p10^40*X1^5-p10^40*X2^5+1 sage: t.factor() (-1) * (p10^8*X2 - 1) * (p10^8*X1 - 1) * (p10^18*X1*X2 - 1) * (p10^32*X2^4 + p10^24*X2^3 + p10^16*X2^2 + p10^8*X2 + 1) * (p10^32*X1^4 + p10^24*X1^3 + p10^16*X1^2 + p10^8*X1 + 1) * (p10^72*X1^4*X2^4 + p10^54*X1^3*X2^3 + p10^36*X1^2*X2^2 + p10^18*X1*X2 + 1) """ ################################################################################ #       Copyright (C) 2007 Joel B. Mohler #                          William Stein # #  Distributed under the terms of the GNU General Public License (GPL) # #                  http://www.gnu.org/licenses/ ################################################################################ import sage.rings.polynomial.polydict as polydict from sage.rings.integer_ring import ZZ from sage.rings.integer import Integer from sage.structure.factorization import Factorization from sage.rings.arith import gcd,next_prime from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing def pAdicExpansion(c,prime): assert(c>=0) l = [] while c > 0: digit = c % prime c = c//prime if 2*digit > prime: digit -= prime c += 1 l.append( digit ) return l def poly_factory(ring,pd): from sage.rings.polynomial.multi_polynomial_element import MPolynomial_polydict return MPolynomial_polydict(ring,polydict.PolyDict(pd,zero=ZZ(0),force_int_exponents=False,force_etuples=False)) def eval_coeff_list_at_gen(mpoly_Ring, gen, coeffList): ngens = mpoly_Ring.ngens() pd = {} t = [0] * ngens for i in range(len(coeffList)): if not coeffList[i].is_zero(): t[gen] = i pd[polydict.ETuple(t)] = coeffList[i] return poly_factory(mpoly_Ring,pd) def _factor_single_variable(p,i): coeffs = [0]*(p.total_degree()+1) for m,c in p.dict().iteritems(): coeffs[m[i]] = c uni_results = PolynomialRing(ZZ,'x')(coeffs).factor() #print uni_results return Factorization([(eval_coeff_list_at_gen(p.parent(),i,fac.list()),e) for fac,e in uni_results],unit=p.parent()(ZZ(uni_results.unit())),sort=False) def _factor_kronecker_core(p,coeffSize,to_spec): ngens = p.parent().ngens() #print p if len(to_spec) == 1: return _factor_single_variable(p,to_spec[0][0]) else: index_to_specialize = to_spec[0][0] #print index_to_specialize gens = p.parent().gens() carry_clearance = 3 while True: prime = next_prime(carry_clearance*(coeffSize+1)) pd = {} for m,c in p.dict().iteritems(): if m[index_to_specialize] > 0: m = list(m) c = c*prime**m[index_to_specialize] coeffSize = max(coeffSize,c.abs()) m[index_to_specialize] = 0 try: pd[polydict.ETuple(m)] += c except KeyError: pd[polydict.ETuple(m)] = c specialized_p = poly_factory(p.parent(),pd) results = _factor_kronecker_core(specialized_p,coeffSize,to_spec[1:]) unit = results.unit() #print results # read back from specialized factors fac_results = [] for fac, e in results: fac_result = 0 for m,c in fac.dict().iteritems(): sign=1 if c < 0: sign = -1 c *= -1 c = eval_coeff_list_at_gen(p.parent(),index_to_specialize,pAdicExpansion(c,prime)) fac_result += c*poly_factory(p.parent(),{m:ZZ(sign)}) if fac.is_constant(): # fac was an integer, it may have been unfactored #  Thus, fac_result is a poly in one variable which may need further factoring. more_details = _factor_single_variable(fac_result, index_to_specialize) fac_results += list(more_details) unit *= more_details.unit() else: fac_results.append((fac_result,e)) f=Factorization(fac_results,unit=unit,sort=False) # test cases indicate that a carry_clearance of 3 is always sufficient # However, I can't find a proof of this so I'm iterating to make sure. # That is, specializing a variable always results in at least as many factors as the unspecialized version -- however, sometimes it results in more factors! if f.value() == p: break else: carry_clearance += 1 return f def factor_ZZx(p): if p.parent().base() is not ZZ: raise TypeError, "I can only factor polys over ZZ this way so far" if p == 0: raise ArithmeticError, "Prime factorization of 0 not defined." one = p.parent()(1) ngens = p.parent().ngens() # select a variable to specialize pdict = p.dict() min_exp = None max_exp = None gcd_exp = None for m,c in pdict.iteritems(): if min_exp is None: min_exp = list(m) max_exp = list(m) gcd_exp = list(m) else: for i in range(ngens): min_exp[i] = min(min_exp[i],m[i]) max_exp[i] = max(max_exp[i],m[i]) gcd_exp[i] = gcd(gcd_exp[i],m[i]) results = [] if max(min_exp) > 0: for i in range(ngens): if min_exp[i] > 0: results.append((p.parent().gen(i),min_exp[i])) s = [(i,max_exp[i]-min_exp[i]) for i in range(ngens) if max_exp[i] > min_exp[i]] if len(s) == 0: # constant times a monomial if p.dict().values()[0] != 1: results.append((p.parent()(p.dict().values()[0]),1)) return Factorization(results,unit=one) elif len(s) == 1: # a single variable, coerce to ZZ[x] relabeled_p = p if max(min_exp) > 0: pd = {} for m,c in pdict.iteritems(): t = [(m[i]-min_exp[i]) for i in range(ngens)] pd[polydict.ETuple(t)] = c relabeled_p = poly_factory(p.parent(),pd) uni_results = _factor_single_variable(relabeled_p,s[0][0]) #print uni_results return Factorization(results + list(uni_results),unit=uni_results.unit()) else: # lots of variables involved (i.e. > 1) # see if it is a good idea to relabel something like x^n -> x gcd_exp = None max_coeff = 0 for m,c in pdict.iteritems(): max_coeff = max(max_coeff,c.abs()) if gcd_exp is None: gcd_exp = [m[i]-min_exp[i] for i in range(ngens)] else: for i in range(ngens): gcd_exp[i] = gcd(gcd_exp[i],m[i]-min_exp[i]) gcd_exp = [g if g>0 else 1 for g in gcd_exp] if len([i for i in range(ngens) if gcd_exp[i] > 1]) > 0 or max(min_exp) > 0: pd = {} for m,c in pdict.iteritems(): t = [(m[i]-min_exp[i])//gcd_exp[i] for i in range(ngens)] pd[polydict.ETuple(t)] = c relabeled_p = poly_factory(p.parent(),pd) else: relabeled_p = p max_exp = [(max_exp[i]-min_exp[i])//(gcd_exp[i]) for i in range(ngens)] min_exp = [0] * ngens # this was the point after all -- i.e. to factor out obvious common factors to_specialize = [(i,max_exp[i]) for i in range(ngens) if max_exp[i] > 0] to_specialize.sort(key=lambda x:x[1]) kron_results = _factor_kronecker_core(relabeled_p,max_coeff,to_specialize) # re-expand x->x^n if compressed units = kron_results.unit() if len([i for i in range(ngens) if gcd_exp[i] > 1]) > 0: for fac,e in kron_results: pd = {} fac_min_exp = None fac_max_exp = None fac_maxcoeff = 0 must_refactor = False for m,c in fac.dict().iteritems(): fac_maxcoeff = max(fac_maxcoeff,c.abs()) if not must_refactor: for i in range(ngens): if m[i] > 0 and gcd_exp[i] != 1: must_refactor = True t = [m[i]*gcd_exp[i] for i in range(ngens)] if fac_min_exp is None: fac_min_exp = list(m) fac_max_exp = list(m) else: for i in range(ngens): fac_min_exp[i] = min(fac_min_exp[i],m[i]) fac_max_exp[i] = max(fac_max_exp[i],m[i]) pd[polydict.ETuple(t)] = c if must_refactor: to_specialize = [(i,fac_max_exp[i]) for i in range(ngens) if fac_max_exp[i] > 0] to_specialize.sort(key=lambda x:x[1]) fac_results = _factor_kronecker_core(poly_factory(p.parent(),pd),fac_maxcoeff,to_specialize) units *= fac_results.unit() results += list(fac_results) else: results.append( (fac,e) ) else: results += list(kron_results) return Factorization(results,unit=units)