Ticket #12417: trac_12417-pfd.2.patch

File trac_12417-pfd.2.patch, 29.9 KB (added by araichev, 10 years ago)

Tested on Sage 4.8

  • sage/categories/quotient_fields.py

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