Ticket #10519: trac_10519-v7.patch

File trac_10519-v7.patch, 144.0 KB (added by araichev, 4 years ago)
  • new file asymptotics_multivariate_generating_functions.py

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