Ticket #10519: trac_10519-asymptotic_multiple-v3.patch

File trac_10519-asymptotic_multiple-v3.patch, 136.0 KB (added by chapoton, 5 years ago)
  • sage/combinat/all.py

    # HG changeset patch
    # User Alex Raichev <alex.raichev@gmail.com>
    # Date 1344471313 -43200
    # Node ID 9fb72a90633451c8903259b6d5cfe38252902521
    # Parent  d06cf4b2215d37d3a87a58f65ac53234502dd471
    #10519: asymptotics for multivariate generating functions
    
    diff --git a/sage/combinat/all.py b/sage/combinat/all.py
    a b from degree_sequences import DegreeSeque 
    118118from cyclic_sieving_phenomenon import CyclicSievingPolynomial, CyclicSievingCheck
    119119
    120120from sidon_sets import sidon_sets
     121
     122from amgf import FFPD
  • new file sage/combinat/amgf.py

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