Ticket #12417: trac_12417-pfd.patch

File trac_12417-pfd.patch, 29.0 KB (added by malb, 10 years ago)
  • sage/categories/quotient_fields.py

    # HG changeset patch
    # User Alex Raichev <tortoise.said@gmail.com>
    # Date 1328736876 -46800
    # Node ID f8971d29e69f37b2d56e728a3ece3dea0af53ff0
    # Parent  e9fd586e05a45dbf67a1eefe8a1094e8be6fb8a1
    initial version; fails
    
    diff --git a/sage/categories/quotient_fields.py b/sage/categories/quotient_fields.py
    a b  
    244244            return (self.numerator().factor(*args, **kwds) /
    245245                    self.denominator().factor(*args, **kwds))
    246246
    247         def partial_fraction_decomposition(self):
     247        def partial_fraction_decomposition_univariate(self):
    248248            """
    249249            Decomposes fraction field element into a whole part and a list of
    250250            fraction field elements over prime power denominators.
     
    334334                parts.append(n/d)
    335335            return whole, parts
    336336
     337        def partial_fraction_decomposition(self, factored=False):
     338            r"""
     339            Returns a partial fraction decomposition of ``self``.
     340            See the note section below for more details.
     341
     342            INPUT:
     343
     344            - ``factored`` - (Optional; default=False) If True, then the
     345                denominators of the partial fractions will be factored.
     346                Note: factoring happens anyway as part of the algorithm.
     347
     348            EXAMPLES::
     349
     350            Extend Sage's old partial fraction decomposition function::
     351
     352                sage: f = 26/15
     353                sage: pfd = f.partial_fraction_decomposition()
     354                sage: print pfd
     355                (1, [1/3, 2/5])
     356                sage: print pfd[0] + sum(pfd[1]) == f
     357                True
     358
     359            Extend Sage's old partial fraction decomposition function::
     360
     361                sage: R.<x> = PolynomialRing(QQ)
     362                sage: f = x + 1/(x^4 - 4)
     363                sage: pfd = f.partial_fraction_decomposition()
     364                sage: print pfd
     365                (x, [1/4/(x^2 - 2), -1/4/(x^2 + 2)])
     366                sage: print pfd[0] + sum(pfd[1]) == f
     367                True
     368
     369            New stuff::
     370
     371                sage: R.<x,y>= PolynomialRing(QQ)
     372                sage: f = x + 1/(x*y*(x*y - 1))
     373                sage: pfd1 = f.partial_fraction_decomposition()
     374                sage: pfd2 = f.partial_fraction_decomposition(factored=True)
     375                sage: print pfd1
     376                (x, [1/(x*y - 1), (-1)/(x*y)])
     377                sage: print pfd2
     378                (x, [[1, [(x*y - 1, 1)]], [-1, [(y, 1), (x, 1)]]])
     379                sage: print pfd1[0] + sum(pfd1[1]) == f
     380                True
     381
     382            ::
     383
     384                sage: R.<x,y,z>= PolynomialRing(QQ)
     385                sage: f = x + (x^2 + y)/(x*y^2*(x*z + 1)^2*(y*z - 1))
     386                sage: pfd1 = f.partial_fraction_decomposition()
     387                sage: pfd2 = f.partial_fraction_decomposition(factored=True)
     388                sage: print pfd2
     389                (x, [[x^2*z^2 + y*z^2, [(x, 1), (y*z - 1, 1)]], [-x^2*y*z - y^2*z - x^2 - y, [(x, 1), (y, 2)]], [-x^3*z^4 - x*y*z^4 - 2*x^2*z^3 - 2*y*z^3, [(y*z - 1, 1), (x*z + 1, 2)]], [x^3*y*z^3 + x*y^2*z^3 + x^3*z^2 + 2*x^2*y*z^2 + x*y*z^2 + 2*y^2*z^2 + 2*x^2*z + 2*y*z, [(y, 2), (x*z + 1, 2)]]])
     390                sage: print pfd1[0] + sum(pfd1[1]) == f
     391                True
     392
     393
     394            TESTS:
     395
     396            Fails over inexact fields, since the call to quo_rem()
     397            fails over inexact fields::
     398
     399                sage: R.<x,y>= PolynomialRing(RealField(20))
     400                sage: f = x + 1/(x*y*(x*y - 1))
     401                sage: print f
     402                (x^3*y^2 - x^2*y + 1.0000)/(x^2*y^2 - x*y)
     403                sage: pfd = f.partial_fraction_decomposition()
     404                Traceback (most recent call last):
     405                ...
     406                KeyError: '1.000e+00'
     407
     408            Fails over rings whose coefficient rings are not fields::
     409
     410                sage: R.<x,y>= PolynomialRing(ZZ)
     411                sage: f = x + 1/(x*y*(x*y - 1))
     412                sage: pfd1 = f.partial_fraction_decomposition()
     413                sage: print pfd1
     414                sage: print pfd1[0] + sum(pfd1[1]) == f
     415                Traceback (most recent call last):
     416                ...
     417                NotImplementedError: Factorization of multivariate polynomials over non-fields is not implemented.
     418
     419            Works over finite fields::
     420
     421                sage: R.<x,y,z>= PolynomialRing(GF(2))
     422                sage: f = x + (x^2 + y)/(x*y^2*(x*z + 1)^2*(y*z - 1))
     423                sage: pfd1 = f.partial_fraction_decomposition()
     424                sage: pfd2 = f.partial_fraction_decomposition(factored=True)
     425                sage: print pfd2
     426                (x, [[x^2*z^2 + y*z^2, [(x, 1), (y*z + 1, 1)]], [x^2*y*z + y^2*z + x^2 + y, [(x, 1), (y, 2)]], [x^3*z^4 + x*y*z^4, [(y*z + 1, 1), (x*z + 1, 2)]], [x^3*y*z^3 + x*y^2*z^3 + x^3*z^2 + x*y*z^2, [(y, 2), (x*z + 1, 2)]]])
     427                sage: print pfd1[0] + sum(pfd1[1]) == f
     428                True
     429
     430            NOTE::
     431
     432            This function extends partial_fraction_decomposition_univariate() to
     433            work over any element of a fraction field of a multivariate
     434            polynomial ring over a field.
     435           
     436            What is a partial fraction decomposition in this case?
     437             
     438            Let `f = P/Q` be a rational expression where `P` and `Q` lie in a
     439            multivariate polynomial ring `R` over a field.
     440            Let `Q_1^{e_1} \cdots Q_m^{e_m}` be the
     441            unique factorization of `Q` in `R` into irreducible factors,
     442            and let `d` be the number of indeterminates of `R`.
     443            Then by Theorem 1 of [Lein1978]_, `f` can be written as a sum of
     444            rational expressions in the form `\sum P_A/\prod_{j\in A} Q_j^{b_j}`,
     445            where the `b_j \le e_j` are positive integers, the `P_A` are in `R`, and
     446            the sum is taken over all subsets `A \subseteq \{ 1, \ldots, d \}`
     447            such that `S:= \{ Q_j : j\in A \}` is algebraically independent and
     448            the ideal generated by `S` is not all of `R`
     449            (that is, by the Nullstellensatz, the polynomials of `S` have a
     450            common zero in `K^d`, where `K` is the algebraic closure of the
     451            coefficient field of `R`).
     452
     453            Following [Lein1978]_, i call any such decomposition of `f` a
     454            `\emph{partial fraction decomposition}`.
     455            Such a decomposition is not unique.
     456
     457            The algorithm herein comes from [Lein1978]_. 
     458            By the way, [Lein1978]_ has a typo in equation (c) on the
     459            second page and equation (b) on the third page. 
     460            The right sides of (c) and (b) should be multiplied by `P`.
     461
     462            REFERENCES:
     463
     464            .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
     465            several variables into partial fractions", Soviet Math. (Iz. VUZ)
     466            22 (1978), no. 10, 35--38.
     467
     468            AUTHORS:
     469
     470            - Alex Raichev (2012-02-02)
     471            """
     472            # Use old Sage partial fraction decomposition in univariate case.
     473            R = self.denominator().parent()
     474            if R.ngens() < 2:
     475                return self.partial_fraction_decomposition_univariate()
     476
     477            # Use this new Sage partial fraction decomposition in multivariate case.     
     478            Q = self.denominator()
     479            whole, P = self.numerator().quo_rem(Q)
     480            fractions = (P/Q)._to_refd().partial_fraction_decomposition()
     481            if factored:
     482                return whole, [t.list() for t in fractions]
     483            else:
     484                return whole, [t.quotient() for t in fractions]
     485       
     486        def _to_refd(self):
     487            r"""
     488            Converts ``self`` into a REFD object.
     489            Assumes that the total degree of the numerator of ``f`` is
     490            less than the total degree of the denominator of ``f``.
     491
     492            Used to format ``f`` in preparation for a partial fraction
     493            decomposition.
     494
     495            EXAMPLES::
     496
     497                sage: R.<x,y>= PolynomialRing(QQ)
     498                sage: f = (x + 1)/(x*y*(x*y+1)^2)
     499                sage: print f._to_refd()
     500                [x + 1, [(y, 1), (x, 1), (x*y + 1, 2)]]
     501
     502            AUTHORS:
     503
     504            - Alex Raichev (2012-02-02)
     505            """   
     506            P = self.numerator()
     507            Q = self.denominator()
     508            R = P.parent()
     509
     510            if Q in R.base_ring():
     511                return REFD(P/Q, [])
     512
     513            Qs= Q.factor(proof=False)   # Singular needs `proof=False` for finite fields.
     514            if Qs.unit() != 1:
     515                P = P/Qs.unit()
     516            Qs= [(q,e) for q,e in Qs] 
     517            return REFD(P, Qs)
     518                   
    337519        def derivative(self, *args):
    338520            r"""
    339521            The derivative of this rational function, with respect to variables
     
    469651            return self.__class__(R, num, den,
    470652                coerce=False, reduce=False)
    471653
     654class REFD(object):
     655    r"""
     656    Represents a rational expression with factored denominator (REFD)
     657    `P/(Q_1^e_1 \cdots Q_m^e_m)` by storing the
     658    parts `P` and `[(Q_1,e_1), \ldots, (Q_m,e_m)]`, where
     659    `P, Q_1, \ldots, Q_m` are elements of a common polynomial ring `R`,
     660    `Q_1, \ldots, Q_m` are irreducible elements of R, and
     661    `e_1, \ldots, e_m` are positive integers.
     662    A polynomial `P` is represented as `[P, (,)]`.
     663
     664    Used in partial_fraction_decomposition()
     665
     666    AUTHORS:
     667
     668    - Alex Raichev (2012-02-02)
     669    """
     670
     671    def __init__(self, numerator, denominators, reduce=True):
     672        r"""
     673        EXAMPLES::
     674
     675            sage: R.<x,y>= PolynomialRing(QQ)
     676            sage: Qs = [x,1], [y,1], [x*y+1,1]
     677            sage: f = REFD(x, Qs)
     678            sage: print f
     679            [1, [(y, 1), (x*y + 1, 1)]]
     680            sage: ff = REFD(x, Qs, reduce=False)
     681            sage: print ff
     682            [x, [(y, 1), (x, 1), (x*y + 1, 1)]]
     683
     684        ::
     685
     686            sage: R.<x,y>= PolynomialRing(QQ)
     687            sage: f = REFD(x + y, [(x + y,1)])
     688            sage: print f
     689            [1, []]
     690        """
     691        if denominators:
     692            R = denominators[0][0].parent()
     693        else:
     694            R = numerator.parent()
     695
     696        self._numerator = numerator
     697        self._denominators = sorted([tuple(t) for t in denominators])
     698        self._base_ring = R
     699
     700        if not reduce:
     701            return
     702
     703        # Reduce rational expression
     704        P = numerator
     705        Qs = denominators
     706        newP = P
     707        newQs = []
     708        for Q,e in Qs:
     709            a = 0
     710            while Q.divides(newP) and a < e:
     711                a += 1
     712                newP = R(newP/Q)
     713            if e - a > 0:
     714                newQs.append((Q,e-a))
     715
     716        self._numerator = newP
     717        self._denominators = newQs
     718        self._base_ring = R
     719
     720    def __str__(self):
     721        r"""
     722        EXAMPLES::
     723
     724            sage: R.<x,y>= PolynomialRing(QQ)
     725            sage: Qs = [x,1], [y,1], [x*y+1,1]
     726            sage: f = REFD(x, Qs)
     727            sage: print f
     728            [1, [(y, 1), (x*y + 1, 1)]]
     729        """
     730        return str(self.list())
     731
     732    def __eq__(self, other):
     733        r"""
     734        EXAMPLES::
     735
     736            sage: R.<x,y>= PolynomialRing(QQ)
     737            sage: Qs = [x,1], [y,1], [x*y+1,1]
     738            sage: f = REFD(x, Qs)
     739            sage: print f == ff
     740            True
     741            sage: g = REFD(y, Qs)
     742            sage: print g == f
     743            False
     744        """
     745        if isinstance(other, REFD) \
     746        and self.quotient() == other.quotient():
     747            return True
     748        else:
     749            return False
     750
     751    def __ne__(self, other):
     752        r"""
     753        EXAMPLES::
     754
     755            sage: R.<x,y>= PolynomialRing(QQ)
     756            sage: Qs = [x,1], [y,1], [x*y+1,1]
     757            sage: f = REFD(x, Qs)
     758            sage: print f != ff
     759            False
     760            sage: g = REFD(y, Qs)
     761            sage: print g != f
     762            True
     763        """
     764        if not self == other:
     765            return True
     766        else:
     767            return False
     768
     769    def __lt__(self, other):
     770        r"""
     771        EXAMPLES::
     772
     773            sage: R.<x,y>= PolynomialRing(QQ)
     774            sage: Qs = [x,1], [y,1], [x*y+1,1]
     775            sage: f = REFD(x, Qs)
     776            sage: ff = REFD(x, Qs, reduce=False)
     777            sage: g = REFD(x + 1, Qs)
     778            sage: h = REFD(x + 2, Qs)
     779            sage: print f < ff
     780            False
     781            sage: print f < g
     782            True
     783            sage: print g < f
     784            False
     785            sage: print g < h
     786            True
     787        """
     788        snum = self.numerator()
     789        onum = other.numerator()
     790        sdenoms = self.denominators()
     791        odenoms = other.denominators()
     792        sdenom = self.denominator()
     793        odenom = other.denominator()
     794
     795        if self == other:
     796            return False
     797        elif len(sdenoms) < len(odenoms) or \
     798        (len(sdenoms) == len(odenoms) and sdenom < odenom) or \
     799        (len(sdenoms) == len(odenoms) and sdenom == odenom and snum < onum):
     800            return True
     801        else:
     802            return False
     803
     804    def __gt__(self, other):
     805        r"""
     806        EXAMPLES::
     807
     808            sage: R.<x,y>= PolynomialRing(QQ)
     809            sage: Qs = [x,1], [y,1], [x*y+1,1]
     810            sage: f = REFD(x, Qs)
     811            sage: ff = REFD(x, Qs, reduce=False)
     812            sage: g = REFD(x + 1, Qs)
     813            sage: h = REFD(x + 2, Qs)
     814            sage: print f > ff
     815            False
     816            sage: print f > g
     817            False
     818            sage: print g > f
     819            True
     820            sage: print g > h
     821            False
     822        """
     823        if (not self == other) and (not self < other):
     824            return True
     825        else:
     826            return False
     827
     828    def __add__(self, other):
     829        r"""
     830        EXAMPLES::
     831
     832            sage: R.<x,y>= PolynomialRing(QQ)
     833            sage: Qs = [x,1], [y,1]
     834            sage: f = REFD(x, Qs)
     835            sage: ff = REFD(x, Qs, reduce=False)
     836            sage: print f + ff
     837            [2, [(y, 1), (x, 1)]]
     838            sage: print f + 2
     839            [2*y + 1, [(y, 1)]]
     840        """
     841        R = self.base_ring()
     842        P = self.numerator()
     843        Q = self.denominator()
     844        Qs = self.denominators()
     845
     846        if not isinstance(other, REFD):
     847            try:
     848                t = R(other)
     849            except:
     850                print "Fail! %s is not an element of %s" % (other, R)
     851            else:
     852                return REFD(P + R(other)*Q, Qs)
     853
     854        # Compute the sum's numerator
     855        num = (self.quotient() + other.quotient()).numerator()
     856
     857        # Compute the sum's denominators.       
     858        denoms = self.denominators() + other.denominators()
     859        if not denoms:
     860            # Then just added two polynomials. Done.
     861            return REFD(num, denoms, reduce=False)
     862
     863        denoms.sort()
     864        common_denoms = []
     865        current = list(denoms[0])
     866        for i in range(1,len(denoms)):
     867            if current[0] == denoms[i][0]:
     868                current[1] = max(current[1],denoms[i][1])
     869            else:
     870                # Append current to common_denoms and update current
     871                common_denoms.append(tuple(current))
     872                current = list(denoms[i])
     873        common_denoms.append(tuple(current))         
     874        return REFD(num, common_denoms, reduce=False)   
     875
     876    def __radd__(self, x):
     877        r"""
     878        EXAMPLES::
     879
     880            sage: R.<x,y>= PolynomialRing(QQ)
     881            sage: Qs = [x,1], [y,1]
     882            sage: f = REFD(x, Qs)
     883            sage: print 2 + f
     884            [2*y + 1, [(y, 1)]]       
     885        """
     886        return self + x
     887
     888    def list(self):
     889        return [self.numerator(), self.denominators()]
     890
     891    def numerator(self):
     892        return self._numerator
     893
     894    def denominators(self):
     895        return self._denominators
     896
     897    def denominator(self):
     898        from sage.misc.misc import prod
     899   
     900        return prod([q**e for q,e in self._denominators])
     901
     902    def base_ring(self):
     903        return self._base_ring
     904
     905    def quotient(self):
     906        return self.numerator()/self.denominator()
     907
     908    def partial_fraction_decomposition(self):
     909        r"""
     910        Return a partial fraction decomposition of ``self`` as a REFDSum object.
     911
     912        EXAMPLES::
     913
     914            sage: R.<x,y>= PolynomialRing(QQ)
     915            sage: f = (x^2 + y)/(x*y*(x*y + 1)^2)
     916            sage: pfd = f._to_refd().partial_fraction_decomposition()
     917            sage: print pfd
     918            [[-x^3*y - x*y^2 - 2*x^2 - 2*y, [(x*y + 1, 2)]], [x^2 + y, [(y, 1), (x, 1)]]]
     919            sage: print sum([t.quotient() for t in pfd]) == f
     920            True
     921
     922        ::
     923
     924            sage: R.<x,y,z>= PolynomialRing(QQ)
     925            sage: f = (x^2 + y)/(x*y^2*(x*z + 1)^2*(y*z - 1))
     926            sage: pfd = f._to_refd().partial_fraction_decomposition()
     927            sage: print pfd
     928            [[x^2*z^2 + y*z^2, [(x, 1), (y*z - 1, 1)]], [-x^2*y*z - y^2*z - x^2 - y, [(x, 1), (y, 2)]], [-x^3*z^4 - x*y*z^4 - 2*x^2*z^3 - 2*y*z^3, [(y*z - 1, 1), (x*z + 1, 2)]], [x^3*y*z^3 + x*y^2*z^3 + x^3*z^2 + 2*x^2*y*z^2 + x*y*z^2 + 2*y^2*z^2 + 2*x^2*z + 2*y*z, [(y, 2), (x*z + 1, 2)]]]
     929            sage: print sum([t.quotient() for t in pfd]) == f
     930            True
     931
     932        NOTE::
     933
     934        See the note section of the __main__ function
     935        partial_fraction_decomposition().
     936        """
     937        return REFDSum([self]).decompose_via_nullstellensatz()\
     938        .decompose_via_algebraic_dependence()
     939
     940class REFDSum(list):
     941    r"""
     942    A list representing the sum of REFD objects with distinct
     943    denominator lists.
     944
     945    Used in partial_fraction_decomposition()
     946
     947    AUTHORS:
     948
     949    - Alex Raichev (2012-02-02)
     950    """
     951
     952    def simplify(self):
     953        r"""
     954        Adds terms in ``self`` with the same denominator.
     955
     956        EXAMPLES::
     957
     958            sage: R.<x,y>= PolynomialRing(QQ)
     959            sage: f = REFD(1, [(x,1), (y,1), (x*y + 1,1)])
     960            sage: g = REFD(x, [(x,1), (y,1), (x*y + 1,1)])
     961            sage: s = REFDSum([f, g, f])
     962            sage: t = s.simplify()
     963            sage: print s
     964            [[1, [(x, 1), (y, 1), (x*y + 1, 1)]], [1, [(y, 1), (x*y + 1, 1)]], [1, [(x, 1), (y, 1), (x*y + 1, 1)]]]
     965            sage: print t
     966            [[1, [(y, 1), (x*y + 1, 1)]], [2, [(y, 1), (x, 1), (x*y + 1, 1)]]]
     967            sage: print sum(t) == sum(s)
     968            True
     969        """
     970        # Combine like terms.
     971        terms = sorted(self)
     972        new_terms = []
     973        total = terms[0]
     974        for i in range(len(terms) - 1):
     975            if  terms[i].denominators() == terms[i + 1].denominators():
     976                total += terms[i + 1]
     977            else:
     978                new_terms.append(total)
     979                total = terms[i + 1]
     980        new_terms.append(total)
     981
     982        return REFDSum(new_terms)
     983
     984    def __str__(self):
     985        return str(self.list())
     986
     987    def list(self):
     988        return [t.list() for t in self]
     989
     990    def base_ring(self):
     991        return self[0].base_ring()
     992
     993    def decompose_via_algebraic_dependence(self):
     994        r"""
     995        Returns a REFDSum decomposition of ``self`` according to
     996        Lemma 2 of [Lein1978]_.
     997        Recursive.
     998
     999        EXAMPLES::
     1000
     1001            sage: R.<x,y>= PolynomialRing(QQ)
     1002            sage: f = (x^2 + y)/(x*y*(x*y + 1)^2)
     1003            sage: pfd = REFDSum([f._to_refd()]).decompose_via_algebraic_dependence()
     1004            sage: print pfd
     1005            [[-x^3*y - x*y^2 - 2*x^2 - 2*y, [(x*y + 1, 2)]], [x^2 + y, [(y, 1), (x, 1)]]]
     1006            sage: print sum([t.quotient() for t in pfd]) == f
     1007            True
     1008
     1009        ::
     1010
     1011            sage: R.<x,y,z>= PolynomialRing(QQ)
     1012            sage: f = (x^2 + y)/(x*y^2*(x*z + 1)^2*(y*z - 1))
     1013            sage: pfd = REFDSum([f._to_refd()]).decompose_via_algebraic_dependence()
     1014            sage: print pfd
     1015            [[2*x^3*y*z + 2*x*y^2*z + 2*x^3 + 2*x*y, [(y, 4)]], [2*x^3 + 2*x*y, [(y, 4), (y*z - 1, 1)]], [2*x^3*y*z + 2*x*y^2*z + 2*x^3 + 2*x*y, [(x*z + 1, 2), (y, 4)]], [-x^5*y^3*z^3 - x^3*y^4*z^3 - x^5*y^2*z^2 - x^3*y^3*z^2 - x^5*y*z - x^3*y^2*z - x^5 - x^3*y, [(x*z + 1, 2), (y, 6)]], [-x^4*z^2 - x^2*y*z^2 - 2*x^3*z - 2*x*y*z + x^2 + y, [(y*z - 1, 1), (y, 2), (x, 1)]], [2*x^3 + 2*x*y, [(y, 4), (y*z - 1, 1), (x*z + 1, 2)]], [-x^5 - x^3*y, [(y, 6), (y*z - 1, 1), (x*z + 1, 2)]]]
     1016            sage: print sum([t.quotient() for t in pfd]) == f
     1017            True
     1018
     1019        NOTE::
     1020
     1021        Let `f = P/Q` be a rational expression where `P` and `Q` lie in a
     1022        (possibly multivariate) polynomial ring `R`.
     1023        Let `Q_1^{e_1} \cdots Q_m^{e_m}` be the
     1024        unique factorization of `Q` in `R` into irreducible factors,
     1025        and let `d` be the number of indeterminates of `R`.
     1026        If `Q_1, \ldots, Q_m` are algebraically dependent, then
     1027        by Lemma 2 of [Lein1978]_, `f` can be written as a sum of
     1028        rational expressions in the form `\sum P_A/\prod_{j\in A} Q_j^{b_j}`,
     1029        where the `b_j` are positive integers, the `P_A` are in `R`, and
     1030        the sum is taken over all size `< m` subsets
     1031        `A \subseteq \{ 1, \ldots, d \}` such that `S:= \{ Q_j : j\in A \}` is
     1032        algebraically independent.
     1033        In other words, `f` can be decomposed into a sum of rational
     1034        expressions, each of which has a denominator whose set of irreducible
     1035        factors is an algebraically independent subset of
     1036        `\{Q_1, \ldots, Q_m\}`.
     1037
     1038        This function decomposes every rational expression `f` in ``self``
     1039        as described above.
     1040
     1041        By the way, [Lein1978]_ has a typo in equation (c) on the
     1042        second page and equation (b) on the third page. 
     1043        The right sides of (c) and (b) should be multiplied by `P`.
     1044
     1045        REFERENCES:
     1046
     1047        .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
     1048        several variables into partial fractions", Soviet Math. (Iz. VUZ)
     1049        22 (1978), no. 10, 35--38.
     1050
     1051        AUTHORS:
     1052
     1053        - Alex Raichev (2011-01-10)   
     1054        """
     1055        from sage.structure.sequence import Sequence
     1056        from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     1057       
     1058        parts = REFDSum([])
     1059        R = self.base_ring()
     1060        done_decomposing = True
     1061        for t in self:
     1062            decomp = REFDSum([t])
     1063            P = t.numerator()
     1064            Qs = t.denominators()
     1065            m = len(Qs)
     1066            # The polynomials [q**e for q,e in Qs] are algebraically dependent
     1067            # iff the polynomials [q for q,e in Qs] are algebraically dependent.
     1068            G = Sequence([q for q,e in Qs], R).algebraic_dependence().gens()[0]
     1069            if G:
     1070                # Then [q**e for q,e in Qs] are algebraically dependent and
     1071                # hence we can decompose t.
     1072                done_decomposing = False
     1073
     1074                # Todo: speed up step below by using
     1075                # G to calculate F.  See [Lein1978]_ Lemma 1.
     1076                F = Sequence([q**e for q,e in Qs], R).algebraic_dependence().gens()[0]
     1077                new_vars = F.parent().gens()
     1078
     1079                # Note that each new_vars[j] corresponds to Qs[j] such that
     1080                # F([q**e for q,e in Qs]) = 0.
     1081                # Assume that F.parent() has negdeglex term order
     1082                # so that F.lt() is indeed the monomial we want below.
     1083                # Use F to rewrite P/Qs into a sum of REFDs,
     1084                # each with < m denominator factors.
     1085                FF = (F.lt() - F)/(F.lc()) 
     1086                numers = map(mul,zip(FF.coefficients(),FF.monomials()))
     1087                e = list(F.lt().exponents())[0:m]
     1088                denoms = [(new_vars[j], e[0][j] + 1) for j in range(m)]
     1089                decomp_temp = REFDSum([REFD(a, denoms) for a in numers]).simplify()
     1090                decomp = REFDSum()
     1091                # Substitute in Qs.
     1092                Qpowsub = dict([(new_vars[j],Qs[j][0]**Qs[j][1]) for j in range(m)])
     1093                for f in decomp_temp:
     1094                    gnum = P*F.parent()(f.numerator()).subs(Qpowsub)
     1095                    gdenoms = []
     1096                    for q,e in f.denominators():
     1097                        j = new_vars.index(q)
     1098                        gdenoms.append((Qs[j][0], Qs[j][1]*e))
     1099                        # if q in new_vars:
     1100                        #     j = new_vars.index(q)
     1101                        #     gdenoms.append((Qs[j][0], Qs[j][1]*e))
     1102                        # else:
     1103                        #     # Occasionally q is an integer.
     1104                        #     gdenoms.append((q,e))
     1105                    decomp.append(REFD(gnum,gdenoms))
     1106            parts.extend(decomp)
     1107        parts = parts.simplify()       
     1108        if done_decomposing:
     1109            return parts
     1110        else:     
     1111            return parts.decompose_via_algebraic_dependence()
     1112
     1113    def decompose_via_nullstellensatz(self):
     1114        r"""
     1115        Returns a REFDSum decomposition of ``self`` according to
     1116        Lemma 3 of [Lein1978]_.
     1117        Recursive.
     1118
     1119        EXAMPLES::
     1120
     1121            sage: R.<x,y>= PolynomialRing(QQ)
     1122            sage: f = (x^2 + y)/(x*y*(x*y + 1)^2)
     1123            sage: pfd = REFDSum([f._to_refd()]).decompose_via_nullstellensatz()
     1124            sage: print pfd
     1125            [[-x^3*y - x*y^2 - 2*x^2 - 2*y, [(x*y + 1, 2)]], [x^2 + y, [(y, 1), (x, 1)]]]
     1126            sage: print sum([t.quotient() for t in pfd]) == f
     1127            True
     1128
     1129        ::
     1130
     1131            sage: R.<x,y,z>= PolynomialRing(QQ)
     1132            sage: f = (x^2 + y)/(x*y^2*(x*z + 1)^2*(y*z - 1))
     1133            sage: pfd = REFDSum([f._to_refd()]).decompose_via_nullstellensatz()
     1134            sage: print pfd
     1135            [[x^2*z^2 + y*z^2, [(x, 1), (y*z - 1, 1)]], [-x^2*y*z - y^2*z - x^2 - y, [(x, 1), (y, 2)]], [-x^3*z^4 - x*y*z^4 - 2*x^2*z^3 - 2*y*z^3, [(y*z - 1, 1), (x*z + 1, 2)]], [x^3*y*z^3 + x*y^2*z^3 + x^3*z^2 + 2*x^2*y*z^2 + x*y*z^2 + 2*y^2*z^2 + 2*x^2*z + 2*y*z, [(y, 2), (x*z + 1, 2)]]]
     1136            sage: print sum([t.quotient() for t in pfd]) == f
     1137            True
     1138
     1139        NOTE::
     1140
     1141        Let `f = P/Q` be a rational expression where `P` and `Q` lie in a
     1142        (possibly multivariate) polynomial ring `R`.
     1143        Let `Q_1^{e_1} \cdots Q_m^{e_m}` be the
     1144        unique factorization of `Q` in `R` into irreducible factors,
     1145        and let `d` be the number of indeterminates of `R`.
     1146        If the ideal generated by `Q_1, \ldots, Q_m` is not all of `R`
     1147        (that is, by the Nullstellensatz, the polynomials have a
     1148        common zerozero in `K^d`, where `K` is the algebraic closure of the
     1149        coefficient field of `R`),
     1150        then by Lemma 3 of [Lein1978]_, `f` can be written as a sum of
     1151        rational expressions in the form `\sum P_A/\prod_{j\in A} Q_j^{b_j}`,
     1152        where the `b_j \le e_j` are positive integers, the `P_A` are in `R`, and
     1153        the sum is taken over all size `< m` subsets
     1154        `A \subseteq \{ 1, \ldots, d \}`.
     1155        In other words, `f` can be decomposed into a sum of rational
     1156        expressions, each of which has a denominator whose set of irreducible
     1157        factors is a size `< m` subset of `\{Q_1, \ldots, Q_m\}`.
     1158
     1159        This function decomposes every rational expression `f` in ``self``
     1160        as described above.
     1161
     1162        By the way, [Lein1978]_ has a typo in equation (c) on the
     1163        second page and equation (b) on the third page. 
     1164        The right sides of (c) and (b) should be multiplied by `P`.
     1165
     1166        REFERENCES:
     1167
     1168        .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
     1169        several variables into partial fractions", Soviet Math. (Iz. VUZ)
     1170        22 (1978), no. 10, 35--38.
     1171
     1172        AUTHORS:
     1173
     1174        - Alex Raichev (2011-01-10)
     1175        """
     1176       
     1177        parts = REFDSum([])
     1178        R= self.base_ring()
     1179        done_decomposing = True
     1180        for t in self:
     1181            decomp = REFDSum([t])
     1182            P= t.numerator()
     1183            Qs= t.denominators()
     1184            m= len(Qs)
     1185            if R(1) in R.ideal([q for q,e in Qs]):
     1186                # Then we can decompose t.
     1187                done_decomposing = False
     1188                L = R(1).lift(R.ideal([q**e for q,e in Qs]))
     1189                decomp = REFDSum([
     1190                REFD(P*L[i],[(Qs[j][0], Qs[j][1]) for j in range(m) if j != i])\
     1191                for i in range(m) if L[i] != 0])
     1192                # The procedure above yields no like terms to combine.
     1193            parts.extend(decomp)
     1194        parts = parts.simplify() 
     1195        if done_decomposing:
     1196            return parts
     1197        else:     
     1198            return parts.decompose_via_nullstellensatz()
  • sage/rings/polynomial/multi_polynomial_sequence.py

    diff --git a/sage/rings/polynomial/multi_polynomial_sequence.py b/sage/rings/polynomial/multi_polynomial_sequence.py
    a b  
    158158from sage.rings.quotient_ring_element import QuotientRingElement
    159159from sage.rings.polynomial.multi_polynomial_ideal import MPolynomialIdeal
    160160from sage.rings.polynomial.multi_polynomial import is_MPolynomial
     161from sage.structure.sequence import Sequence_generic
    161162from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    162163
    163164from sage.interfaces.singular import singular