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 b import polydict 
    5454import polydict
    5555
    5656from sage.structure.factorization import Factorization
     57from multi_polynomial_factor import factor_ZZx
    5758
    5859from sage.rings.polynomial.polynomial_singular_interface import Polynomial_singular_repr
    5960
    class MPolynomial_polydict(Polynomial_si 
    10731074            sage: k.factor()
    10741075            (s^2 + 2/3) * (x + s*y)^2 * (x + (-s)*y)^5 * (x^2 + s*x*y + s^2*y^2)^5
    10751076        """
     1077        if self.parent().base() is ZZ:
     1078            return factor_ZZx(self)
    10761079        # I do not think this applied anymore.  Or at least it's
    10771080        # more relevant to optimizing the NTL build.
    10781081        #\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
    - +  
     1"""
     2Factorization algorithm for multivariate polynomials over ZZ.
     3The algorithm uses Kronecker's Trick -- that is specialize a polynomial
     4in n variables by substituting a large (prime) integer to get a polynomial
     5in n-1 variables.  Iterate until we have a bi-variate or uni-variate and
     6factor directly.
     7
     8AUTHORS:
     9    -- Joel B. Mohler: Initial implementation
     10
     11TESTS:
     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################################################################################
     30import sage.rings.polynomial.polydict as polydict
     31from sage.rings.integer_ring import ZZ
     32from sage.rings.integer import Integer
     33from sage.structure.factorization import Factorization
     34from sage.rings.arith import gcd,next_prime
     35from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     36
     37def 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
     49def 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
     53def 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
     63def _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
     71def _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
     127def 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)