| 1 | """ |
| 2 | Factorization algorithm for multivariate polynomials over ZZ. |
| 3 | The algorithm uses Kronecker's Trick -- that is specialize a polynomial |
| 4 | in n variables by substituting a large (prime) integer to get a polynomial |
| 5 | in n-1 variables. Iterate until we have a bi-variate or uni-variate and |
| 6 | factor directly. |
| 7 | |
| 8 | AUTHORS: |
| 9 | -- Joel B. Mohler: Initial implementation |
| 10 | |
| 11 | TESTS: |
| 12 | sage: R.<x,y,z> = ZZ[] |
| 13 | sage: f = (1-x)*(z-y) |
| 14 | sage: f.factor() |
| 15 | (-1) * (-1*y + z) * (x - 1) |
| 16 | sage: R.<p10,g0,g1,g2,g3,g4,X1,X2>=ZZ[] |
| 17 | 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 |
| 18 | sage: t.factor() |
| 19 | (-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) |
| 20 | """ |
| 21 | |
| 22 | ################################################################################ |
| 23 | # Copyright (C) 2007 Joel B. Mohler <joel@kiwistrawberry.us> |
| 24 | # William Stein <wstein@gmail.com> |
| 25 | # |
| 26 | # Distributed under the terms of the GNU General Public License (GPL) |
| 27 | # |
| 28 | # http://www.gnu.org/licenses/ |
| 29 | ################################################################################ |
| 30 | import sage.rings.polynomial.polydict as polydict |
| 31 | from sage.rings.integer_ring import ZZ |
| 32 | from sage.rings.integer import Integer |
| 33 | from sage.structure.factorization import Factorization |
| 34 | from sage.rings.arith import gcd,next_prime |
| 35 | from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing |
| 36 | |
| 37 | def pAdicExpansion(c,prime): |
| 38 | assert(c>=0) |
| 39 | l = [] |
| 40 | while c > 0: |
| 41 | digit = c % prime |
| 42 | c = c//prime |
| 43 | if 2*digit > prime: |
| 44 | digit -= prime |
| 45 | c += 1 |
| 46 | l.append( digit ) |
| 47 | return l |
| 48 | |
| 49 | def poly_factory(ring,pd): |
| 50 | from sage.rings.polynomial.multi_polynomial_element import MPolynomial_polydict |
| 51 | return MPolynomial_polydict(ring,polydict.PolyDict(pd,zero=ZZ(0),force_int_exponents=False,force_etuples=False)) |
| 52 | |
| 53 | def eval_coeff_list_at_gen(mpoly_Ring, gen, coeffList): |
| 54 | ngens = mpoly_Ring.ngens() |
| 55 | pd = {} |
| 56 | t = [0] * ngens |
| 57 | for i in range(len(coeffList)): |
| 58 | if not coeffList[i].is_zero(): |
| 59 | t[gen] = i |
| 60 | pd[polydict.ETuple(t)] = coeffList[i] |
| 61 | return poly_factory(mpoly_Ring,pd) |
| 62 | |
| 63 | def _factor_single_variable(p,i): |
| 64 | coeffs = [0]*(p.total_degree()+1) |
| 65 | for m,c in p.dict().iteritems(): |
| 66 | coeffs[m[i]] = c |
| 67 | uni_results = PolynomialRing(ZZ,'x')(coeffs).factor() |
| 68 | #print uni_results |
| 69 | 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) |
| 70 | |
| 71 | def _factor_kronecker_core(p,coeffSize,to_spec): |
| 72 | ngens = p.parent().ngens() |
| 73 | #print p |
| 74 | if len(to_spec) == 1: |
| 75 | return _factor_single_variable(p,to_spec[0][0]) |
| 76 | else: |
| 77 | index_to_specialize = to_spec[0][0] |
| 78 | #print index_to_specialize |
| 79 | gens = p.parent().gens() |
| 80 | carry_clearance = 3 |
| 81 | while True: |
| 82 | prime = next_prime(carry_clearance*(coeffSize+1)) |
| 83 | pd = {} |
| 84 | for m,c in p.dict().iteritems(): |
| 85 | if m[index_to_specialize] > 0: |
| 86 | m = list(m) |
| 87 | c = c*prime**m[index_to_specialize] |
| 88 | coeffSize = max(coeffSize,c.abs()) |
| 89 | m[index_to_specialize] = 0 |
| 90 | try: |
| 91 | pd[polydict.ETuple(m)] += c |
| 92 | except KeyError: |
| 93 | pd[polydict.ETuple(m)] = c |
| 94 | specialized_p = poly_factory(p.parent(),pd) |
| 95 | results = _factor_kronecker_core(specialized_p,coeffSize,to_spec[1:]) |
| 96 | unit = results.unit() |
| 97 | #print results |
| 98 | # read back from specialized factors |
| 99 | fac_results = [] |
| 100 | for fac, e in results: |
| 101 | fac_result = 0 |
| 102 | for m,c in fac.dict().iteritems(): |
| 103 | sign=1 |
| 104 | if c < 0: |
| 105 | sign = -1 |
| 106 | c *= -1 |
| 107 | c = eval_coeff_list_at_gen(p.parent(),index_to_specialize,pAdicExpansion(c,prime)) |
| 108 | fac_result += c*poly_factory(p.parent(),{m:ZZ(sign)}) |
| 109 | if fac.is_constant(): |
| 110 | # fac was an integer, it may have been unfactored |
| 111 | # Thus, fac_result is a poly in one variable which may need further factoring. |
| 112 | more_details = _factor_single_variable(fac_result, index_to_specialize) |
| 113 | fac_results += list(more_details) |
| 114 | unit *= more_details.unit() |
| 115 | else: |
| 116 | fac_results.append((fac_result,e)) |
| 117 | f=Factorization(fac_results,unit=unit,sort=False) |
| 118 | # test cases indicate that a carry_clearance of 3 is always sufficient |
| 119 | # However, I can't find a proof of this so I'm iterating to make sure. |
| 120 | # That is, specializing a variable always results in at least as many factors as the unspecialized version -- however, sometimes it results in more factors! |
| 121 | if f.value() == p: |
| 122 | break |
| 123 | else: |
| 124 | carry_clearance += 1 |
| 125 | return f |
| 126 | |
| 127 | def factor_ZZx(p): |
| 128 | if p.parent().base() is not ZZ: |
| 129 | raise TypeError, "I can only factor polys over ZZ this way so far" |
| 130 | |
| 131 | if p == 0: |
| 132 | raise ArithmeticError, "Prime factorization of 0 not defined." |
| 133 | |
| 134 | one = p.parent()(1) |
| 135 | ngens = p.parent().ngens() |
| 136 | # select a variable to specialize |
| 137 | pdict = p.dict() |
| 138 | min_exp = None |
| 139 | max_exp = None |
| 140 | gcd_exp = None |
| 141 | for m,c in pdict.iteritems(): |
| 142 | if min_exp is None: |
| 143 | min_exp = list(m) |
| 144 | max_exp = list(m) |
| 145 | gcd_exp = list(m) |
| 146 | else: |
| 147 | for i in range(ngens): |
| 148 | min_exp[i] = min(min_exp[i],m[i]) |
| 149 | max_exp[i] = max(max_exp[i],m[i]) |
| 150 | gcd_exp[i] = gcd(gcd_exp[i],m[i]) |
| 151 | |
| 152 | results = [] |
| 153 | if max(min_exp) > 0: |
| 154 | for i in range(ngens): |
| 155 | if min_exp[i] > 0: |
| 156 | results.append((p.parent().gen(i),min_exp[i])) |
| 157 | |
| 158 | s = [(i,max_exp[i]-min_exp[i]) for i in range(ngens) if max_exp[i] > min_exp[i]] |
| 159 | if len(s) == 0: |
| 160 | # constant times a monomial |
| 161 | if p.dict().values()[0] != 1: |
| 162 | results.append((p.parent()(p.dict().values()[0]),1)) |
| 163 | return Factorization(results,unit=one) |
| 164 | elif len(s) == 1: |
| 165 | # a single variable, coerce to ZZ[x] |
| 166 | relabeled_p = p |
| 167 | if max(min_exp) > 0: |
| 168 | pd = {} |
| 169 | for m,c in pdict.iteritems(): |
| 170 | t = [(m[i]-min_exp[i]) for i in range(ngens)] |
| 171 | pd[polydict.ETuple(t)] = c |
| 172 | relabeled_p = poly_factory(p.parent(),pd) |
| 173 | uni_results = _factor_single_variable(relabeled_p,s[0][0]) |
| 174 | #print uni_results |
| 175 | return Factorization(results + list(uni_results),unit=uni_results.unit()) |
| 176 | else: |
| 177 | # lots of variables involved (i.e. > 1) |
| 178 | # see if it is a good idea to relabel something like x^n -> x |
| 179 | gcd_exp = None |
| 180 | max_coeff = 0 |
| 181 | for m,c in pdict.iteritems(): |
| 182 | max_coeff = max(max_coeff,c.abs()) |
| 183 | if gcd_exp is None: |
| 184 | gcd_exp = [m[i]-min_exp[i] for i in range(ngens)] |
| 185 | else: |
| 186 | for i in range(ngens): |
| 187 | gcd_exp[i] = gcd(gcd_exp[i],m[i]-min_exp[i]) |
| 188 | |
| 189 | gcd_exp = [g if g>0 else 1 for g in gcd_exp] |
| 190 | |
| 191 | if len([i for i in range(ngens) if gcd_exp[i] > 1]) > 0 or max(min_exp) > 0: |
| 192 | pd = {} |
| 193 | for m,c in pdict.iteritems(): |
| 194 | t = [(m[i]-min_exp[i])//gcd_exp[i] for i in range(ngens)] |
| 195 | pd[polydict.ETuple(t)] = c |
| 196 | relabeled_p = poly_factory(p.parent(),pd) |
| 197 | else: |
| 198 | relabeled_p = p |
| 199 | |
| 200 | max_exp = [(max_exp[i]-min_exp[i])//(gcd_exp[i]) for i in range(ngens)] |
| 201 | min_exp = [0] * ngens # this was the point after all -- i.e. to factor out obvious common factors |
| 202 | |
| 203 | to_specialize = [(i,max_exp[i]) for i in range(ngens) if max_exp[i] > 0] |
| 204 | to_specialize.sort(key=lambda x:x[1]) |
| 205 | |
| 206 | kron_results = _factor_kronecker_core(relabeled_p,max_coeff,to_specialize) |
| 207 | # re-expand x->x^n if compressed |
| 208 | units = kron_results.unit() |
| 209 | if len([i for i in range(ngens) if gcd_exp[i] > 1]) > 0: |
| 210 | for fac,e in kron_results: |
| 211 | pd = {} |
| 212 | fac_min_exp = None |
| 213 | fac_max_exp = None |
| 214 | fac_maxcoeff = 0 |
| 215 | must_refactor = False |
| 216 | for m,c in fac.dict().iteritems(): |
| 217 | fac_maxcoeff = max(fac_maxcoeff,c.abs()) |
| 218 | if not must_refactor: |
| 219 | for i in range(ngens): |
| 220 | if m[i] > 0 and gcd_exp[i] != 1: |
| 221 | must_refactor = True |
| 222 | t = [m[i]*gcd_exp[i] for i in range(ngens)] |
| 223 | if fac_min_exp is None: |
| 224 | fac_min_exp = list(m) |
| 225 | fac_max_exp = list(m) |
| 226 | else: |
| 227 | for i in range(ngens): |
| 228 | fac_min_exp[i] = min(fac_min_exp[i],m[i]) |
| 229 | fac_max_exp[i] = max(fac_max_exp[i],m[i]) |
| 230 | pd[polydict.ETuple(t)] = c |
| 231 | if must_refactor: |
| 232 | to_specialize = [(i,fac_max_exp[i]) for i in range(ngens) if fac_max_exp[i] > 0] |
| 233 | to_specialize.sort(key=lambda x:x[1]) |
| 234 | fac_results = _factor_kronecker_core(poly_factory(p.parent(),pd),fac_maxcoeff,to_specialize) |
| 235 | units *= fac_results.unit() |
| 236 | results += list(fac_results) |
| 237 | else: |
| 238 | results.append( (fac,e) ) |
| 239 | else: |
| 240 | results += list(kron_results) |
| 241 | return Factorization(results,unit=units) |