Ticket #12417: trac_12417-pfd.3.patch

File trac_12417-pfd.3.patch, 41.8 KB (added by araichev, 9 years ago)

Tested on Sage 5.0.1

  • sage/categories/quotient_fields.py

    # HG changeset patch
    # User Alex Raichev <alex.raichev@gmail.com>
    # Date 1340772096 -43200
    # Node ID 8b83d9f8c0a1dcfeb427f74e3ec20c8b7d17aa3e
    # Parent  01e0779db32182ef38b6bbeec105ba5f821a40a3
    #12417: Refactored quotient_fields.py so that all partial fraction decompositions go through the REFD class
    
    diff --git a/sage/categories/quotient_fields.py b/sage/categories/quotient_fields.py
    a b  
    11r"""
    22Quotient fields
     3
     4REFERENCES:
     5   
     6    .. [Raic2012]  Alexander Raichev. "Leinartas's partial fraction decomposition". `<http://arxiv.org/abs/1206.4740>`_.
    37"""
    48#*****************************************************************************
    59#  Copyright (C) 2008 Teresa Gomez-Diaz (CNRS) <Teresa.Gomez-Diaz@univ-mlv.fr>
     
    1216from sage.categories.category_singleton import Category_singleton
    1317from sage.misc.cachefunc import cached_method
    1418from sage.misc.abstract_method import abstract_method
     19from functools import total_ordering
    1520
    1621class QuotientFields(Category_singleton):
    1722    """
     
    5762            """
    5863            Greatest common divisor
    5964
    60             NOTE:
     65            NOTE::
    6166
    6267            In a field, the greatest common divisor is not very
    6368            informative, as it is only determined up to a unit. But in
     
    136141            """
    137142            Least common multiple
    138143
    139             NOTE:
     144            NOTE::
    140145
    141146            In a field, the least common multiple is not very
    142147            informative, as it is only determined up to a unit. But in
     
    245250            return (self.numerator().factor(*args, **kwds) /
    246251                    self.denominator().factor(*args, **kwds))
    247252
    248         def partial_fraction_decomposition(self):
    249             """
    250             Decomposes fraction field element into a whole part and a list of
    251             fraction field elements over prime power denominators.
    252            
    253             The sum will be equal to the original fraction.
    254            
    255             AUTHORS:
    256            
    257             - Robert Bradshaw (2007-05-31)
    258            
    259             EXAMPLES::
    260            
    261                 sage: S.<t> = QQ[]
    262                 sage: q = 1/(t+1) + 2/(t+2) + 3/(t-3); q
    263                 (6*t^2 + 4*t - 6)/(t^3 - 7*t - 6)
    264                 sage: whole, parts = q.partial_fraction_decomposition(); parts
    265                 [3/(t - 3), 1/(t + 1), 2/(t + 2)]
    266                 sage: sum(parts) == q
    267                 True
    268                 sage: q = 1/(t^3+1) + 2/(t^2+2) + 3/(t-3)^5
    269                 sage: whole, parts = q.partial_fraction_decomposition(); parts
    270                 [1/3/(t + 1), 3/(t^5 - 15*t^4 + 90*t^3 - 270*t^2 + 405*t - 243), (-1/3*t + 2/3)/(t^2 - t + 1), 2/(t^2 + 2)]
    271                 sage: sum(parts) == q
    272                 True
    273            
    274             We do the best we can over inexact fields::
    275            
    276                 sage: R.<x> = RealField(20)[]
    277                 sage: q = 1/(x^2 + x + 2)^2 + 1/(x-1); q
    278                 (x^4 + 2.0000*x^3 + 5.0000*x^2 + 5.0000*x + 3.0000)/(x^5 + x^4 + 3.0000*x^3 - x^2 - 4.0000)
    279                 sage: whole, parts = q.partial_fraction_decomposition(); parts
    280                 [1.0000/(x - 1.0000), 1.0000/(x^4 + 2.0000*x^3 + 5.0000*x^2 + 4.0000*x + 4.0000)]
    281                 sage: sum(parts)
    282                 (x^4 + 2.0000*x^3 + 5.0000*x^2 + 5.0000*x + 3.0000)/(x^5 + x^4 + 3.0000*x^3 - x^2 - 4.0000)
    283    
    284             TESTS:
    285 
    286             We test partial fraction for irreducible denominators::
    287            
    288                 sage: R.<x> = ZZ[]
    289                 sage: q = x^2/(x-1)
    290                 sage: q.partial_fraction_decomposition()
    291                 (x + 1, [1/(x - 1)])
    292                 sage: q = x^10/(x-1)^5
    293                 sage: whole, parts = q.partial_fraction_decomposition()
    294                 sage: whole + sum(parts) == q
    295                 True
    296 
    297             And also over finite fields (see trac #6052, #9945)::
    298            
    299                 sage: R.<x> = GF(2)[]
    300                 sage: q = (x+1)/(x^3+x+1)
    301                 sage: q.partial_fraction_decomposition()
    302                 (0, [(x + 1)/(x^3 + x + 1)])
    303 
    304                 sage: R.<x> = GF(11)[]
    305                 sage: q = x + 1 + 1/(x+1) + x^2/(x^3 + 2*x + 9)
    306                 sage: q.partial_fraction_decomposition()
    307                 (x + 1, [1/(x + 1), x^2/(x^3 + 2*x + 9)])
    308            
    309             And even the rationals::
    310            
    311                 sage: (26/15).partial_fraction_decomposition()
    312                 (1, [1/3, 2/5])
    313             """
    314             from sage.misc.misc import prod
    315             denom = self.denominator()
    316             whole, numer = self.numerator().quo_rem(denom)
    317             factors = denom.factor()
    318             if factors.unit() != 1:
    319                 numer *= ~factors.unit()
    320             if len(factors) == 1:
    321                 return whole, [numer/r**e for r,e in factors]
    322             if not self.parent().is_exact():
    323                 # factors not grouped in this case
    324                 all = {}
    325                 for r in factors: all[r[0]] = 0
    326                 for r in factors: all[r[0]] += r[1]
    327                 factors = all.items()
    328                 factors.sort() # for doctest consistency
    329             factors = [r**e for r,e in factors]
    330             parts = []
    331             for d in factors:
    332                 # note that the product below is non-empty, since the case
    333                 # of only one factor has been dealt with above
    334                 n = numer * prod([r for r in factors if r != d]).inverse_mod(d) % d # we know the inverse exists as the two are relatively prime
    335                 parts.append(n/d)
    336             return whole, parts
    337 
    338253        def derivative(self, *args):
    339254            r"""
    340255            The derivative of this rational function, with respect to variables
     
    470385            return self.__class__(R, num, den,
    471386                coerce=False, reduce=False)
    472387
     388        def partial_fraction_decomposition(self, kind=None, factor=True):
     389            r"""
     390            Return a partial fraction decomposition of ``self`` of kind
     391            ``kind``.
     392
     393            INPUT:
     394
     395            - ``kind`` - (Optional; default=None) The string
     396              'nullstellensatz' or 'algebraic_dependence', indicating what kind
     397              of multivariate decomposition to return in case ``self`` is
     398              multivariate.
     399              The default multivariate decomposition is a Leinartas
     400              decomposition.
     401            - ``factor`` - (Optional; default=True) If True, then the
     402              output will be a REFDSum instance (which features factored
     403              denominators).
     404              Note that factoring the denominator of ``self`` happens anyway as
     405              part of the algorithm.
     406
     407            OUTPUT:
     408
     409            A sum of fraction field elements or, if factor=True,
     410            a REFDSum instance that is a partial fraction decomposition of
     411            ``self``.
     412
     413            EXAMPLES::
     414
     415            No variables::
     416
     417                sage: f = 64/45
     418                sage: print f.partial_fraction_decomposition()
     419                [(1, []), (2, [(3, 2)]), (1, [(5, 1)])]
     420                sage: print f.partial_fraction_decomposition(factor=False)
     421                [1, 2/9, 1/5]
     422                sage: print sum(f.partial_fraction_decomposition(factor=False)) == f
     423                True
     424
     425            One variable::
     426
     427                sage: R.<x> = PolynomialRing(QQ)
     428                sage: f = 1 + 2/x + 2*x/(x^2 - 4*x + 8)
     429                sage: print f
     430                (x^3 + 16)/(x^3 - 4*x^2 + 8*x)
     431                sage: d1 = f.partial_fraction_decomposition()
     432                sage: print d1
     433                [(1, []), (2, [(x, 1)]), (2*x, [(x^2 - 4*x + 8, 1)])]
     434                sage: d1.sum().quotient() == f
     435                True
     436                sage: print f.partial_fraction_decomposition(factor=False)
     437                [1, 2/x, 2*x/(x^2 - 4*x + 8)]
     438
     439            Two variables::
     440
     441                sage: R.<x,y>= PolynomialRing(QQ)
     442                sage: f = 1/x + 1/y + 1/(x + y)
     443                sage: print f
     444                (x^2 + 3*x*y + y^2)/(x^2*y + x*y^2)
     445                sage: print f.partial_fraction_decomposition(kind='nullstellensatz')
     446                [(0, []), (x^2 + 3*x*y + y^2, [(y, 1), (x, 1), (x + y, 1)])]
     447                sage: print f.partial_fraction_decomposition(kind='algebraic_dependence')
     448                [(0, []), (x^2 + 3*x*y + y^2, [(y, 2), (x, 1)]), (-x^2 - 3*x*y - y^2, [(y, 2), (x + y, 1)])]
     449                sage: print f.partial_fraction_decomposition()
     450                [(0, []), (x^2 + 3*x*y + y^2, [(y, 2), (x, 1)]), (-x^2 - 3*x*y - y^2, [(y, 2), (x + y, 1)])]
     451
     452            Three variables over a finite field::
     453
     454                sage: R.<x,y,z>= PolynomialRing(GF(2))
     455                sage: f = x + 1/x + 1/y + 1/z + 1/(x*y + z)
     456                sage: print f
     457                (x^3*y^2*z + x^2*y*z^2 + x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2)/(x^2*y^2*z + x*y*z^2)
     458                sage: d1 = f.partial_fraction_decomposition()
     459                sage: print d1
     460                [(x, []), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (x*y + z, 1)]), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (y, 1), (x, 1)])]
     461                sage: d2 = f.partial_fraction_decomposition()
     462                sage: print d2
     463                [(x, []), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (x*y + z, 1)]), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (y, 1), (x, 1)])]
     464                sage: print d2.sum().quotient() == f
     465                True
     466
     467            NOTE::
     468
     469            See the docstrings of the REFD class methods
     470            leinartas_decomposition(),
     471            nullstellensatz_decomposition(), and
     472            algebraic_dependence_decomposition()
     473            for definitions of the different kinds of partial fraction
     474            decomposition options available.
     475
     476            AUTHORS:
     477
     478            - Robert Bradshaw (2007-05-31)
     479            - Alexander Raichev (2012-06-25)
     480            """
     481            f = REFD(quotient=self)
     482            if f.ring().ngens() <= 1:
     483                decomp = f.univariate_decomposition()
     484            else:
     485                if kind == 'nullstellensatz':
     486                    decomp = f.nullstellensatz_decomposition()
     487                elif kind == 'algebraic_dependence':
     488                    decomp = f.algebraic_dependence_decomposition()
     489                else:
     490                    decomp = f.leinartas_decomposition()
     491            if factor:
     492                return decomp
     493            else:
     494                return [r.quotient() for r in decomp]
     495
     496
     497@total_ordering
     498class REFD(object):
     499    r"""
     500    Represents a rational expression with factored denominator (REFD)
     501    `p/(q_1^{e_1} \cdots q_m^{e_m})` by storing the parts `p` and
     502    `[(q_1,e_1), \ldots, (q_m,e_m)]`.
     503    Here `p, q_1, \ldots, q_m` are elements of a 0- or 1-variate factorial
     504    polynomial ring `R`, `q_1, \ldots, q_m` are distinct irreducible elements
     505    of `R`, and `e_1, \ldots, e_m` are positive integers.
     506    An element `r` of `R` is represented as `[r, (,)]`.
     507
     508    AUTHORS:
     509
     510    - Alexander Raichev (2012-06-25)
     511    """
     512
     513    def __init__(self, numerator=None, denominator_list=None, quotient=None,
     514                 reduce_=True):
     515        r"""
     516        Create a REFD instance.
     517
     518        INPUT:
     519
     520        - ``numerator`` - (Optional; default=None) An element `p` of a
     521          0- or 1-variate factorial polynomial ring `R`.
     522        - ``denominator_list`` - (Optional; default=None) A list of the form
     523          `[(q_1,e_1), \ldots, (q_m,e_m)]` where the `q_1, \ldots, q_m` are
     524          distinct irreducible elements of `R` and the `e_i` are positive
     525          integers.
     526        - ``quotient`` - (Optional; default=None) An element of a field of
     527          fractions of a factorial ring.
     528        - ``reduce_`` - (Optional; default=True) If True, then represent
     529          `p/(q_1^{e_1} \cdots q_m^{e_m})` in lowest terms.
     530          If False, then won't attempt to divide `p` by any of the `q_i`.
     531
     532        OUTPUT:
     533
     534        A REFD instance representing the rational expression
     535        `p/(q_1^{e_1} \cdots q_m^{e_m})`.
     536        To get a non-None output, one of ``numerator`` or ``quotient`` must be
     537        non-None.
     538
     539        EXAMPLES::
     540
     541            sage: from sage.categories.quotient_fields import REFD
     542            sage: f = REFD(64, [(3, 2), (5, 1)])
     543            sage: print f
     544            (64, [(3, 2), (5, 1)])
     545
     546        ::
     547       
     548            sage: R.<x,y> = PolynomialRing(QQ)
     549            sage: qs = [x, 1], [y, 1], [x*y+1, 1]
     550            sage: f = REFD(x, qs)
     551            sage: print f
     552            (1, [(y, 1), (x*y + 1, 1)])
     553            sage: ff = REFD(x, qs, reduce_=False)
     554            sage: print ff
     555            (x, [(y, 1), (x, 1), (x*y + 1, 1)])
     556
     557        ::
     558       
     559            sage: f = REFD(x + y, [(x + y, 1)])
     560            sage: print f
     561            (1, [])
     562
     563        ::
     564       
     565            sage: f = 64/45
     566            sage: print REFD(quotient=f)
     567            (64, [(3, 2), (5, 1)])
     568
     569        ::
     570       
     571            sage: R.<x> = PolynomialRing(QQ)
     572            sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
     573            sage: print REFD(quotient=f)
     574            (5*x^7 - 5*x^6 + 5/3*x^5 - 5/3*x^4 + 2*x^3 - 2/3*x^2 + 1/3*x - 1/3, [(x - 1, 1), (x, 1), (x^2 + 1/3, 1)])
     575
     576        ::
     577
     578            sage: R.<x,y> = PolynomialRing(QQ)
     579            sage: f = 2*y/(5*(x^3 - 1)*(y + 1))
     580            sage: print REFD(quotient=f)
     581            (2/5*y, [(y + 1, 1), (x - 1, 1), (x^2 + x + 1, 1)])
     582
     583        Singular throws a 'not implemented' error when trying to factor in
     584        a multivariate polynomial ring over an inexact field ::
     585
     586            sage: R.<x,y>= PolynomialRing(CC)
     587            sage: f = (x + 1)/(x*y*(x*y + 1)^2)
     588            sage: REFD(quotient=f)
     589            Traceback (most recent call last):
     590            ...
     591            TypeError: Singular error:
     592               ? not implemented
     593               ? error occurred in or before STDIN line 17: `def sage9=factorize(sage8);`
     594
     595        """
     596        from sage.rings.polynomial.polynomial_ring import is_PolynomialRing
     597        from sage.rings.polynomial.multi_polynomial import is_MPolynomial
     598       
     599        if quotient is not None:
     600            p = quotient.numerator()
     601            q = quotient.denominator()
     602            R = p.parent()
     603            self._ring = R
     604            if (is_PolynomialRing(R) or is_MPolynomial(R)) and \
     605               (q in R.base_ring()):
     606                # R is a polynomial ring and q is in its base ring.
     607                # No factoring needed.
     608                self._numerator = quotient
     609                self._denominator_list = [] 
     610            else:
     611                # Otherwise factor q.
     612                try:
     613                    qs = q.factor()
     614                except NotImplementedError:
     615                    # Singular's factor() needs `proof=False`.
     616                    qs = q.factor(proof=False)   
     617                if qs.unit() != 1:
     618                    p = p/qs.unit()
     619                qs = sorted([tuple(t) for t in qs]) # Sort for consitency.
     620                self._numerator = p
     621                self._denominator_list = qs 
     622            # Done. No reducing needed as Sage automatically reduces fractions.
     623            return
     624
     625        if numerator is not None:
     626            R = numerator.parent()     
     627        else:
     628            R = None
     629        self._numerator = numerator
     630        self._denominator_list = sorted([tuple(t) for t in denominator_list])
     631        self._ring = R
     632
     633        if (numerator is not None) and reduce_:
     634            # Reduce rational expression if possible by canceling
     635            # irreducible factors in denominator_list.
     636            p = numerator
     637            qs = denominator_list
     638            p_new = p
     639            qs_new = []
     640            for (q, e) in qs:
     641                a = 0
     642                while q.divides(p_new) and a < e:
     643                    a += 1
     644                    p_new = R(p_new/q)
     645                if e - a > 0:
     646                    qs_new.append((q,e-a))
     647            self._numerator = p_new
     648            self._denominator_list = qs_new
     649
     650    def numerator(self):
     651        r"""
     652        Return the numerator of ``self``.
     653        """
     654        return self._numerator
     655
     656    def denominator_list(self):
     657        r"""
     658        Return the denominator list of ``self``.
     659        """
     660        return self._denominator_list
     661
     662    def denominator(self):
     663        r"""
     664        Return the denominator of ``self``.
     665        """
     666        from sage.misc.misc import prod
     667        return prod([q**e for q,e in self.denominator_list()])
     668
     669    def ring(self):
     670        r"""
     671        Return the ring of ``self``.
     672        """
     673        return self._ring
     674
     675    def list(self):
     676        r"""
     677        Convert ``self`` into a list.
     678        """
     679        return (self.numerator(), self.denominator_list())
     680
     681    def quotient(self):
     682        r"""
     683        Convert ``self`` into a field of fractions element.
     684        """
     685        return self.numerator()/self.denominator()
     686
     687    def __str__(self):
     688        return str(self.list())
     689
     690    def __eq__(self, other):
     691        r"""
     692        Two REFD instances are equal iff they represent the same
     693        rational expression.
     694
     695        EXAMPLES::
     696
     697            sage: from sage.categories.quotient_fields import REFD
     698            sage: R.<x,y>= PolynomialRing(QQ)
     699            sage: qs = [x,1], [y,1], [x*y+1,1]
     700            sage: f = REFD(x, qs)
     701            sage: ff = REFD(x, qs, reduce_=False)
     702            sage: f == ff
     703            True
     704            sage: g = REFD(y, qs)
     705            sage: g == f
     706            False
     707        """
     708        return self.quotient() == other.quotient()
     709
     710    def __ne__(self, other):
     711        r"""
     712        EXAMPLES::
     713
     714            sage: from sage.categories.quotient_fields import REFD
     715            sage: R.<x,y>= PolynomialRing(QQ)
     716            sage: qs = [x,1], [y,1], [x*y+1,1]
     717            sage: f = REFD(x, qs)
     718            sage: ff = REFD(x, qs, reduce_=False)
     719            sage: f != ff
     720            False
     721            sage: g = REFD(y, qs)
     722            sage: g != f
     723            True
     724        """
     725        return not (self == other)
     726
     727    def __lt__(self, other):
     728        r"""
     729        REFD A is less than REFD B iff the denominator list of A
     730        is shorter than that of B, or
     731        the denominator list lengths are equal and
     732        the denominator of A is less than that of B in the ambient ring, or
     733        the denominator list lengths are equal and the denominators are equal
     734        and the the numerator of A is less than that of B in the ambient ring.
     735
     736        EXAMPLES::
     737
     738            sage: from sage.categories.quotient_fields import REFD
     739            sage: R.<x,y>= PolynomialRing(QQ)
     740            sage: qs = [x,1], [y,1], [x*y+1,1]
     741            sage: f = REFD(x, qs)
     742            sage: ff = REFD(x, qs, reduce_=False)
     743            sage: g = REFD(x + 1, qs)
     744            sage: h = REFD(x + 2, qs)
     745            sage: print f < ff
     746            False
     747            sage: print f < g
     748            True
     749            sage: print g < f
     750            False
     751            sage: print g < h
     752            True
     753        """
     754        sn = self.numerator()
     755        on = other.numerator()
     756        sdl = self.denominator_list()
     757        odl = other.denominator_list()
     758        sd = self.denominator()
     759        od = other.denominator()
     760
     761        if self == other:
     762            return False
     763        elif len(sdl) < len(odl) or\
     764            (len(sdl) == len(odl) and sd < od) or\
     765            (len(sdl) == len(odl) and sd == od and sn < on):
     766            return True
     767        else:
     768            return False     
     769
     770    def univariate_decomposition(self):
     771        r"""
     772        Return the usual 0- or 1-variate partial fraction decomposition
     773        of ``self`` as a REFDSum instance.
     774        Assume that ``self`` lies in the field of fractions
     775        of a 0- or 1-variate factorial polynomial ring.
     776
     777        EXAMPLES::
     778
     779        No variables::
     780
     781            sage: from sage.categories.quotient_fields import REFD
     782            sage: f = 64/45
     783            sage: decomp = REFD(quotient=f).univariate_decomposition()
     784            sage: print decomp
     785            [(1, []), (2, [(3, 2)]), (1, [(5, 1)])]
     786            sage: print decomp.sum().quotient() == f
     787            True
     788
     789        One variable::
     790
     791            sage: R.<x> = PolynomialRing(QQ)
     792            sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
     793            sage: print f
     794            (15*x^7 - 15*x^6 + 5*x^5 - 5*x^4 + 6*x^3 - 2*x^2 + x - 1)/(3*x^4 - 3*x^3 + x^2 - x)
     795            sage: decomp = REFD(quotient=f).univariate_decomposition()
     796            sage: print decomp
     797            [(5*x^3, []), (1, [(x - 1, 1)]), (1, [(x, 1)]), (1/3, [(x^2 + 1/3, 1)])]
     798            sage: print decomp.sum().quotient() == f
     799            True
     800
     801        One variable over a finite field::
     802
     803            sage: R.<x> = PolynomialRing(GF(2))
     804            sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
     805            sage: print f
     806            (x^6 + x^4 + 1)/(x^3 + x)
     807            sage: decomp = REFD(quotient=f).univariate_decomposition()
     808            sage: print decomp
     809            [(x^3, []), (1, [(x, 1)]), (x, [(x + 1, 2)])]
     810            sage: print decomp.sum().quotient() == f
     811            True
     812
     813        One variable over an inexact field::
     814
     815            sage: R.<x> = PolynomialRing(CC)
     816            sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
     817            sage: print f
     818            (15.0000000000000*x^7 - 15.0000000000000*x^6 + 5.00000000000000*x^5 - 5.00000000000000*x^4 + 6.00000000000000*x^3 - 2.00000000000000*x^2 + x - 1.00000000000000)/(3.00000000000000*x^4 - 3.00000000000000*x^3 + x^2 - x)
     819            sage: decomp = REFD(quotient=f).univariate_decomposition()
     820            sage: print decomp
     821            [(5.00000000000000*x^3, []), (1.00000000000000, [(x - 1.00000000000000, 1)]), (-0.288675134594813*I, [(x - 0.577350269189626*I, 1)]), (1.00000000000000, [(x, 1)]), (0.288675134594813*I, [(x + 0.577350269189626*I, 1)])]
     822            sage: print decomp.sum().quotient() == f # Rounding error coming
     823            False
     824
     825        NOTE::
     826
     827        Let `f = p/q` be a rational expression where `p` and `q` lie in a
     828        0- or 1-variate factorial polynomial ring `R`.
     829        Let `q_1^{e_1} \cdots q_m^{e_m}` be the
     830        unique factorization of `q` in `R` into irreducible factors.
     831        Then `f` can be written uniquely as
     832
     833        (*)   `p_0 + \sum_{i=1}^{m} p_i/ q_i^{e_i}`,
     834
     835        for some `p_j \in R`.
     836        I call (*) the *usual partial fraction decomposition* of f.
     837
     838        AUTHORS:
     839
     840        - Robert Bradshaw (2007-05-31)
     841        - Alexander Raichev (2012-06-25)
     842        """
     843        from sage.misc.misc import prod
     844       
     845        p = self.numerator()
     846        q = self.denominator()
     847        whole, p = p.quo_rem(q)
     848        qs = self.denominator_list()
     849        decomp = [REFD(whole, [])]
     850        for (a, m) in qs:
     851            numer = p * prod([b**n for (b, n) in qs if b != a])\
     852                    .inverse_mod(a**m) % (a**m) 
     853            # The inverse exists because the product and a**m
     854            # are relatively prime.
     855            decomp.append(REFD(numer, [(a, m)]))
     856        return REFDSum(decomp)
     857
     858    def nullstellensatz_certificate(self):
     859        r"""
     860        Let `[(q_1, e_1), \ldots, (q_m,e_m)]` be the denominator list of   
     861        ``self``.
     862        Return a list of polynomials `h_1, \ldots, h_m` in ``self.ring()``
     863        that satisfies `h_1 q_1 + \cdots + h_m q_m = 1` if it exists.
     864        Otherwise return None.
     865        Only works for multivariate ``self``.
     866
     867        EXAMPLES::
     868
     869            sage: from sage.categories.quotient_fields import REFD
     870            sage: R.<x,y> = PolynomialRing(QQ)
     871            sage: f = 1/(x^2 * (x*y + 1))
     872            sage: print f
     873            1/(x^3*y + x^2)
     874            sage: ff = REFD(quotient=f)           
     875            sage: L = ff.nullstellensatz_certificate()
     876            sage: print L
     877            [y^2, -x*y + 1]
     878            sage: qs = ff.denominator_list()
     879            sage: sum([L[i]*qs[i][0]**qs[i][1] for i in range(len(qs))]) == 1
     880            True
     881
     882        ::
     883
     884            sage: f = 1/(x*y)
     885            sage: L = REFD(quotient=f).nullstellensatz_certificate()
     886            sage: L is None
     887            True
     888
     889        """
     890
     891        R = self.ring()
     892        qs = self.denominator_list()
     893        J = R.ideal([q**e for q, e in qs])
     894        if R(1) in J:
     895            return R(1).lift(J)
     896        else:
     897            return None
     898
     899    def nullstellensatz_decomposition(self):
     900        r"""
     901        Return a Nullstellensatz decomposition of ``self`` as a REFDSum
     902        instance.
     903        Recursive.
     904        Only works for multivariate ``self``.
     905
     906        EXAMPLES::
     907
     908            sage: from sage.categories.quotient_fields import REFD
     909            sage: R.<x,y> = PolynomialRing(QQ)
     910            sage: f = 1/(x*(x*y + 1))
     911            sage: decomp = REFD(quotient=f).nullstellensatz_decomposition()
     912            sage: print decomp
     913            [(0, []), (1, [(x, 1)]), (-y, [(x*y + 1, 1)])]
     914            sage: decomp.sum().quotient() == f
     915            True
     916            sage: for r in decomp:
     917            ...       L = r.nullstellensatz_certificate()
     918            ...       L is None
     919            True
     920            True
     921            True
     922
     923        NOTE::
     924
     925        Let `f = p/q` be a rational expression where `p` and `q` lie in a
     926        `d`-variate polynomial ring `K[X]` for some field `K` and `d \ge 1`.
     927        Let `q_1^{e_1} \cdots q_m^{e_m}` be the
     928        unique factorization of `q` in `K[X]` into irreducible factors and
     929        let `V_i` be the algebraic variety `\{x\in L^d: q_i(x) = 0\}` of `q_i`
     930         over the algebraic closure `L` of `K`.
     931        By [Raic2012]_, `f` can be written as
     932
     933        (*)   `\sum p_A/\prod_{i \in A} q_i^{e_i}`,
     934
     935        where the `p_A` are in `K[X]` and the sum is taken over all subsets
     936        `A \subseteq \{1, \ldots, m\}` such that
     937        `\cap_{i\in A} T_i \neq \emptyset`.
     938
     939        I call (*) a *Nullstellensatz decomposition* of `f`.
     940        Nullstellensatz decompositions are not unique.
     941
     942        The algorithm used comes from [Raic2012]_.       
     943        """       
     944        decomp = REFDSum()
     945        L = self.nullstellensatz_certificate()
     946        if L is None:
     947            # No decomposing possible.
     948            decomp = REFDSum([self])
     949        else:
     950            # Decompose recursively.
     951            p = self.numerator()
     952            qs = self.denominator_list()           
     953            m = len(qs)
     954            iteration1 = REFDSum([REFD(p*L[i],[(qs[j][0], qs[j][1])
     955                                  for j in range(m) if j != i])
     956                                  for i in range(m) if L[i] != 0])
     957
     958            # Now decompose each REFD of iteration1.
     959            for r in iteration1:
     960                decomp.extend(r.nullstellensatz_decomposition())
     961
     962        # Simplify and return result.
     963        return decomp.combine_like_terms().whole_and_parts()
     964
     965    def algebraic_dependence_certificate(self):
     966        r"""
     967        Return the ideal `J` of annihilating polynomials for the set
     968        of polynomials ``[q**e for (q,e) in self.denominator_list()]``,
     969        which could be the zero ideal.
     970        The ideal `J` lies in a polynomial ring over the field
     971        ``self.ring().base_ring()`` that has
     972        ``m = len(self.denominator_list())`` indeterminates.
     973
     974        EXAMPLES:
     975
     976            sage: from sage.categories.quotient_fields import REFD
     977            sage: R.<x,y> = PolynomialRing(QQ)
     978            sage: f = 1/(x^2 * (x*y + 1) * y^3)
     979            sage: ff = REFD(quotient=f)
     980            sage: J = ff.algebraic_dependence_certificate()
     981            sage: print J
     982            Ideal (1 - 6*T2 + 15*T2^2 - 20*T2^3 + 15*T2^4 - T0^2*T1^3 - 6*T2^5 + T2^6) of Multivariate Polynomial Ring in T0, T1, T2 over Rational Field
     983            sage: g = J.gens()[0]
     984            sage: qs = ff.denominator_list()
     985            sage: g(*(q**e for q, e in qs)) == 0
     986            True
     987
     988        ::
     989
     990            sage: f = 1/(x^3 * y^2)
     991            sage: J = REFD(quotient=f).algebraic_dependence_certificate()
     992            sage: print J
     993            Ideal (0) of Multivariate Polynomial Ring in T0, T1 over Rational Field
     994
     995        """
     996        from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     997       
     998        R = self.ring()
     999        qs = self.denominator_list()       
     1000        if not qs:
     1001            return R.ideal()    # The zero ideal.
     1002        m = len(qs)
     1003        F = R.base_ring()
     1004        Xs = list(R.gens())
     1005        d = len(Xs)
     1006
     1007        # Expand R by 2*m new variables.
     1008        S = 'S'
     1009        while S in [str(x) for x in Xs]:
     1010            S = S + 'S'
     1011        Ss = [S + str(i) for i in range(m)]
     1012        T = 'T'
     1013        while T in [str(x) for x in Xs]:
     1014            T = T + 'T'
     1015        Ts = [T + str(i) for i in range(m)]
     1016
     1017        Vs = [str(x) for x in Xs] + Ss + Ts
     1018        RR = PolynomialRing(F, Vs)
     1019        Xs = RR.gens()[:d]
     1020        Ss = RR.gens()[d: d + m]
     1021        Ts = RR.gens()[d + m: d + 2 * m]
     1022
     1023        # Compute the appropriate elimination ideal.
     1024        J = RR.ideal([ Ss[j] - RR(qs[j][0]) for j in range(m)] +\
     1025                 [ Ss[j]**qs[j][1] - Ts[j] for j in range(m)])
     1026        J = J.elimination_ideal(Xs + Ss)
     1027
     1028        # Coerce J into the polynomial ring in the indeteminates Ts[m:].
     1029        # I choose the negdeglex order because i find it useful in my work.
     1030        RRR = PolynomialRing(F, [str(t) for t in Ts], order='negdeglex')
     1031        return RRR.ideal(J)
     1032
     1033    def algebraic_dependence_decomposition(self, whole_and_parts=True):
     1034        r"""
     1035        Return an algebraic dependence decomposition of ``self`` as a REFDSum
     1036        instance.
     1037        Recursive.
     1038
     1039        EXAMPLES::
     1040
     1041            sage: from sage.categories.quotient_fields import REFD
     1042            sage: R.<x,y> = PolynomialRing(QQ)
     1043            sage: f = 1/(x^2 * (x*y + 1) * y^3)
     1044            sage: decomp = REFD(quotient=f).algebraic_dependence_decomposition()
     1045            sage: print decomp
     1046            [(0, []), (-x, [(x*y + 1, 1)]), (x^2*y^2 - x*y + 1, [(y, 3), (x, 2)])]
     1047            sage: print decomp.sum().quotient() == f
     1048            True
     1049            sage: for r in decomp:
     1050            ...       J = r.algebraic_dependence_certificate()
     1051            ...       Z = J.ring().ideal()    # The zero ideal
     1052            ...       J == Z
     1053            True
     1054            True
     1055            True
     1056
     1057        NOTE::
     1058
     1059        Let `f = p/q` be a rational expression where `p` and `q` lie in a
     1060        `d`-variate polynomial ring `K[X]` for some field `K`.
     1061        Let `q_1^{e_1} \cdots q_m^{e_m}` be the
     1062        unique factorization of `q` in `K[X]` into irreducible factors and
     1063        let `V_i` be the algebraic variety `\{x\in L^d: q_i(x) = 0\}` of `q_i`
     1064         over the algebraic closure `L` of `K`.
     1065        By [Raic2012]_, `f` can be written as
     1066
     1067        (*)   `\sum p_A/\prod_{i \in A} q_i^{b_i}`,
     1068
     1069        where the `b_i` are positive integers, the `p_A` are in `K[X]`,
     1070        and the sum is taken over all subsets `A \subseteq \{1, \ldots, m\}`
     1071        such that `|A| \le d` and `\{q_i : i\in A\}` is algebraically
     1072        independent.
     1073
     1074        I call (*) an *algebraic dependence decomposition* of `f`.
     1075        Algebraic dependence decompositions are not unique.
     1076
     1077        The algorithm used comes from [Raic2012]_.
     1078        """
     1079        from sage.misc.misc import prod
     1080       
     1081        decomp = REFDSum()
     1082        R = self.ring()
     1083        J = self.algebraic_dependence_certificate()
     1084        if not J:
     1085            # No decomposing possible.
     1086            decomp = REFDSum([self])
     1087        else:
     1088            # Decompose recursively.
     1089            p = self.numerator()
     1090            qs = self.denominator_list()
     1091            m = len(qs)
     1092            g = J.gens()[0]     # An annihilating polynomial for the qs.
     1093            new_vars = J.ring().gens()
     1094            # Note that each new_vars[j] corresponds to qs[j] such that
     1095            # g([q**e for q,e in qs]) = 0.
     1096            # Assuming here that g.parent() has negdeglex term order
     1097            # so that g.lt() is indeed the monomial we want below.
     1098            # Use g to rewrite r into a sum of REFDs,
     1099            # each with < m distinct denominator factors.
     1100            gg = (g.lt() - g)/(g.lc()) 
     1101            numers = map(prod, zip(gg.coefficients(), gg.monomials()))
     1102            e = list(g.lt().exponents())[0:m]
     1103            denoms = [(new_vars[j], e[0][j] + 1) for j in range(m)]
     1104            # Write r in terms of new_vars,
     1105            # cancel factors in the denominator, and combine like terms.
     1106            iteration1_temp = REFDSum([REFD(a, denoms) for a in numers]).\
     1107                              combine_like_terms()
     1108            # Substitute in qs.
     1109            qpowsub = dict([(new_vars[j],qs[j][0]**qs[j][1])
     1110                            for j in range(m)])
     1111            iteration1 = REFDSum()
     1112            for r in iteration1_temp:
     1113                num1 = p*g.parent()(r.numerator()).subs(qpowsub)
     1114                denoms1 = []
     1115                for q,e in r.denominator_list():
     1116                    j = new_vars.index(q)
     1117                    denoms1.append((qs[j][0], qs[j][1]*e))
     1118                iteration1.append(REFD(num1, denoms1))
     1119            # Now decompose each REFD of iteration1.
     1120            for r in iteration1:
     1121                decomp.extend(r.algebraic_dependence_decomposition())
     1122
     1123        # Simplify and return result.
     1124        return decomp.combine_like_terms().whole_and_parts()
     1125
     1126    def leinartas_decomposition(self):
     1127        r"""
     1128        Return a Leinartas decomposition of ``self`` as a REFDSum instance.
     1129
     1130        EXAMPLES::
     1131
     1132            sage: from sage.categories.quotient_fields import REFD
     1133            sage: R.<x,y> = PolynomialRing(QQ)
     1134            sage: f = 1/x + 1/y + 1/(x*y + 1)
     1135            sage: decomp = REFD(quotient=f).leinartas_decomposition()
     1136            sage: print decomp
     1137            [(0, []), (1, [(x*y + 1, 1)]), (x + y, [(y, 1), (x, 1)])]
     1138            sage: print decomp.sum().quotient() == f
     1139            True
     1140            sage: for r in decomp:
     1141            ...       L = r.nullstellensatz_certificate()
     1142            ...       print L is None
     1143            ...       J = r.algebraic_dependence_certificate()
     1144            ...       Z = J.ring().ideal()    # The zero ideal
     1145            ...       print J == Z
     1146            True
     1147            True
     1148            True
     1149            True
     1150            True
     1151            True
     1152
     1153        ::
     1154       
     1155            sage: R.<x,y,z>= PolynomialRing(GF(2, 'a'))
     1156            sage: f = 1/(x * y * z * (x*y + z))
     1157            sage: decomp = REFD(quotient=f).leinartas_decomposition()
     1158            sage: print decomp
     1159            [(0, []), (1, [(z, 2), (x*y + z, 1)]), (1, [(z, 2), (y, 1), (x, 1)])]
     1160            sage: print decomp.sum().quotient() == f
     1161            True
     1162
     1163        NOTE::
     1164
     1165        Let `f = p/q` be a rational expression where `p` and `q` lie in a
     1166        `d`-variate polynomial ring `K[X]` for some field `K`.
     1167        Let `q_1^{e_1} \cdots q_m^{e_m}` be the
     1168        unique factorization of `q` in `K[X]` into irreducible factors and
     1169        let `V_i` be the algebraic variety `\{x\in L^d: q_i(x) = 0\}` of `q_i`
     1170         over the algebraic closure `L` of `K`.
     1171        By [Raic2012]_, `f` can be written as
     1172
     1173        (*)   `\sum p_A/\prod_{i \in A} q_i^{b_i}`,
     1174
     1175        where the `b_i` are positive integers, the `p_A` are in `K[X]`,
     1176        and the sum is taken over all subsets `A \subseteq \{1, \ldots, m\}`
     1177        such that
     1178        (1) `|A| \le d`,
     1179        (2) `\cap_{i\in A} T_i \neq \emptyset`, and
     1180        (3) `\{q_i : i\in A\}` is algebraically independent.
     1181
     1182        In particular, any rational expression in `d` variables
     1183        can be represented as a sum of rational expressions
     1184        whose denominators each contain at most `d` distinct irreducible
     1185        factors.
     1186
     1187        I call (*) a *Leinartas decomposition* of `f`.
     1188        Leinartas decompositions are not unique.
     1189
     1190        The algorithm used comes from [Raic2012]_.
     1191        """
     1192        d = len(self.ring().gens())
     1193        if d == 1:
     1194            # Sage's lift() function doesn't work in
     1195            # univariate polynomial rings.
     1196            # So nullstellensatz_decomposition() won't work.
     1197            # So use algebraic_dependence_decomposition(),
     1198            # which is sufficient.
     1199            temp = REFDSum([self])
     1200        else:
     1201            temp = self.nullstellensatz_decomposition()
     1202        decomp = REFDSum()
     1203        for r in temp:
     1204            decomp.extend(r.algebraic_dependence_decomposition())
     1205
     1206        # Simplify and return result.
     1207        return decomp.combine_like_terms().whole_and_parts()
     1208
     1209class REFDSum(list):
     1210    r"""
     1211    A list representing the sum of REFD objects with distinct
     1212    denominator lists.
     1213
     1214    AUTHORS:
     1215
     1216    - Alexander Raichev (2012-06-25)
     1217    """
     1218    def __str__(self):
     1219        return str([r.list() for r in self])
     1220
     1221    def __eq__(self, other):
     1222        return sorted(self) == sorted(other)
     1223
     1224    def __ne__(self, other):
     1225        return not (self == other)
     1226
     1227    def ring(self):
     1228        r"""
     1229        Return the ring of ``self``.
     1230        """
     1231        if self:
     1232            return self[0].ring()
     1233        else:
     1234            return None
     1235
     1236    def whole_and_parts(self):
     1237        r"""
     1238        Rewrite ``self`` as a REFDSum of a (possibly zero) polynomial
     1239        REFD followed by reduced rational expression REFDs.
     1240        Only useful for multivariate decompositions.
     1241
     1242        EXAMPLES::
     1243
     1244            sage: from sage.categories.quotient_fields import REFD
     1245            sage: from sage.categories.quotient_fields import REFDSum
     1246            sage: R = PolynomialRing(QQ, 'x, y')
     1247            sage: x, y = R.gens()
     1248            sage: f = x**2 + 3*y + 1/x + 1/y
     1249            sage: f = REFD(quotient=f)
     1250            sage: print f
     1251            (x^3*y + 3*x*y^2 + x + y, [(y, 1), (x, 1)])
     1252            sage: print REFDSum([f]).whole_and_parts()
     1253            [(x^2 + 3*y, []), (x + y, [(y, 1), (x, 1)])]
     1254
     1255        """
     1256        whole = 0
     1257        parts = []
     1258        R = self.ring()
     1259        for r in self:
     1260            # Since r has already passed through REFD.__init__()'s reducing
     1261            # procedure, r is already in lowest terms.
     1262            # Check if can write r as a mixed fraction: whole + fraction.
     1263            p = r.numerator()
     1264            q = r.denominator()
     1265            if q == 1:
     1266                # r is already a whole polynomial
     1267                whole += p
     1268            else:
     1269                # r is a fraction
     1270                num = R(r.numerator())  # Coerce into R in case its a constant.
     1271                denom = r.denominator()
     1272                a, b = num.quo_rem(denom)
     1273                whole += a
     1274                parts.append(REFD(b, r.denominator_list(), reduce_=False))
     1275        return REFDSum([REFD(whole, ())] + parts)
     1276
     1277    def combine_like_terms(self):
     1278        r"""
     1279        Combine terms in ``self`` with the same denominator.
     1280        Only useful for multivariate decompositions.
     1281
     1282        EXAMPLES::
     1283
     1284            sage: from sage.categories.quotient_fields import REFD
     1285            sage: from sage.categories.quotient_fields import REFDSum
     1286            sage: R.<x,y>= PolynomialRing(QQ)
     1287            sage: f = REFD(quotient=1/(x * y * (x*y + 1)))
     1288            sage: g = REFD(quotient=x/(x * y * (x*y + 1)))
     1289            sage: s = REFDSum([f, g, f])
     1290            sage: t = s.combine_like_terms()
     1291            sage: print s
     1292            [(1, [(y, 1), (x, 1), (x*y + 1, 1)]), (1, [(y, 1), (x*y + 1, 1)]), (1, [(y, 1), (x, 1), (x*y + 1, 1)])]
     1293            sage: print t
     1294            [(1, [(y, 1), (x*y + 1, 1)]), (2, [(y, 1), (x, 1), (x*y + 1, 1)])]
     1295
     1296        """
     1297        if not self:
     1298            return self
     1299
     1300        # Combine like terms.
     1301        refds = sorted(self)
     1302        new_refds = []
     1303        temp = refds[0]
     1304        for f in refds[1:]:
     1305            if  temp.denominator_list() == f.denominator_list():
     1306                # Add f to temp.
     1307                num = temp.numerator() + f.numerator()
     1308                temp = REFD(num, temp.denominator_list())
     1309            else:
     1310                # Append temp to new_refds and update temp.
     1311                new_refds.append(temp)
     1312                temp = f
     1313        new_refds.append(temp)
     1314        return REFDSum(new_refds)
     1315
     1316    def sum(self):
     1317        r"""
     1318        Return the sum of the REFDs in ``self`` as a REFD.
     1319
     1320        EXAMPLES::
     1321
     1322            sage: from sage.categories.quotient_fields import REFD
     1323            sage: from sage.categories.quotient_fields import REFDSum
     1324            sage: R.<x,y> = PolynomialRing(QQ)
     1325            sage: qs = (x, 1), (y, 1), (x*y + 1, 1)
     1326            sage: f = REFD(quotient=1/(x * y * (x*y + 1)))
     1327            sage: g = REFD(quotient=x*y/(x * y * (x*y + 1)))
     1328            sage: print REFDSum([f, g]).sum()
     1329            (1, [(y, 1), (x, 1)]) 
     1330
     1331        """
     1332        if not self:
     1333            return self
     1334
     1335        # Compute the sum's numerator and denominator.
     1336        summy = sum((f.quotient() for f in self))
     1337        num = summy.numerator()
     1338        denom = summy.denominator()
     1339        if denom == 1:
     1340            # Done
     1341            return REFD(num, [])
     1342
     1343        # Compute the irreducible factors of denom.       
     1344        # Could use the factor() command here, but that might be slow.
     1345        factors = []
     1346        for f in self:
     1347            factors.extend([q for q,e in f.denominator_list()])
     1348
     1349        # Eliminate repeats from factors and sort
     1350        factors = sorted(list(set(factors)))
     1351
     1352        # The irreducible factors of denom lie in factors.
     1353        # Use this fact to build a denominator list for denom.
     1354        denom_list = []
     1355        for q in factors:
     1356            e = 0
     1357            quo, rem = denom.quo_rem(q)
     1358            while rem == 0:
     1359                e += 1
     1360                denom = quo
     1361                quo, rem = denom.quo_rem(q)
     1362            if e > 0:
     1363                denom_list.append((q, e))
     1364        return REFD(num, denom_list, reduce_=False)
     1365 No newline at end of file