Ticket #10519: trac_10519-v5.patch

File trac_10519-v5.patch, 144.6 KB (added by chapoton, 5 years ago)
  • new file sage/combinat/asymptotics_multivariate_generating_functions.py

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