Ticket #10519: trac_10519.2.patch

File trac_10519.2.patch, 137.6 KB (added by araichev, 5 years ago)

Major rewrite.

  • sage/combinat/all.py

    # HG changeset patch
    # User Alex Raichev <alex.raichev@gmail.com>
    # Date 1344471313 -43200
    # Node ID e8f5867568b9d729afbaed0ea3329d02aaa84626
    # Parent  90cdc6a8ad4e4bed83142b4ca0e6c551acd73204
    #10519: fixed
    
    diff --git a/sage/combinat/all.py b/sage/combinat/all.py
    a b  
    118118from cyclic_sieving_phenomenon import CyclicSievingPolynomial, CyclicSievingCheck
    119119
    120120from sidon_sets import sidon_sets
     121
     122from amgf import FFPD
     123 No newline at end of file
  • new file sage/combinat/amgf.py

    diff --git a/sage/combinat/amgf.py b/sage/combinat/amgf.py
    new file mode 100644
    - +  
     1r"""
     2Let $F(x) = \sum_{\nu \in \NN^d} F_{\nu} x^\nu$ be a multivariate power series with complex coefficients that converges in a neighborhood of the origin. Assume that $F = G/H$ for some functions $G$ and $H$ holomorphic in a neighborhood of the origin.
     3Assume also that $H$ is a polynomial.
     4
     5This Python module for use within `Sage <http://www.sagemath.org>`_ computes asymptotics for the coefficients $F_{r \alpha}$ as $r \to \infty$ with $r \alpha \in \NN^d$ for $\alpha$ in a permissible subset of $d$-tuples of positive reals.
     6More specifically, it computes arbitrary terms of the asymptotic expansion for $F_{r \alpha}$ when the asymptotics are controlled by a multiple point of the alegbraic variety $H = 0$.
     7
     8The algorithms and formulas implemented here come from [RaWi2008a]_
     9and [RaWi2012]_.
     10
     11.. [AiYu1983] I.A. Aizenberg and A.P. Yuzhakov, "Integral representations and residues in multidimensional complex analysis", Translations of Mathematical Monographs, 58. American Mathematical Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X.
     12
     13.. [Raic2012] Alexander Raichev, "Leinartas's partial fraction decomposition", `<http://arxiv.org/abs/1206.4740>`_.
     14
     15.. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", Electronic Journal of Combinatorics, Vol. 15 (2008), R89, `<http://arxiv.org/pdf/0803.2914.pdf>`_.
     16
     17.. [RaWi2012] Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", To appear in 2012 in the Online Journal of Analytic Combinatorics, `<http://arxiv.org/pdf/1009.5715.pdf>`_.
     18
     19AUTHORS:
     20
     21- Alexander Raichev (2008-10-01): Initial version
     22- Alexander Raichev (2010-09-28): Corrected many functions
     23- Alexander Raichev (2010-12-15): Updated documentation
     24- Alexander Raichev (2011-03-09): Fixed a division by zero bug in relative_error()
     25- Alexander Raichev (2011-04-26): Rewrote in object-oreinted style
     26- Alexander Raichev (2011-05-06): Fixed bug in cohomologous_integrand() and fixed _crit_cone_combo() to work in SR
     27- Alexander Raichev (2012-08-06): Major rewrite. Created class FFPD and moved functions around.
     28
     29EXAMPLES::
     30
     31A smooth point example (Example 5.4 of [RaWi2008a]_)::
     32
     33    sage: from sage.combinat.amgf import *
     34    sage: R.<x,y> = PolynomialRing(QQ)
     35    sage: q = 1/2
     36    sage: qq = q.denominator()
     37    sage: H = 1 - q*x + q*x*y - x^2*y
     38    sage: Hfac = H.factor()
     39    sage: G = (1 - q*x)/Hfac.unit()
     40    sage: F = FFPD(G, Hfac)
     41    sage: alpha = list(qq*vector([2, 1 - q]))
     42    sage: print alpha
     43    [4, 1]
     44    sage: I = F.smooth_critical_ideal(alpha)
     45    sage: print I
     46    Ideal (y^2 - 2*y + 1, x + 1/4*y - 5/4) of Multivariate Polynomial Ring
     47    in x, y over Rational Field
     48    sage: s = solve(I.gens(), [SR(x) for x in R.gens()], solution_dict=true)
     49    sage: print s
     50    [{y: 1, x: 1}]
     51    sage: p = s[0]
     52    sage: asy = F.asymptotics(p, alpha, 1) # long time
     53    Creating auxiliary functions...
     54    Computing derivatives of auxiallary functions...
     55    Computing derivatives of more auxiliary functions...
     56    Computing second order differential operator actions...
     57    sage: print asy # long time
     58    (1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*r^(1/3)), 1,
     59    1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*r^(1/3)))
     60    sage: print F.relative_error(asy[0], alpha, [1, 2, 4, 8, 16], asy[1]) # long time
     61    Calculating errors table in the form
     62    exponent, scaled Maclaurin coefficient, scaled asymptotic values,
     63    relative errors...
     64    [((4, 1), 0.1875000000, [0.1953794675], [-0.04202382689]), ((8, 2),
     65    0.1523437500, [0.1550727862], [-0.01791367323]), ((16, 4), 0.1221771240,
     66    [0.1230813519], [-0.007400959228]), ((32, 8), 0.09739671811,
     67    [0.09768973377], [-0.003008475766]), ((64, 16), 0.07744253816,
     68    [0.07753639308], [-0.001211929722])]
     69   
     70A multiple point example (Example 6.5 of [RaWi2012]_)::
     71
     72    sage: R.<x,y>= PolynomialRing(QQ)
     73    sage: H = (1 - 2*x - y)**2 * (1 - x - 2*y)**2
     74    sage: Hfac = H.factor()
     75    sage: G = 1/Hfac.unit()
     76    sage: F = FFPD(G, Hfac)
     77    sage: print F
     78    (1, [(x + 2*y - 1, 2), (2*x + y - 1, 2)])
     79    sage: I = F.singular_ideal()
     80    sage: print I
     81    Ideal (x - 1/3, y - 1/3) of Multivariate Polynomial Ring in x, y over
     82    Rational Field
     83    sage: p = {x: 1/3, y: 1/3}
     84    sage: print F.is_convenient_multiple_point(p)
     85    (True, 'convenient in variables [x, y]')
     86    sage: alpha = (var('a'), var('b'))
     87    sage: print F.asymptotic_decomposition(alpha) # long time
     88    [(0, []), (-1/9*(2*a^2*y^2 - 5*a*b*x*y + 2*b^2*x^2)*r^2/(x^2*y^2) +
     89    1/9*(5*(a + b)*x*y - 6*a*y^2 - 6*b*x^2)*r/(x^2*y^2) - 1/9*(4*x^2 - 5*x*y
     90    + 4*y^2)/(x^2*y^2), [(x + 2*y - 1, 1), (2*x + y - 1, 1)])]
     91    sage: print F.asymptotics(p, alpha, 2) # long time
     92    (-3*((2*a^2 - 5*a*b + 2*b^2)*r^2 + (a + b)*r +
     93    3)*((1/3)^(-b)*(1/3)^(-a))^r, (1/3)^(-b)*(1/3)^(-a), -3*(2*a^2 - 5*a*b +
     94    2*b^2)*r^2 - 3*(a + b)*r - 9)
     95    sage: alpha = [4, 3]
     96    sage: asy = F.asymptotics(p, alpha, 2) # long time
     97    sage: print asy # long time
     98    (3*(10*r^2 - 7*r - 3)*2187^r, 2187, 30*r^2 - 21*r - 9)
     99    sage: print F.relative_error(asy[0], alpha, [1, 2, 4, 8], asy[1]) # long time
     100    Calculating errors table in the form
     101    exponent, scaled Maclaurin coefficient, scaled asymptotic values,
     102    relative errors...
     103    [((4, 3), 30.72702332, [0.0000000000], [1.000000000]), ((8, 6),
     104    111.9315678, [69.00000000], [0.3835519207]), ((16, 12), 442.7813138,
     105    [387.0000000], [0.1259793763]), ((32, 24), 1799.879232, [1743.000000],
     106    [0.03160169385])]
     107"""
     108#*****************************************************************************
     109#       Copyright (C) 2008 Alexander Raichev <tortoise.said@gmail.com>
     110#
     111#  Distributed under the terms of the GNU General Public License (GPL)
     112#                  http://www.gnu.org/licenses/
     113#*****************************************************************************
     114
     115from functools import total_ordering
     116
     117# Sage libraries
     118from sage.categories.unique_factorization_domains   import UniqueFactorizationDomains
     119from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     120from sage.rings.polynomial.polynomial_ring  import is_PolynomialRing
     121from sage.rings.polynomial.multi_polynomial_ring_generic    import is_MPolynomialRing
     122from sage.symbolic.ring import SR
     123from sage.geometry.cone import Cone
     124from sage.calculus.functional   import diff
     125from sage.calculus.functions    import jacobian
     126from sage.calculus.var  import function, var
     127from sage.combinat.cartesian_product    import CartesianProduct
     128from sage.combinat.combinat     import stirling_number1
     129from sage.combinat.permutation  import Permutation
     130from sage.combinat.tuple    import UnorderedTuples
     131from sage.functions.log     import exp, log
     132from sage.functions.other   import factorial, gamma, sqrt
     133from sage.matrix.constructor    import matrix
     134from sage.misc.misc import add
     135from sage.misc.misc_c   import prod
     136from sage.misc.mrange   import cartesian_product_iterator, mrange
     137from sage.modules.free_module_element   import vector
     138from sage.rings.arith   import binomial
     139from sage.rings.all import CC
     140from sage.rings.fraction_field  import FractionField
     141from sage.rings.integer import Integer
     142from sage.rings.integer_ring    import ZZ
     143from sage.rings.polynomial.polynomial_ring_constructor  import PolynomialRing
     144from sage.rings.rational_field  import QQ
     145from sage.sets.set  import Set
     146from sage.structure.sage_object import SageObject
     147from sage.symbolic.constants    import pi
     148from sage.symbolic.relation import solve
     149
     150@total_ordering
     151class FFPD(object):
     152    r"""
     153    Represents a fraction with factored polynomial denominator (FFPD)
     154    $p/(q_1^{e_1} \cdots q_n^{e_n})$ by storing the parts $p$ and
     155    $[(q_1, e_1), \ldots, (q_n, e_n)]$.
     156    Here $q_1, \ldots, q_n$ are elements of a 0- or multi-variate factorial
     157    polynomial ring $R$ , $q_1, \ldots, q_n$ are distinct irreducible elements
     158    of $R$ , $e_1, \ldots, e_n$ are positive integers, and $p$ is a function
     159    of the indeterminates of $R$ (a Sage Symbolic Expression).
     160    An element $r$ with no polynomial denominator is represented as $[r, (,)]$.
     161
     162    AUTHORS:
     163
     164    - Alexander Raichev (2012-07-26)
     165    """
     166
     167    def __init__(self, numerator=None, denominator_factored=None,
     168                 quotient=None, reduce_=True):
     169        r"""
     170        Create a FFPD instance.
     171
     172        INPUT:
     173
     174        - ``numerator`` - (Optional; default=None) An element $p$ of a
     175          0- or 1-variate factorial polynomial ring $R$.
     176        - ``denominator_factored`` - (Optional; default=None)
     177          A list of the form
     178          $[(q_1, e_1), \ldots, (q_n, e_n)]$ where the $q_1, \ldots, q_n$ are
     179          distinct irreducible elements of $R$ and the $e_i$ are positive
     180          integers.
     181        - ``quotient`` - (Optional; default=None) An element of a field of
     182          fractions of a factorial ring.
     183        - ``reduce_`` - (Optional; default=True) If True, then represent
     184          $p/(q_1^{e_1} \cdots q_n^{e_n})$ in lowest terms.
     185          If False, then won't attempt to divide $p$ by any of the $q_i$.
     186
     187        OUTPUT:
     188
     189        A FFPD instance representing the rational expression
     190        $p/(q_1^{e_1} \cdots q_n^{e_n})$.
     191        To get a non-None output, one of ``numerator`` or ``quotient`` must be
     192        non-None.
     193
     194        EXAMPLES::
     195
     196            sage: from sage.combinat.amgf import *
     197
     198            sage: R.<x, y> = PolynomialRing(QQ)
     199            sage: df = [x, 1], [y, 1], [x*y+1, 1]
     200            sage: f = FFPD(x, df)
     201            sage: print f
     202            (1, [(y, 1), (x*y + 1, 1)])
     203            sage: ff = FFPD(x, df, reduce_=False)
     204            sage: print ff
     205            (x, [(y, 1), (x, 1), (x*y + 1, 1)])
     206
     207        ::
     208       
     209            sage: f = FFPD(x + y, [(x + y, 1)])
     210            sage: print f
     211            (1, [])
     212
     213        ::
     214       
     215            sage: R.<x> = PolynomialRing(QQ)
     216            sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
     217            sage: print FFPD(quotient=f)
     218            (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,
     219            [(x - 1, 1), (x, 1), (x^2 + 1/3, 1)])
     220
     221        ::
     222
     223            sage: R.<x, y> = PolynomialRing(QQ)
     224            sage: f = 2*y/(5*(x^3 - 1)*(y + 1))
     225            sage: print FFPD(quotient=f)
     226            (2/5*y, [(y + 1, 1), (x - 1, 1), (x^2 + x + 1, 1)])
     227           
     228        ::
     229       
     230            sage: R.<x, y>= PolynomialRing(QQ)
     231            sage: p = 1/x^2
     232            sage: q = 3*x**2*y
     233            sage: qs = q.factor()
     234            sage: f = FFPD(p/qs.unit(), qs)
     235            sage: print f
     236            (1/(3*x^2), [(y, 1), (x, 2)])
     237           
     238        ::
     239       
     240            sage: R.<x, y> = PolynomialRing(QQ)
     241            sage: f = FFPD(cos(x)*x*y^2, [(x, 2), (y, 1)])
     242            sage: print f
     243            (x*y^2*cos(x), [(y, 1), (x, 2)])
     244
     245        ::
     246           
     247            sage: R.<x, y> = PolynomialRing(QQ)
     248            sage: G = exp(x + y)
     249            sage: H = (1 - 2*x - y) * (1 - x - 2*y)
     250            sage: a = FFPD(quotient=G/H)
     251            sage: print a
     252            (e^(x + y)/(2*x^2 + 5*x*y + 2*y^2 - 3*x - 3*y + 1), [])
     253            sage: print a._ring
     254            None
     255            sage: b = FFPD(G, H.factor())
     256            sage: print b
     257            (e^(x + y), [(x + 2*y - 1, 1), (2*x + y - 1, 1)])
     258            sage: print b._ring
     259            Multivariate Polynomial Ring in x, y over Rational Field
     260           
     261        Singular throws a 'not implemented' error when trying to factor in
     262        a multivariate polynomial ring over an inexact field ::
     263
     264            sage: R.<x, y>= PolynomialRing(CC)
     265            sage: f = (x + 1)/(x*y*(x*y + 1)^2)
     266            sage: FFPD(quotient=f)
     267            Traceback (most recent call last):
     268            ...
     269            TypeError: Singular error:
     270               ? not implemented
     271               ? error occurred in or before STDIN line 17:
     272               `def sage9=factorize(sage8);`
     273
     274        """       
     275        # Attributes are
     276        # self._numerator           
     277        # self._denominator_factored   
     278        # self._ring
     279        if quotient is not None:
     280            p = quotient.numerator()
     281            q = quotient.denominator()
     282            R = q.parent()
     283            self._numerator = quotient
     284            self._denominator_factored = [] 
     285            if is_PolynomialRing(R) or is_MPolynomialRing(R):
     286                self._ring = R
     287                if not R(q).is_unit():
     288                    # Factor q
     289                    try:
     290                        df = q.factor()
     291                    except NotImplementedError:
     292                        # Singular's factor() needs 'proof=False'.
     293                        df = q.factor(proof=False)   
     294                    self._numerator = p/df.unit()
     295                    df = sorted([tuple(t) for t in df]) # Sort for consitency.
     296                    self._denominator_factored = df
     297            else:
     298                self._ring = None
     299            # Done. No reducing needed, as Sage reduced quotient beforehand.
     300            return
     301           
     302        self._numerator = numerator
     303        if denominator_factored:
     304            self._denominator_factored = sorted([tuple(t) for t in
     305                                             denominator_factored])
     306            self._ring = denominator_factored[0][0].parent()     
     307        else:
     308            self._denominator_factored = []
     309            self._ring = None
     310        R = self._ring
     311        if R is not None and numerator in R and reduce_:
     312            # Reduce fraction if possible.
     313            numer = R(self._numerator)
     314            df = self._denominator_factored
     315            new_df = []
     316            for (q, e) in df:
     317                ee = e
     318                quo, rem = numer.quo_rem(q)
     319                while rem == 0 and ee > 0:
     320                    ee -= 1
     321                    numer = quo
     322                    quo, rem = numer.quo_rem(q)
     323                if ee > 0:
     324                    new_df.append((q, ee))
     325            self._numerator = numer
     326            self._denominator_factored = new_df
     327
     328    def numerator(self):
     329        r"""
     330        Return the numerator of ``self``.
     331       
     332        EXAMPLES::
     333       
     334            sage: from sage.combinat.amgf import *
     335            sage: R.<x,y>= PolynomialRing(QQ)
     336            sage: H = (1 - x - y - x*y)**2*(1-x)
     337            sage: Hfac = H.factor()
     338            sage: G = exp(y)/Hfac.unit()
     339            sage: F = FFPD(G, Hfac)
     340            sage: print F.numerator()
     341            -e^y
     342        """
     343        return self._numerator
     344
     345    def denominator(self):
     346        r"""
     347        Return the denominator of ``self``.
     348       
     349        EXAMPLES::
     350       
     351            sage: from sage.combinat.amgf import *
     352            sage: R.<x,y>= PolynomialRing(QQ)
     353            sage: H = (1 - x - y - x*y)**2*(1-x)
     354            sage: Hfac = H.factor()
     355            sage: G = exp(y)/Hfac.unit()
     356            sage: F = FFPD(G, Hfac)
     357            sage: print F.denominator()
     358            x^3*y^2 + 2*x^3*y + x^2*y^2 + x^3 - 2*x^2*y - x*y^2 - 3*x^2 - 2*x*y
     359            - y^2 + 3*x + 2*y - 1
     360        """
     361        return prod([q**e for q, e in self.denominator_factored()])
     362
     363    def denominator_factored(self):
     364        r"""
     365        Return the factorization in ``self.ring()`` of the denominator of
     366        ``self`` but without the unit part.
     367       
     368        EXAMPLES::
     369       
     370            sage: from sage.combinat.amgf import *
     371            sage: R.<x,y>= PolynomialRing(QQ)
     372            sage: H = (1 - x - y - x*y)**2*(1-x)
     373            sage: Hfac = H.factor()
     374            sage: G = exp(y)/Hfac.unit()
     375            sage: F = FFPD(G, Hfac)
     376            sage: print F.denominator_factored()
     377            [(x - 1, 1), (x*y + x + y - 1, 2)]
     378        """
     379        return self._denominator_factored
     380
     381    def ring(self):
     382        r"""
     383        Return the ring of the denominator of ``self``, which is
     384        None in the case where ``self`` doesn't have a denominator.
     385       
     386        EXAMPLES::
     387       
     388            sage: from sage.combinat.amgf import *
     389            sage: R.<x,y>= PolynomialRing(QQ)
     390            sage: H = (1 - x - y - x*y)**2*(1-x)
     391            sage: Hfac = H.factor()
     392            sage: G = exp(y)/Hfac.unit()
     393            sage: F = FFPD(G, Hfac)
     394            sage: print F.ring()
     395            Multivariate Polynomial Ring in x, y over Rational Field
     396            sage: F = FFPD(quotient=G/H)
     397            sage: print F
     398            (e^y/(x^3*y^2 + 2*x^3*y + x^2*y^2 + x^3 - 2*x^2*y - x*y^2 - 3*x^2 -
     399            2*x*y - y^2 + 3*x + 2*y - 1), [])
     400            sage: print F.ring()
     401            None
     402        """
     403        return self._ring
     404       
     405    def dimension(self):
     406        r"""
     407        Return the number of indeterminates of ``self.ring()``.
     408       
     409        EXAMPLES::
     410       
     411            sage: from sage.combinat.amgf import *
     412            sage: R.<x,y>= PolynomialRing(QQ)
     413            sage: H = (1 - x - y - x*y)**2*(1-x)
     414            sage: Hfac = H.factor()
     415            sage: G = exp(y)/Hfac.unit()
     416            sage: F = FFPD(G, Hfac)
     417            sage: print F.dimension()
     418            2
     419        """       
     420        R = self.ring()
     421        if is_PolynomialRing(R) or is_MPolynomialRing(R):
     422            return R.ngens()
     423        else:
     424            return None
     425       
     426    def list(self):
     427        r"""
     428        Convert ``self`` into a list for printing.
     429        """
     430        return (self.numerator(), self.denominator_factored())
     431
     432    def quotient(self):
     433        r"""
     434        Convert ``self`` into a quotient.
     435       
     436        EXAMPLES::
     437
     438            sage: from sage.combinat.amgf import *
     439
     440            sage: R.<x,y>= PolynomialRing(QQ)
     441            sage: H = (1 - x - y - x*y)**2*(1-x)
     442            sage: Hfac = H.factor()
     443            sage: G = exp(y)/Hfac.unit()
     444            sage: F = FFPD(G, Hfac)
     445            sage: print F
     446            (-e^y, [(x - 1, 1), (x*y + x + y - 1, 2)])
     447            sage: print F.quotient()
     448            -e^y/(x^3*y^2 + 2*x^3*y + x^2*y^2 + x^3 - 2*x^2*y - x*y^2 - 3*x^2 -
     449            2*x*y - y^2 + 3*x + 2*y - 1)
     450        """
     451        return self.numerator()/self.denominator()
     452
     453    def __str__(self):
     454        return str(self.list())
     455
     456    def __eq__(self, other):
     457        r"""
     458        Two FFPD instances are equal iff they represent the same
     459        fraction.
     460
     461        EXAMPLES::
     462
     463            sage: from sage.combinat.amgf import *
     464
     465            sage: R.<x, y>= PolynomialRing(QQ)
     466            sage: df = [x, 1], [y, 1], [x*y+1, 1]
     467            sage: f = FFPD(x, df)
     468            sage: ff = FFPD(x, df, reduce_=False)
     469            sage: f == ff
     470            True
     471            sage: g = FFPD(y, df)
     472            sage: g == f
     473            False
     474           
     475        ::
     476       
     477            sage: R.<x, y> = PolynomialRing(QQ)
     478            sage: G = exp(x + y)
     479            sage: H = (1 - 2*x - y) * (1 - x - 2*y)
     480            sage: a = FFPD(quotient=G/H)
     481            sage: b = FFPD(G, H.factor())
     482            sage: if (a == b):
     483            ...       print 'Yep'
     484            ...       
     485            Yep
     486        """
     487        return self.quotient() == other.quotient()
     488
     489    def __ne__(self, other):
     490        r"""
     491        EXAMPLES::
     492
     493            sage: from sage.combinat.amgf import *
     494
     495            sage: R.<x, y>= PolynomialRing(QQ)
     496            sage: df = [x, 1], [y, 1], [x*y+1, 1]
     497            sage: f = FFPD(x, df)
     498            sage: ff = FFPD(x, df, reduce_=False)
     499            sage: f != ff
     500            False
     501            sage: g = FFPD(y, df)
     502            sage: g != f
     503            True
     504        """
     505        return not (self == other)
     506
     507    def __lt__(self, other):
     508        r"""
     509        FFPD A is less than FFPD B iff
     510        (the denominator factorization of A is shorter than that of B) or
     511        (the denominator factorization lengths are equal and
     512        the denominator of A is less than that of B in their ring) or
     513        (the denominator factorization lengths are equal and the
     514        denominators are equal and the numerator of A is less than that of B
     515        in their ring).
     516
     517        EXAMPLES::
     518
     519            sage: from sage.combinat.amgf import *
     520
     521            sage: R.<x, y>= PolynomialRing(QQ)
     522            sage: df = [x, 1], [y, 1], [x*y+1, 1]
     523            sage: f = FFPD(x, df)
     524            sage: ff = FFPD(x, df, reduce_=False)
     525            sage: g = FFPD(y, df)
     526            sage: h = FFPD(exp(x), df)
     527            sage: i = FFPD(sin(x + 2), df)
     528            sage: print f, ff
     529            (1, [(y, 1), (x*y + 1, 1)]) (x, [(y, 1), (x, 1), (x*y + 1, 1)])
     530            sage: print f < ff
     531            True
     532            sage: print f < g
     533            True
     534            sage: print g < h
     535            True
     536            sage: print h < i
     537            False
     538        """
     539        sn = self.numerator()
     540        on = other.numerator()
     541        sdf = self.denominator_factored()
     542        odf = other.denominator_factored()
     543        sd = self.denominator()
     544        od = other.denominator()
     545
     546        if len(sdf) < len(odf) or\
     547          (len(sdf) == len(odf) and sd < od) or\
     548          (len(sdf) == len(odf) and sd == od and sn < on):
     549            return True
     550        else:
     551            return False     
     552
     553    def univariate_decomposition(self):
     554        r"""
     555        Return the usual univariate partial fraction decomposition
     556        of ``self`` as a FFPDSum instance.
     557        Assume that ``self`` lies in the field of fractions
     558        of a univariate factorial polynomial ring.
     559
     560        EXAMPLES::
     561
     562            sage: from sage.combinat.amgf import *
     563
     564        One variable::
     565
     566            sage: R.<x> = PolynomialRing(QQ)
     567            sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
     568            sage: print f
     569            (15*x^7 - 15*x^6 + 5*x^5 - 5*x^4 + 6*x^3 - 2*x^2 + x - 1)/(3*x^4 -
     570            3*x^3 + x^2 - x)
     571            sage: decomp = FFPD(quotient=f).univariate_decomposition()
     572            sage: print decomp
     573            [(5*x^3, []), (1, [(x - 1, 1)]), (1, [(x, 1)]),
     574            (1/3, [(x^2 + 1/3, 1)])]
     575            sage: print decomp.sum().quotient() == f
     576            True
     577
     578        One variable with numerator in symbolic ring::
     579
     580            sage: R.<x> = PolynomialRing(QQ)
     581            sage: f = 5*x^3 + 1/x + 1/(x-1) + exp(x)/(3*x^2 + 1)
     582            sage: print f
     583            e^x/(3*x^2 + 1) + ((5*(x - 1)*x^3 + 2)*x - 1)/((x - 1)*x)
     584            sage: decomp = FFPD(quotient=f).univariate_decomposition()
     585            sage: print decomp
     586            [(e^x/(3*x^2 + 1) + ((5*(x - 1)*x^3 + 2)*x - 1)/((x - 1)*x), [])]
     587
     588        One variable over a finite field::
     589
     590            sage: R.<x> = PolynomialRing(GF(2))
     591            sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
     592            sage: print f
     593            (x^6 + x^4 + 1)/(x^3 + x)
     594            sage: decomp = FFPD(quotient=f).univariate_decomposition()
     595            sage: print decomp
     596            [(x^3, []), (1, [(x, 1)]), (x, [(x + 1, 2)])]
     597            sage: print decomp.sum().quotient() == f
     598            True
     599
     600        One variable over an inexact field::
     601
     602            sage: R.<x> = PolynomialRing(CC)
     603            sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
     604            sage: print f
     605            (15.0000000000000*x^7 - 15.0000000000000*x^6 + 5.00000000000000*x^5
     606            - 5.00000000000000*x^4 + 6.00000000000000*x^3 -
     607            2.00000000000000*x^2 + x - 1.00000000000000)/(3.00000000000000*x^4
     608            - 3.00000000000000*x^3 + x^2 - x)
     609            sage: decomp = FFPD(quotient=f).univariate_decomposition()
     610            sage: print decomp
     611            [(5.00000000000000*x^3, []), (1.00000000000000,
     612            [(x - 1.00000000000000, 1)]), (-0.288675134594813*I,
     613            [(x - 0.577350269189626*I, 1)]), (1.00000000000000, [(x, 1)]),
     614            (0.288675134594813*I, [(x + 0.577350269189626*I, 1)])]
     615            sage: print decomp.sum().quotient() == f # Rounding error coming
     616            False
     617
     618        NOTE::
     619
     620        Let $f = p/q$ be a rational expression where $p$ and $q$ lie in a
     621        univariate factorial polynomial ring $R$.
     622        Let $q_1^{e_1} \cdots q_n^{e_n}$ be the
     623        unique factorization of $q$ in $R$ into irreducible factors.
     624        Then $f$ can be written uniquely as
     625
     626        (*)   $p_0 + \sum_{i=1}^{m} p_i/ q_i^{e_i}$,
     627
     628        for some $p_j \in R$.
     629        I call (*) the *usual partial fraction decomposition* of f.
     630
     631        AUTHORS:
     632
     633        - Robert Bradshaw (2007-05-31)
     634        - Alexander Raichev (2012-06-25)
     635        """       
     636        if self.dimension() is None or self.dimension() > 1:
     637            return FFPDSum([self])
     638
     639        R = self.ring()
     640        p = self.numerator()
     641        q = self.denominator()
     642        if p in R:
     643            whole, p = p.quo_rem(q)
     644        else:
     645            whole = p
     646            p = R(1)
     647        df = self.denominator_factored()
     648        decomp = [FFPD(whole, [])]
     649        for (a, m) in df:
     650            numer = p * prod([b**n for (b, n) in df if b != a]).\
     651                    inverse_mod(a**m) % (a**m) 
     652            # The inverse exists because the product and a**m
     653            # are relatively prime.
     654            decomp.append(FFPD(numer, [(a, m)]))
     655        return FFPDSum(decomp)
     656
     657    def nullstellensatz_certificate(self):
     658        r"""
     659        Let $[(q_1, e_1), \ldots, (q_n, e_n)]$ be the denominator factorization
     660        of ``self``.
     661        Return a list of polynomials $h_1, \ldots, h_m$ in ``self.ring()``
     662        that satisfies $h_1 q_1 + \cdots + h_m q_n = 1$ if it exists.
     663        Otherwise return None.
     664        Only works for multivariate ``self``.
     665
     666        EXAMPLES::
     667
     668            sage: from sage.combinat.amgf import *
     669
     670            sage: R.<x, y> = PolynomialRing(QQ)
     671            sage: G = sin(x)
     672            sage: H = x^2 * (x*y + 1)
     673            sage: f = FFPD(G, H.factor())           
     674            sage: L = f.nullstellensatz_certificate()
     675            sage: print L
     676            [y^2, -x*y + 1]
     677            sage: df = f.denominator_factored()
     678            sage: sum([L[i]*df[i][0]**df[i][1] for i in xrange(len(df))]) == 1
     679            True
     680
     681        ::
     682
     683            sage: f = 1/(x*y)
     684            sage: L = FFPD(quotient=f).nullstellensatz_certificate()
     685            sage: L is None
     686            True
     687
     688        """
     689
     690        R = self.ring()
     691        if R is None:
     692            return None
     693           
     694        df = self.denominator_factored()
     695        J = R.ideal([q**e for q, e in df])
     696        if R(1) in J:
     697            return R(1).lift(J)
     698        else:
     699            return None
     700
     701    def nullstellensatz_decomposition(self):
     702        r"""
     703        Return a Nullstellensatz decomposition of ``self`` as a FFPDSum
     704        instance.
     705        Recursive.
     706        Only works for multivariate ``self``.
     707
     708        EXAMPLES::
     709
     710            sage: from sage.combinat.amgf import *
     711
     712            sage: R.<x, y> = PolynomialRing(QQ)
     713            sage: f = 1/(x*(x*y + 1))
     714            sage: decomp = FFPD(quotient=f).nullstellensatz_decomposition()
     715            sage: print decomp
     716            [(0, []), (1, [(x, 1)]), (-y, [(x*y + 1, 1)])]
     717            sage: decomp.sum().quotient() == f
     718            True
     719            sage: for r in decomp:
     720            ...       L = r.nullstellensatz_certificate()
     721            ...       L is None
     722            ...       
     723            True
     724            True
     725            True
     726
     727        ::
     728       
     729            sage: R.<x, y> = PolynomialRing(QQ)
     730            sage: G = sin(y)
     731            sage: H = x*(x*y + 1)
     732            sage: f = FFPD(G, H.factor())
     733            sage: decomp = f.nullstellensatz_decomposition()
     734            sage: print decomp
     735            [(0, []), (sin(y), [(x, 1)]), (-y*sin(y), [(x*y + 1, 1)])]
     736            sage: if decomp.sum().quotient() == G/H:
     737            ...       print 'yep'
     738            ...       
     739            yep
     740            sage: for r in decomp:
     741            ...       L = r.nullstellensatz_certificate()
     742            ...       L is None
     743            ...       
     744            True
     745            True
     746            True
     747
     748        NOTE::
     749
     750        Let $f = p/q$ where $q$ lies in a $d$ -variate polynomial ring $K[X]$ for some field $K$ and $d \ge 1$.
     751        Let $q_1^{e_1} \cdots q_n^{e_n}$ be the
     752        unique factorization of $q$ in $K[X]$ into irreducible factors and
     753        let $V_i$ be the algebraic variety $\{x \in L^d: q_i(x) = 0\}$ of
     754        $q_i$ over the algebraic closure $L$ of $K$.
     755        By [Raic2012]_, $f$ can be written as
     756
     757        (*)   $\sum p_A/\prod_{i \in A} q_i^{e_i}$,
     758
     759        where the $p_A$ are products of $p$ and elements in $K[X]$ and
     760        the sum is taken over all subsets
     761        $A \subseteq \{1, \ldots, m\}$ such that
     762        $\cap_{i\in A} T_i \neq \emptyset$.
     763
     764        I call (*) a *Nullstellensatz decomposition* of $f$.
     765        Nullstellensatz decompositions are not unique.
     766
     767        The algorithm used comes from [Raic2012]_.       
     768        """       
     769        L = self.nullstellensatz_certificate()
     770        if L is None:
     771            # No decomposing possible.
     772            return FFPDSum([self])
     773   
     774        # Otherwise decompose recursively.
     775        decomp = FFPDSum()
     776        p = self.numerator()
     777        df = self.denominator_factored()           
     778        m = len(df)
     779        iteration1 = FFPDSum([FFPD(p*L[i],[df[j]
     780                              for j in xrange(m) if j != i])
     781                              for i in xrange(m) if L[i] != 0])
     782
     783        # Now decompose each FFPD of iteration1.
     784        for r in iteration1:
     785            decomp.extend(r.nullstellensatz_decomposition())
     786
     787        # Simplify and return result.
     788        return decomp.combine_like_terms().whole_and_parts()
     789
     790    def algebraic_dependence_certificate(self):
     791        r"""
     792        Return the ideal $J$ of annihilating polynomials for the set
     793        of polynomials ``[q**e for (q, e) in self.denominator_factored()]``,
     794        which could be the zero ideal.
     795        The ideal $J$ lies in a polynomial ring over the field
     796        ``self.ring().base_ring()`` that has
     797        ``m = len(self.denominator_factored())`` indeterminates.
     798        Return None if ``self.ring()`` is None.
     799
     800        EXAMPLES::
     801
     802            sage: from sage.combinat.amgf import *
     803
     804            sage: R.<x, y> = PolynomialRing(QQ)
     805            sage: f = 1/(x^2 * (x*y + 1) * y^3)
     806            sage: ff = FFPD(quotient=f)
     807            sage: J = ff.algebraic_dependence_certificate()
     808            sage: print J
     809            Ideal (1 - 6*T2 + 15*T2^2 - 20*T2^3 + 15*T2^4 - T0^2*T1^3 -
     810            6*T2^5  + T2^6) of Multivariate Polynomial Ring in
     811            T0, T1, T2 over Rational Field
     812            sage: g = J.gens()[0]
     813            sage: df = ff.denominator_factored()
     814            sage: g(*(q**e for q, e in df)) == 0
     815            True
     816
     817        ::
     818       
     819            sage: R.<x, y> = PolynomialRing(QQ)
     820            sage: G = exp(x + y)
     821            sage: H = x^2 * (x*y + 1) * y^3
     822            sage: ff = FFPD(G, H.factor())
     823            sage: J = ff.algebraic_dependence_certificate()
     824            sage: print J
     825            Ideal (1 - 6*T2 + 15*T2^2 - 20*T2^3 + 15*T2^4 - T0^2*T1^3 -
     826            6*T2^5 + T2^6) of Multivariate Polynomial Ring in
     827            T0, T1, T2 over Rational Field
     828            sage: g = J.gens()[0]
     829            sage: df = ff.denominator_factored()
     830            sage: g(*(q**e for q, e in df)) == 0
     831            True
     832           
     833        ::
     834       
     835            sage: f = 1/(x^3 * y^2)
     836            sage: J = FFPD(quotient=f).algebraic_dependence_certificate()
     837            sage: print J
     838            Ideal (0) of Multivariate Polynomial Ring in T0, T1 over
     839            Rational Field
     840
     841        ::
     842       
     843            sage: f = sin(1)/(x^3 * y^2)
     844            sage: J = FFPD(quotient=f).algebraic_dependence_certificate()
     845            sage: print J
     846            None
     847        """
     848        R = self.ring()
     849        if R is None:
     850            return None
     851           
     852        df = self.denominator_factored()       
     853        if not df:
     854            return R.ideal()    # The zero ideal.
     855        m = len(df)
     856        F = R.base_ring()
     857        Xs = list(R.gens())
     858        d = len(Xs)
     859
     860        # Expand R by 2*m new variables.
     861        S = 'S'
     862        while S in [str(x) for x in Xs]:
     863            S = S + 'S'
     864        Ss = [S + str(i) for i in xrange(m)]
     865        T = 'T'
     866        while T in [str(x) for x in Xs]:
     867            T = T + 'T'
     868        Ts = [T + str(i) for i in xrange(m)]
     869
     870        Vs = [str(x) for x in Xs] + Ss + Ts
     871        RR = PolynomialRing(F, Vs)
     872        Xs = RR.gens()[:d]
     873        Ss = RR.gens()[d: d + m]
     874        Ts = RR.gens()[d + m: d + 2 * m]
     875
     876        # Compute the appropriate elimination ideal.
     877        J = RR.ideal([ Ss[j] - RR(df[j][0]) for j in xrange(m)] +\
     878                 [ Ss[j]**df[j][1] - Ts[j] for j in xrange(m)])
     879        J = J.elimination_ideal(Xs + Ss)
     880
     881        # Coerce J into the polynomial ring in the indeteminates Ts[m:].
     882        # I choose the negdeglex order because i find it useful in my work.
     883        RRR = PolynomialRing(F, [str(t) for t in Ts], order ='negdeglex')
     884        return RRR.ideal(J)
     885
     886    def algebraic_dependence_decomposition(self, whole_and_parts=True):
     887        r"""
     888        Return an algebraic dependence decomposition of ``self`` as a FFPDSum
     889        instance.
     890        Recursive.
     891
     892        EXAMPLES::
     893
     894            sage: from sage.combinat.amgf import *
     895
     896            sage: R.<x, y> = PolynomialRing(QQ)
     897            sage: f = 1/(x^2 * (x*y + 1) * y^3)
     898            sage: ff = FFPD(quotient=f)
     899            sage: decomp = ff.algebraic_dependence_decomposition()
     900            sage: print decomp
     901            [(0, []), (-x, [(x*y + 1, 1)]), (x^2*y^2 - x*y + 1,
     902            [(y, 3), (x, 2)])]
     903            sage: print decomp.sum().quotient() == f
     904            True
     905            sage: for r in decomp:
     906            ...       J = r.algebraic_dependence_certificate()
     907            ...       J is None or J == J.ring().ideal()  # The zero ideal
     908            ...       
     909            True
     910            True
     911            True
     912
     913        ::
     914       
     915            sage: R.<x, y> = PolynomialRing(QQ)
     916            sage: G = sin(x)
     917            sage: H = x^2 * (x*y + 1) * y^3
     918            sage: f = FFPD(G, H.factor())
     919            sage: decomp = f.algebraic_dependence_decomposition()
     920            sage: print decomp
     921            [(0, []), (x^4*y^3*sin(x), [(x*y + 1, 1)]),
     922            (-(x^5*y^5 - x^4*y^4 + x^3*y^3 - x^2*y^2 + x*y - 1)*sin(x),
     923            [(y, 3), (x, 2)])]
     924            sage: if decomp.sum().quotient() == G/H:
     925            ...       print 'yep'
     926            ...       
     927            yep
     928            sage: for r in decomp:
     929            ...       J = r.algebraic_dependence_certificate()
     930            ...       J is None or J == J.ring().ideal()
     931            ...       
     932            True
     933            True
     934            True
     935           
     936        NOTE::
     937
     938        Let $f = p/q$ where $q$ lies in a $d$ -variate polynomial ring
     939        $K[X]$ for some field $K$.
     940        Let $q_1^{e_1} \cdots q_n^{e_n}$ be the
     941        unique factorization of $q$ in $K[X]$ into irreducible factors and
     942        let $V_i$ be the algebraic variety $\{x\in L^d: q_i(x) = 0\}$ of
     943        $q_i$ over the algebraic closure $L$ of $K$.
     944        By [Raic2012]_, $f$ can be written as
     945
     946        (*)   $\sum p_A/\prod_{i \in A} q_i^{b_i}$,
     947
     948        where the $b_i$ are positive integers, each $p_A$ is a products
     949        of $p$ and an element in $K[X]$,
     950        and the sum is taken over all subsets
     951        $A \subseteq \{1, \ldots, m\}$ such that $|A| \le d$ and
     952        $\{q_i : i\in A\}$ is algebraically independent.
     953
     954        I call (*) an *algebraic dependence decomposition* of $f$.
     955        Algebraic dependence decompositions are not unique.
     956
     957        The algorithm used comes from [Raic2012]_.
     958        """       
     959        J = self.algebraic_dependence_certificate()
     960        if not J:
     961            # No decomposing possible.
     962            return FFPDSum([self])
     963           
     964        # Otherwise decompose recursively.
     965        decomp = FFPDSum()
     966        p = self.numerator()
     967        df = self.denominator_factored()
     968        m = len(df)
     969        g = J.gens()[0]     # An annihilating polynomial for df.
     970        new_vars = J.ring().gens()
     971        # Note that each new_vars[j] corresponds to df[j] such that
     972        # g([q**e for q, e in df]) = 0.
     973        # Assuming here that g.parent() has negdeglex term order
     974        # so that g.lt() is indeed the monomial we want below.
     975        # Use g to rewrite r into a sum of FFPDs,
     976        # each with < m distinct denominator factors.
     977        gg = (g.lt() - g)/(g.lc()) 
     978        numers = map(prod, zip(gg.coefficients(), gg.monomials()))
     979        e = list(g.lt().exponents())[0:m]
     980        denoms = [(new_vars[j], e[0][j] + 1) for j in xrange(m)]
     981        # Write r in terms of new_vars,
     982        # cancel factors in the denominator, and combine like terms.
     983        iteration1_temp = FFPDSum([FFPD(a, denoms) for a in numers]).\
     984                          combine_like_terms()
     985        # Substitute in df.
     986        qpowsub = dict([(new_vars[j], df[j][0]**df[j][1])
     987                        for j in xrange(m)])
     988        iteration1 = FFPDSum()
     989        for r in iteration1_temp:
     990            num1 = p*g.parent()(r.numerator()).subs(qpowsub)
     991            denoms1 = []
     992            for q, e in r.denominator_factored():
     993                j = new_vars.index(q)
     994                denoms1.append((df[j][0], df[j][1]*e))
     995            iteration1.append(FFPD(num1, denoms1))
     996        # Now decompose each FFPD of iteration1.
     997        for r in iteration1:
     998            decomp.extend(r.algebraic_dependence_decomposition())
     999
     1000        # Simplify and return result.
     1001        return decomp.combine_like_terms().whole_and_parts()
     1002
     1003    def leinartas_decomposition(self):
     1004        r"""
     1005        Return a Leinartas decomposition of ``self`` as a FFPDSum instance.
     1006
     1007        EXAMPLES::
     1008
     1009            sage: from sage.combinat.amgf import *
     1010
     1011            sage: R.<x, y> = PolynomialRing(QQ)
     1012            sage: f = 1/x + 1/y + 1/(x*y + 1)
     1013            sage: decomp = FFPD(quotient=f).leinartas_decomposition()
     1014            sage: print decomp
     1015            [(0, []), (1, [(x*y + 1, 1)]), (x + y, [(y, 1), (x, 1)])]
     1016            sage: print decomp.sum().quotient() == f
     1017            True
     1018            sage: for r in decomp:
     1019            ...       L = r.nullstellensatz_certificate()
     1020            ...       print L is None
     1021            ...       J = r.algebraic_dependence_certificate()
     1022            ...       print J is None or J == J.ring().ideal()
     1023            ...       
     1024            True
     1025            True
     1026            True
     1027            True
     1028            True
     1029            True
     1030
     1031        ::
     1032       
     1033            sage: R.<x, y> = PolynomialRing(QQ)
     1034            sage: f = sin(x)/x + 1/y + 1/(x*y + 1)
     1035            sage: G = f.numerator()
     1036            sage: H = R(f.denominator())
     1037            sage: ff = FFPD(G, H.factor())
     1038            sage: decomp = ff.leinartas_decomposition()
     1039            sage: print decomp
     1040            [(0, []), (-(x*y^2*sin(x) + x^2*y + x*y + y*sin(x) + x)*y,
     1041            [(y, 1)]), ((x*y^2*sin(x) + x^2*y + x*y + y*sin(x) + x)*x*y,
     1042            [(x*y + 1, 1)]), (x*y^2*sin(x) + x^2*y + x*y + y*sin(x) + x,
     1043            [(y, 1), (x, 1)])]
     1044            sage: if decomp.sum().quotient() == f:
     1045            ...       print 'yep'
     1046            ...       
     1047            yep
     1048            sage: for r in decomp:
     1049            ...       L = r.nullstellensatz_certificate()
     1050            ...       print L is None
     1051            ...       J = r.algebraic_dependence_certificate()
     1052            ...       print J is None or J == J.ring().ideal()
     1053            ...       
     1054            True
     1055            True
     1056            True
     1057            True
     1058            True
     1059            True
     1060            True
     1061            True
     1062
     1063        ::
     1064       
     1065            sage: R.<x, y, z>= PolynomialRing(GF(2, 'a'))
     1066            sage: f = 1/(x * y * z * (x*y + z))
     1067            sage: decomp = FFPD(quotient=f).leinartas_decomposition()
     1068            sage: print decomp
     1069            [(0, []), (1, [(z, 2), (x*y + z, 1)]),
     1070            (1, [(z, 2), (y, 1), (x, 1)])]
     1071            sage: print decomp.sum().quotient() == f
     1072            True
     1073
     1074        NOTE::
     1075
     1076        Let $f = p/q$ where $q$ lies in a $d$ -variate polynomial ring $K[X]$
     1077        for some field $K$.
     1078        Let $q_1^{e_1} \cdots q_n^{e_n}$ be the
     1079        unique factorization of $q$ in $K[X]$ into irreducible factors and
     1080        let $V_i$ be the algebraic variety
     1081        $\{x\in L^d: q_i(x) = 0\}$ of $q_i$ over the algebraic closure
     1082        $L$ of $K$.
     1083        By [Raic2012]_, $f$ can be written as
     1084
     1085        (*)   $\sum p_A/\prod_{i \in A} q_i^{b_i}$,
     1086
     1087        where the $b_i$ are positive integers, each $p_A$ is a product of $p$ and an element of $K[X]$,
     1088        and the sum is taken over all subsets $A \subseteq \{1, \ldots, m\}$ such that
     1089        (1) $|A| \le d$,
     1090        (2) $\cap_{i\in A} T_i \neq \emptyset$, and
     1091        (3) $\{q_i : i\in A\}$ is algebraically independent.
     1092
     1093        In particular, any rational expression in $d$ variables
     1094        can be represented as a sum of rational expressions
     1095        whose denominators each contain at most $d$ distinct irreducible
     1096        factors.
     1097
     1098        I call (*) a *Leinartas decomposition* of $f$.
     1099        Leinartas decompositions are not unique.
     1100
     1101        The algorithm used comes from [Raic2012]_.
     1102        """
     1103        d = self.dimension()
     1104        if d == 1:
     1105            # Sage's lift() function doesn't work in
     1106            # univariate polynomial rings.
     1107            # So nullstellensatz_decomposition() won't work.
     1108            # So use algebraic_dependence_decomposition(),
     1109            # which is sufficient.
     1110            temp = FFPDSum([self])
     1111        else:
     1112            temp = self.nullstellensatz_decomposition()
     1113        decomp = FFPDSum()
     1114        for r in temp:
     1115            decomp.extend(r.algebraic_dependence_decomposition())
     1116
     1117        # Simplify and return result.
     1118        return decomp.combine_like_terms().whole_and_parts()
     1119
     1120    def cohomology_decomposition(self):
     1121        r"""
     1122        Let $p/(q_1^{e_1} \cdots q_n^{e_n})$ be the fraction represented
     1123        by ``self`` and let $K[x_1, \ldots, x_d]$ be the polynomial ring
     1124        in which the $q_i$ lie.
     1125        Assume that $n \le d$ and that the gradients of the $q_i$ are linearly
     1126        independent at all points in the intersection
     1127        $V_1 \cap \ldots \cap V_n$ of the algebraic varieties
     1128        $V_i = \{x \in L^d: q_i(x) = 0 \}$, where $L$ is the algebraic closure
     1129        of the field $K$.
     1130        Return a FFPDSum $f$ such that the differential form
     1131        $f dx_1 \wedge \cdots \wedge dx_d$ is de Rham cohomologous to the
     1132        differential form
     1133        $p/(q_1^{e_1} \cdots q_n^{e_n}) dx_1 \wedge \cdots \wedge dx_d$
     1134        and such that the denominator of each summand of $f$ contains
     1135        no repeated irreducible factors.
     1136 
     1137        EXAMPLES::
     1138
     1139            sage: from sage.combinat.amgf import *
     1140
     1141            sage: R.<x, y>= PolynomialRing(QQ)
     1142            sage: print FFPD(1, [(x, 1), (y, 2)]).cohomology_decomposition()
     1143            [(0, [])]
     1144           
     1145        ::
     1146
     1147            sage: R.<x, y>= PolynomialRing(QQ)
     1148            sage: p = 1
     1149            sage: qs = [(x*y - 1, 1), (x**2 + y**2 - 1, 2)]
     1150            sage: f = FFPD(p, qs)
     1151            sage: print f.cohomology_decomposition()
     1152            [(0, []), (4/3*x*y + 4/3, [(x^2 + y^2 - 1, 1)]),
     1153            (1/3, [(x*y - 1, 1), (x^2 + y^2 - 1, 1)])]
     1154           
     1155        NOTES:
     1156
     1157        The algorithm used here comes from the proof of Theorem 17.4 of
     1158        [AiYu1983]_.         
     1159
     1160        AUTHORS:
     1161
     1162        - Alexander Raichev (2008-10-01, 2011-01-15, 2012-07-31)
     1163        """
     1164        R = self.ring()       
     1165        df = self.denominator_factored()
     1166        n = len(df)
     1167        if R is None or sum([e for (q, e) in df]) <= n:
     1168            # No decomposing possible.
     1169            return FFPDSum([self])
     1170       
     1171        # Otherwise decompose recursively.   
     1172        decomp = FFPDSum()       
     1173        p = self.numerator()
     1174        qs = [q for (q, e) in df]
     1175        X = sorted(R.gens())
     1176        var_sets_n = Set(X).subsets(n)
     1177
     1178        # Compute Jacobian determinants for qs.
     1179        dets = []
     1180        for v in var_sets_n:
     1181            # Sort v according to the term order of R.
     1182            x = sorted(v)
     1183            jac = jacobian(qs, x)
     1184            dets.append(R(jac.determinant()))
     1185                   
     1186        # Get a Nullstellensatz certificate for qs and dets.
     1187        L = R(1).lift(R.ideal(qs + dets))
     1188                 
     1189        # Do first iteration of decomposition.
     1190        iteration1 = FFPDSum()
     1191        # Contributions from qs.
     1192        for i in xrange(n):
     1193            if L[i] == 0:
     1194                continue
     1195            # Cancel one df[i] from denominator.
     1196            new_df = [list(t) for t in df]
     1197            if new_df[i][1] > 1:
     1198                new_df[i][1] -= 1
     1199            else:
     1200                del(new_df[i])
     1201            iteration1.append(FFPD(p*L[i], new_df))
     1202               
     1203        # Contributions from dets.
     1204        # Compute each contribution's cohomologous form using
     1205        # the least index j such that new_df[j][1] > 1.
     1206        # Know such an index exists by first 'if' statement at
     1207        # the top.
     1208        for j in xrange(n):
     1209            if df[j][1] > 1:
     1210                J = j
     1211                break
     1212        new_df = [list(t) for t in df]
     1213        new_df[J][1] -= 1
     1214        for k in xrange(var_sets_n.cardinality()):
     1215            if L[n + k] == 0:
     1216                continue
     1217            # Sort variables according to the term order of R.
     1218            x = sorted(var_sets_n[k])
     1219            # Compute Jacobian in the Symbolic Ring.
     1220            jac = jacobian([SR(p*L[n + k])] +
     1221                           [SR(qs[j]) for j in xrange(n) if j != J],
     1222                           [SR(xx) for xx in x])
     1223            det = jac.determinant()             
     1224            psign = FFPD.permutation_sign(x, X)           
     1225            iteration1.append(FFPD((-1)**J*det/\
     1226                                   (psign*new_df[J][1]),
     1227                                   new_df))
     1228           
     1229        # Now decompose each FFPD of iteration1.
     1230        for r in iteration1:
     1231            decomp.extend(r.cohomology_decomposition())
     1232           
     1233        # Simplify and return result.
     1234        return decomp.combine_like_terms().whole_and_parts()
     1235       
     1236    @staticmethod
     1237    def permutation_sign(s, u):
     1238        r"""
     1239        This function returns the sign of the permutation on
     1240        ``1, ..., len(u)`` that is induced by the sublist
     1241        ``s`` of ``u``.
     1242        For internal use by cohomology_decomposition().
     1243
     1244        INPUT:
     1245
     1246        - ``s`` - A sublist of ``u``.
     1247        - ``u`` - A list.
     1248
     1249        OUTPUT:
     1250
     1251        The sign of the permutation obtained by taking indices
     1252        within ``u`` of the list ``s + sc``, where ``sc`` is ``u``
     1253        with the elements of ``s`` removed.
     1254
     1255        EXAMPLES::
     1256
     1257            sage: from sage.combinat.amgf import *
     1258
     1259            sage: u = ['a','b','c','d','e']
     1260            sage: s = ['b','d']
     1261            sage: FFPD.permutation_sign(s, u)
     1262            -1
     1263            sage: s = ['d','b']
     1264            sage: FFPD.permutation_sign(s, u)
     1265            1
     1266
     1267        AUTHORS:
     1268
     1269        - Alexander Raichev (2008-10-01, 2012-07-31)
     1270        """
     1271        # Convert lists to lists of numbers in {1,..., len(u)}
     1272        A = [i + 1  for i in xrange(len(u))]
     1273        B = [u.index(x) + 1 for x in s]
     1274       
     1275        C = sorted(list(Set(A).difference(Set(B))))
     1276        P = Permutation(B + C)
     1277        return P.signature()
     1278
     1279    def asymptotic_decomposition(self, alpha, asy_var=None):
     1280        r"""
     1281        Return a FFPDSum that has the same asymptotic expansion as ``self``
     1282        in the direction ``alpha`` but each of whose FFPDs has a
     1283        denominator factorization of the form $[(q_1, 1), \ldots, (q_n, 1)]$,
     1284        where ``n`` is at most ``d = self.dimension()``.
     1285        The output results from a Leinartas decomposition followed by a
     1286        cohomology decomposition.
     1287
     1288        INPUT:
     1289       
     1290        - ``alpha`` - A d-tuple of positive integers or symbolic variables.
     1291        - ``asy_var`` - (Optional; default=None) A symbolic variable with
     1292          respect to which to compute asymptotics.
     1293          If None is given the set ``asy_var = var('r')``.
     1294       
     1295        EXAMPLES::
     1296       
     1297            sage: from sage.combinat.amgf import *
     1298            sage: R.<x, y>= PolynomialRing(QQ)
     1299            sage: H = (1 - 2*x -y)*(1 - x -2*y)**2
     1300            sage: Hfac = H.factor()
     1301            sage: G = 1/Hfac.unit()
     1302            sage: F = FFPD(G, Hfac)
     1303            sage: alpha = var('a, b')
     1304            sage: r = var('r')
     1305            sage: print F.asymptotic_decomposition(alpha, r) # long time
     1306            [(0, []), (-1/3*(a*y - 2*b*x)*r/(x*y) + 1/3*(2*x - y)/(x*y),
     1307            [(x + 2*y - 1, 1), (2*x + y - 1, 1)])]
     1308           
     1309        AUTHORS:
     1310       
     1311        - Alexander Raichev (2012-08-01)
     1312        """
     1313        R = self.ring()
     1314        if R is None:
     1315            return None
     1316           
     1317        d = self.dimension()       
     1318        n = len(self.denominator_factored())
     1319        X = [SR(x) for x in R.gens()]
     1320
     1321        # Reduce number of distinct factors in denominator of self
     1322        # down to at most d.
     1323        decomp1 = FFPDSum([self])
     1324        if n > d:
     1325            decomp1 = decomp1[0].leinartas_decomposition()
     1326                     
     1327        # Reduce to no repeated factors in denominator of each element
     1328        # of decomp1.
     1329        # Compute the cohomology decomposition for each
     1330        # Cauchy differential form generated by each element of decomp.
     1331        if asy_var is None:
     1332            asy_var = var('r')
     1333        cauchy_stuff = prod([X[j]**(-alpha[j]*asy_var - 1) for j in xrange(d)])
     1334        decomp2 = FFPDSum()
     1335        for f in decomp1:
     1336            ff = FFPD(f.numerator()*cauchy_stuff,
     1337                      f.denominator_factored())
     1338            decomp2.extend(ff.cohomology_decomposition())           
     1339        decomp2 = decomp2.combine_like_terms()       
     1340       
     1341        # Divide out cauchy_stuff from integrands.
     1342        decomp3 = FFPDSum()
     1343        for f in decomp2:
     1344            ff = FFPD((f.numerator()/cauchy_stuff).\
     1345                      simplify_full().collect(asy_var),
     1346                      f.denominator_factored())
     1347            decomp3.append(ff)
     1348       
     1349        return decomp3
     1350
     1351    def asymptotics(self, p, alpha, N, asy_var=None, numerical=0):
     1352        r"""
     1353        Return the first $N$ terms (some of which could be zero)
     1354        of the asymptotic expansion of the Maclaurin ray coefficients
     1355        $F_{r\alpha}$ of the function $F$ represented by ``self``
     1356        as $r\to\infty$, where $r$ = ``asy_var`` and ``alpha`` is a tuple of
     1357        positive integers of length ``d = self.dimension()``.
     1358        Assume that $F$ is holomorphic in a neighborhood of the origin,
     1359        that the denominator factorization of ``self`` is also the unique
     1360        factorization of the denominator of $F$ in the local analytic ring
     1361        at $p$ (not just in the polynomial ring ``self.ring()``),
     1362        that $p$ is a convenient strictly minimal smooth or multiple point
     1363        with all nonzero coordinates that is critical and nondegenerate
     1364        for ``alpha``.
     1365   
     1366        INPUT:
     1367   
     1368        - ``p`` - A dictionary with keys that can be coerced to equal
     1369          ``self.ring().gens()``.
     1370        - ``alpha`` - A tuple of length ``self.dimension()`` of
     1371          positive integers or, if $p$ is a smooth point,
     1372          possibly of symbolic variables.
     1373        - ``N`` - A positive integer.
     1374        - ``numerical`` - (Optional; default=0) A natural number.
     1375          If numerical > 0, then return a numerical approximation
     1376          of $F_{r \alpha}$ with ``numerical`` digits of precision.
     1377          Otherwise return exact values.
     1378        - ``asy_var`` - (Optional; default=None) A  symbolic variable.
     1379          The variable of the asymptotic expansion.
     1380          If none is given, ``var('r')`` will be assigned.
     1381   
     1382        OUTPUT:
     1383   
     1384        The tuple (asy, exp_scale, subexp_part).
     1385        Here asy is the sum of the first $N$ terms (some of which might be 0)
     1386        of the asymptotic expansion of $F_{r\alpha}$ as $r\to\infty$;   
     1387        exp_scale**r is the exponential factor of asy;
     1388        subexp_part is the subexponential factor of asy. 
     1389   
     1390        EXAMPLES::
     1391   
     1392            sage: from sage.combinat.amgf import *
     1393           
     1394        A smooth point example::
     1395   
     1396            sage: R.<x,y>= PolynomialRing(QQ)
     1397            sage: H = (1 - x - y - x*y)**2
     1398            sage: Hfac = H.factor()
     1399            sage: G = 1/Hfac.unit()
     1400            sage: F = FFPD(G, Hfac)
     1401            sage: print F           
     1402            (1, [(x*y + x + y - 1, 2)])
     1403            sage: alpha = [4, 3]
     1404            sage: p = {y: 1/3, x: 1/2}
     1405            sage: asy = F.asymptotics(p, alpha, 2) # long time
     1406            Creating auxiliary functions...
     1407            Computing derivatives of auxiallary functions...
     1408            Computing derivatives of more auxiliary functions...
     1409            Computing second order differential operator actions...
     1410            sage: print asy # long time
     1411            (1/6000*(3600*sqrt(2)*sqrt(3)*sqrt(5)*sqrt(r)/sqrt(pi) +
     1412            463*sqrt(2)*sqrt(3)*sqrt(5)/(sqrt(pi)*sqrt(r)))*432^r, 432,
     1413            1/6000*(3600*sqrt(5)*r +
     1414            463*sqrt(5))*sqrt(2)*sqrt(3)/(sqrt(pi)*sqrt(r)))
     1415            sage: print F.relative_error(asy[0], alpha, [1, 2, 4, 8, 16], asy[1]) # long time
     1416            Calculating errors table in the form
     1417                exponent, scaled Maclaurin coefficient, scaled asymptotic
     1418                values, relative errors...
     1419                [((4, 3), 2.083333333, [2.092576110], [-0.004436533009]),
     1420                ((8, 6), 2.787374614, [2.790732875], [-0.001204811281]),
     1421                ((16, 12), 3.826259447, [3.827462310], [-0.0003143703383]),
     1422                ((32, 24), 5.328112821, [5.328540787], [-0.00008032229296]),
     1423                ((64, 48), 7.475927885, [7.476079664], [-0.00002030233658])]
     1424   
     1425        A multiple point example::
     1426   
     1427            sage: R.<x,y,z>= PolynomialRing(QQ)
     1428            sage: H = (4 - 2*x - y - z)**2*(4 - x - 2*y - z)
     1429            sage: Hfac = H.factor()
     1430            sage: G = 16/Hfac.unit()
     1431            sage: F = FFPD(G, Hfac)
     1432            sage: print F
     1433            (-16, [(x + 2*y + z - 4, 1), (2*x + y + z - 4, 2)])
     1434            sage: alpha = [3, 3, 2]
     1435            sage: p = {x: 1, y: 1, z: 1}
     1436            sage: asy = F.asymptotics(p, alpha, 2) # long time
     1437            Creating auxiliary functions...
     1438            Computing derivatives of auxiliary functions...
     1439            Computing derivatives of more auxiliary functions...
     1440            Computing second-order differential operator actions...
     1441            sage: print asy # long time
     1442            (4/3*sqrt(3)*sqrt(r)/sqrt(pi) +
     1443            47/216*sqrt(3)/(sqrt(pi)*sqrt(r)), 1,
     1444            1/216*(288*sqrt(3)*r + 47*sqrt(3))/(sqrt(pi)*sqrt(r)))
     1445            sage: print F.relative_error(asy[0], alpha, [1, 2, 4, 8], asy[1]) # long time
     1446            Calculating errors table in the form
     1447            exponent, scaled Maclaurin coefficient, scaled asymptotic values,
     1448            relative errors...
     1449            [((3, 3, 2), 0.9812164307, [1.515572606], [-0.5445854340]),
     1450            ((6, 6, 4),
     1451            1.576181132, [1.992989399], [-0.2644418580]),
     1452            ((12, 12, 8), 2.485286378,
     1453            [2.712196351], [-0.09130133851]), ((24, 24, 16), 3.700576827,
     1454            [3.760447895], [-0.01617884750])]
     1455       
     1456        NOTES:
     1457   
     1458        The algorithms used here come from [RaWi2008a]_ and [RaWi2012]_.
     1459   
     1460        AUTHORS:
     1461   
     1462        - Alexander Raichev (2008-10-01, 2010-09-28, 2011-04-27, 2012-08-03)
     1463        """
     1464        R = self.ring()
     1465        if R is None:
     1466            return None
     1467           
     1468        # Coerce keys of p into R.
     1469        p = FFPD.coerce_point(R, p)
     1470       
     1471        if asy_var is None:
     1472            asy_var = var('r')   
     1473        d = self.dimension()
     1474        X = list(R.gens())   
     1475        alpha = list(alpha)
     1476        df = self.denominator_factored()
     1477        # Find greatest i such that X[i] is a convenient coordinate,
     1478        # that is, such that for all (h, e) in df, we have
     1479        # (X[i]*diff(h, X[i])).subs(p) != 0. 
     1480        # Assuming such an i exists.
     1481        i = d - 1
     1482        while 0 in [(X[i]*diff(h, X[i])).subs(p) for (h, e) in df]:
     1483            i -= 1
     1484        coordinate = i
     1485           
     1486        # Decompose self into a sum of simpler pieces.
     1487        # Where each piece has a denominator with at most d distinct
     1488        # factors and no repeated factors.
     1489        decomp = self.asymptotic_decomposition(alpha, asy_var)
     1490               
     1491        # Sum asymptotic expansions of terms of decomp.
     1492        asy_expansions = []
     1493        for f in decomp:
     1494            if f.quotient() == 0:
     1495                continue
     1496            critical_cone_p = f.critical_cone(p, coordinate)   
     1497            if [alpha[i] for i in xrange(d)] == [True for i in xrange(d)] and\
     1498              alpha not in critical_cone_p[1]:
     1499                # f does not contribute to asymptotics of self
     1500                # in direction alpha.
     1501                continue   
     1502            n = len(f.denominator_factored())
     1503            if n == 1:
     1504                # Smooth point.
     1505                asy_expansions.append(
     1506                f.asymptotics_smooth(p, alpha, N, asy_var,
     1507                                     coordinate, numerical))
     1508            else:
     1509                # Multiple point with n <= d.
     1510                asy_expansions.append(
     1511                f.asymptotics_multiple(p, alpha, N, asy_var,
     1512                                       coordinate, numerical))
     1513        if asy_expansions:
     1514            asy = sum([a[0] for a in asy_expansions])
     1515            exp_scale = asy_expansions[0][1] # Same for all.
     1516            subexp_part = sum([a[2] for a in asy_expansions]).\
     1517                               simplify_full()
     1518            return (asy, exp_scale, subexp_part)
     1519        else:
     1520            return None
     1521           
     1522    def asymptotics_smooth(self, p, alpha, N, asy_var, coordinate=None,
     1523                           numerical=0):
     1524        r"""
     1525        Does what asymptotics() does but only in the case of a
     1526        convenient smooth point and assuming that the denominator of ``self``
     1527        contains no repeated factors.
     1528
     1529        INPUT:
     1530   
     1531        - ``p`` - A dictionary with keys that can be coerced to equal
     1532          ``self.ring().gens()``.
     1533        - ``alpha`` - A tuple of length ``d = self.dimension()`` of
     1534          positive integers or, if $p$ is a smooth point,
     1535          possibly of symbolic variables.
     1536        - ``N`` - A positive integer.
     1537        - ``asy_var`` - (Optional; default=None) A symbolic variable.
     1538          The variable of the asymptotic expansion.
     1539          If none is given, ``var('r')`` will be assigned.
     1540        - ``coordinate``- (Optional; default=None) An integer in
     1541          {0, ..., d-1} indicating a convenient coordinate to base
     1542          the asymptotic calculations on.
     1543          If None is assigned, then choose ``coordinate=d-1``.
     1544        - ``numerical`` - (Optional; default=0) A natural number.
     1545          If numerical > 0, then return a numerical approximation of the
     1546          Maclaurin ray coefficients of ``self`` with ``numerical`` digits
     1547          of precision. 
     1548          Otherwise return exact values.
     1549       
     1550        NOTES:
     1551
     1552        The formulas used for computing the asymptotic expansions are
     1553        Theorems 3.2 and 3.3 [RaWi2008a]_ with the exponent of H equal to 1.
     1554        Theorem 3.2 is a specialization of Theorem 3.4 of [RaWi2012]_
     1555        with $n=1$.
     1556
     1557        EXAMPLES::
     1558       
     1559            sage: from sage.combinat.amgf import *
     1560            sage: R.<x, y>= PolynomialRing(QQ)
     1561            sage: H = 1-x-y-x*y
     1562            sage: Hfac = H.factor()
     1563            sage: G = 1/Hfac.unit()
     1564            sage: F = FFPD(G, Hfac)
     1565            sage: alpha = [3, 2]
     1566            sage: p = {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3}
     1567            sage: print F.asymptotics_smooth(p, alpha, 2, var('r'), numerical=3) # long time
     1568            Creating auxiliary functions...
     1569            Computing derivatives of auxiallary functions...
     1570            Computing derivatives of more auxiliary functions...
     1571            Computing second order differential operator actions...
     1572            ((0.369/sqrt(r) - 0.0186/r^(3/2))*71.2^r, 71.2,
     1573            0.369/sqrt(r) - 0.0186/r^(3/2))
     1574           
     1575        ::
     1576       
     1577            sage: R.<x, y> = PolynomialRing(QQ)
     1578            sage: q = 1/2
     1579            sage: qq = q.denominator()
     1580            sage: H = 1 - q*x + q*x*y - x^2*y
     1581            sage: Hfac = H.factor()
     1582            sage: G = (1 - q*x)/Hfac.unit()
     1583            sage: F = FFPD(G, Hfac)
     1584            sage: alpha = list(qq*vector([2, 1 - q]))
     1585            sage: print alpha
     1586            [4, 1]
     1587            sage: p = {x: 1, y: 1}
     1588            sage: print F.asymptotics_smooth(p, alpha, 5, var('r')) # long time
     1589            Creating auxiliary functions...
     1590            Computing derivatives of auxiallary functions...
     1591            Computing derivatives of more auxiliary functions...
     1592            Computing second order differential operator actions...
     1593            (1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*r^(1/3)) -
     1594            1/96*2^(1/3)*sqrt(3)*gamma(2/3)/(pi*r^(5/3)), 1,
     1595            1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*r^(1/3)) -
     1596            1/96*2^(1/3)*sqrt(3)*gamma(2/3)/(pi*r^(5/3)))
     1597           
     1598        AUTHORS:
     1599
     1600        - Alexander Raichev (2008-10-01, 2010-09-28, 2012-08-01)
     1601        """
     1602        R = self.ring()
     1603        if R is None:
     1604            return None
     1605               
     1606        d = self.dimension()
     1607        I = sqrt(-Integer(1))
     1608        # Coerce everything into the Symbolic Ring.
     1609        X = [SR(x) for x in R.gens()]
     1610        G = SR(self.numerator())
     1611        H = SR(self.denominator())
     1612        p = dict([(SR(x), p[x]) for x in R.gens()])
     1613        alpha = [SR(a) for a in alpha]
     1614       
     1615        # Put given convenient coordinate at end of variable list.
     1616        if coordinate is not None:
     1617            x = X.pop(coordinate)
     1618            X.append(x)
     1619            a = alpha.pop(coordinate)
     1620            alpha.append(a)
     1621       
     1622        # If p is a tuple of rationals, then compute with it directly.
     1623        # Otherwise, compute symbolically and plug in p at the end.
     1624        if vector(p.values()) in QQ**d:
     1625            P = p
     1626        else:
     1627            sP = [var('p' + str(j)) for j in xrange(d)]
     1628            P = dict( [(X[j], sP[j]) for j in xrange(d)] )
     1629            p = dict( [(sP[j], p[X[j]]) for j in xrange(d)] )
     1630
     1631        # Setup.
     1632        print "Creating auxiliary functions..."
     1633        # Implicit functions.
     1634        h = function('h', *tuple(X[:d - 1]))
     1635        U = function('U', *tuple(X))       
     1636        # All other functions are defined in terms of h, U, and
     1637        # explicit functions.
     1638        Gcheck = -G/U * (h/X[d - 1])
     1639        A = Gcheck.subs({X[d - 1]: Integer(1)/h})/h
     1640        t = 't'
     1641        while t in [str(x) for x in X]:
     1642            t = t + 't'
     1643        T = [var(t + str(i)) for i in xrange(d - 1)]
     1644        e = dict([(X[i], P[X[i]]*exp(I*T[i])) for i in xrange(d - 1)])
     1645        ht = h.subs(e)
     1646        Ht = H.subs(e).subs({X[d - 1]: Integer(1)/ht})
     1647        At = A.subs(e) 
     1648        Phit = -log(P[X[d - 1]]*ht) + \
     1649               I * add([alpha[i]/alpha[d - 1]*T[i] for i in xrange(d - 1)])
     1650        Tstar = dict([(t, Integer(0)) for t in T])
     1651        # Store h and U and all their derivatives evaluated at P.
     1652        atP = P.copy()
     1653        atP.update({h.subs(P): Integer(1)/P[X[d - 1]]})
     1654
     1655        # Compute the derivatives of h up to order 2*N, evaluate at P, 
     1656        # and store in atP. 
     1657        # Keep a copy of unevaluated h derivatives for use in the case
     1658        # d = 2 and v > 2 below.
     1659        hderivs1 = {}   # First derivatives of h.
     1660        for i in xrange(d - 1):
     1661            s = solve( diff(H.subs({X[d - 1]: Integer(1)/h}), X[i]),
     1662                      diff(h, X[i]))[0].rhs().simplify()
     1663            hderivs1.update({diff(h, X[i]): s})
     1664            atP.update({diff(h, X[i]).subs(P): s.subs(P).subs(atP)})
     1665        hderivs = FFPD.diff_all(h, X[0: d - 1], 2*N, sub=hderivs1, rekey=h)
     1666        for k in hderivs.keys():
     1667            atP.update({k.subs(P):hderivs[k].subs(atP)})
     1668   
     1669        # Compute the derivatives of U up to order 2*N and evaluate at P. 
     1670        # To do this, differentiate H = U*Hcheck over and over, evaluate at P,
     1671        # and solve for the derivatives of U at P.
     1672        # Need the derivatives of H with short keys to pass on
     1673        # to diff_prod later.
     1674        Hderivs = FFPD.diff_all(H, X, 2*N, ending=[X[d - 1]], sub_final=P)
     1675        print "Computing derivatives of auxiallary functions..."
     1676        # For convenience in checking if all the nontrivial derivatives of U
     1677        # at p are zero a few line below, store the value of U(p) in atP
     1678        # instead of in Uderivs.
     1679        Uderivs ={}
     1680        atP.update({U.subs(P): diff(H, X[d - 1]).subs(P)})
     1681        end = [X[d - 1]]
     1682        Hcheck = X[d - 1] - Integer(1)/h
     1683        k = H.polynomial(CC).degree() - 1
     1684        if k == 0:
     1685            # Then we can conclude that all higher derivatives of U are zero.
     1686            for l in xrange(1, 2*N + 1):
     1687                for s in UnorderedTuples(X, l):
     1688                    Uderivs[diff(U, s).subs(P)] = Integer(0)
     1689        elif k > 0 and k < 2*N:
     1690            all_zero = True
     1691            Uderivs =  FFPD.diff_prod(Hderivs, U, Hcheck, X,
     1692                                     range(1, k + 1), end, Uderivs, atP)
     1693            # Check for a nonzero U derivative.
     1694            if Uderivs.values() != [Integer(0)  for i in xrange(len(Uderivs))]:
     1695                all_zero = False
     1696            if all_zero:
     1697                # Then, using a proposition at the end of [RaWi2012], we can
     1698                # conclude that all higher derivatives of U are zero.
     1699                for l in xrange(k + 1, 2*N +1):
     1700                    for s in UnorderedTuples(X, l):
     1701                        Uderivs.update({diff(U, s).subs(P): Integer(0)})
     1702            else:
     1703                # Have to compute the rest of the derivatives.
     1704                Uderivs = FFPD.diff_prod(Hderivs, U, Hcheck, X,
     1705                                         range(k + 1, 2*N + 1), end, Uderivs,
     1706                                         atP)
     1707        else:
     1708            Uderivs = FFPD.diff_prod(Hderivs, U, Hcheck, X,
     1709                                     range(1, 2*N + 1), end, Uderivs, atP)
     1710        atP.update(Uderivs)
     1711
     1712        # In general, this algorithm is not designed to handle the case of a
     1713        # singular Phit''(Tstar). 
     1714        # However, when d = 2 the algorithm can cope.
     1715        if d == 2:
     1716            # Compute v, the order of vanishing at Tstar of Phit.
     1717            # It is at least 2.
     1718            v = Integer(2)
     1719            Phitderiv = diff(Phit, T[0], 2)
     1720            splat = Phitderiv.subs(Tstar).subs(atP).subs(p).simplify()
     1721            while splat == 0:
     1722                v += 1
     1723                if v > 2*N: 
     1724                    # Then need to compute more derivatives of h for atP.
     1725                    hderivs.update({diff(h, X[0], v):
     1726                                    diff(hderivs[diff(h, X[0], v - 1)],
     1727                                    X[0]).subs(hderivs1)})
     1728                    atP.update({diff(h, X[0], v).subs(P):       
     1729                                hderivs[diff(h, X[0], v)].subs(atP)})
     1730                Phitderiv = diff(Phitderiv, T[0])
     1731                splat = Phitderiv.subs(Tstar).subs(atP).subs(p).simplify()   
     1732        if d == 2 and v > 2:     
     1733            t = T[0]  # Simplify variable names.
     1734            a = splat/factorial(v)
     1735            Phitu = Phit - a*t**v
     1736   
     1737            # Compute all partial derivatives of At and Phitu
     1738            # up to orders 2*(N - 1) and 2*(N - 1) + v, respectively,
     1739            # in case v is even.
     1740            # Otherwise, compute up to orders N - 1 and N - 1 + v,
     1741            # respectively. 
     1742            # To speed up later computations,
     1743            # create symbolic functions AA and BB
     1744            # to stand in for the expressions At and Phitu, respectively.
     1745            print "Computing derivatives of more auxiliary functions..."
     1746            AA = function('AA', t)
     1747            BB = function('BB', t)
     1748            if v.mod(2) == 0:
     1749                At_derivs = FFPD.diff_all(At, T, 2*N - 2,
     1750                                          sub=hderivs1, sub_final=[Tstar, atP],
     1751                                          rekey=AA)
     1752                Phitu_derivs = FFPD.diff_all(Phitu, T, 2*N - 2 +v,
     1753                                             sub=hderivs1,
     1754                                             sub_final=[Tstar, atP],
     1755                                             zero_order=v + 1, rekey=BB)
     1756            else:
     1757                At_derivs = FFPD.diff_all(At, T, N - 1, sub=hderivs1,
     1758                                          sub_final=[Tstar, atP], rekey=AA)
     1759                Phitu_derivs = FFPD.diff_all(Phitu, T, N - 1 + v,
     1760                                             sub=hderivs1,
     1761                                             sub_final=[Tstar, atP],
     1762                                             zero_order=v + 1 , rekey=BB)
     1763            AABB_derivs = At_derivs
     1764            AABB_derivs.update(Phitu_derivs)
     1765            AABB_derivs[AA] = At.subs(Tstar).subs(atP)
     1766            AABB_derivs[BB] = Phitu.subs(Tstar).subs(atP)
     1767            print "Computing second order differential operator actions..."
     1768            DD = FFPD.diff_op_simple(AA, BB, AABB_derivs, t, v, a, N)
     1769           
     1770            # Plug above into asymptotic formula.
     1771            L = []
     1772            if v.mod(2) == 0:
     1773                for k in xrange(N):
     1774                    L.append(add([
     1775                    (-1)**l * gamma((2*k + v*l + 1)/v)/\
     1776                    (factorial(l) * factorial(2*k + v*l))*\
     1777                    DD[(k, l)] for l in xrange(0, 2*k + 1) ]))
     1778                chunk = a**(-1/v)/(pi*v)*add([
     1779                        alpha[d - 1 ]**(-(2*k + 1)/v)*\
     1780                        L[k]*asy_var**(-(2*k + 1)/v) for k in xrange(N) ])
     1781            else:
     1782                zeta = exp(I*pi/(2*v))
     1783                for k in xrange(N):
     1784                    L.append(add([
     1785                    (-1)**l*gamma((k + v*l + 1)/v)/\
     1786                    (factorial(l)*factorial(k + v*l))*\
     1787                    (zeta**(k + v*l + 1) +\
     1788                    (-1)**(k + v*l)*zeta**(-(k + v*l + 1)))*\
     1789                    DD[(k, l)] for l in xrange(0, k + 1) ]))
     1790                chunk = abs(a)**(-1/v)/(2*pi*v)*add([
     1791                        alpha[d - 1]**(-(k + 1)/v)*\
     1792                        L[k] *asy_var**(-(k + 1)/v) for k in xrange(N) ])
     1793                       
     1794        # Asymptotics for d >= 2 case. 
     1795        # A singular Phit''(Tstar) will cause a crash in this case.
     1796        else:
     1797            Phit1 = jacobian(Phit, T).subs(hderivs1)
     1798            a = jacobian(Phit1, T).subs(hderivs1).subs(Tstar).subs(atP)
     1799            a_inv = a.inverse()
     1800            Phitu = Phit - (Integer(1)/Integer(2))*matrix([T])*\
     1801                    a*matrix([T]).transpose()
     1802            Phitu = Phitu[0][0]
     1803            # Compute all partial derivatives of At and Phitu up to
     1804            # orders 2*N-2 and 2*N, respectively. 
     1805            # Take advantage of the fact that At and Phitu
     1806            # are sufficiently differentiable functions so that mixed partials
     1807            # are equal.  Thus only need to compute representative partials.
     1808            # Choose nondecreasing sequences as representative differentiation-
     1809            # order sequences.
     1810            # To speed up later computations,
     1811            # create symbolic functions AA and BB
     1812            # to stand in for the expressions At and Phitu, respectively.
     1813            print "Computing derivatives of more auxiliary functions..."
     1814            AA = function('AA', *tuple(T))
     1815            At_derivs = FFPD.diff_all(At, T, 2*N - 2, sub=hderivs1,
     1816                                      sub_final =[Tstar, atP], rekey=AA)
     1817            BB = function('BB', *tuple(T))
     1818            Phitu_derivs = FFPD.diff_all(Phitu, T, 2*N, sub=hderivs1,
     1819                                         sub_final =[Tstar, atP], rekey=BB,
     1820                                         zero_order=3)
     1821            AABB_derivs = At_derivs
     1822            AABB_derivs.update(Phitu_derivs)
     1823            AABB_derivs[AA] = At.subs(Tstar).subs(atP)
     1824            AABB_derivs[BB] = Phitu.subs(Tstar).subs(atP)
     1825            print "Computing second order differential operator actions..."
     1826            DD = FFPD.diff_op(AA, BB, AABB_derivs, T, a_inv, 1 , N)
     1827   
     1828            # Plug above into asymptotic formula.
     1829            L =[]
     1830            for k in xrange(N):
     1831                L.append(add([
     1832                         DD[(0, k, l)]/((-1)**k*2**(l+k)*\
     1833                         factorial(l)*factorial(l+k))
     1834                         for l in xrange(0, 2*k + 1) ]))
     1835            chunk = add([ (2*pi)**((1 - d)/Integer(2))*\
     1836                    a.determinant()**(-Integer(1)/Integer(2))*\
     1837                    alpha[d - 1]**((Integer(1) - d)/Integer(2) - k)*L[k]*\
     1838                    asy_var**((Integer(1) - d)/Integer(2) - k)
     1839                    for k in xrange(N) ])
     1840         
     1841        chunk = chunk.subs(p).simplify()
     1842        coeffs = chunk.coefficients(asy_var)
     1843        coeffs.reverse()
     1844        coeffs = coeffs[:N]
     1845        if numerical:
     1846            subexp_part = add([co[0].subs(p).n(digits=numerical)*\
     1847                          asy_var**co[1] for co in coeffs])
     1848            exp_scale = prod([(P[X[i]]**(-alpha[i])).subs(p)
     1849                              for i in xrange(d)]).n(digits=numerical)
     1850        else:
     1851            subexp_part = add([co[0].subs(p)*asy_var**co[1] for co in coeffs])
     1852            exp_scale = prod([(P[X[i]]**(-alpha[i])).subs(p)
     1853                             for i in xrange(d)])
     1854        return (exp_scale**asy_var*subexp_part, exp_scale, subexp_part)
     1855
     1856    def asymptotics_multiple(self, p, alpha, N, asy_var, coordinate=None,
     1857                             numerical=0):
     1858        r"""
     1859        Does what asymptotics() does but only in the case of a
     1860        convenient multiple point nondegenerate for alpha and assuming that
     1861        that the number of distinct irreducible factors of the denominator of
     1862        ``self`` is at most ``self.dimension()`` and that no factors
     1863        are repeated.
     1864        Assume also that ``p.values()`` are not symbolic variables.
     1865 
     1866        INPUT:
     1867       
     1868        - ``p`` - A dictionary with keys that can be coerced to equal
     1869          ``self.ring().gens()``.
     1870        - ``alpha`` - A tuple of length ``d = self.dimension()`` of
     1871          positive integers or, if $p$ is a smooth point,
     1872          possibly of symbolic variables.
     1873        - ``N`` - A positive integer.
     1874        - ``asy_var`` - (Optional; default=None) A symbolic variable.
     1875          The variable of the asymptotic expansion.
     1876          If none is given, ``var('r')`` will be assigned.
     1877        - ``coordinate``- (Optional; default=None) An integer in
     1878          {0, ..., d-1} indicating a convenient coordinate to base
     1879          the asymptotic calculations on.
     1880          If None is assigned, then choose ``coordinate=d-1``.
     1881        - ``numerical`` - (Optional; default=0) A natural number.
     1882          If numerical > 0, then return a numerical approximation of the
     1883          Maclaurin ray coefficients of ``self`` with ``numerical`` digits
     1884          of precision. 
     1885          Otherwise return exact values.
     1886           
     1887        NOTES:
     1888
     1889        The formulas used for computing the asymptotic expansion are
     1890        Theorem 3.4 and Theorem 3.7 of [RaWi2012]_.
     1891
     1892        EXAMPLES::
     1893       
     1894            sage: from sage.combinat.amgf import *
     1895            sage: R.<x, y, z>= PolynomialRing(QQ)
     1896            sage: H = (4 - 2*x - y - z)*(4 - x -2*y - z)
     1897            sage: Hfac = H.factor()
     1898            sage: G = 16/Hfac.unit()
     1899            sage: F = FFPD(G, Hfac)
     1900            sage: print F
     1901            (16, [(x + 2*y + z - 4, 1), (2*x + y + z - 4, 1)])
     1902            sage: p = {x: 1, y: 1, z: 1}
     1903            sage: alpha = [3, 3, 2]
     1904            sage: print F.asymptotics_multiple(p, alpha, 2, var('r')) # long time
     1905            Creating auxiliary functions...
     1906            Computing derivatives of auxiliary functions...
     1907            Computing derivatives of more auxiliary functions...
     1908            Computing second-order differential operator actions...
     1909            (4/3*sqrt(3)/(sqrt(pi)*sqrt(r)) -
     1910            25/216*sqrt(3)/(sqrt(pi)*r^(3/2)), 1,
     1911            4/3*sqrt(3)/(sqrt(pi)*sqrt(r)) - 25/216*sqrt(3)/(sqrt(pi)*r^(3/2)))
     1912           
     1913        ::
     1914       
     1915            sage: R.<x, y, z>= PolynomialRing(QQ)
     1916            sage: H = (1 - x*(1 + y))*(1 - z*x**2*(1 + 2*y))
     1917            sage: Hfac = H.factor()
     1918            sage: G = 1/Hfac.unit()
     1919            sage: F = FFPD(G, Hfac)
     1920            sage: print F
     1921            (1, [(x*y + x - 1, 1), (2*x^2*y*z + x^2*z - 1, 1)])
     1922            sage: p = {x: 1/2, z: 4/3, y: 1}
     1923            sage: alpha = [8, 3, 3]
     1924            sage: print F.asymptotics_multiple(p, alpha, 2, var('r'), coordinate=1) # long time
     1925            Creating auxiliary functions...
     1926            Computing derivatives of auxiliary functions...
     1927            Computing derivatives of more auxiliary functions...
     1928            Computing second-order differential operator actions...
     1929            (1/172872*(24696*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(r)) -
     1930            1231*sqrt(3)*sqrt(7)/(sqrt(pi)*r^(3/2)))*108^r, 108,
     1931            1/7*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(r)) -
     1932            1231/172872*sqrt(3)*sqrt(7)/(sqrt(pi)*r^(3/2)))
     1933           
     1934        ::
     1935       
     1936            sage: R.<x, y>= PolynomialRing(QQ)
     1937            sage: H = (1 - 2*x - y) * (1 - x - 2*y)
     1938            sage: Hfac = H.factor()
     1939            sage: G = exp(x + y)/Hfac.unit()
     1940            sage: F = FFPD(G, Hfac)
     1941            sage: print F
     1942            (e^(x + y), [(x + 2*y - 1, 1), (2*x + y - 1, 1)])
     1943            sage: p = {x: 1/3, y: 1/3}
     1944            sage: alpha = (var('a'), var('b'))
     1945            sage: print F.asymptotics_multiple(p, alpha, 2, var('r')) # long time
     1946            (3*((1/3)^(-b)*(1/3)^(-a))^r*e^(2/3), (1/3)^(-b)*(1/3)^(-a),
     1947            3*e^(2/3))
     1948           
     1949        AUTHORS:
     1950
     1951        - Alexander Raichev (2008-10-01, 2010-09-28, 2012-08-02)       
     1952        """
     1953        from itertools import product
     1954       
     1955        R = self.ring()
     1956        if R is None:
     1957            return None
     1958       
     1959        # Coerce keys of p into R.
     1960        p = FFPD.coerce_point(R, p)
     1961           
     1962        d = self.dimension()
     1963        I = sqrt(-Integer(1))
     1964        # Coerce everything into the Symbolic Ring.
     1965        X = [SR(x) for x in R.gens()]
     1966        G = SR(self.numerator())
     1967        H = [SR(h) for (h, e) in self.denominator_factored()] 
     1968        Hprod = prod(H)
     1969        n = len(H)
     1970        P = dict([(SR(x), p[x]) for x in R.gens()])
     1971        Sstar = self.crit_cone_combo(p, alpha, coordinate)
     1972       
     1973        # Put the given convenient variable at end of variable list.
     1974        if coordinate is not None:
     1975            x = X.pop(coordinate)
     1976            X.append(x)
     1977            a = alpha.pop(coordinate)
     1978            alpha.append(a)
     1979       
     1980       
     1981        # Case n = d.
     1982        if n == d:
     1983            det = jacobian(H, X).subs(P).determinant().abs()
     1984            exp_scale = prod([(P[X[i]]**(-alpha[i])).subs(P)
     1985                              for i in xrange(d)] )
     1986            subexp_part = G.subs(P)/(det*prod(P.values()))
     1987            if numerical:
     1988                exp_scale = exp_scale.n(digits=numerical)
     1989                subexp_part = subexp_part.n(digits=numerical)
     1990            return (exp_scale**asy_var*subexp_part, exp_scale, subexp_part)
     1991       
     1992        # Case n < d.   
     1993        # If P is a tuple of rationals, then compute with it directly.
     1994        # Otherwise, compute symbolically and plug in P at the end.
     1995        if vector(P.values()) not in QQ**d:
     1996            sP = [var('p' + str(j)) for j in xrange(d)]
     1997            P = dict( [(X[j], sP[j]) for j in xrange(d)] )
     1998            p = dict( [(sP[j], p[X[j]]) for j in xrange(d)] )
     1999
     2000        # Setup.
     2001        print "Creating auxiliary functions..."       
     2002        # Create T and S variables.
     2003        t = 't'
     2004        while t in [str(x) for x in X]:
     2005            t = t + 't'
     2006        T = [var(t + str(i)) for i in xrange(d - 1)]
     2007        s = 's'
     2008        while s in [str(x) for x in X]:
     2009            s = s + 't'
     2010        S = [var(s + str(i)) for i in xrange(n - 1)]
     2011        Sstar = dict([(S[j], Sstar[j]) for j in xrange(n - 1)])
     2012        thetastar = dict([(t, Integer(0)) for t in T])
     2013        thetastar.update(Sstar)
     2014        # Create implicit functions.
     2015        h = [function('h' + str(j), *tuple(X[:d - 1])) for j in xrange(n)]
     2016        U = function('U', *tuple(X))         
     2017        # All other functions are defined in terms of h, U, and
     2018        # explicit functions.
     2019        Hcheck =  prod([X[d - 1] - Integer(1)/h[j] for j in xrange(n)])
     2020        Gcheck = -G/U * prod([-h[j]/X[d - 1] for j in xrange(n)])
     2021        A = [(-1)**(n - 1)*X[d - 1]**(-n + j)*\
     2022             diff(Gcheck.subs({X[d - 1]: Integer(1)/X[d - 1]}), X[d - 1], j)
     2023             for j in xrange(n)]
     2024        e = dict([(X[i], P[X[i]]*exp(I*T[i])) for i in xrange(d - 1)])
     2025        ht = [hh.subs(e) for hh in h]
     2026        Ht = [H[j].subs(e).subs({X[d - 1]: Integer(1)/ht[j]})
     2027              for j in xrange(n)]
     2028        hsumt = add([S[j]*ht[j] for j in xrange(n - 1)]) +\
     2029                (Integer(1) - add(S))*ht[n - 1]
     2030        At = [AA.subs(e).subs({X[d - 1]: hsumt}) for AA in A]
     2031        Phit = -log(P[X[d - 1]]*hsumt) +\
     2032               I*add([alpha[i]/alpha[d - 1]*T[i] for i in xrange(d - 1)])
     2033        # atP Stores h and U and all their derivatives evaluated at C.
     2034        atP = P.copy()
     2035        atP.update(dict([(hh.subs(P), Integer(1)/P[X[d - 1]]) for hh in h]))
     2036
     2037        # Compute the derivatives of h up to order 2*N and evaluate at P. 
     2038        hderivs1 = {}   # First derivatives of h.
     2039        for (i, j) in mrange([d - 1, n]):
     2040            s = solve(diff(H[j].subs({X[d - 1]: Integer(1)/h[j]}), X[i]),
     2041                      diff(h[j], X[i]))[0].rhs().simplify()
     2042            hderivs1.update({diff(h[j], X[i]): s})
     2043            atP.update({diff(h[j], X[i]).subs(P): s.subs(P).subs(atP)})
     2044        hderivs = FFPD.diff_all(h, X[0:d - 1], 2*N, sub=hderivs1, rekey=h)
     2045        for k in hderivs.keys():
     2046            atP.update({k.subs(P): hderivs[k].subs(atP)})
     2047
     2048        # Compute the derivatives of U up to order 2*N - 2 + min{n, N} - 1 and
     2049        # evaluate at P. 
     2050        # To do this, differentiate H = U*Hcheck over and over, evaluate at P,
     2051        # and solve for the derivatives of U at P.
     2052        # Need the derivatives of H with short keys to pass on to
     2053        # diff_prod later.
     2054        print "Computing derivatives of auxiliary functions..."
     2055        m = min(n, N)
     2056        end = [X[d-1] for j in xrange(n)]
     2057        Hprodderivs = FFPD.diff_all(Hprod, X, 2*N - 2 + n, ending=end,
     2058                                    sub_final=P)
     2059        atP.update({U.subs(P): diff(Hprod, X[d - 1], n).subs(P)/factorial(n)})
     2060        Uderivs ={}
     2061        k = Hprod.polynomial(CC).degree() - n
     2062        if k == 0:
     2063            # Then we can conclude that all higher derivatives of U are zero.
     2064            for l in xrange(1, 2*N - 2 + m):
     2065                for s in UnorderedTuples(X, l):
     2066                    Uderivs[diff(U, s).subs(P)] = Integer(0)
     2067        elif k > 0 and k < 2*N - 2 + m - 1:
     2068            all_zero = True
     2069            Uderivs = FFPD.diff_prod(Hprodderivs, U, Hcheck, X,
     2070                                     range(1, k + 1), end, Uderivs, atP)
     2071            # Check for a nonzero U derivative.
     2072            if Uderivs.values() != [Integer(0)  for i in xrange(len(Uderivs))]:
     2073                all_zero = False
     2074            if all_zero:
     2075                # Then all higher derivatives of U are zero.
     2076                for l in xrange(k + 1, 2*N - 2 + m):
     2077                    for s in UnorderedTuples(X, l):
     2078                        Uderivs.update({diff(U, s).subs(P): Integer(0)})
     2079            else:
     2080                # Have to compute the rest of the derivatives.
     2081                Uderivs = FFPD.diff_prod(Hprodderivs, U, Hcheck, X,
     2082                                         range(k + 1, 2*N - 2 + m), end,
     2083                                         Uderivs, atP)
     2084        else:
     2085            Uderivs = FFPD.diff_prod(Hprodderivs, U, Hcheck, X,
     2086                                     range(1, 2*N - 2 + m), end, Uderivs, atP)
     2087        atP.update(Uderivs)
     2088        Phit1 = jacobian(Phit, T + S).subs(hderivs1)
     2089        a = jacobian(Phit1, T + S).subs(hderivs1).subs(thetastar).subs(atP)
     2090        a_inv = a.inverse()
     2091        Phitu = Phit - (Integer(1)/Integer(2))*matrix([T + S])*a*\
     2092                matrix([T + S]).transpose()
     2093        Phitu = Phitu[0][0]
     2094
     2095        # Compute all partial derivatives of At and Phitu up to orders 2*N - 2
     2096        # and 2*N, respectively.  Take advantage of the fact that At and Phitu
     2097        # are sufficiently differentiable functions so that mixed partials
     2098        # are equal.  Thus only need to compute representative partials.
     2099        # Choose nondecreasing sequences as representative differentiation-
     2100        # order sequences.
     2101        # To speed up later computations, create symbolic functions AA and BB
     2102        # to stand in for the expressions At and Phitu respectively.
     2103        print "Computing derivatives of more auxiliary functions..."
     2104        AA = [function('A' + str(j), *tuple(T + S)) for j in xrange(n)]
     2105        At_derivs = FFPD.diff_all(At, T + S, 2*N - 2, sub=hderivs1,
     2106                                  sub_final =[thetastar, atP], rekey=AA)
     2107        BB = function('BB', *tuple(T + S))
     2108        Phitu_derivs = FFPD.diff_all(Phitu, T + S, 2*N, sub=hderivs1,
     2109                                     sub_final =[thetastar, atP], rekey=BB,
     2110                                     zero_order=3)
     2111        AABB_derivs = At_derivs
     2112        AABB_derivs.update(Phitu_derivs)
     2113        for j in xrange(n):
     2114            AABB_derivs[AA[j]] = At[j].subs(thetastar).subs(atP)
     2115        AABB_derivs[BB] = Phitu.subs(thetastar).subs(atP)
     2116       
     2117        print "Computing second-order differential operator actions..."
     2118        DD = FFPD.diff_op(AA, BB, AABB_derivs, T + S, a_inv, n, N)
     2119        L = {}
     2120        for (j, k) in product(xrange(min(n, N)), xrange(max(0, N - 1 - n), N)):
     2121            if j + k <= N - 1:
     2122                L[(j, k)] = add([DD[(j, k, l)]/((-1)**k*2**(k + l)*\
     2123                                 factorial(l)*factorial(k + l))
     2124                                 for l in xrange(2*k + 1)])
     2125        det = a.determinant()**(-Integer(1)/Integer(2))*\
     2126              (2*pi)**((n - d)/Integer(2))   
     2127        chunk = det*add([
     2128                (alpha[d - 1]*asy_var)**((n - d)/Integer(2) - q)*\
     2129                add([L[(j, k)]*binomial(n - 1, j)*\
     2130                     stirling_number1(n - j, n + k - q)*(-1)**(q - j - k)
     2131                     for (j, k) in product(xrange(min(n - 1, q) + 1),               
     2132                                           xrange(max(0, q - n), q + 1))
     2133                     if j + k <= q])
     2134                for q in xrange(N)])
     2135        chunk = chunk.subs(P).simplify()
     2136        coeffs = chunk.coefficients(asy_var)
     2137        coeffs.reverse()
     2138        coeffs = coeffs[:N]
     2139        if numerical:
     2140            subexp_part = add([co[0].subs(p).n(digits=numerical)*asy_var**co[1]
     2141                               for co in coeffs])
     2142            exp_scale = prod([(P[X[i]]**(-alpha[i])).subs(p)
     2143                              for i in xrange(d)]).n(digits=numerical)
     2144        else:
     2145            subexp_part = add([co[0].subs(p)*asy_var**co[1] for co in coeffs])
     2146            exp_scale = prod([(P[X[i]]**(-alpha[i])).subs(p)
     2147                              for i in xrange(d)] )   
     2148        return (exp_scale**asy_var*subexp_part, exp_scale, subexp_part)
     2149
     2150    @staticmethod
     2151    def subs_all(f, sub, simplify=False):
     2152        r"""
     2153        Return the items of $f$ substituted by the dictionaries
     2154        of $sub$ in order of their appearance in $sub$.
     2155
     2156        INPUT:
     2157
     2158        - ``f`` - An individual or list of symbolic expressions or dictionaries
     2159        - ``sub`` - An individual or list of dictionaries.
     2160        - ``simplify`` - Boolean (default: False).
     2161
     2162        OUTPUT:
     2163
     2164        The items of $f$ substituted by the dictionaries of $sub$ in order of
     2165        their appearance in $sub$. 
     2166        The subs() command is used.
     2167        If simplify is True, then simplify() is used after substitution.
     2168
     2169        EXAMPLES::
     2170
     2171            sage: from sage.combinat.amgf import *
     2172
     2173            sage: var('x, y, z')
     2174            (x, y, z)
     2175            sage: a = {x:1}
     2176            sage: b = {y:2}
     2177            sage: c = {z:3}
     2178            sage: FFPD.subs_all(x + y + z, a)
     2179            y + z + 1
     2180            sage: FFPD.subs_all(x + y + z, [c, a])
     2181            y + 4
     2182            sage: FFPD.subs_all([x + y + z, y^2], b)
     2183            [x + z + 2, 4]
     2184            sage: FFPD.subs_all([x + y + z, y^2], [b, c])
     2185            [x + 5, 4]
     2186
     2187        ::
     2188
     2189            sage: var('x, y')
     2190            (x, y)
     2191            sage: a = {'foo': x**2 + y**2, 'bar': x - y}
     2192            sage: b = {x: 1 , y: 2}
     2193            sage: FFPD.subs_all(a, b)
     2194            {'foo': 5, 'bar': -1}
     2195
     2196        AUTHORS:
     2197
     2198        - Alexander Raichev (2009-05-05)
     2199        """
     2200        singleton = False
     2201        if not isinstance(f, list):
     2202            f = [f]
     2203            singleton = True
     2204        if not isinstance(sub, list):
     2205            sub = [sub]
     2206        g = []
     2207        for ff in f:
     2208            for D in sub:
     2209                if isinstance(ff, dict):
     2210                    ff = dict( [(k, ff[k].subs(D)) for k in ff.keys()] )
     2211                else:
     2212                    ff = ff.subs(D)
     2213            g.append(ff)
     2214        if singleton and simplify:
     2215            if isinstance(g[Integer(0) ], dict):
     2216                return g[Integer(0) ]
     2217            else:
     2218                return g[Integer(0) ].simplify()
     2219        elif singleton and not simplify:
     2220            return g[Integer(0) ]
     2221        elif not singleton and simplify:
     2222            G = []
     2223            for gg in g:
     2224                if isinstance(gg, dict):
     2225                    G.append(gg)
     2226                else:
     2227                    G.append(gg.simplify())
     2228            return G
     2229        else:
     2230            return g
     2231
     2232    @staticmethod
     2233    def diff_all(f, V, n, ending=[], sub=None, sub_final=None,
     2234                 zero_order=0, rekey=None):
     2235        r"""
     2236        Return a dictionary of representative mixed partial
     2237        derivatives of $f$ from order 1 up to order $n$ with respect to the
     2238        variables in $V$.
     2239        The default is to key the dictionary by all nondecreasing sequences
     2240        in $V$ of length 1 up to length $n$.
     2241        For internal use.
     2242
     2243        INPUT:
     2244
     2245        - ``f`` - An individual or list of $\mathcal{C}^{n+1}$ functions.
     2246        - ``V`` - A list of variables occurring in $f$.
     2247        - ``n`` - A natural number.
     2248        - ``ending`` - A list of variables in $V$.
     2249        - ``sub`` - An individual or list of dictionaries.
     2250        - ``sub_final`` - An individual or list of dictionaries.
     2251        - ``rekey`` - A callable symbolic function in $V$ or list thereof.
     2252        - ``zero_order`` - A natural number.
     2253
     2254        OUTPUT:
     2255
     2256        The dictionary ${s_1:deriv_1,..., s_r:deriv_r}$.
     2257        Here $s_1,\ldots, s_r$ is a listing of
     2258        all nondecreasing sequences of length 1 up to length $n$ over the
     2259        alphabet $V$, where $w > v$ in $X$ iff $str(w) > str(v)$, and
     2260        $deriv_j$ is the derivative of $f$ with respect to the derivative
     2261        sequence $s_j$ and simplified with respect to the substitutions in
     2262        $sub$ and evaluated at $sub_final$.
     2263        Moreover, all derivatives with respect to sequences of length less than
     2264        $zero_order$ (derivatives of order less than $zero_order$ )
     2265        will be made zero.
     2266
     2267        If $rekey$ is nonempty, then $s_1,\ldots, s_r$ will be replaced by the
     2268        symbolic derivatives of the functions in $rekey$.
     2269
     2270        If $ending$ is nonempty, then every derivative sequence $s_j$ will be
     2271        suffixed by $ending$.
     2272
     2273        EXAMPLES::
     2274
     2275        I'd like to print the entire dictionaries, but that doesn't yield
     2276        consistent output order for doctesting.
     2277        Order of keys changes.::
     2278
     2279            sage: from sage.combinat.amgf import *
     2280            sage: f = function('f', x)
     2281            sage: dd = FFPD.diff_all(f, [x], 3)
     2282            sage: dd[(x, x, x)]
     2283            D[0, 0, 0](f)(x)
     2284
     2285        ::
     2286
     2287            sage: d1 = {diff(f, x): 4*x^3}
     2288            sage: dd = FFPD.diff_all(f,[x], 3, sub=d1)
     2289            sage: dd[(x, x, x)]
     2290            24*x
     2291
     2292        ::
     2293
     2294            sage: dd = FFPD.diff_all(f,[x], 3, sub=d1, rekey=f)
     2295            sage: dd[diff(f, x, 3)]
     2296            24*x
     2297
     2298        ::
     2299
     2300            sage: a = {x:1}
     2301            sage: dd = FFPD.diff_all(f,[x], 3, sub=d1, rekey=f, sub_final=a)
     2302            sage: dd[diff(f, x, 3)]
     2303            24
     2304
     2305        ::
     2306
     2307            sage: X = var('x, y, z')
     2308            sage: f = function('f',*X)
     2309            sage: dd = FFPD.diff_all(f, X, 2, ending=[y, y, y])
     2310            sage: dd[(z, y, y, y)]
     2311            D[1, 1, 1, 2](f)(x, y, z)
     2312
     2313        ::
     2314
     2315            sage: g = function('g',*X)
     2316            sage: dd = FFPD.diff_all([f, g], X, 2)
     2317            sage: dd[(0, y, z)]
     2318            D[1, 2](f)(x, y, z)
     2319
     2320        ::
     2321
     2322            sage: dd[(1, z, z)]
     2323            D[2, 2](g)(x, y, z)
     2324
     2325        ::
     2326
     2327            sage: f = exp(x*y*z)
     2328            sage: ff = function('ff',*X)
     2329            sage: dd = FFPD.diff_all(f, X, 2, rekey=ff)
     2330            sage: dd[diff(ff, x, z)]
     2331            x*y^2*z*e^(x*y*z) + y*e^(x*y*z)
     2332
     2333        AUTHORS:
     2334
     2335        - Alexander Raichev (2008-10-01, 2009-04-01, 2010-02-01)
     2336        """
     2337        singleton=False
     2338        if not isinstance(f, list):
     2339            f = [f]
     2340            singleton=True
     2341
     2342        # Build the dictionary of derivatives iteratively from a list
     2343        # of nondecreasing derivative-order sequences.
     2344        derivs = {}
     2345        r = len(f)
     2346        if ending:
     2347            seeds = [ending]
     2348            start = Integer(1)
     2349        else:
     2350            seeds = [[v] for v in V]
     2351            start = Integer(2)
     2352        if singleton:
     2353            for s in seeds:
     2354                derivs[tuple(s)] = FFPD.subs_all(diff(f[0], s), sub)
     2355            for l in xrange(start, n + 1):
     2356                for t in UnorderedTuples(V, l):
     2357                    s = tuple(t + ending)
     2358                    derivs[s] = FFPD.subs_all(diff(derivs[s[1:]], s[0]), sub)
     2359        else:
     2360            # Make the dictionary keys of the form (j, sequence of variables),
     2361            # where j in range(r).
     2362            for s in seeds:
     2363                value = FFPD.subs_all([diff(f[j], s) for j in xrange(r)], sub)
     2364                derivs.update(dict([(tuple([j]+s), value[j])
     2365                                    for j in xrange(r)]))
     2366            for l in xrange(start, n + 1):
     2367                for t in UnorderedTuples(V, l):
     2368                    s = tuple(t + ending)
     2369                    value = FFPD.subs_all([diff(derivs[(j,) + s[1:]],
     2370                                          s[0]) for j in xrange(r)], sub)
     2371                    derivs.update(dict([((j,) + s, value[j])
     2372                                        for j in xrange(r)]))
     2373        if zero_order:
     2374            # Zero out all the derivatives of order < zero_order
     2375            if singleton:
     2376                for k in derivs.keys():
     2377                    if len(k) < zero_order:
     2378                        derivs[k] = Integer(0)
     2379            else:
     2380                # Ignore the first of element of k, which is an index.
     2381                for k in derivs.keys():
     2382                    if len(k) - 1  < zero_order:
     2383                        derivs[k] = Integer(0)
     2384        if sub_final:
     2385            # Substitute sub_final into the values of derivs.
     2386            for k in derivs.keys():
     2387                derivs[k] = FFPD.subs_all(derivs[k], sub_final) 
     2388        if rekey:
     2389            # Rekey the derivs dictionary by the value of rekey.
     2390            F = rekey
     2391            if singleton:
     2392                # F must be a singleton.
     2393                derivs = dict( [(diff(F, list(k)), derivs[k])
     2394                                for k in derivs.keys()] )
     2395            else:
     2396                # F must be a list.
     2397                derivs = dict( [(diff(F[k[0]], list(k)[1:]), derivs[k])
     2398                                for k in derivs.keys()] )
     2399        return derivs
     2400
     2401    @staticmethod
     2402    def diff_op(A, B, AB_derivs, V, M, r, N):
     2403        r"""
     2404        Return the derivatives $DD^(l+k)(A[j] B^l)$ evaluated at a point
     2405        $p$ for various natural numbers $j, k, l$ which depend on $r$ and $N$. 
     2406        Here $DD$ is a specific second-order linear differential operator
     2407        that depends on $M$ , $A$ is a list of symbolic functions,
     2408        $B$ is symbolic function, and $AB_derivs$ contains all the derivatives
     2409        of $A$ and $B$ evaluated at $p$ that are necessary for the computation.
     2410        For internal use by the functions asymptotics_smooth() and
     2411        asymptotics_multiple().
     2412
     2413        INPUT:
     2414
     2415        - ``A`` - A single or length ``r`` list of symbolic functions in the
     2416          variables ``V``.
     2417        - ``B`` - A symbolic function in the variables ``V``.
     2418        - ``AB_derivs`` - A dictionary whose keys are the (symbolic)
     2419          derivatives of ``A[0], ..., A[r-1]`` up to order ``2*N-2`` and
     2420          the (symbolic) derivatives of ``B`` up to order ``2*N``. 
     2421          The values of the dictionary are complex numbers that are
     2422          the keys evaluated at a common point $p$.
     2423        - ``V`` - The variables of the ``A[j]`` and ``B``.
     2424        - ``M`` - A symmetric $l \times l$ matrix, where $l$ is the
     2425          length of ``V``.
     2426        - ``r, N`` - Natural numbers.
     2427
     2428        OUTPUT:
     2429
     2430        A dictionary whose keys are natural number tuples of the form
     2431        $(j, k, l)$, where $l \le 2k$, $j \le r-1$, and $j+k \le N-1$,
     2432        and whose values are $DD^(l+k)(A[j] B^l)$ evaluated at a point
     2433        $p$, where $DD$ is the linear second-order differential operator 
     2434        $-\sum_{i=0}^{l-1} \sum_{j=0}^{l-1} M[i][j]
     2435        \partial^2 /(\partial V[j] \partial V[i])$.
     2436
     2437        EXAMPLES::
     2438
     2439            sage: from sage.combinat.amgf import *
     2440
     2441            sage: T = var('x, y')
     2442            sage: A = function('A',*tuple(T))
     2443            sage: B = function('B',*tuple(T))
     2444            sage: AB_derivs = {}
     2445            sage: M = matrix([[1, 2],[2, 1]])
     2446            sage: DD = FFPD.diff_op(A, B, AB_derivs, T, M, 1, 2)
     2447            sage: DD.keys()
     2448            [(0, 1, 2), (0, 1, 1), (0, 1, 0), (0, 0, 0)]
     2449            sage: len(DD[(0, 1, 2)])
     2450            246
     2451
     2452        AUTHORS:
     2453
     2454        - Alexander Raichev (2008-10-01, 2010-01-12)
     2455        """
     2456        if not isinstance(A, list):
     2457            A = [A]
     2458
     2459        # First, compute the necessary product derivatives of A and B.
     2460        product_derivs = {}
     2461        for (j, k) in mrange([r, N]):
     2462            if j + k < N:
     2463                for l in xrange(2*k + 1):
     2464                    for s in UnorderedTuples(V, 2*(k + l)):
     2465                        product_derivs[tuple([j, k, l] + s)] = \
     2466                        diff(A[j]*B**l, s).subs(AB_derivs)
     2467
     2468        # Second, compute DD^(k+l)(A[j]*B^l)(p) and store values in dictionary.
     2469        DD = {}
     2470        rows = M.nrows()
     2471        for (j, k) in mrange([r, N]):
     2472            if j + k < N:
     2473                for l in xrange(2*k + 1):
     2474                    # Take advantage of the symmetry of M by ignoring
     2475                    # the upper-diagonal entries of M and multiplying by
     2476                    # appropriate powers of 2.
     2477                    if k + l == 0 :
     2478                        DD[(j, k, l)] = product_derivs[(j, k, l)]
     2479                        continue
     2480                    S = [(a, b) for (a, b) in mrange([rows, rows]) if b <= a]
     2481                    P =  cartesian_product_iterator([S for i in range(k+l)])
     2482                    diffo = Integer(0)
     2483                    for t in P:
     2484                        if product_derivs[(j, k, l) + FFPD.diff_seq(V, t)] !=\
     2485                          Integer(0):
     2486                            MM = Integer(1)
     2487                            for (a, b) in t:
     2488                                MM = MM * M[a][b]
     2489                                if a != b:
     2490                                    MM = Integer(2) *MM
     2491                            diffo = diffo + MM*product_derivs[(j, k, l) +\
     2492                                           FFPD.diff_seq(V, t)]
     2493                    DD[(j, k, l)] = (-Integer(1) )**(k+l)*diffo
     2494        return DD
     2495       
     2496    @staticmethod
     2497    def diff_seq(V, s):
     2498        r"""
     2499        Given a list ``s`` of tuples of natural numbers, return the
     2500        list of elements of ``V`` with indices the elements of the elements
     2501        of ``s``.
     2502        This function is for internal use by the function diff_op().
     2503
     2504        INPUT:
     2505
     2506        - ``V`` - A list.
     2507        - ``s`` - A list of tuples of natural numbers in the interval
     2508          ``range(len(V))``.
     2509
     2510        OUTPUT:
     2511
     2512        The tuple ``tuple([V[tt] for tt in sorted(t)])``, where ``t`` is the
     2513        list of elements of the elements of ``s``.
     2514
     2515        EXAMPLES::
     2516
     2517            sage: from sage.combinat.amgf import *
     2518
     2519            sage: V = list(var('x, t, z'))
     2520            sage: FFPD.diff_seq(V,([0, 1],[0, 2, 1],[0, 0]))
     2521            (x, x, x, x, t, t, z)
     2522
     2523        AUTHORS:
     2524
     2525        - Alexander Raichev (2009-05-19)
     2526        """
     2527        t = []
     2528        for ss in s:
     2529            t.extend(ss)
     2530        return tuple([V[tt] for tt in sorted(t)])
     2531
     2532    @staticmethod
     2533    def diff_op_simple(A, B, AB_derivs, x, v, a, N):
     2534        r"""
     2535        Return $DD^(e k + v l)(A B^l)$ evaluated at a point $p$ for
     2536        various natural numbers $e, k, l$ that depend on $v$ and $N$. 
     2537        Here $DD$ is a specific linear differential operator that depends
     2538        on $a$ and $v$ , $A$ and $B$ are symbolic functions, and $AB_derivs$
     2539        contains all the derivatives of $A$ and $B$ evaluated at $p$ that are
     2540        necessary for the computation.
     2541        For internal use by the function asymptotics_smooth().
     2542
     2543        INPUT:
     2544
     2545        - ``A, B`` - Symbolic functions in the variable ``x``.
     2546        - ``AB_derivs`` - A dictionary whose keys are the (symbolic)
     2547          derivatives of ``A`` up to order ``2*N`` if ``v`` is even or
     2548          ``N`` if ``v`` is odd and the (symbolic) derivatives of ``B``
     2549          up to order ``2*N + v`` if ``v`` is even or ``N + v``
     2550          if ``v`` is odd. 
     2551          The values of the dictionary are complex numbers that are
     2552          the keys evaluated at a common point $p$.
     2553        - ``x`` - Symbolic variable.
     2554        - ``a`` - A complex number.
     2555        - ``v, N`` - Natural numbers.
     2556
     2557        OUTPUT:
     2558
     2559        A dictionary whose keys are natural number pairs of the form $(k, l)$,
     2560        where $k < N$ and $l \le 2k$ and whose values are
     2561        $DD^(e k + v l)(A B^l)$ evaluated at a point $p$. 
     2562        Here $e=2$ if $v$ is even, $e=1$ if $v$ is odd, and $DD$ is the
     2563        linear differential operator
     2564        $(a^{-1/v} d/dt)$ if $v$ is even and
     2565        $(|a|^{-1/v} i \text{sgn}(a) d/dt)$ if $v$ is odd.
     2566
     2567        EXAMPLES::
     2568
     2569            sage: from sage.combinat.amgf import *
     2570
     2571            sage: A = function('A', x)
     2572            sage: B = function('B', x)
     2573            sage: AB_derivs = {}
     2574            sage: FFPD.diff_op_simple(A, B, AB_derivs, x, 3, 2, 2)
     2575            {(1, 0): 1/2*I*2^(2/3)*D[0](A)(x), (0, 0): A(x), (1, 1):
     2576            1/4*(A(x)*D[0, 0, 0, 0](B)(x) + B(x)*D[0, 0, 0, 0](A)(x) +
     2577            4*D[0](A)(x)*D[0, 0, 0](B)(x) + 4*D[0](B)(x)*D[0, 0, 0](A)(x) +
     2578            6*D[0, 0](A)(x)*D[0, 0](B)(x))*2^(2/3)}
     2579
     2580        AUTHORS:
     2581
     2582        - Alexander Raichev (2010-01-15) 
     2583        """
     2584        I = sqrt(-Integer(1))
     2585        DD = {}
     2586        if v.mod(Integer(2)) == Integer(0) :
     2587            for k in xrange(N):
     2588                for l in xrange(2*k + 1):
     2589                    DD[(k, l)] = (a**(-Integer(1)/v))**(2*k + v*l)*\
     2590                                 diff(A*B**l, x, 2*k + v*l).subs(AB_derivs)
     2591        else:
     2592            for k in xrange(N):
     2593                for l in xrange(k + 1):
     2594                    DD[(k, l)] = (abs(a)**(-Integer(1)/v)*I*\
     2595                                 a/abs(a))**(k+v*l)*\
     2596                                 diff(A*B**l, x, k + v*l).subs(AB_derivs)
     2597        return DD
     2598
     2599    @staticmethod
     2600    def diff_prod(f_derivs, u, g, X, interval, end, uderivs, atc):
     2601        r"""
     2602        Take various derivatives of the equation $f = ug$,
     2603        evaluate them at a point $c$, and solve for the derivatives of $u$.
     2604        For internal use by the function asymptotics_multiple().
     2605
     2606        INPUT:
     2607
     2608        - ``f_derivs`` - A dictionary whose keys are all tuples of the form
     2609          ``s + end``, where ``s`` is a sequence of variables from ``X`` whose
     2610          length lies in ``interval``, and whose values are the derivatives
     2611          of a function $f$ evaluated at $c$.
     2612        - ``u`` - A callable symbolic function.
     2613        - ``g`` - An expression or callable symbolic function.
     2614        - ``X`` - A list of symbolic variables.
     2615        - ``interval`` - A list of positive integers. 
     2616          Call the first and last values $n$ and $nn$, respectively.
     2617        - ``end`` - A possibly empty list of repetitions of the
     2618          variable ``z``, where ``z`` is the last element of ``X``.
     2619        - ``uderivs`` - A dictionary whose keys are the symbolic
     2620          derivatives of order 0 to order $n-1$ of ``u`` evaluated at $c$
     2621          and whose values are the corresponding derivatives evaluated
     2622          at $c$.
     2623        - ``atc`` - A dictionary whose keys are the keys of $c$ and all
     2624          the symbolic derivatives of order 0 to order $nn$ of ``g``
     2625          evaluated $c$ and whose values are the corresponding
     2626          derivatives evaluated at $c$.
     2627
     2628        OUTPUT:
     2629
     2630        A dictionary whose keys are the derivatives of ``u`` up to order
     2631        $nn$ and whose values are those derivatives evaluated at $c$.
     2632
     2633        EXAMPLES::
     2634
     2635        I'd like to print out the entire dictionary, but that does not give
     2636        consistent output for doctesting. 
     2637        Order of keys changes ::
     2638
     2639            sage: from sage.combinat.amgf import *
     2640            sage: u = function('u', x)
     2641            sage: g = function('g', x)
     2642            sage: fd = {(x,):1,(x, x):1}
     2643            sage: ud = {u(x=2): 1}
     2644            sage: atc = {x: 2, g(x=2): 3, diff(g, x)(x=2): 5}
     2645            sage: atc[diff(g, x, x)(x=2)] = 7
     2646            sage: dd = FFPD.diff_prod(fd, u, g, [x], [1, 2], [], ud, atc)
     2647            sage: dd[diff(u, x, 2)(x=2)]
     2648            22/9
     2649
     2650        NOTES:
     2651
     2652        This function works by differentiating the equation $f = ug$
     2653        with respect to the variable sequence ``s + end``,
     2654        for all tuples ``s`` of ``X`` of lengths in ``interval``,
     2655        evaluating at the point $c$ ,
     2656        and solving for the remaining derivatives of ``u``. 
     2657        This function assumes that ``u`` never appears in the
     2658        differentiations of $f = ug$ after evaluating at $c$.
     2659
     2660        AUTHORS:
     2661
     2662        - Alexander Raichev (2009-05-14, 2010-01-21)
     2663        """
     2664        for l in interval:
     2665            D = {}
     2666            rhs = []
     2667            lhs = []
     2668            for t in UnorderedTuples(X, l):
     2669                s = t + end
     2670                lhs.append(f_derivs[tuple(s)])
     2671                rhs.append(diff(u*g, s).subs(atc).subs(uderivs))
     2672                # Since Sage's solve command can't take derivatives as variable
     2673                # names, make new variables based on t to stand in for
     2674                # diff(u, t) and store them in D.
     2675                D[diff(u, t).subs(atc)] = var('zing' +\
     2676                                              ''.join([str(x) for x in t]))
     2677            eqns = [lhs[i] == rhs[i].subs(uderivs).subs(D)
     2678                    for i in xrange(len(lhs))]
     2679            variables = D.values()
     2680            sol = solve(eqns,*variables, solution_dict=True)
     2681            uderivs.update(FFPD.subs_all(D, sol[Integer(0) ]))
     2682        return uderivs   
     2683
     2684    def crit_cone_combo(self, p, alpha, coordinate=None):
     2685        r"""
     2686        Return an auxiliary point associated to the multiple
     2687        point ``p`` of the factors ``self``. 
     2688        For internal use by asymptotics_multiple().
     2689
     2690        INPUT:
     2691
     2692        - ``p`` - A dictionary with keys that can be coerced to equal
     2693          ``self.ring().gens()``.
     2694        - ``alpha`` - A list of rationals.
     2695
     2696        OUTPUT:
     2697
     2698        A solution of the matrix equation $y \Gamma = \alpha'$ for $y$ ,
     2699        where $\Gamma$ is the matrix given by
     2700        ``[FFPD.direction(v) for v in self.log_grads(p)]`` and $\alpha'$
     2701        is ``FFPD.direction(alpha)``
     2702
     2703        EXAMPLES::
     2704
     2705            sage: from sage.combinat.amgf import *
     2706
     2707            sage: R.<x, y>= PolynomialRing(QQ)
     2708            sage: p = exp(x)
     2709            sage: df = [(1 - 2*x - y, 1), (1 - x - 2*y, 1)]
     2710            sage: f = FFPD(p, df)
     2711            sage: p = {x: 1/3, y: 1/3}
     2712            sage: alpha = (var('a'), var('b'))
     2713            sage: print f.crit_cone_combo(p, alpha)
     2714            [1/3*(2*a - b)/b, -2/3*(a - 2*b)/b]
     2715           
     2716        NOTES:
     2717
     2718        Use this function only when $\Gamma$ is well-defined and
     2719        there is a unique solution to the matrix equation
     2720        $y \Gamma = \alpha'$. 
     2721        Fails otherwise.
     2722
     2723        AUTHORS:
     2724
     2725        - Alexander Raichev (2008-10-01, 2008-11-25, 2009-03-04, 2010-09-08,
     2726                        2010-12-02, 2012-08-02)
     2727        """
     2728        # Assuming here that each log_grads(f) has nonzero final component.
     2729        # Then 'direction' will not throw a division by zero error.
     2730        R = self.ring()
     2731        if R is None:
     2732            return None
     2733       
     2734        # Coerce keys of p into R.
     2735        p = FFPD.coerce_point(R, p)
     2736           
     2737        d = self.dimension()
     2738        n = len(self.denominator_factored())
     2739        Gamma = matrix([FFPD.direction(v, coordinate)
     2740                        for v in self.log_grads(p)])
     2741        beta = FFPD.direction(alpha, coordinate)
     2742        # solve_left() fails when working in SR :-(.
     2743        # So use solve() instead.
     2744        # Gamma.solve_left(vector(beta))
     2745        V = [var('sss'+str(i)) for i in range(n)]
     2746        M = matrix(V)*Gamma
     2747        eqns = [M[0][j] == beta[j] for j in range(d)]
     2748        s = solve(eqns, V, solution_dict=True)[0]  # Assume a unique solution.
     2749        return [s[v] for v in V]
     2750           
     2751    @staticmethod
     2752    def direction(v, coordinate=None):
     2753        r"""
     2754        Returns ``[vv/v[coordinate] for vv in v]`` where
     2755        ``coordinate`` is the last index of v if not specified otherwise.
     2756
     2757        INPUT:
     2758
     2759        - ``v`` - A vector.
     2760        - ``coordinate`` - (Optional; default=None) An index for ``v``.
     2761
     2762        EXAMPLES::
     2763
     2764            sage: from sage.combinat.amgf import *
     2765
     2766            sage: FFPD.direction([2, 3, 5])
     2767            (2/5, 3/5, 1)
     2768            sage: FFPD.direction([2, 3, 5], 0)
     2769            (1, 3/2, 5/2)
     2770
     2771        AUTHORS:
     2772
     2773        - Alexander Raichev (2010-08-25)
     2774        """
     2775        if coordinate is None:
     2776            coordinate = len(v) - 1     
     2777        return tuple([vv/v[coordinate] for vv in v])
     2778
     2779    def grads(self, p):
     2780        r"""
     2781        Return a list of the gradients of the polynomials
     2782        ``[q for (q, e) in self.denominator_factored()]`` evalutated at ``p``.
     2783
     2784        INPUT:
     2785
     2786        - ``p`` - (Optional: default=None) A dictionary whose keys are the
     2787          generators of ``self.ring()``.
     2788
     2789
     2790        EXAMPLES::
     2791
     2792            sage: from sage.combinat.amgf import *
     2793
     2794            sage: R.<x, y>= PolynomialRing(QQ)
     2795            sage: p = exp(x)
     2796            sage: df = [(x**3 + 3*y^2, 5), (x*y, 2), (y, 1)]
     2797            sage: f = FFPD(p, df)
     2798            sage: print f
     2799            (e^x, [(y, 1), (x*y, 2), (x^3 + 3*y^2, 5)])
     2800            sage: print R.gens()
     2801            (x, y)
     2802            sage: p = None
     2803            sage: print f.grads(p)
     2804            [(0, 1), (y, x), (3*x^2, 6*y)]
     2805       
     2806        ::
     2807       
     2808            sage: p = {x: sqrt(2), y: var('a')}
     2809            sage: print f.grads(p)
     2810            [(0, 1), (a, sqrt(2)), (6, 6*a)]
     2811           
     2812        AUTHORS:
     2813
     2814        - Alexander Raichev (2009-03-06)
     2815        """
     2816        R = self.ring()
     2817        if R is None:
     2818            return
     2819           
     2820        # Coerce keys of p into R.
     2821        p = FFPD.coerce_point(R, p)
     2822       
     2823        X = R.gens()
     2824        d = self.dimension()
     2825        H = [h for (h, e) in self.denominator_factored()]
     2826        n = len(H)
     2827        return [tuple([diff(H[i], X[j]).subs(p) for j in xrange(d)])
     2828                for i in xrange(n)]
     2829
     2830    def log_grads(self, p):
     2831        r"""
     2832        Return a list of the logarithmic gradients of the polynomials
     2833        ``[q for (q, e) in self.denominator_factored()]`` evalutated at ``p``.
     2834
     2835        INPUT:
     2836
     2837        - ``p`` - (Optional: default=None) A dictionary whose keys are the
     2838          generators of ``self.ring()``.
     2839           
     2840        NOTE:
     2841       
     2842        The logarithmic gradient of a function $f$ at point $p$ is the vector
     2843        $(x_1 \partial_1 f(x), \ldots, x_d \partial_d f(x) )$ evaluated at
     2844        $p$.
     2845
     2846
     2847        EXAMPLES::
     2848
     2849            sage: from sage.combinat.amgf import *
     2850
     2851            sage: R.<x, y>= PolynomialRing(QQ)
     2852            sage: p = exp(x)
     2853            sage: df = [(x**3 + 3*y^2, 5), (x*y, 2), (y, 1)]
     2854            sage: f = FFPD(p, df)
     2855            sage: print f
     2856            (e^x, [(y, 1), (x*y, 2), (x^3 + 3*y^2, 5)])
     2857            sage: print R.gens()
     2858            (x, y)
     2859            sage: p = None
     2860            sage: print f.log_grads(p)
     2861            [(0, y), (x*y, x*y), (3*x^3, 6*y^2)]
     2862           
     2863        ::
     2864       
     2865            sage: p = {x: sqrt(2), y: var('a')}
     2866            sage: print f.log_grads(p)
     2867            [(0, a), (sqrt(2)*a, sqrt(2)*a), (6*sqrt(2), 6*a^2)]
     2868
     2869        AUTHORS:
     2870
     2871        - Alexander Raichev (2009-03-06)
     2872        """
     2873        R = self.ring()
     2874        if R is None:
     2875            return None
     2876           
     2877        # Coerce keys of p into R.
     2878        p = FFPD.coerce_point(R, p)
     2879       
     2880        X = R.gens()
     2881        d = self.dimension()
     2882        H = [h for (h, e) in self.denominator_factored()]
     2883        n = len(H)
     2884        return [tuple([(X[j]*diff(H[i], X[j])).subs(p) for j in xrange(d)])
     2885                for i in xrange(n)]
     2886
     2887    def critical_cone(self, p, coordinate=None):
     2888        r"""
     2889        Return the critical cone of the convenient multiple point ``p``.
     2890       
     2891        INPUT:
     2892
     2893        - ``p`` - A dictionary with keys that can be coerced to equal
     2894          ``self.ring().gens()`` and values in a field.
     2895        - ``coordinate`` - (Optional; default=None) A natural number.
     2896
     2897        OUTPUT:
     2898
     2899        A list of vectors that generate the critical cone of ``p`` and
     2900        the cone itself, which is None if the values of ``p`` don't lie in QQ.
     2901        Divide logarithmic gradients by their component ``coordinate`` entries.
     2902        If ``coordinate=None``, then search from d-1 down to 0 for the
     2903        first index j such that for all i ``self.log_grads()[i][j] != 0``
     2904        and set ``coordinate=j``.
     2905
     2906        EXAMPLES::
     2907
     2908            sage: from sage.combinat.amgf import *
     2909
     2910            sage: R.<x, y, z>= PolynomialRing(QQ)
     2911            sage: G = 1
     2912            sage: H = (1 - x*(1 + y))*(1 - z*x**2*(1 + 2*y))
     2913            sage: Hfac = H.factor()
     2914            sage: G = 1/Hfac.unit()
     2915            sage: F = FFPD(G, Hfac)
     2916            sage: p = {x: 1/2, y: 1, z: 4/3}
     2917            sage: print F.critical_cone(p)
     2918            ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N)
     2919           
     2920        AUTHORS:
     2921
     2922        - Alexander Raichev (2010-08-25, 2012-08-02)
     2923        """
     2924        R = self.ring()
     2925        if R is None:
     2926            return
     2927         
     2928        # Coerce keys of p into R.
     2929        p = FFPD.coerce_point(R, p)
     2930 
     2931        X = R.gens()
     2932        d = self.dimension()
     2933        lg = self.log_grads(p)
     2934        n = len(lg)
     2935        if coordinate not in xrange(d):
     2936            # Search from d-1 down to 0 for a coordinate j such that
     2937            # for all i we have lg[i][j] != 0.
     2938            # One is guaranteed to exist in the case of a convenient multiple
     2939            # point.
     2940            for j in reversed(xrange(d)):
     2941                if 0 not in [lg[i][j] for i in xrange(n)]:
     2942                    coordinate = j
     2943                    break
     2944        Gamma = [FFPD.direction(v, coordinate) for v in lg]
     2945        try:
     2946            cone = Cone(Gamma)
     2947        except TypeError:
     2948            cone = None
     2949        return (Gamma, cone)
     2950   
     2951    def is_convenient_multiple_point(self, p):
     2952        r"""
     2953        Return True if ``p`` is a convenient multiple point of ``self`` and
     2954        False otherwise.
     2955        Also return a short comment.
     2956
     2957        INPUT:
     2958
     2959        - ``p`` - A dictionary with keys that can be coerced to equal
     2960          ``self.ring().gens()``.
     2961
     2962        OUTPUT:
     2963
     2964        A pair (verdict, comment).
     2965        In case ``p`` is a convenient multiple point, verdict=True and
     2966        comment ='No problem'.
     2967        In case ``p`` is not, verdict=False and comment is string explaining
     2968        why ``p`` fails to be a convenient multiple point.
     2969       
     2970        EXAMPLES::
     2971
     2972            sage: from sage.combinat.amgf import *
     2973
     2974            sage: R.<x, y, z>= PolynomialRing(QQ)
     2975            sage: H = (1 - x*(1 + y))*(1 - z*x**2*(1 + 2*y))
     2976            sage: df = H.factor()
     2977            sage: G = 1/df.unit()
     2978            sage: F = FFPD(G, df)
     2979            sage: p1 = {x: 1/2, y: 1, z: 4/3}
     2980            sage: p2 = {x: 1, y: 2, z: 1/2}
     2981            sage: print F.is_convenient_multiple_point(p1)
     2982            (True, 'convenient in variables [x, y]')
     2983            sage: print F.is_convenient_multiple_point(p2)
     2984            (False, 'not a singular point')
     2985
     2986        NOTES:
     2987
     2988        See [RaWi2012]_ for more details.
     2989
     2990        AUTHORS:
     2991
     2992        - Alexander Raichev (2011-04-18, 2012-08-02)
     2993        """
     2994        R = self.ring()
     2995        if R is None:
     2996            return
     2997           
     2998        # Coerce keys of p into R.
     2999        p = FFPD.coerce_point(R, p)
     3000       
     3001        H = [h for (h, e) in self.denominator_factored()]
     3002        n = len(H)
     3003        d = self.dimension()
     3004
     3005        # Test 1: Are the factors in H zero at p?
     3006        if [h.subs(p) for h in H] != [0 for h in H]:
     3007            # Failed test 1.  Move on to next point.
     3008            return (False, 'not a singular point')
     3009
     3010        # Test 2: Are the factors in H smooth at p?
     3011        grads = self.grads(p)
     3012        for v in grads:
     3013            if v == [0 for i in xrange(d)]:
     3014                return (False, 'not smooth point of factors')
     3015         
     3016        # Test 3: Do the factors in H intersect transversely at p?
     3017        if n <= d:
     3018            M = matrix(grads)
     3019            if M.rank() != n:
     3020                return (False, 'not a transverse intersection')
     3021        else:
     3022            # Check all sub-multisets of grads of size d.
     3023            for S in Subsets(grads, d, submultiset=True):
     3024                M = matrix(S)
     3025                if M.rank() != d:
     3026                    return (False, 'not a transverse intersection')
     3027           
     3028        # Test 4: Is p convenient?
     3029        M = matrix(self.log_grads(p))
     3030        convenient_coordinates = []
     3031        for j in xrange(d):
     3032            if 0 not in M.columns()[j]:
     3033                convenient_coordinates.append(j)
     3034        if not convenient_coordinates:
     3035            return (False, 'multiple point but not convenient')
     3036
     3037        # Tests all passed
     3038        X = R.gens()
     3039        return (True, 'convenient in variables %s' %\
     3040                [X[i] for i in convenient_coordinates])
     3041
     3042    def singular_ideal(self):
     3043        r"""
     3044        Let $R$ be the ring of ``self`` and $H$ its denominator.
     3045        Let $Hred$ be the reduction (square-free part) of $H$.
     3046        Return the ideal in $R$ generated by $Hred$ and
     3047        its partial derivatives.
     3048        If the coefficient field of $R$ is algebraically closed,
     3049        then the output is the ideal of the singular locus (which is a variety)
     3050        of the variety of $H$.
     3051       
     3052        EXAMPLES::
     3053
     3054            sage: from sage.combinat.amgf import *
     3055
     3056            sage: R.<x, y, z>= PolynomialRing(QQ)
     3057            sage: H = (1 - x*(1 + y))**3*(1 - z*x**2*(1 + 2*y))
     3058            sage: df = H.factor()
     3059            sage: G = 1/df.unit()
     3060            sage: F = FFPD(G, df)
     3061            sage: F.singular_ideal()
     3062            Ideal (x*y + x - 1, y^2 - 2*y*z + 2*y - z + 1, x*z + y - 2*z + 1)
     3063            of Multivariate Polynomial Ring in x, y, z over Rational Field
     3064
     3065        AUTHORS:
     3066
     3067        - Alexander Raichev (2008-10-01, 2008-11-20, 2010-12-03, 2011-04-18,
     3068                        2012-08-03)
     3069        """
     3070        R = self.ring()
     3071        if R is None:
     3072            return
     3073           
     3074        Hred = prod([h for (h, e) in self.denominator_factored()])
     3075        J = R.ideal([Hred] + Hred.gradient())
     3076        return R.ideal(J.groebner_basis())
     3077
     3078    def smooth_critical_ideal(self, alpha):
     3079        r"""
     3080        Let $R$ be the ring of ``self`` and $H$ its denominator.
     3081        Return the ideal in $R$ of smooth critical points of the variety
     3082        of $H$ for the direction ``alpha``.
     3083        If the variety $V$ of $H$ has no smooth points, then return the ideal
     3084        in $R$ of $V$.
     3085
     3086        INPUT:
     3087
     3088        - ``alpha`` - A d-tuple of positive integers and/or symbolic entries,
     3089          where d = ``self.ring().ngens()``.
     3090
     3091        EXAMPLES::
     3092
     3093            sage: from sage.combinat.amgf import *
     3094
     3095            sage: R.<x, y> = PolynomialRing(QQ)
     3096            sage: H = (1-x-y-x*y)^2
     3097            sage: Hfac = H.factor()
     3098            sage: G = 1/Hfac.unit()
     3099            sage: F = FFPD(G, Hfac)
     3100            sage: alpha = var('a1, a2')
     3101            sage: F.smooth_critical_ideal(alpha)
     3102            Ideal (y^2 + 2*a1/a2*y - 1, x + (a2/(-a1))*y + (-a2 + a1)/(-a1))
     3103            of Multivariate Polynomial Ring in x, y over Fraction Field of
     3104            Multivariate Polynomial Ring in a2, a1 over Rational Field
     3105           
     3106        ::
     3107       
     3108            sage: R.<x, y> = PolynomialRing(QQ)
     3109            sage: H = (1-x-y-x*y)^2
     3110            sage: Hfac = H.factor()
     3111            sage: G = 1/Hfac.unit()
     3112            sage: F = FFPD(G, Hfac)
     3113            sage: alpha = [7/3, var('a')]
     3114            sage: F.smooth_critical_ideal(alpha)
     3115            Ideal (y^2 + (-14/(-3*a))*y - 1, x + (-3/7*a)*y + 3/7*a - 1) of
     3116            Multivariate Polynomial Ring in x, y over Fraction Field of
     3117            Univariate Polynomial Ring in a over Rational Field
     3118           
     3119        NOTES:
     3120
     3121        See [RaWi2012]_ for more details.
     3122
     3123        AUTHORS:
     3124
     3125        - Alexander Raichev (2008-10-01, 2008-11-20, 2009-03-09, 2010-12-02,
     3126                        2011-04-18, 2012-08-03)
     3127        """
     3128        R = self.ring()
     3129        if R is None:
     3130            return
     3131           
     3132        Hred = prod([h for (h, e) in self.denominator_factored()])
     3133        K = R.base_ring()
     3134        d = self.dimension()
     3135       
     3136        # Expand K by the variables of alpha if there are any.
     3137        indets = []
     3138        for a in alpha:
     3139            if a not in K and a in SR:
     3140                indets.append(a)
     3141        indets = list(Set(indets))   # Delete duplicates in indets.
     3142        if indets:
     3143            L = FractionField(PolynomialRing(K, indets))
     3144            S = R.change_ring(L)
     3145            # Coerce alpha into L.
     3146            alpha = [L(a) for a in alpha]
     3147        else:
     3148            S = R
     3149       
     3150        # Find smooth, critical points for alpha.
     3151        X = S.gens()
     3152        Hred = S(Hred)
     3153        J = S.ideal([Hred] +\
     3154            [alpha[d - 1]*X[i]*diff(Hred, X[i]) -\
     3155             alpha[i]*X[d - 1]*diff(Hred, X[d - 1])
     3156             for i in xrange(d - 1)])
     3157        return S.ideal(J.groebner_basis())
     3158
     3159    def maclaurin_coefficients(self, multi_indices, numerical=0):
     3160        r"""
     3161        Returns the Maclaurin coefficients of self that have multi-indices
     3162        ``alpha``, ``2*alpha``, ..., ``r*alpha``.
     3163
     3164        INPUT:
     3165
     3166        - ``multi_indices`` - A list of tuples of positive integers, where
     3167          each tuple has length ``self.ring().ngens()``.
     3168        - ``numerical`` - (Optional; default=0) A natural number.
     3169          If positive, return numerical
     3170          approximations of coefficients with ``numerical`` digits of
     3171          accuracy.
     3172           
     3173        OUTPUT:
     3174
     3175        A dictionary of the form
     3176        (nu, Maclaurin coefficient of index nu of self).
     3177
     3178        EXAMPLES::
     3179       
     3180            sage: from sage.combinat.amgf import *
     3181            sage: R.<x, y, z> = PolynomialRing(QQ)
     3182            sage: H = (4 - 2*x - y - z) * (4 - x - 2*y - z)
     3183            sage: Hfac = H.factor()
     3184            sage: G = 16/Hfac.unit()
     3185            sage: F = FFPD(G, Hfac)
     3186            sage: alpha = vector([3, 3, 2])
     3187            sage: interval = [1, 2, 4]
     3188            sage: S = [r*alpha for r in interval]
     3189            sage: print F.maclaurin_coefficients(S, numerical=10)
     3190            {(6, 6, 4): 0.7005249476, (12, 12, 8): 0.5847732654,
     3191            (3, 3, 2): 0.7849731445}
     3192           
     3193        NOTES:
     3194       
     3195        Uses iterated univariate Maclaurin expansions.
     3196        Slow.
     3197       
     3198        AUTHORS:
     3199
     3200        - Alexander Raichev (2011-04-08, 2012-08-03)
     3201        """
     3202        R = self.ring()
     3203        if R is None:
     3204           return
     3205
     3206        d = self.dimension()
     3207
     3208        # Create biggest multi-index needed.
     3209        alpha = []
     3210        for i in xrange(d):
     3211           alpha.append(max((nu[i] for nu in multi_indices)))
     3212
     3213        # Compute Maclaurin expansion of self up to index alpha.
     3214        # Use iterated univariate expansions.
     3215        # Slow!
     3216        f = SR(self.quotient())
     3217        X = [SR(x) for x in R.gens()]
     3218        for i in xrange(d):
     3219           f = f.taylor(X[i], 0, alpha[i])
     3220        F = R(f)
     3221
     3222        # Collect coefficients.
     3223        coeffs = {}
     3224        X = R.gens()
     3225        for nu in multi_indices:
     3226           monomial = prod([X[i]**nu[i] for i in xrange(d)])
     3227           coeffs[tuple(nu)] = F.monomial_coefficient(monomial)
     3228           if numerical:
     3229               coeffs[tuple(nu)] = coeffs[tuple(nu)].n(digits=numerical)
     3230        return coeffs
     3231
     3232    def relative_error(self, approx, alpha, interval, exp_scale=Integer(1),
     3233                       digits=10):
     3234        r"""
     3235        Returns the relative error between the values of the Maclaurin
     3236        coefficients of ``self`` with multi-indices ``r alpha`` for ``r`` in
     3237        ``interval`` and the values of the functions (of the variable ``r``)
     3238        in ``approx``.
     3239
     3240        INPUT:
     3241
     3242        - ``approx`` - An individual or list of symbolic expressions in
     3243          one variable.
     3244        - ``alpha`` - A list of positive integers of length
     3245          ``self.ring().ngens()``
     3246        - ``interval`` - A list of positive integers.
     3247        - ``exp_scale`` - (Optional; default=1) A number.
     3248
     3249        OUTPUT:
     3250
     3251        A list whose entries are of the form
     3252        ``[r*alpha, a_r, b_r, err_r]`` for ``r`` in ``interval``.
     3253        Here ``r*alpha`` is a tuple; ``a_r`` is the ``r*alpha`` (multi-index)
     3254        coefficient of the Maclaurin series for ``self`` divided by
     3255        ``exp_scale**r``;
     3256        ``b_r`` is a list of the values of the functions in ``approx``
     3257        evaluated at ``r`` and divided by ``exp_scale**m``;
     3258        ``err_r`` is the list of relative errors
     3259        ``(a_r - f)/a_r`` for ``f`` in ``b_r``.
     3260        All outputs are decimal approximations.
     3261
     3262        EXAMPLES::
     3263
     3264            sage: from sage.combinat.amgf import *
     3265
     3266            sage: R.<x, y>= PolynomialRing(QQ)
     3267            sage: H = 1 - x - y - x*y
     3268            sage: Hfac = H.factor()
     3269            sage: G = 1/Hfac.unit()
     3270            sage: F = FFPD(G, Hfac)
     3271            sage: alpha = [1, 1]
     3272            sage: r = var('r')
     3273            sage: a1 = (0.573/sqrt(r))*5.83^r
     3274            sage: a2 = (0.573/sqrt(r) - 0.0674/r^(3/2))*5.83^r
     3275            sage: es = 5.83
     3276            sage: F.relative_error([a1, a2], alpha, [1, 2, 4, 8], es) # long time
     3277            Calculating errors table in the form
     3278            exponent, scaled Maclaurin coefficient, scaled asymptotic values,
     3279            relative errors...
     3280            [((1, 1), 0.5145797599, [0.5730000000, 0.5056000000],
     3281            [-0.1135300000, 0.01745066667]), ((2, 2), 0.3824778089,
     3282            [0.4051721856, 0.3813426871], [-0.05933514614, 0.002967810973]),
     3283            ((4, 4), 0.2778630595, [0.2865000000, 0.2780750000],
     3284            [-0.03108344267, -0.0007627515584]), ((8, 8), 0.1991088276,
     3285            [0.2025860928, 0.1996074055], [-0.01746414394, -0.002504047242])]
     3286           
     3287           
     3288        AUTHORS:
     3289
     3290        - Alexander Raichev (2009-05-18, 2011-04-18, 2012-08-03)
     3291        """
     3292
     3293        if not isinstance(approx, list):
     3294            approx = [approx]
     3295        av = approx[0].variables()[0]
     3296
     3297        print "Calculating errors table in the form"
     3298        print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..."           
     3299
     3300        # Get Maclaurin coefficients of self.
     3301        multi_indices = [r*vector(alpha) for r in interval]
     3302        mac = self.maclaurin_coefficients(multi_indices, numerical=digits)
     3303        #mac = self.old_maclaurin_coefficients(alpha, max(interval))
     3304        mac_approx = {}
     3305        stats = []
     3306        for r in interval:
     3307            beta = tuple(r*vector(alpha))
     3308            mac[beta] = (mac[beta]/exp_scale**r).n(digits=digits)
     3309            mac_approx[beta] = [(f.subs({av:r})/exp_scale**r).n(digits=digits)
     3310                                for f in approx]
     3311            stats_row = [beta, mac[beta], mac_approx[beta]]
     3312            if mac[beta] == 0:
     3313                 stats_row.extend([None for a in mac_approx[beta]])
     3314            else:
     3315                 stats_row.append([(mac[beta] - a)/mac[beta]
     3316                                   for a in mac_approx[beta]])
     3317            stats.append(tuple(stats_row))
     3318        return stats
     3319
     3320    @staticmethod
     3321    def coerce_point(R, p):
     3322        r"""
     3323        Coerce the keys of the dictionary ``p`` into the ring ``R``.
     3324        Assume this is possible.
     3325       
     3326        AUTHORS:
     3327
     3328        - Alexander Raichev (2009-05-18, 2011-04-18, 2012-08-03)
     3329        """
     3330        result = p
     3331        if p is not None and p.keys() and p.keys()[0].parent() != R:
     3332            try:
     3333                result = dict([(x, p[SR(x)]) for x in R.gens()])
     3334            except TypeError:
     3335                pass
     3336        return result
     3337           
     3338           
     3339class FFPDSum(list):
     3340    r"""
     3341    A list representing the sum of FFPD objects with distinct
     3342    denominator factorizations.
     3343
     3344    AUTHORS:
     3345
     3346    - Alexander Raichev (2012-06-25)
     3347    """
     3348    def __str__(self):
     3349        return str([r.list() for r in self])
     3350
     3351    def __eq__(self, other):
     3352        return sorted(self) == sorted(other)
     3353
     3354    def __ne__(self, other):
     3355        return not (self == other)
     3356
     3357    def ring(self):
     3358        r"""
     3359        Return the polynomial ring of the denominators of ``self``.
     3360        If ``self`` doesn't have any denominators, then return None.
     3361        """
     3362        for r in self:
     3363            R = r.ring()
     3364            if R is not None:
     3365                return R
     3366        return None
     3367
     3368    def whole_and_parts(self):
     3369        r"""
     3370        Rewrite ``self`` as a FFPDSum of a (possibly zero) polynomial
     3371        FFPD followed by reduced rational expression FFPDs.
     3372        Only useful for multivariate decompositions.
     3373
     3374        EXAMPLES::
     3375
     3376            sage: from sage.combinat.amgf import *
     3377
     3378            sage: R.<x, y> = PolynomialRing(QQ, 'x, y')
     3379            sage: f = x**2 + 3*y + 1/x + 1/y
     3380            sage: f = FFPD(quotient=f)
     3381            sage: print f
     3382            (x^3*y + 3*x*y^2 + x + y, [(y, 1), (x, 1)])
     3383            sage: print FFPDSum([f]).whole_and_parts()
     3384            [(x^2 + 3*y, []), (x + y, [(y, 1), (x, 1)])]
     3385
     3386        ::
     3387       
     3388            sage: R.<x, y> = PolynomialRing(QQ)
     3389            sage: f = cos(x)**2 + 3*y + 1/x + 1/y
     3390            sage: print f
     3391            1/x + 1/y + cos(x)^2 + 3*y
     3392            sage: G = f.numerator()
     3393            sage: H = R(f.denominator())
     3394            sage: f = FFPD(G, H.factor())
     3395            sage: print f
     3396            (x*y*cos(x)^2 + 3*x*y^2 + x + y, [(y, 1), (x, 1)])
     3397            sage: print FFPDSum([f]).whole_and_parts()
     3398            [(0, []), (x*y*cos(x)^2 + 3*x*y^2 + x + y, [(y, 1), (x, 1)])]
     3399        """
     3400        whole = 0
     3401        parts = []
     3402        R = self.ring()
     3403        for r in self:
     3404            # Since r has already passed through FFPD.__init__()'s reducing
     3405            # procedure, r is already in lowest terms.
     3406            # Check if can write r as a mixed fraction: whole + fraction.
     3407            p = r.numerator()
     3408            q = r.denominator()
     3409            if q == 1:
     3410                # r is already whole
     3411                whole += p
     3412            else:
     3413                try:
     3414                    # Coerce p into R and divide p by q
     3415                    p = R(p)
     3416                    a, b = p.quo_rem(q)
     3417                except TypeError:
     3418                    # p is not in R and so can't divide p by q
     3419                    a = 0
     3420                    b = p
     3421                whole += a
     3422                parts.append(FFPD(b, r.denominator_factored(), reduce_=False))
     3423        return FFPDSum([FFPD(whole, ())] + parts)
     3424
     3425    def combine_like_terms(self):
     3426        r"""
     3427        Combine terms in ``self`` with the same denominator.
     3428        Only useful for multivariate decompositions.
     3429
     3430        EXAMPLES::
     3431
     3432            sage: from sage.combinat.amgf import *
     3433
     3434            sage: R.<x, y>= PolynomialRing(QQ)
     3435            sage: f = FFPD(quotient=1/(x * y * (x*y + 1)))
     3436            sage: g = FFPD(quotient=x/(x * y * (x*y + 1)))
     3437            sage: s = FFPDSum([f, g, f])
     3438            sage: t = s.combine_like_terms()
     3439            sage: print s
     3440            [(1, [(y, 1), (x, 1), (x*y + 1, 1)]), (1, [(y, 1), (x*y + 1, 1)]),
     3441            (1, [(y, 1), (x, 1), (x*y + 1, 1)])]
     3442            sage: print t
     3443            [(1, [(y, 1), (x*y + 1, 1)]), (2, [(y, 1), (x, 1), (x*y + 1, 1)])]
     3444
     3445        ::
     3446       
     3447            sage: R.<x, y>= PolynomialRing(QQ)
     3448            sage: H = x * y * (x*y + 1)
     3449            sage: f = FFPD(1, H.factor())
     3450            sage: g = FFPD(exp(x + y), H.factor())
     3451            sage: s = FFPDSum([f, g])
     3452            sage: print s
     3453            [(1, [(y, 1), (x, 1), (x*y + 1, 1)]), (e^(x + y), [(y, 1), (x, 1),
     3454            (x*y + 1, 1)])]
     3455            sage: t = s.combine_like_terms()
     3456            sage: print t
     3457            [(e^(x + y) + 1, [(y, 1), (x, 1), (x*y + 1, 1)])]
     3458        """
     3459        if not self:
     3460            return self
     3461
     3462        # Combine like terms.
     3463        FFPDs = sorted(self)
     3464        new_FFPDs = []
     3465        temp = FFPDs[0]
     3466        for f in FFPDs[1:]:
     3467            if  temp.denominator_factored() == f.denominator_factored():
     3468                # Add f to temp.
     3469                num = temp.numerator() + f.numerator()
     3470                temp = FFPD(num, temp.denominator_factored())
     3471            else:
     3472                # Append temp to new_FFPDs and update temp.
     3473                new_FFPDs.append(temp)
     3474                temp = f
     3475        new_FFPDs.append(temp)
     3476        return FFPDSum(new_FFPDs)
     3477
     3478    def sum(self):
     3479        r"""
     3480        Return the sum of the FFPDs in ``self`` as a FFPD.
     3481
     3482        EXAMPLES::
     3483
     3484            sage: from sage.combinat.amgf import *
     3485
     3486            sage: R.<x, y> = PolynomialRing(QQ)
     3487            sage: df = (x, 1), (y, 1), (x*y + 1, 1)
     3488            sage: f = FFPD(2, df)
     3489            sage: g = FFPD(2*x*y, df)
     3490            sage: print FFPDSum([f, g])
     3491            [(2, [(y, 1), (x, 1), (x*y + 1, 1)]), (2, [(x*y + 1, 1)])]
     3492            sage: print FFPDSum([f, g]).sum()
     3493            (2, [(y, 1), (x, 1)]) 
     3494
     3495        ::
     3496           
     3497            sage: R.<x, y> = PolynomialRing(QQ)
     3498            sage: f = FFPD(cos(x), [(x, 2)])
     3499            sage: g = FFPD(cos(y), [(x, 1), (y, 2)])
     3500            sage: print FFPDSum([f, g])
     3501            [(cos(x), [(x, 2)]), (cos(y), [(y, 2), (x, 1)])]
     3502            sage: print FFPDSum([f, g]).sum()
     3503            (y^2*cos(x) + x*cos(y), [(y, 2), (x, 2)])
     3504
     3505
     3506        """       
     3507        if not self:
     3508            return self
     3509
     3510        # Compute the sum's numerator and denominator.
     3511        R = self.ring()
     3512        summy = sum((f.quotient() for f in self))
     3513        numer = summy.numerator()
     3514        denom = R(summy.denominator())
     3515
     3516        # Compute the sum's denominator factorization.       
     3517        # Could use the factor() command, but it's probably faster to use
     3518        # the irreducible factors of the denominators of self.
     3519        df = [] # The denominator factorization for the sum.
     3520        if denom == 1:
     3521            # Done
     3522            return FFPD(numer, df, reduce_=False)
     3523
     3524        factors = []
     3525        for f in self:
     3526            factors.extend([q for (q, e) in f.denominator_factored()])
     3527
     3528        # Eliminate repeats from factors and sort.
     3529        factors = sorted(list(set(factors)))
     3530
     3531        # The irreducible factors of denom lie in factors.
     3532        # Use this fact to build df.
     3533        for q in factors:
     3534            e = 0
     3535            quo, rem = denom.quo_rem(q)
     3536            while rem == 0:
     3537                e += 1
     3538                denom = quo
     3539                quo, rem = denom.quo_rem(q)
     3540            if e > 0:
     3541                df.append((q, e))
     3542        return FFPD(numer, df, reduce_=False)
     3543 No newline at end of file