Ticket #10519: trac_10519-fixed.patch

File trac_10519-fixed.patch, 114.0 KB (added by davidloeffler, 5 years ago)

Patch against 5.0

  • sage/combinat/all.py

    # HG changeset patch
    # User Alex Raichev <tortoise.said@gmail.com>
    # Date 1331399002 0
    # Node ID b7819b98b3c6ed2078a83c4e08e0c59fe4b91357
    # Parent  68e89260148df131fec1708338ceba3ea964b2bb
    #10519: initial version
    
    diff --git a/sage/combinat/all.py b/sage/combinat/all.py
    a b  
    116116from cyclic_sieving_phenomenon import CyclicSievingPolynomial, CyclicSievingCheck
    117117
    118118from sidon_sets import sidon_sets
     119
     120from amgf import QuasiRationalExpression
  • new file sage/combinat/amgf.py

    diff --git a/sage/combinat/amgf.py b/sage/combinat/amgf.py
    new file mode 100644
    - +  
     1r"""
     2Asymptotics of Maclaurin coefficients of generating functions
     3
     4This code relates to analytic combinatorics. More specifically, it is a
     5collection of functions designed to compute asymptotics of Maclaurin
     6coefficients of certain classes of multivariate generating functions.
     7
     8The main function asymptotics() returns the first `N` terms of
     9the asymptotic expansion of the Maclaurin coefficients `F_{n\alpha}`
     10of the multivariate meromorphic function `F=G/H` as `n\to\infty`.
     11It assumes that `F` is holomorphic in a neighborhood of the origin,
     12that `H` is a polynomial, and that asymptotics in the direction of
     13`\alpha` (a tuple of positive integers) are controlled by convenient
     14smooth or multiple points.
     15
     16For an introduction to this subject, see [PeWi2008]_.
     17The main algorithms and formulas implemented here come from [RaWi2008a]_
     18and [RaWi2011]_.
     19
     20REFERENCES:
     21
     22.. [AiYu1983]  I. A. A\u\izenberg and A. P. Yuzhakov, "Integral
     23representations and residues in multidimensional complex analysis",
     24Translations of Mathematical Monographs, 58. American Mathematical
     25Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X.
     26
     27.. [DeLo2006]  Wolfram Decker and Christoph Lossen, "Computing in
     28algebraic geometry", Chapter 7.1, Springer-Verlag, 2006.
     29
     30.. [DiEm2005]  Alicia Dickenstein and Ioannis Z. Emiris (editors),
     31"Solving polynomial equations", Chapter 9.0, Springer-Verlag, 2005.
     32
     33.. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
     34several variables into partial fractions", Soviet Math. (Iz. VUZ)
     3522 (1978), no. 10, 35--38.
     36
     37.. [PeWi2008]  Robin Pemantle and Mark C. Wilson, "Twenty combinatorial
     38examples of asymptotics derived from multivariate generating
     39functions", SIAM Rev. (2008) vol. 50 (2) pp. 199-272.
     40
     41.. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
     42coefficients of multivariate generating functions: improvements
     43for smooth points", Electronic Journal of Combinatorics, Vol. 15
     44(2008), R89.
     45
     46.. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
     47coefficients of multivariate generating functions: improvements
     48for smooth points", To appear.
     49
     50AUTHORS:
     51
     52- Alex Raichev (2008-10-01) : Initial version
     53- Alex Raichev (2010-09-28) : Corrected many functions
     54- Alex Raichev (2010-12-15) : Updated documentation
     55- Alex Raichev (2011-03-09) : Fixed a division by zero bug in relative_error()
     56- Alex Raichev (2011-04-26) : Rewrote in object-oreinted style
     57- Alex Raichev (2011-05-06) : Fixed bug in cohomologous_integrand() and
     58    fixed _crit_cone_combo() to work in SR
     59
     60EXAMPLES:
     61
     62These come from [RaWi2008a]_ and [RaWi2011]_.
     63A smooth point example::
     64
     65    sage: R.<x,y>= PolynomialRing(QQ)
     66    sage: G= 1
     67    sage: H= 1-x-y-x*y
     68    sage: F= QuasiRationalExpression(G,H)
     69    sage: alpha= [3,2]
     70    sage: F.maclaurin_coefficients(alpha,8)
     71    {(21, 14): 1279919346549, (6, 4): 1289, (3, 2): 25, (18, 12): 19403906105, (24, 16): 85275509086721, (15, 10): 298199265, (9, 6): 75517, (12, 8): 4673345}
     72    sage: I= F.smooth_critical(alpha); I
     73    Ideal (y^2 + 3*y - 1, x - 2/3*y - 1/3) of Multivariate Polynomial Ring in x, y over Rational Field
     74    sage: s= solve(I.gens(),F.variables(),solution_dict=true); s
     75    [{y: -1/2*sqrt(13) - 3/2, x: -1/3*sqrt(13) - 2/3}, {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3}]
     76    sage: c= s[1]
     77    sage: asys= [F.asymptotics(c,alpha,n,numerical=3) for n in [1..2]]; asys
     78    Initializing auxiliary functions...
     79    Calculating derivatives of auxiallary functions...
     80    Calculating derivatives of more auxiliary functions...
     81    Calculating actions of the second order differential operator...
     82    Initializing auxiliary functions...
     83    Calculating derivatives of auxiallary functions...
     84    Calculating derivatives of more auxiliary functions...
     85    Calculating actions of the second order differential operator...
     86    [(0.369*71.2^n/sqrt(n), 71.2, 0.369/sqrt(n)), ((0.369/sqrt(n) - 0.0186/n^(3/2))*71.2^n, 71.2, 0.369/sqrt(n) - 0.0186/n^(3/2))]
     87    sage: F.relative_error([f[0] for f in asys],alpha,[1,2,4,8,16],asys[0][1])
     88    Calculating errors table in the form
     89    exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors...
     90    [[(3, 2), 0.351196289062500, [0.369140625000000, 0.351000000000000], [-0.0510940551757812, 0.00174000000000000]], [(6, 4), 0.254364013671875, [0.261021839148940, 0.254558441227157], [-0.0262197902254473, 0.000151011402221735]], [(12, 8), 0.181976318359375, [0.184570312500000, 0.182000000000000], [-0.0142545700073242, -0.00151000000000000]], [(24, 16), 0.129287719726562, [0.130510919574470, 0.129683383669613], [-0.00947434154988702, -0.00267741572252445]], [(48, 32), 0.0914535522460938, [0.0922851562500000, 0.0920000000000000], [-0.00909328460693359, -0.00592000000000000]]]
     91
     92Another smooth point example. Here it turns out that the terms 2--4 of the
     93asymptotic expansion are 0.
     94
     95::
     96
     97    sage: R.<x,y> = PolynomialRing(QQ)
     98    sage: q= 1/2
     99    sage: qq= q.denominator()
     100    sage: G= 1-q*x
     101    sage: H= 1-q*x +q*x*y -x^2*y
     102    sage: F= QuasiRationalExpression(G,H)
     103    sage: alpha= list(qq*vector([2,1-q])); alpha
     104    [4, 1]
     105    sage: I= F.smooth_critical(alpha); I
     106    Ideal (y^2 - 2*y + 1, x + 1/4*y - 5/4) of Multivariate Polynomial Ring in x, y over Rational Field
     107    sage: s= solve(I.gens(),F.variables(),solution_dict=true); s
     108    [{y: 1, x: 1}]
     109    sage: c= s[0]
     110    sage: asy= F.asymptotics(c,alpha,5); asy
     111    Initializing auxiliary functions...
     112    Calculating derivatives of auxiallary functions...
     113    Calculating derivatives of more auxiliary functions...
     114    Calculating actions of the second order differential operator...
     115    (1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*n^(1/3)) - 1/96*2^(1/3)*sqrt(3)*gamma(2/3)/(pi*n^(5/3)), 1, 1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*n^(1/3)) - 1/96*2^(1/3)*sqrt(3)*gamma(2/3)/(pi*n^(5/3)))
     116    sage: F.relative_error(asy[0],alpha,[1,2],asy[1])
     117    Calculating errors table in the form
     118    exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors...
     119    [[(4, 1), 0.187500000000000, [0.185581424449526], [0.0102324029358597]], [(8, 2), 0.152343750000000, [0.151986595969759], [0.00234439568568301]]]
     120
     121A multiple point example.
     122
     123::
     124
     125    sage: R.<x,y,z>= PolynomialRing(QQ)
     126    sage: G= 1
     127    sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y))
     128    sage: F= QuasiRationalExpression(G,H)
     129    sage: I= F.singular_points(); I
     130    Ideal (x*y + x - 1, y^2 - 2*y*z + 2*y - z + 1, x*z + y - 2*z + 1) of Multivariate Polynomial Ring in x, y, z over Rational Field
     131    sage: c= {x: 1/2, y: 1, z: 4/3}
     132    sage: F.is_cmp(c)
     133    [({y: 1, z: 4/3, x: 1/2}, True, 'all good')]
     134    sage: cc= F.critical_cone(c); cc
     135    ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N)
     136    sage: alpha= [8,3,3]
     137    sage: alpha in cc[1]
     138    True
     139    sage: asy= F.asymptotics(c,alpha,2); asy
     140    Initializing auxiliary functions...
     141    Calculating derivatives of auxiliary functions...
     142    Calculating derivatives of more auxiliary functions...
     143    Calculating second-order differential operator actions...
     144    (1/172872*(24696*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n)) - 1231*sqrt(3)*sqrt(7)/(sqrt(pi)*n^(3/2)))*108^n, 108, 1/7*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n)) - 1231/172872*sqrt(3)*sqrt(7)/(sqrt(pi)*n^(3/2)))
     145
     146Another multiple point example.
     147
     148::
     149
     150    sage: R.<x,y,z>= PolynomialRing(QQ)
     151    sage: G= 16
     152    sage: H= (4-2*x-y-z)^2*(4-x-2*y-z)
     153    sage: F= QuasiRationalExpression(G,H)
     154    sage: var('a1,a2,a3')
     155    (a1, a2, a3)
     156    sage: F.cohomologous_integrand([a1,a2,a3])
     157    [[-16*(a2*z - 2*a3*y)*n/(y*z) + 16*(2*y - z)/(y*z), [x + 2*y + z - 4, 1], [2*x + y + z - 4, 1]]]
     158    sage: I= F.singular_points(); I
     159    Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field
     160    sage: c= {x:1, y:1, z:1}
     161    sage: F.is_cmp(c)
     162    [({y: 1, z: 1, x: 1}, True, 'all good')]
     163    sage: alpha= [3,3,2]
     164    sage: cc= F.critical_cone(c); cc
     165    ([(1, 2, 1), (2, 1, 1)], 2-d cone in 3-d lattice N)
     166    sage: alpha in cc[1]
     167    True
     168    sage: asy= F.asymptotics(c,alpha,2); asy
     169    Initializing auxiliary functions...
     170    Calculating derivatives of auxiliary functions...
     171    Calculating derivatives of more auxiliary functions...
     172    Calculating second-order differential operator actions...
     173    (4/3*sqrt(3)*sqrt(n)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(n)), 1, 4/3*sqrt(3)*sqrt(n)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(n)))
     174
     175"""
     176#*****************************************************************************
     177#       Copyright (C) 2008 Alexander Raichev <tortoise.said@gmail.com>
     178#
     179#  Distributed under the terms of the GNU General Public License (GPL)
     180#                  http://www.gnu.org/licenses/
     181#*****************************************************************************
     182
     183from copy import copy, deepcopy
     184
     185# Sage libraries
     186from sage.calculus.functional                           import diff
     187from sage.calculus.functions                            import jacobian
     188from sage.calculus.var                                  import function, var
     189from sage.combinat.cartesian_product                    import CartesianProduct
     190from sage.combinat.combinat                             import stirling_number1
     191from sage.combinat.permutation                          import Permutation
     192from sage.combinat.tuple                                import UnorderedTuples
     193from sage.functions.log                                 import exp, log
     194from sage.functions.other                               import factorial, gamma, sqrt
     195from sage.geometry.cone                                 import Cone
     196from sage.matrix.constructor                            import matrix
     197from sage.misc.misc                                     import add, ellipsis_range, mul
     198from sage.misc.misc_c                                   import prod
     199from sage.misc.mrange                                   import cartesian_product_iterator, mrange
     200from sage.modules.free_module_element                   import vector
     201from sage.rings.arith                                   import binomial
     202from sage.rings.all                                     import CC
     203from sage.rings.fraction_field                          import FractionField
     204from sage.rings.integer                                 import Integer
     205from sage.rings.integer_ring                            import ZZ
     206from sage.rings.polynomial.polynomial_ring_constructor  import PolynomialRing
     207from sage.rings.rational_field                          import QQ
     208from sage.sets.set                                      import Set
     209from sage.structure.sage_object                         import SageObject
     210from sage.symbolic.constants                            import pi
     211from sage.symbolic.relation                             import solve
     212from sage.symbolic.ring                                 import SR
     213
     214
     215# Functions to incorporate later into existing Sage classes.
     216def algebraic_dependence(fs):
     217    r"""
     218    This function returns an irreducible annihilating polynomial for the
     219    polynomials in `fs`, if those polynomials are algebraically dependent.
     220    Otherwise it returns 0.
     221
     222    INPUT:
     223
     224    - ``fs`` - A list of polynomials `f_1,\ldots,f_r' from a common polynomial
     225        ring.
     226
     227    OUTPUT:
     228
     229    If the polynomials in `fs` are algebraically dependent, then the
     230    output is an irreducible polynomial `g` in `K[T_1,\ldots,T_r]` such that
     231    `g(f_1,\ldots,f_r) = 0`.
     232    Here `K` is the coefficient ring of self and `T_1,\ldots,T_r` are
     233    indeterminates different from those of self.
     234    If the polynomials in `fs` are algebraically independent, then the output
     235    is the zero polynomial.
     236
     237    EXAMPLES::
     238
     239        sage: from sage.combinat.amgf import *
     240        sage: R.<x>= PolynomialRing(QQ)
     241        sage: fs= [x-3, x^2+1]
     242        sage: g= algebraic_dependence(fs); g
     243        10 + 6*T0 - T1 + T0^2
     244        sage: g(fs)
     245        0
     246
     247    ::
     248
     249        sage: R.<w,z>= PolynomialRing(QQ)
     250        sage: fs= [w,  (w^2 + z^2 - 1)^2, w*z - 2]
     251        sage: g= algebraic_dependence(fs); g
     252        16 + 32*T2 - 8*T0^2 + 24*T2^2 - 8*T0^2*T2 + 8*T2^3 + 9*T0^4 - 2*T0^2*T2^2 + T2^4 - T0^4*T1 + 8*T0^4*T2 - 2*T0^6 + 2*T0^4*T2^2 + T0^8
     253        sage: g(fs)
     254        0
     255        sage: g.parent()
     256        Multivariate Polynomial Ring in T0, T1, T2 over Rational Field
     257
     258    ::
     259
     260        sage: R.<x,y,z>= PolynomialRing(GF(7))
     261        sage: fs= [x*y,y*z]
     262        sage: g= algebraic_dependence(fs); g
     263        0
     264
     265    AUTHORS:
     266
     267    - Alex Raichev (2011-01-06)
     268    """
     269    r= len(fs)
     270    R= fs[Integer(0) ].parent()
     271    B= R.base_ring()
     272    Xs= list(R.gens())
     273    d= len(Xs)
     274
     275    # Expand R by r new variables.
     276    T= 'T'
     277    while T in [str(x) for x in Xs]:
     278        T= T+'T'
     279    Ts= [T + str(j) for j in range(r)]
     280    RR= PolynomialRing(B,d+r,tuple(Xs+Ts))
     281    Vs= list(RR.gens())
     282    Xs= Vs[Integer(0) :d]
     283    Ts= Vs[d:]
     284
     285    # Find an irreducible annihilating polynomial g for the fs if there is one.
     286    J= RR.ideal([ Ts[j] -RR(fs[j]) for j in range(r)])
     287    JJ= J.elimination_ideal(Xs)
     288    g= JJ.gens()[Integer(0) ]
     289
     290    # Shrink the ambient ring of g to include only Ts.
     291    # I chose the negdeglex order simply because i find it useful in my work.
     292    RRR= PolynomialRing(B,r,tuple(Ts),order='negdeglex')
     293    return RRR(g)
     294#-------------------------------------------------------------------------------
     295def partial_fraction_decomposition(f):
     296    r"""
     297    Return a partial fraction decomposition of the rational expression `f`.
     298    Works for univariate and mulitivariate `f`.
     299
     300    INPUT:
     301
     302    - ``f`` - An element of the field of fractions `F` of a polynomial ring
     303        `R` whose coefficients ring is a field.
     304        In the univariate case, the coefficient ring doesn't have to be a field.
     305
     306    OUTPUT:
     307
     308    A tuple `whole,parts` where `whole \in R` and `parts` is the list of
     309    terms (in `F`) of a partial fraction decomposition of `f - whole`.
     310    See the notes below for more details.
     311
     312    EXAMPLES::
     313
     314        sage: from sage.combinat.amgf import *
     315        sage: S.<t> = QQ[]
     316        sage: q = 1/(t+1) + 2/(t+2) + 3/(t-3)
     317        sage: whole, parts = partial_fraction_decomposition(q); parts
     318        [3/(t - 3), 1/(t + 1), 2/(t + 2)]
     319        sage: whole +sum(parts) == q
     320        True
     321
     322    Notice that there is a whole part below despite the appearance of q ::
     323
     324        sage: R.<x,y>= PolynomialRing(QQ)
     325        sage: q= 1/( x *y *(x*y-1) )
     326        sage: whole,parts= partial_fraction_decomposition(q)
     327        sage: whole, parts
     328        (-1, [x*y/(x*y - 1), (-1)/(x*y)])
     329        sage: q == whole +sum(parts)
     330        True
     331
     332    ::
     333
     334        sage: R.<x,y>= PolynomialRing(QQ)
     335        sage: q= x +1/x +1/(x*y-2) +1/(x^2+y^2-1)
     336        sage: whole,parts= partial_fraction_decomposition(q)
     337        sage: whole, parts
     338        (x, [(1/2*x^3*y^2 + 1/2*x*y^4 + 1/2*x^3*y + 1/2*x^2*y^2 + 1/2*x*y^3 - x^2*y - 1/2*x*y^2 - y^3 - 3/2*x*y + y)/(x^3*y + x*y^3 - 2*x^2 - x*y - 2*y^2 + 2), (-1/2*x^3*y - 1/2*x*y^3 - 1/2*x^3 - 1/2*x^2*y - 1/2*x*y^2 + x^2 + 1/2*x*y + y^2 + 3/2*x - 1)/(x^3 + x*y^2 - x)])
     339        sage: whole +sum(parts)==q
     340        True
     341
     342    ::
     343
     344        sage: R.<x,y,z>= PolynomialRing(QQ)
     345        sage: q= 1/x +1/(x*y-z-2)^2 +1/(x^2+y^2 +z^2-1)^3 +1/(x*y-3)
     346        sage: whole,parts= partial_fraction_decomposition(q)
     347        sage: whole +sum(parts)==q
     348        True
     349
     350    NOTES:
     351
     352    In the case of univariate `f` this function calls the old univariate
     353    partial fraction decomposition function.
     354    In the multivariate case, it uses a different notion of and algorithm for
     355    partial fraction decompositions.
     356
     357    Let `f= P/Q` where `P,Q \in R`, let `Q_1^{e_1} \cdots Q_m^{e_m}` be the
     358    unique factorization of `Q` in `R` into irreducible factors, and let `d` be
     359    the number of indeterminates of `R`.
     360    Then `f` can be written as a sum `\sum P_A/\prod_{j\in A} Q_j^{b_j}`,
     361    where the `b_j \le e_j` are positive integers, the `P_A` are in `R`, and
     362    the sum is taken over all size `\le m` subsets `A \subseteq \{ 1, \ldots, d \}`
     363    such that `S:= \{ Q_j : j\in A \}` is algebraically independent and the
     364    ideal generated by `S` is not all of `R` (that is, by the Nullstellensatz,
     365    the polynomials of `S` have no common roots in the algebraic closure of the
     366    coefficient field of `R`).
     367    Following [Lein1978]_, i call any such decomposition of `f` a
     368    `\emph{partial fraction decomposition}`.
     369
     370    The algorithm used below comes from Theorem 1, Lemma 2, and Lemma 3 of
     371    [Lein1978]_.
     372    By the way, that paper has a typo in equation (c) on the
     373    second page and equation (b) on the third page.
     374    The right sides of (c) and (b) should be multiplied by `P`.
     375
     376    REFERENCES:
     377
     378    .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
     379    several variables into partial fractions", Soviet Math. (Iz. VUZ)
     380    22 (1978), no. 10, 35--38.
     381
     382    AUTHORS:
     383
     384    - Alex Raichev (2011-01-10)
     385    """
     386    R= f.denominator().parent()
     387    d= len(R.gens())
     388    if d==Integer(1) :
     389        return f.partial_fraction_decomposition()
     390    Q= f.denominator()
     391    whole,P= f.numerator().quo_rem(Q)
     392    parts= [format_quotient(Integer(1) /Q)]
     393    # Decompose via nullstellensatz trick
     394    # (which is faster than the algebraic dependence trick)
     395    parts= decompose_via_nullstellensatz(parts)
     396    # Decompose via algebraic dependence trick
     397    parts= decompose_via_algebraic_dependence(parts)
     398    # Rewrite parts back in terms of rational expressions
     399    new_parts=[]
     400    for p in parts:
     401        f= unformat_quotient(p)
     402        if f.denominator() == Integer(1) :
     403            whole= whole +f
     404        else:
     405            new_parts.append(f)
     406    # Put P back in
     407    new_parts= [P*f for f in new_parts]
     408    return whole, new_parts
     409#-------------------------------------------------------------------------------
     410def combine_like_terms(parts,rationomial=True):
     411    r"""
     412    Combines like terms of the fractions represented by parts.
     413    For use by partial_fraction_decomposition() above.
     414
     415    EXAMPLES:
     416
     417    sage: from sage.combinat.amgf import *
     418    sage: R.<x,y>= PolynomialRing(QQ)
     419    sage: parts =[[1, [x*y, 1]], [x, [x*y, 1]]]
     420    sage: combine_like_terms(parts)
     421    [[x + 1, [y, 1], [x, 1]]]
     422
     423    ::
     424
     425    sage: R.<x>= PolynomialRing(QQ)
     426    sage: parts =[[1, [x, 1]], [x-1, [x, 1]]]
     427    sage: combine_like_terms(parts)
     428    [[1]]
     429
     430    AUTHORS:
     431
     432    - Alex Raichev (2011-01-10)
     433    """
     434
     435    if parts == []: return parts
     436    # Sort parts by denominators
     437    new_parts= sorted(parts, key= lambda p:p[Integer(1) :])
     438    # Combine items of parts with same denominators.
     439    newnew_parts=[]
     440    left,right=Integer(0) ,Integer(1)
     441    glom= new_parts[left][Integer(0) ]
     442    while right <= len(new_parts)-Integer(1) :
     443        if new_parts[left][Integer(1) :] == new_parts[right][Integer(1) :]:
     444            glom= glom +new_parts[right][Integer(0) ]
     445        else:
     446            newnew_parts.append([glom]+new_parts[left][Integer(1) :])
     447            left= right
     448            glom= new_parts[left][Integer(0) ]
     449        right= right +Integer(1)
     450    if glom != Integer(0) :
     451        newnew_parts.append([glom]+new_parts[left][Integer(1) :])
     452    if rationomial:
     453        # Reduce fractions in newnew_parts.
     454        # Todo: speed up below by working directly with newnew_parts and
     455        # thereby make fewer calls to format_quotient() which in turn
     456        # calls factor().
     457        newnew_parts= [format_quotient(unformat_quotient(part)) for part in newnew_parts]
     458    return newnew_parts
     459#-------------------------------------------------------------------------------
     460def decompose_via_algebraic_dependence(parts):
     461    r"""
     462    Returns a decomposition of parts.
     463    Used by partial_fraction_decomposition() above.
     464    Implements Lemma 2 of [Lein1978]_.
     465    Recursive.
     466
     467    REFERENCES:
     468
     469    .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
     470    several variables into partial fractions", Soviet Math. (Iz. VUZ)
     471    22 (1978), no. 10, 35--38.
     472
     473    AUTHORS:
     474
     475    - Alex Raichev (2011-01-10)
     476    """
     477    decomposing_done= True
     478    new_parts= []
     479    for p in parts:
     480        p_parts= [p]
     481        P= p[Integer(0) ]
     482        Qs= p[Integer(1) :]
     483        m= len(Qs)
     484        G= algebraic_dependence([q for q,e in Qs])
     485        if G:
     486            # Then the denominator factors are algebraically dependent
     487            # and so we can decompose p.
     488            decomposing_done= False
     489            P= p[Integer(0) ]
     490            Qs= p[Integer(1) :]
     491
     492            # Todo: speed up step below by using
     493            # G to calculate F.  See [Lein1978]_ Lemma 1.
     494            F= algebraic_dependence([q**e for q,e in Qs])
     495            new_vars= F.parent().gens()
     496
     497            # Note that new_vars[j] corresponds to Qs[j] so that
     498            # F([q^e for q,e in Qs]) = 0.
     499            # Assuming below that F.parent() has negdeglex term order
     500            # so that F.lt() is indeed the monomial we want.
     501            FF= (F.lt() -F)/(F.lc())
     502            numers= map(mul,zip(FF.coefficients(),FF.monomials()))
     503            e= list(F.lt().exponents())[Integer(0) :m]
     504            denom= [[new_vars[j], e[Integer(0) ][j]+Integer(1) ] for j in range(m)]
     505
     506            # Before making things messy by substituting in Qs,
     507            # reduce terms and combine like terms.
     508            p_parts_temp= [format_quotient(unformat_quotient([a]+denom)) for a in numers]
     509            p_parts_temp= combine_like_terms(p_parts_temp)
     510
     511            # Substitute Qs into new_p.
     512            Qpowsub= dict([(new_vars[j],Qs[j][Integer(0) ]**Qs[j][Integer(1) ]) for j in range(m)])
     513            p_parts=[]
     514            for x in p_parts_temp:
     515                y= P*F.parent()(x[Integer(0) ]).subs(Qpowsub)
     516                yy=[]
     517                for xx in x[Integer(1) :]:
     518                    if xx[Integer(0) ] in new_vars:
     519                        j= new_vars.index(xx[Integer(0) ])
     520                        yy.append([Qs[j][Integer(0) ],Qs[j][Integer(1) ]*xx[Integer(1) ]])
     521                    else:
     522                        # Occasionally xx[0] is an integer.
     523                        yy.append(xx)
     524                p_parts.append([y]+yy)
     525        # Done one step of decomposing p.  Add it to new_parts.
     526        new_parts.extend(p_parts)
     527    if decomposing_done:
     528        return new_parts
     529    else:
     530        return decompose_via_algebraic_dependence(new_parts)
     531#-------------------------------------------------------------------------------
     532def decompose_via_cohomology(parts):
     533    r"""
     534    Given each nice (described below) differential form
     535    `(P/Q) dx_1 \wedge\cdots\wedge dx_d` in `parts`,
     536    this function returns a differential form equivalent in De Rham
     537    cohomology that has no repeated factors in the denominator.
     538
     539    INPUT:
     540
     541    - ``parts`` - A list of the form `[chunk_1,\ldots,chunk_r]`, where each
     542        `chunk_j` has the form `[P,[Q_1,e_1],\ldots,[Q_m,e_m]]`,
     543        `Q_1,\ldots,Q_m` are irreducible elements of a common polynomial
     544        ring `R` such that their corresponding algebraic varieties
     545        `\{x\in F^d : B_j(x)=0\}` intersect transversely (where `F` is the
     546        algebraic closure of the field of coefficients of `R`),
     547        `e_1,\ldots,e_m` are positive integers, `m \le d`, and
     548        `P` is a symbolic expression in some of the indeterminates of `R`.
     549        Here `[P,[Q_1,e_1],\ldots,[Q_m,e_m]]` represents the fraction
     550        `P/(Q_1^e_1 \cdots Q_m^e_m)`.
     551
     552    OUTPUT:
     553
     554    A list of the form `[chunky_1,\ldots,chunky_s]`, where each
     555    `chunky_j` has the form `[P,[Q_1,1],\ldots,[Q_m,1]]`.
     556
     557    EXAMPLES::
     558
     559        sage: from sage.combinat.amgf import *
     560        sage: R.<x,y>= PolynomialRing(QQ)
     561        sage: decompose_via_cohomology([[ 1, [x,3] ]])
     562        []
     563
     564    ::
     565
     566        sage: R.<x,y>= PolynomialRing(QQ)
     567        sage: decompose_via_cohomology([[ 1, [x,3], [x*y-1,2] ]])
     568        [[-3*y^2, [x, 1], [x*y - 1, 1]], [-5*y^3, [x*y - 1, 1]]]
     569
     570    NOTES:
     571
     572    This is a recursive function thats stops calling itself when all the
     573    `e_j` equal 1 or `parts == []`.
     574    The algorithm used here is that of Theorem 17.4 of
     575    [AiYu1983]_.
     576    The algebraic varieties `\{x\in F^d : Q_j(x)=0\}`
     577    (where `F` is the algebraic closure of the field of coefficients of `R`)
     578    corresponding to the `Q_j` __intersect transversely__ iff for each
     579    point `c` of their intersection and for all `k \le m`,
     580    the Jacobian matrix of any `k` polynomials from
     581    `\{Q_1,\ldots,Q_m\}` has rank equal to `\min\{k,d\}` when evaluated at
     582    `c`.
     583
     584    REFERENCES:
     585
     586    .. [AiYu1983]  I. A. A\u\izenberg and A. P. Yuzhakov, "Integral
     587    representations and residues in multidimensional complex analysis",
     588    Translations of Mathematical Monographs, 58. American Mathematical
     589    Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X.
     590
     591    AUTHORS:
     592
     593    - Alex Raichev (2008-10-01, 2011-01-15)
     594    """
     595    if parts == []: return parts
     596    decomposing_done= True
     597    new_parts= []
     598    R= parts[Integer(0) ][Integer(1) ][Integer(0) ].parent()
     599    V= list(R.gens())
     600    for p in parts:
     601        p_parts= [p]
     602        P= p[Integer(0) ]
     603        Qs= p[Integer(1) :]
     604        m= len(Qs)
     605        if sum([e for q,e in Qs]) > m:
     606            # Then we can decompose p
     607            p_parts= []
     608            decomposing_done=False
     609            dets= []
     610            vars_subsets= Set(V).subsets(m)
     611            for v in vars_subsets:
     612                # Sort variables so that first polynomial ring indeterminate
     613                # declared is first in vars list.
     614                v= sorted(v,reverse=True)
     615                jac= jacobian([q for q,e in Qs],v)
     616                dets.append(jac.determinant())
     617            # Get a Nullstellensatz certificate.
     618            L= R(Integer(1) ).lift(R.ideal([q for q,e in Qs] +dets))
     619            for j in range(m):
     620                if L[j] != Integer(0) :
     621                    # Make a copy of (and not a reference to) the nested list Qs.
     622                    new_Qs = deepcopy(Qs)
     623                    if new_Qs[j][Integer(1) ] > Integer(1) :
     624                        new_Qs[j][Integer(1) ]= new_Qs[j][Integer(1) ] -Integer(1)
     625                    else:
     626                        del new_Qs[j]
     627                    p_parts.append( [P*L[j]] +new_Qs )
     628            for k in range(vars_subsets.cardinality()):
     629                if L[m+k] != Integer(0) :
     630                    new_Qs = deepcopy(Qs)
     631                    for j in range(m):
     632                        if new_Qs[j][Integer(1) ] > Integer(1) :
     633                            new_Qs[j][Integer(1) ]= new_Qs[j][Integer(1) ] -Integer(1)
     634                            v=sorted(vars_subsets[k],reverse=True)
     635                            jac= jacobian([SR(P*L[m+k])] +[ SR(Qs[jj][Integer(0) ]) for \
     636                            jj in range(m) if jj !=j], [SR(vv) for vv in v])
     637                            det= jac.determinant()
     638                            if det != Integer(0) :
     639                                p_parts.append([permutation_sign(v,V) \
     640                                  *(-Integer(1) )**j/new_Qs[j][Integer(1) ] *det] +new_Qs)
     641                            break
     642        # Done one step of decomposing p.  Add it to new_parts.
     643        new_parts.extend(p_parts)
     644        new_parts= combine_like_terms(new_parts,rationomial=False)
     645    if decomposing_done:
     646        return new_parts
     647    else:
     648        return decompose_via_cohomology(new_parts)
     649#-------------------------------------------------------------------------------
     650def decompose_via_nullstellensatz(parts):
     651    r"""
     652    Returns a decomposition of parts.
     653    Used by partial_fraction_decomposition() above.
     654    Implements Lemma 3 of [Lein1978]_.
     655    Recursive.
     656
     657    REFERENCES:
     658
     659    .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
     660    several variables into partial fractions", Soviet Math. (Iz. VUZ)
     661    22 (1978), no. 10, 35--38.
     662
     663    AUTHORS:
     664
     665    - Alex Raichev (2011-01-10)
     666    """
     667    decomposing_done= True
     668    new_parts= []
     669    R= parts[Integer(0) ][Integer(0) ].parent()
     670    for p in parts:
     671        p_parts= [p]
     672        P= p[Integer(0) ]
     673        Qs= p[Integer(1) :]
     674        m= len(Qs)
     675        if R(Integer(1) ) in R.ideal([q for q,e in Qs]):
     676            # Then we can decompose p.
     677            decomposing_done= False
     678            L= R(Integer(1) ).lift(R.ideal([q**e for q,e in Qs]))
     679            p_parts= [ [P*L[i]] + \
     680            [[Qs[j][Integer(0) ],Qs[j][Integer(1) ]] for j in range(m) if j != i] \
     681            for i in range(m) if L[i]!=Integer(0) ]
     682            # The procedure above yields no like terms to combine.
     683        # Done one step of decomposing p.  Add it to new_parts.
     684        new_parts.extend(p_parts)
     685    if decomposing_done:
     686        return new_parts
     687    else:
     688        return decompose_via_nullstellensatz(new_parts)
     689#-------------------------------------------------------------------------------
     690def format_quotient(f):
     691    r"""
     692    Formats `f` for use by partial_fraction_decomposition() above.
     693
     694    AUTHORS:
     695
     696    - Alex Raichev (2011-01-10)
     697    """
     698    P= f.numerator()
     699    Q= f.denominator()
     700    Qs= Q.factor()
     701    if Qs.unit() != Integer(1) :
     702        P= P/Qs.unit()
     703    Qs= sorted([[q,e] for q,e in Qs])  # sorting for future bookkeeping
     704    return [P]+Qs
     705#-------------------------------------------------------------------------------
     706def permutation_sign(v,vars):
     707    r"""
     708    This function returns the sign of the permutation on `1,\ldots,len(vars)`
     709    that is induced by the sublist `v` of `vars`.
     710    For internal use by _cohom_equiv_main().
     711
     712    INPUT:
     713
     714    - ``v`` - A sublist of `vars`.
     715    - ``vars`` - A list of predefined variables or numbers.
     716
     717    OUTPUT:
     718
     719    The sign of the permutation obtained by taking indices (and adding 1)
     720    within `vars` of the list `v,w`, where `w` is the list `vars` with the
     721    elements of `v` removed.
     722
     723    EXAMPLES::
     724
     725        sage: from sage.combinat.amgf import *
     726        sage: vars= ['a','b','c','d','e']
     727        sage: v= ['b','d']
     728        sage: permutation_sign(v,vars)
     729        -1
     730        sage: v= ['d','b']
     731        sage: permutation_sign(v,vars)
     732        1
     733
     734    AUTHORS:
     735
     736    - Alex Raichev (2008-10-01)
     737    """
     738    # Convert variable lists to lists of numbers in {1,...,len(vars)}
     739    A= [x+Integer(1)  for x in range(len(vars))]
     740    B= [vars.index(x)+Integer(1)  for x in v]
     741    C= list(Set(A).difference(Set(B)))
     742    C.sort()
     743    P= Permutation(B+C)
     744    return P.signature()
     745#-------------------------------------------------------------------------------
     746def unformat_quotient(part):
     747    r"""
     748    Unformats `f` for use by partial_fraction_decomposition() above.
     749    Inverse of format_quotient() above.
     750
     751    AUTHORS:
     752
     753    - Alex Raichev (2011-01-10)
     754    """
     755    P= part[Integer(0) ]
     756    Qs= part[Integer(1) :]
     757    Q= prod([q**e for q,e in Qs])
     758    return P/Q
     759#===============================================================================
     760# Class for calculation of asymptotics of multivariate generating functions.
     761class QuasiRationalExpression(SageObject):
     762    "Store an expression G/H, where H comes from a polynomial ring R and \
     763    G comes from R or the Symbolic Ring."
     764    def __init__(self, G, H):
     765        # Store important information about object as attributes of self.
     766        # G, H, H's ring, ring dimension, H's factorization.
     767        self._G = G
     768        self._H = H
     769        R= H.parent()
     770        self._R = R
     771        self._d = len(R.gens())
     772        self._Hfac= list(H.factor())
     773
     774        # Variables of self as elements of the SR.
     775        # Remember that G might be in SR and not in R.
     776        try:
     777            # This fails if G is a constant, for example.
     778            Gv= Set([R(x) for x in G.variables()])
     779        except:
     780            Gv= Set([])
     781        try:
     782            # This fails if H is a constant, for example.
     783            Hv= Set(H.variables())
     784        except:
     785            Hv= Set([])
     786        # Preserve the ring ordering of the variables which some methods below
     787        # depends upon.
     788        V= sorted(list(Gv.union(Hv)),reverse=True)
     789        self._variables = tuple([SR(x) for x in V])
     790#-------------------------------------------------------------------------------
     791# Keeping methods in alphabetical order (ignoring initial single underscores)
     792    def asymptotics(self,c,alpha,N,numerical=Integer(0) ,asy_var=None):
     793        r"""
     794        This function returns the first `N` terms of the asymptotic expansion
     795        of the Maclaurin coefficients `F_{n\alpha}` of the
     796        multivariate meromorphic function `F=G/H` as `n\to\infty`,
     797        where `F = self`.
     798        It assumes that `F` is holomorphic in a neighborhood of the origin,
     799        that `H` is a polynomial, and that `c` is a convenient strictly minimal
     800        smooth or multiple of `F` that is critical for `\alpha`.
     801
     802        INPUT:
     803
     804        - ``alpha`` - A `d`-tuple of positive integers or, if `c` is a smooth
     805            point, possibly of symbolic variables.
     806        - ``c`` - A dictionary with keys `self._variables` and values from a
     807            superfield of the field of `self._R.base_ring()`.
     808        - ``N`` - A positive integer.
     809        - ``numerical`` - A natural number (default: 0).
     810            If k=numerical > 0, then a numerical approximation of the coefficients
     811            of `F_{n\alpha}` with k digits of precision will be returned.
     812            Otherwise exact values will be returned.
     813        - ``asy_var`` - A  symbolic variable (default: None).
     814            The variable of the asymptotic expansion.
     815            If none is given, `var('n')` will be assigned.
     816
     817        OUTPUT:
     818
     819        The tuple `(asy,exp_part,subexp_part)`, where `asy` is first `N` terms
     820        of the asymptotic expansion of the Maclaurin coefficients `F_{n\alpha}`
     821        of the function `F=self` as `n\to\infty`, `exp_part^n` is the exponential
     822        factor of `asy`, and `subexp_part` is the subexponential factor of
     823        `asy`.
     824
     825        EXAMPLES::
     826
     827        A smooth point example ::
     828
     829            sage: R.<x,y>= PolynomialRing(QQ)
     830            sage: G= 1
     831            sage: H= 1-x-y-x*y
     832            sage: F= QuasiRationalExpression(G,H)
     833            sage: alpha= [3,2]
     834            sage: c= {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3}
     835            sage: F.asymptotics(c,alpha,2,numerical=3)
     836            Initializing auxiliary functions...
     837            Calculating derivatives of auxiallary functions...
     838            Calculating derivatives of more auxiliary functions...
     839            Calculating actions of the second order differential operator...
     840            ((0.369/sqrt(n) - 0.0186/n^(3/2))*71.2^n, 71.2, 0.369/sqrt(n) - 0.0186/n^(3/2))
     841
     842        A multiple point example ::
     843
     844            sage: R.<x,y,z>= PolynomialRing(QQ)
     845            sage: G= 1
     846            sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y))
     847            sage: F= QuasiRationalExpression(G,H)
     848            sage: c= {z: 4/3, y: 1, x: 1/2}
     849            sage: alpha= [8,3,3]
     850            sage: F.asymptotics(c,alpha,1)
     851            Initializing auxiliary functions...
     852            Calculating derivatives of auxiliary functions...
     853            Calculating derivatives of more auxiliary functions...
     854            Calculating second-order differential operator actions...
     855            (1/7*sqrt(3)*sqrt(7)*108^n/(sqrt(pi)*sqrt(n)), 108, 1/7*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n)))
     856
     857        NOTES:
     858
     859        A zero `c` of `H` is __strictly minimal__ if there is no zero `x` of `H`
     860        such that `x_j < c_j` for all `0 \le j < d`.
     861        For definitions of the terms "smooth critical point for `\alpha`" and
     862        "multiple critical point for `\alpha`",
     863        see the documentation for _asymptotics_main_smooth() and
     864        _asymptotics_main_multiple(), which are the functions that do most of the
     865        work.
     866
     867        ALGORITHM:
     868
     869        The algorithm used here comes from [RaWi2011]_.
     870
     871        (1) Use Cauchy's multivariate integral formula to write `F_{n\alpha}` as
     872        an integral around a polycirle centered at the origin of the
     873        differential form `\frac{G(x) dx[0] \wedge\cdots\wedge
     874        dx[d-1]}{H(x)x^\alpha}`.
     875
     876        (2) Decompose `G/H` into a sum of partial fractions `P[0] +\cdots+ P[r]`
     877        so that each term of the sum has at most `d` irreducible factors of `H`
     878        in the denominator.
     879
     880        (3) For each differential form `P[j] dx[0] \wedge\cdots\wedge dx[d-1]`,
     881        find an equivalent form `\omega[j]` in de Rham cohomology with no
     882        repeated irreducible factors of `H` in its denominator.
     883
     884        (4) Compute an asymptotic expansion for each integral `\omega[j]`.
     885
     886        (5) Add the expansions.
     887
     888        REFERENCES:
     889
     890        .. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
     891        coefficients of multivariate generating functions: improvements
     892        for smooth points", Electronic Journal of Combinatorics, Vol. 15
     893        (2008), R89.
     894
     895        .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
     896        coefficients of multivariate generating functions: improvements
     897        for smooth points", To appear.
     898
     899        AUTHORS:
     900
     901        - Alex Raichev (2008-10-01, 2010-09-28, 2011-04-27)
     902        """
     903        # The variable for asymptotic expansions.
     904        if not asy_var:
     905            asy_var= var('n')
     906
     907        # Create symbolic (non-ring) variables.
     908        R= self._R
     909        d= self._d
     910        X= list(self._variables)
     911
     912        # Do steps (1)--(3)
     913        new_integrands= self.cohomologous_integrand(alpha,asy_var)
     914
     915        # Coerce everything into the Symbolic Ring, as polynomial ring
     916        # calculations are no longer needed.
     917        # Calculate asymptotics.
     918        cc={}
     919        for k in c.keys():
     920            cc[SR(k)] = SR(c[k])
     921        c= cc
     922        for i in range(len(alpha)):
     923            alpha[i] = SR(alpha[i])
     924        subexp_parts= []
     925        for chunk in new_integrands:
     926            # Convert chunk into Symbolic Ring
     927            GG= SR(chunk[Integer(0) ])
     928            HH= [SR(f) for (f,e) in chunk[Integer(1) :]]
     929            asy= self._asymptotics_main(GG,HH,X,c,asy_var,alpha,N,numerical)
     930            subexp_parts.append(asy[Integer(2) ])
     931        exp_scale= asy[Integer(1) ]  # same for all chunk in new_integrands
     932
     933        # Do step (5).
     934        subexp_part= add(subexp_parts)
     935        return exp_scale**asy_var *subexp_part, exp_scale, subexp_part
     936#--------------------------------------------------------------------------------
     937    def _asymptotics_main(self,G,H,X,c,n,alpha,N,numerical):
     938        r"""
     939        This function is for internal use by asymptotics().
     940        It finds a variable in `X` to use to calculate asymptotics and decides
     941        whether to call _asymptotics_main_smooth() or
     942        _asymptotics_main_multiple().
     943
     944        Does not use `self`.
     945
     946        INPUT:
     947
     948        - ``G`` - A symbolic expression.
     949        - ``H`` - A list of symbolic expressions.
     950        - ``X`` - The list of variables occurring in `G` and `H`.
     951            Call its length `d`.
     952        - ``c`` - A dictionary with `X` as keys and numbers as values.
     953        - ``n`` - The variable of the asymptotic expansion.
     954        - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic
     955            variables if `c` is a smooth point.
     956        - ``N`` - A positive integer.
     957        - ``numerical`` - Natural number.
     958            If k=numerical > 0, then a numerical approximation of the coefficients
     959            of `F_{n\alpha}` with k digits of precision will be returned.
     960            Otherwise exact values will be returned.
     961
     962        OUTPUT:
     963
     964        The same as the function asymptotics().
     965
     966        AUTHORS:
     967
     968        - Alex Raichev (2008-10-01, 2010-09-28)
     969        """
     970        d= len(X)
     971        r= len(H)  # We know 1 <= r <= d.
     972
     973        # Find a good variable x in X to do asymptotics calculations with, that is,
     974        # a variable x in X such that for all h in H, diff(h,x).subs(c) != 0.
     975        # A good variable is guaranteed to exist since we are dealing with
     976        # convenient smooth or multiple points.
     977        # Search for good x in X from back to front (to be consistent with
     978        # [RaWi2008a]_ [RaWi2011]_ which uses X[d-1] as a good variable).
     979        # Put each not good x found at the beginning of X and reshuffle alpha
     980        # in the same way.
     981        x= X[d-Integer(1) ]
     982        beta= copy(alpha)
     983        while Integer(0)  in [diff(h,x).subs(c) for h in H]:
     984            X.pop()
     985            X.insert(Integer(0) ,x)
     986            x= X[d-Integer(1) ]
     987            a= beta.pop()
     988            beta.insert(Integer(0) ,a)
     989        if d==r:
     990            # This is the case of a 'simple' multiple point.
     991            A= G.subs(c) / jacobian(H,X).subs(c).determinant().abs()
     992            return A,Integer(1) ,A
     993        elif r==Integer(1) :
     994            # So 1 = r < d, and we have a smooth point.
     995            return self._asymptotics_main_smooth(G,H[Integer(0) ],X,c,n,beta,N,numerical)
     996        else:
     997            # So 1 < r < d, and we have a non-smooth multiple point.
     998            return self._asymptotics_main_multiple(G,H,X,c,n,beta,N,numerical)
     999#--------------------------------------------------------------------------------
     1000    def _asymptotics_main_multiple(self,G,H,X,c,n,alpha,N,numerical):
     1001        r"""
     1002        This function is for internal use by _asymptotics_main().
     1003        It calculates asymptotics in case they depend upon multiple points.
     1004
     1005        Does not use `self`.
     1006
     1007        INPUT:
     1008
     1009        - ``G`` - A symbolic expression.
     1010        - ``H`` - A list of symbolic expressions.
     1011        - ``X`` - The list of variables occurring in `G` and `H`.
     1012            Call its length `d`.
     1013        - ``c`` - A dictionary with `X` as keys and numbers as values.
     1014        - ``n`` - The variable of the asymptotic expansion.
     1015        - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic
     1016            variables if `c` is a smooth point.
     1017        - ``N`` - A positive integer.
     1018        - ``numerical`` - Natural number.
     1019            If k=numerical > 0, then a numerical approximation of the coefficients
     1020            of `F_{n\alpha}` with k digits of precision will be returned.
     1021            Otherwise exact values will be returned.
     1022
     1023        OUTPUT:
     1024
     1025        The same as the function asymptotics().
     1026
     1027        NOTES:
     1028
     1029        The formula used for computing the asymptotic expansion is
     1030        Theorem 3.4 of [RaWi2011]_.
     1031
     1032        Currently this function cannot handle `c` with symbolic variable keys,
     1033        because _crit_cone_combo() crashes.
     1034
     1035        REFERENCES:
     1036
     1037        .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
     1038        coefficients of multivariate generating functions: improvements
     1039        for smooth points", To appear.
     1040
     1041        AUTHORS:
     1042
     1043        - Alex Raichev (2008-10-01, 2010-09-28)
     1044        """
     1045        I= sqrt(-Integer(1) )
     1046        d= len(X)   # > r > 1
     1047        r= len(H)   # > 1
     1048        C= copy(c)
     1049
     1050        S= [var(self._new_var_name('s',X) + str(j)) for j in range(r-Integer(1) )]
     1051        T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-Integer(1) )]
     1052        Sstar= {}
     1053        temp= self._crit_cone_combo(H,X,c,alpha)
     1054        for j in range(r-Integer(1) ):
     1055            Sstar[S[j]]= temp[j]
     1056        thetastar= dict([(t,Integer(0) ) for t in T])
     1057        thetastar.update(Sstar)
     1058
     1059        # Setup.
     1060        print "Initializing auxiliary functions..."
     1061        Hmul= mul(H)
     1062        h= [function('h'+str(j),*tuple(X[:d-Integer(1) ])) for j in range(r)]    # Implicit functions
     1063        U = function('U',*tuple(X))
     1064        # All other functions are defined in terms of h, U, and explicit functions.
     1065        Hcheck=  mul([ X[d-Integer(1) ] -Integer(1) /h[j] for j in range(r)] )
     1066        Gcheck= -G/U *mul( [-h[j]/X[d-Integer(1) ] for j in range(r)] )
     1067        A= [(-Integer(1) )**(r-Integer(1) ) *X[d-Integer(1) ]**(-r+j)*diff(Gcheck.subs({X[d-Integer(1) ]:Integer(1) /X[d-Integer(1) ]}),X[d-Integer(1) ],j) for j in range(r)]
     1068        e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-Integer(1) )])
     1069        ht= [hh.subs(e) for hh in h]
     1070        Ht= [H[j].subs(e).subs({X[d-Integer(1) ]:Integer(1) /ht[j]}) for j in range(r)]
     1071        hsumt= add([S[j]*ht[j] for j in range(r-Integer(1) )]) +(Integer(1) -add(S))*ht[r-Integer(1) ]
     1072        At= [AA.subs(e).subs({X[d-Integer(1) ]:hsumt}) for AA in A]
     1073        Phit = -log(C[X[d-Integer(1) ]] *hsumt)+ I*add([alpha[i]/alpha[d-Integer(1) ] *T[i] for i in range(d-Integer(1) )])
     1074        # atC Stores h and U and all their derivatives evaluated at C.
     1075        atC = copy(C)
     1076        atC.update(dict( [(hh.subs(C),Integer(1) /C[X[d-Integer(1) ]]) for hh in h ]))
     1077
     1078        # Compute the derivatives of h up to order 2*N and evaluate at C.
     1079        hderivs1= {}   # First derivatives of h.
     1080        for (i,j) in mrange([d-Integer(1) ,r]):
     1081            s= solve( diff(H[j].subs({X[d-Integer(1) ]:Integer(1) /h[j]}),X[i]), diff(h[j],X[i]) )[Integer(0) ].rhs()\
     1082            .simplify()
     1083            hderivs1.update({diff(h[j],X[i]):s})
     1084            atC.update({diff(h[j],X[i]).subs(C):s.subs(C).subs(atC)})
     1085        hderivs = self._diff_all(h,X[Integer(0) :d-Integer(1) ],Integer(2) *N,sub=hderivs1,rekey=h)
     1086        for k in hderivs.keys():
     1087            atC.update({k.subs(C):hderivs[k].subs(atC)})
     1088
     1089        # Compute the derivatives of U up to order 2*N-2+ min{r,N}-1 and evaluate at C.
     1090        # To do this, differentiate H = U*Hcheck over and over, evaluate at C,
     1091        # and solve for the derivatives of U at C.
     1092        # Need the derivatives of H with short keys to pass on to diff_prod later.
     1093        m= min(r,N)
     1094        end= [X[d-Integer(1) ] for j in range(r)]
     1095        Hmulderivs= self._diff_all(Hmul,X,Integer(2) *N-Integer(2) +r,ending=end,sub_final=C)
     1096        print "Calculating derivatives of auxiliary functions..."
     1097        atC.update({U.subs(C):diff(Hmul,X[d-Integer(1) ],r).subs(C)/factorial(r)})
     1098        Uderivs={}
     1099        p= Hmul.polynomial(CC).degree()-r
     1100        if p == Integer(0) :
     1101            # Then, using a proposition at the end of [RaWi2011], we can
     1102            # conclude that all higher derivatives of U are zero.
     1103            for l in (ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )):
     1104                for s in UnorderedTuples(X,l):
     1105                    Uderivs[diff(U,s).subs(C)] = Integer(0)
     1106        elif p > Integer(0)  and p < Integer(2) *N-Integer(2) +m-Integer(1) :
     1107            all_zero= True
     1108            Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,p)),end,Uderivs,atC)
     1109            # Check for a nonzero U derivative.
     1110            if Uderivs.values() != [Integer(0)  for i in range(len(Uderivs))]:
     1111                all_zero= False
     1112            if all_zero:
     1113                # Then, using a proposition at the end of [RaWi2011], we can
     1114                # conclude that all higher derivatives of U are zero.
     1115                for l in (ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )):
     1116                    for s in UnorderedTuples(X,l):
     1117                        Uderivs.update({diff(U,s).subs(C):Integer(0) })
     1118            else:
     1119                # Have to compute the rest of the derivatives.
     1120                Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )),end,Uderivs,atC)
     1121        else:
     1122            Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )),end,Uderivs,atC)
     1123        atC.update(Uderivs)
     1124        Phit1= jacobian(Phit,T+S).subs(hderivs1)
     1125        a= jacobian(Phit1,T+S).subs(hderivs1).subs(thetastar).subs(atC)
     1126        a_inv= a.inverse()
     1127        Phitu= Phit -(Integer(1) /Integer(2) ) *matrix([T+S]) *a *matrix([T+S]).transpose()
     1128        Phitu= Phitu[Integer(0) ][Integer(0) ]
     1129
     1130        # Compute all partial derivatives of At and Phitu up to orders 2*N-2
     1131        # and 2*N, respectively.  Take advantage of the fact that At and Phitu
     1132        # are sufficiently differentiable functions so that mixed partials
     1133        # are equal.  Thus only need to compute representative partials.
     1134        # Choose nondecreasing sequences as representative differentiation-
     1135        # order sequences.
     1136        # To speed up later computations, create symbolic functions AA and BB
     1137        # to stand in for the expressions At and Phitu respectively.
     1138        print "Calculating derivatives of more auxiliary functions..."
     1139        AA= [function('A'+str(j),*tuple(T+S)) for j in range(r)]
     1140        At_derivs= self._diff_all(At,T+S,Integer(2) *N-Integer(2) ,sub=hderivs1,sub_final=[thetastar,atC],rekey=AA)
     1141        BB= function('BB',*tuple(T+S))
     1142        Phitu_derivs= self._diff_all(Phitu,T+S,Integer(2) *N,sub=hderivs1,sub_final=[thetastar,atC],rekey=BB,zero_order=Integer(3) )
     1143        AABB_derivs= At_derivs
     1144        AABB_derivs.update(Phitu_derivs)
     1145        for j in range(r):
     1146            AABB_derivs[AA[j]] = At[j].subs(thetastar).subs(atC)
     1147        AABB_derivs[BB] = Phitu.subs(thetastar).subs(atC)
     1148        print "Calculating second-order differential operator actions..."
     1149        DD= self._diff_op(AA,BB,AABB_derivs,T+S,a_inv,r,N)
     1150
     1151        L={}
     1152        for (j,k) in  CartesianProduct((ellipsis_range(Integer(0) ,Ellipsis,min(r-Integer(1) ,N-Integer(1) ))), (ellipsis_range(max(Integer(0) ,N-Integer(1) -r),Ellipsis,N-Integer(1) ))):
     1153            if j+k <= N-Integer(1) :
     1154                L[(j,k)] = add([ \
     1155                DD[(j,k,l)] /( (-Integer(1) )**k *Integer(2) **(k+l) *factorial(l) *factorial(k+l) ) \
     1156                for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k))] )
     1157        # The next line's QQ coercion is a workaround for the Sage 4.6 bug reported
     1158        # on http://trac.sagemath.org/sage_trac/ticket/8659.
     1159        # Once the bug is fixed, the QQ can be removed.
     1160        det= QQ(a.determinant())**(-Integer(1) /Integer(2) ) *(Integer(2) *pi)**((r-d)/Integer(2) )
     1161        chunk= det *add([ (alpha[d-Integer(1) ]*n)**((r-d)/Integer(2) -q) *add([ \
     1162        L[(j,k)] *binomial(r-Integer(1) ,j) *stirling_number1(r-j,r+k-q) *(-Integer(1) )**(q-j-k) \
     1163        for (j,k) in CartesianProduct((ellipsis_range(Integer(0) ,Ellipsis,min(r-Integer(1) ,q))), (ellipsis_range(max(Integer(0) ,q-r),Ellipsis,q))) if j+k <= q ]) \
     1164        for q in range(N)])
     1165
     1166        chunk= chunk.subs(C).simplify()
     1167        coeffs= chunk.coefficients(n)
     1168        coeffs.reverse()
     1169        coeffs= coeffs[:N]
     1170        if numerical: # If a numerical approximation is desired...
     1171            subexp_part = add( [co[Integer(0) ].subs(c).n(digits=numerical)*n**co[Integer(1) ] \
     1172            for co in coeffs] )
     1173            exp_scale= (Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )) \
     1174            .n(digits=numerical)
     1175        else:
     1176            subexp_part = add( [co[Integer(0) ].subs(c)*n**co[Integer(1) ] for co in coeffs] )
     1177            exp_scale= Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )
     1178        return exp_scale**n*subexp_part, exp_scale, subexp_part
     1179#--------------------------------------------------------------------------------
     1180    def _asymptotics_main_smooth(self,G,H,X,c,n,alpha,N,numerical):
     1181        r"""
     1182        This function is for internal use by _asymptotics_main().
     1183        It calculates asymptotics in case they depend upon smooth points.
     1184
     1185        Does not use `self`.
     1186
     1187        INPUT:
     1188
     1189        - ``G`` - A symbolic expression.
     1190        - ``H`` - A list of symbolic expressions.
     1191        - ``X`` - The list of variables occurring in `G` and `H`.
     1192            Call its length `d`.
     1193        - ``c`` - A dictionary with `X` as keys and numbers as values.
     1194        - ``n`` - The variable of the asymptotic expansion.
     1195        - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic
     1196            variables if `c` is a smooth point.
     1197        - ``N`` - A positive integer.
     1198        - ``numerical`` - Natural number.
     1199            If k=numerical > 0, then a numerical approximation of the coefficients
     1200            of `F_{n\alpha}` with k digits of precision will be returned.
     1201            Otherwise exact values will be returned.
     1202
     1203        OUTPUT:
     1204
     1205        The same as the function asymptotics().
     1206
     1207        NOTES:
     1208
     1209        The formulas used for computing the asymptotic expansions are
     1210        Theorems 3.2 and 3.3 [RaWi2008a]_ with `p=1`.
     1211        Theorem 3.2 is a specialization of Theorem 3.4 of [RaWi2011]_
     1212        with `r=1`.
     1213
     1214        REFERENCES:
     1215
     1216        .. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
     1217        coefficients of multivariate generating functions: improvements
     1218        for smooth points", Electronic Journal of Combinatorics, Vol. 15
     1219        (2008), R89.
     1220
     1221        .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
     1222        coefficients of multivariate generating functions: improvements
     1223        for smooth points", To appear.
     1224
     1225        AUTHORS:
     1226
     1227        - Alex Raichev (2008-10-01, 2010-09-28)
     1228        """
     1229        I= sqrt(-Integer(1) )
     1230        d= len(X) # > 1
     1231
     1232        # If c is a tuple of rationals, then compute with it directly.
     1233        # Otherwise, compute symbolically and plug in c at the end.
     1234        if vector(c.values()) in QQ**d:
     1235            C= c
     1236        else:
     1237            Cs= [var('cs' +str(j)) for j in range(d)]
     1238            C= dict( [(X[j],Cs[j]) for j in range(d)] )
     1239            c= dict( [(Cs[j],c[X[j]]) for j in range(d)] )
     1240
     1241        # Setup.
     1242        print "Initializing auxiliary functions..."
     1243        h= function('h',*tuple(X[:d-Integer(1) ]))    # Implicit functions
     1244        U = function('U',*tuple(X))         #
     1245        # All other functions are defined in terms of h, U, and explicit functions.
     1246        Gcheck = -G/U *(h/X[d-Integer(1) ])
     1247        A= Gcheck.subs({X[d-Integer(1) ]:Integer(1) /h})/h
     1248        T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-Integer(1) )]
     1249        e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-Integer(1) )])
     1250        ht= h.subs(e)
     1251        Ht= H.subs(e).subs({X[d-Integer(1) ]:Integer(1) /ht})
     1252        At= A.subs(e)
     1253        Phit = -log(C[X[d-Integer(1) ]] *ht)\
     1254        + I* add([alpha[i]/alpha[d-Integer(1) ] *T[i] for i in range(d-Integer(1) )])
     1255        Tstar= dict([(t,Integer(0) ) for t in T])
     1256        atC = copy(C)
     1257        atC.update({h.subs(C):Integer(1) /C[X[d-Integer(1) ]]}) # Stores h and U and all their derivatives
     1258                                            # evaluated at C.
     1259
     1260        # Compute the derivatives of h up to order 2*N and evaluate at C and store
     1261        # in atC.  Keep a copy of unevaluated h derivatives for use in the case
     1262        # d=2 and v > 2 below.
     1263        hderivs1= {}   # First derivatives of h.
     1264        for i in range(d-Integer(1) ):
     1265            s= solve( diff(H.subs({X[d-Integer(1) ]:Integer(1) /h}),X[i]), diff(h,X[i]) )[Integer(0) ].rhs()\
     1266            .simplify()
     1267            hderivs1.update({diff(h,X[i]):s})
     1268            atC.update({diff(h,X[i]).subs(C):s.subs(C).subs(atC)})
     1269        hderivs = self._diff_all(h,X[Integer(0) :d-Integer(1) ],Integer(2) *N,sub=hderivs1,rekey=h)
     1270        for k in hderivs.keys():
     1271            atC.update({k.subs(C):hderivs[k].subs(atC)})
     1272
     1273        # Compute the derivatives of U up to order 2*N and evaluate at C.
     1274        # To do this, differentiate H = U*Hcheck over and over, evaluate at C,
     1275        # and solve for the derivatives of U at C.
     1276        # Need the derivatives of H with short keys to pass on to diff_prod later.
     1277        Hderivs= self._diff_all(H,X,Integer(2) *N,ending=[X[d-Integer(1) ]],sub_final=C)
     1278        print "Calculating derivatives of auxiallary functions..."
     1279        # For convenience in checking if all the nontrivial derivatives of U at c
     1280        # are zero a few line below, store the value of U(c) in atC instead of in
     1281        # Uderivs.
     1282        Uderivs={}
     1283        atC.update({U.subs(C):diff(H,X[d-Integer(1) ]).subs(C)})
     1284        end= [X[d-Integer(1) ]]
     1285        Hcheck= X[d-Integer(1) ] - Integer(1) /h
     1286        p= H.polynomial(CC).degree()-Integer(1)
     1287        if p == Integer(0) :
     1288            # Then, using a proposition at the end of [RaWi2011], we can
     1289            # conclude that all higher derivatives of U are zero.
     1290            for l in (ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N)):
     1291                for s in UnorderedTuples(X,l):
     1292                    Uderivs[diff(U,s).subs(C)] = Integer(0)
     1293        elif p > Integer(0)  and p < Integer(2) *N:
     1294            all_zero= True
     1295            Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,p)),end,Uderivs,atC)
     1296            # Check for a nonzero U derivative.
     1297            if Uderivs.values() != [Integer(0)  for i in range(len(Uderivs))]:
     1298                all_zero= False
     1299            if all_zero:
     1300                # Then, using a proposition at the end of [RaWi2011], we can
     1301                # conclude that all higher derivatives of U are zero.
     1302                for l in (ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N)):
     1303                    for s in UnorderedTuples(X,l):
     1304                        Uderivs.update({diff(U,s).subs(C):Integer(0) })
     1305            else:
     1306                # Have to compute the rest of the derivatives.
     1307                Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N)),end,Uderivs,atC)
     1308        else:
     1309            Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N)),end,Uderivs,atC)
     1310        atC.update(Uderivs)
     1311
     1312        # In general, this algorithm is not designed to handle the case of a
     1313        # singular Phit''(Tstar).  However, when d=2 the algorithm can cope.
     1314        if d==Integer(2) :
     1315            # Compute v, the order of vanishing at Tstar of Phit. It is at least 2.
     1316            v=Integer(2)
     1317            Phitderiv= diff(Phit,T[Integer(0) ],Integer(2) )
     1318            splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify()
     1319            while splat==Integer(0) :
     1320                v= v+Integer(1)
     1321                if v > Integer(2) *N:  # Then need to compute more derivatives of h for atC.
     1322                    hderivs.update({diff(h,X[Integer(0) ],v) \
     1323                    :diff(hderivs[diff(h,X[Integer(0) ],v-Integer(1) )],X[Integer(0) ]).subs(hderivs1)})
     1324                    atC.update({diff(h,X[Integer(0) ],v).subs(C) \
     1325                    :hderivs[diff(h,X[Integer(0) ],v)].subs(atC)})
     1326                Phitderiv= diff(Phitderiv,T[Integer(0) ])
     1327                splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify()
     1328        if d==Integer(2)  and v>Integer(2) :
     1329            t= T[Integer(0) ]  # Simplify variable names.
     1330            a= splat/factorial(v)
     1331            Phitu= Phit -a*t**v
     1332
     1333            # Compute all partial derivatives of At and Phitu up to orders 2*(N-1)
     1334            # and 2*(N-1)+v, respectively, in case v is even.
     1335            # Otherwise, compute up to orders N-1 and N-1+v, respectively.
     1336            # To speed up later computations, create symbolic functions AA and BB
     1337            # to stand in for the expressions At and Phitu, respectively.
     1338            print "Calculating derivatives of more auxiliary functions..."
     1339            AA= function('AA',t)
     1340            BB= function('BB',t)
     1341            if v.mod(Integer(2) )==Integer(0) :
     1342                At_derivs= self._diff_all(At,T,Integer(2) *N-Integer(2) , \
     1343                sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
     1344                Phitu_derivs= self._diff_all(Phitu,T,Integer(2) *N-Integer(2) +v, \
     1345                sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+Integer(1) ,rekey=BB)
     1346            else:
     1347                At_derivs= self._diff_all(At,T,N-Integer(1) ,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
     1348                Phitu_derivs= self._diff_all(Phitu,T,N-Integer(1) +v,sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+Integer(1) ,rekey=BB)
     1349            AABB_derivs= At_derivs
     1350            AABB_derivs.update(Phitu_derivs)
     1351            AABB_derivs[AA] = At.subs(Tstar).subs(atC)
     1352            AABB_derivs[BB] = Phitu.subs(Tstar).subs(atC)
     1353            print "Calculating actions of the second order differential operator..."
     1354            DD= self._diff_op_simple(AA,BB,AABB_derivs,t,v,a,N)
     1355            # Plug above into asymptotic formula.
     1356            L = []
     1357            if v.mod(Integer(2) ) == Integer(0) :
     1358                for k in range(N):
     1359                    L.append( add([ \
     1360                    (-Integer(1) )**l *gamma((Integer(2) *k+v*l+Integer(1) )/v) \
     1361                    / (factorial(l) *factorial(Integer(2) *k+v*l)) \
     1362                    * DD[(k,l)] for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)) ]) )
     1363                chunk= a**(-Integer(1) /v) /(pi*v) *add([ alpha[d-Integer(1) ]**(-(Integer(2) *k+Integer(1) )/v) \
     1364                * L[k] *n**(-(Integer(2) *k+Integer(1) )/v) for k in range(N) ])
     1365            else:
     1366                zeta= exp(I*pi/(Integer(2) *v))
     1367                for k in range(N):
     1368                    L.append( add([ \
     1369                    (-Integer(1) )**l *gamma((k+v*l+Integer(1) )/v) \
     1370                    / (factorial(l) *factorial(k+v*l)) \
     1371                    * (zeta**(k+v*l+Integer(1) ) +(-Integer(1) )**(k+v*l)*zeta**(-(k+v*l+Integer(1) ))) \
     1372                    * DD[(k,l)] for l in (ellipsis_range(Integer(0) ,Ellipsis,k)) ]) )
     1373                chunk= abs(a)**(-Integer(1) /v) /(Integer(2) *pi*v) *add([ alpha[d-Integer(1) ]**(-(k+Integer(1) )/v) \
     1374                * L[k] *n**(-(k+Integer(1) )/v) for k in range(N) ])
     1375        # Asymptotics for d>=2 case.  A singular Phit''(Tstar) will cause a crash
     1376        # in this case.
     1377        else:
     1378            Phit1= jacobian(Phit,T).subs(hderivs1)
     1379            a= jacobian(Phit1,T).subs(hderivs1).subs(Tstar).subs(atC)
     1380            a_inv= a.inverse()
     1381            Phitu= Phit -(Integer(1) /Integer(2) ) *matrix([T]) *a *matrix([T]).transpose()
     1382            Phitu= Phitu[Integer(0) ][Integer(0) ]
     1383            # Compute all partial derivatives of At and Phitu up to orders 2*N-2
     1384            # and 2*N, respectively.  Take advantage of the fact that At and Phitu
     1385            # are sufficiently differentiable functions so that mixed partials
     1386            # are equal.  Thus only need to compute representative partials.
     1387            # Choose nondecreasing sequences as representative differentiation-
     1388            # order sequences.
     1389            # To speed up later computations, create symbolic functions AA and BB
     1390            # to stand in for the expressions At and Phitu respectively.
     1391            print "Calculating derivatives of more auxiliary functions..."
     1392            AA= function('AA',*tuple(T))
     1393            At_derivs= self._diff_all(At,T,Integer(2) *N-Integer(2) ,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
     1394            BB= function('BB',*tuple(T))
     1395            Phitu_derivs= self._diff_all(Phitu,T,Integer(2) *N,sub=hderivs1,sub_final=[Tstar,atC],rekey=BB,zero_order=Integer(3) )
     1396            AABB_derivs= At_derivs
     1397            AABB_derivs.update(Phitu_derivs)
     1398            AABB_derivs[AA] = At.subs(Tstar).subs(atC)
     1399            AABB_derivs[BB] = Phitu.subs(Tstar).subs(atC)
     1400            print "Calculating actions of the second order differential operator..."
     1401            DD= self._diff_op(AA,BB,AABB_derivs,T,a_inv,Integer(1) ,N)
     1402
     1403            # Plug above into asymptotic formula.
     1404            L=[]
     1405            for k in range(N):
     1406                L.append( add([ \
     1407                  DD[(Integer(0) ,k,l)] / ( (-Integer(1) )**k *Integer(2) **(l+k) *factorial(l) *factorial(l+k) ) \
     1408                  for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k))]) )
     1409            chunk= add([ (Integer(2) *pi)**((Integer(1) -d)/Integer(2) ) *a.determinant()**(-Integer(1) /Integer(2) ) \
     1410                  *alpha[d-Integer(1) ]**((Integer(1) -d)/Integer(2)  -k) *L[k] \
     1411                  *n**((Integer(1) -d)/Integer(2) -k) for k in range(N) ])
     1412
     1413        chunk= chunk.subs(c).simplify()
     1414        coeffs= chunk.coefficients(n)
     1415        coeffs.reverse()
     1416        coeffs= coeffs[:N]
     1417        if numerical: # If a numerical approximation is desired...
     1418            subexp_part = add( [co[Integer(0) ].subs(c).n(digits=numerical)*n**co[Integer(1) ] for co in coeffs] )
     1419            exp_scale=(Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )) \
     1420            .n(digits=numerical)
     1421        else:
     1422            subexp_part = add( [co[Integer(0) ].subs(c)*n**co[Integer(1) ] for co in coeffs] )
     1423            exp_scale= Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )
     1424        return exp_scale**n*subexp_part, exp_scale, subexp_part
     1425#-----------------------------------------------------------------------------
     1426    def _crit_cone_combo(self,fs,X,c,alpha):
     1427        r"""
     1428        This function returns an auxiliary point associated to the multiple
     1429        point `c` of the factors `fs`.
     1430        It is for internal use by _asymptotics_main_multiple().
     1431
     1432        INPUT:
     1433
     1434        - ``fs`` - A list of expressions in the variables of `X`.
     1435        - ``X`` - A list of variables.
     1436        - ``c`` - A dictionary with keys `X` and values in some field.
     1437        - ``alpha`` - A list of rationals.
     1438
     1439        OUTPUT:
     1440
     1441        A solution of the matrix equation `y Gamma = a` for `y`,
     1442        where `Gamma` is the matrix whose `j`th row is
     1443        _direction(_log_grad(fj,X,c)) where `fj`
     1444        is the `j`th item in `fs` and where `a` is _direction(alpha).
     1445
     1446        EXAMPLES::
     1447
     1448            sage: R.<x,y,z>= PolynomialRing(QQ)
     1449            sage: G,H= 1,1
     1450            sage: F= QuasiRationalExpression(G,H)
     1451            sage: fs= [x + 2*y + z - 4, 2*x + y + z - 4]
     1452            sage: c= {x:1,y:1,z:1}
     1453            sage: alpha= [2,1,1]
     1454            sage: F._crit_cone_combo(fs,R.gens(),c,alpha)
     1455            [0, 1]
     1456
     1457        NOTES:
     1458
     1459        Use this function only when `Gamma` is well-defined and
     1460        there is a unique solution to the matrix equation `y Gamma = a`.
     1461        Fails otherwise.
     1462
     1463        AUTHORS:
     1464
     1465        - Alex Raichev (2008-10-01, 2008-11-25, 2009-03-04, 2010-09-08,
     1466            2010-12-02)
     1467        """
     1468        # Assuming here that each _log_grad(f) has nonzero final component.
     1469        # Then 'direction' will not throw a division by zero error.
     1470        d= len(X)
     1471        r= len(fs)
     1472        Gamma= matrix([self._direction(self._log_grad(f,X,c)) for f in fs])
     1473        # solve_left() fails when working in SR :-(. So use solve() instead.
     1474        #s= Gamma.solve_left(vector(alpha)/alpha[d-1])
     1475        V= [var('sss'+str(i)) for i in range(r)]
     1476        M= matrix(V)*Gamma
     1477        eqns= [M[Integer(0) ][i]== alpha[i]/alpha[d-Integer(1) ] for i in range(d)]
     1478        s= solve(eqns,V,solution_dict=True)[Integer(0) ]  # Assuming a unique solution.
     1479        return [s[v] for v in V]
     1480# B ==========================================================================
     1481# C ==========================================================================
     1482    def cohomologous_integrand(self,alpha,asy_var=None):
     1483        r"""
     1484        This function takes a multivariate Cauchy type integral
     1485        `\int F / x^{asy_var\alpha+1} dx`, where `F=self`, and breaks it up
     1486        into a list of nicer Cauchy type integrals for the purposes of
     1487        computing asymptotics of the original integral as `asy_var\to\infty`.
     1488        The sum of the nicer integrals is de Rham cohomologous to the original
     1489        integral.
     1490        It assumes that algebraic varieties corresponding to the irreducible
     1491        factors of `self._H` intersect transversely (see notes below).
     1492
     1493        INPUT:
     1494
     1495        - ``alpha`` - A list of positive integers or symbolic variables.
     1496        - ``asy_var`` - A symbolic variable (default: None).
     1497            Eventually set to `var('n')` if None is given.
     1498
     1499        OUTPUT:
     1500
     1501        A list of the form `[chunk_1,\ldots,chunk_r]`, where each
     1502        `chunk_j` has the form `[P,[B_1,1],\ldots,[B_l,1]]`.
     1503        Here `l \le d`, `P` is a symbolic expression in the indeterminates of
     1504        `R` and `asy_var`, `\{B_1,\ldots,B_l\} \subseteq \{Q_1,\ldots,Q_m\}`,
     1505        and `[P,[B_1,1],\ldots,[B_l,1]]` represents the integral
     1506        `\int P/(B_1 \cdots B_l) dx`.
     1507
     1508        EXAMPLES::
     1509
     1510        sage: R.<x,y>= PolynomialRing(QQ)
     1511        sage: G= 9*exp(x+y)
     1512        sage: H= (3-2*x-y)*(3-x-2*y)
     1513        sage: F= QuasiRationalExpression(G,H)
     1514        sage: alpha= [4,3]
     1515        sage: F.cohomologous_integrand(alpha)
     1516        [[9*e^(x + y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]]
     1517
     1518        sage: R.<x,y>= PolynomialRing(QQ)
     1519        sage: G= 9*exp(x+y)
     1520        sage: H= (3-2*x-y)^2*(3-x-2*y)
     1521        sage: F= QuasiRationalExpression(G,H)
     1522        sage: alpha= [4,3]
     1523        sage: F.cohomologous_integrand(alpha)
     1524        [[-3*(3*x*e^x - 8*y*e^x)*n*e^y/(x*y) - 3*((x - 2)*y*e^x + x*e^x)*e^y/(x*y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]]
     1525
     1526        sage: R.<x,y,z>= PolynomialRing(QQ)
     1527        sage: G= 16
     1528        sage: H= (4-2*x-y-z)^2*(4-x-2*y-z)
     1529        sage: F= QuasiRationalExpression(G,H)
     1530        sage: alpha= [3,3,2]
     1531        sage: F.cohomologous_integrand(alpha)
     1532        [[16*(4*y - 3*z)*n/(y*z) + 16*(2*y - z)/(y*z), [x + 2*y + z - 4, 1], [2*x + y + z - 4, 1]]]
     1533
     1534        NOTES:
     1535
     1536        The varieties  corresponding to `Q_1,\ldots,Q_m`
     1537        __intersect transversely__ iff for each point `c` of their intersection
     1538        and for all `k \le l`, the Jacobian matrix of any `k` polynomials from
     1539        `\{Q_1,\ldots,Q_m\}` has rank equal to `\min\{k,d\}` when evaluated at
     1540        `c`.
     1541
     1542        ALGORITHM:
     1543
     1544        Let `asy_var= n` and consider the integral around a polycirle centered
     1545        at the origin of the `d`-variate differential form
     1546        `\frac{G(x) dx_1 \wedge\cdots\wedge dx_d}{H(x) x^{n\alpha+1}}`, where
     1547        `G=self._G` and `H=self._H`.
     1548
     1549        (1) Decompose `G/H` into a sum of partial fractions `P_1 +\cdots+ P_r`
     1550        so that each term of the sum has at most `d` irreducible factors of `H`
     1551        in the denominator.
     1552
     1553        (2) For each differential form `P_j dx_1 \wedge\cdots\wedge dx_d`,
     1554        find an equivalent form `\omega_j` in de Rham cohomology with no
     1555        repeated irreducible factors of `H` in its denominator.
     1556
     1557        AUTHOR:
     1558
     1559        - Alex Raichev (2010-09-22)
     1560        """
     1561        if not asy_var:
     1562            asy_var = var('n')
     1563
     1564        # Create symbolic (non-ring) variables.
     1565        G= self._G
     1566        H= self._H
     1567        R= self._R
     1568        d= self._d
     1569        X= self._variables
     1570        # Prepare input for partial_fraction_decomposition() which only works
     1571        # for functions in the field of fractions of R.
     1572        if G in R:
     1573            numer=Integer(1)
     1574            F= G/H
     1575        else:
     1576            numer= G
     1577            F= R(Integer(1) )/H
     1578        nstuff= Integer(1) /mul([X[j]**(alpha[j]*asy_var+Integer(1) ) for j in range(d)])
     1579
     1580        # Do steps (2) and (1).
     1581        integrands= []
     1582        whole,parts= partial_fraction_decomposition(F)
     1583        for f in parts:
     1584            a= format_quotient(f)
     1585            integrands.append( [a[Integer(0) ]*numer*nstuff] + a[Integer(1) :] )
     1586
     1587        # Do step (3).
     1588        ce= decompose_via_cohomology(integrands)
     1589        ce_new= []
     1590        for a in ce:
     1591            ce_new.append( [(a[Integer(0) ]/nstuff).simplify_full().collect(n)] + a[Integer(1) :] )
     1592        return ce_new
     1593#-------------------------------------------------------------------------------
     1594    def critical_cone(self,c,coordinate=None):
     1595        r"""
     1596        Returns the critical cone of a convenient multiple point `c`.
     1597
     1598        INPUT:
     1599
     1600        - ``c`` - A dictionary with keys `self.variables()` and values
     1601            in a field.
     1602        - ``coordinate`` - (optional) A natural number.
     1603
     1604        OUTPUT:
     1605
     1606        A list of vectors that generate the critical cone of `c` and
     1607        the cone itself if the values of `c` lie in QQ.
     1608
     1609        EXAMPLES::
     1610
     1611            sage: R.<x,y,z>= PolynomialRing(QQ)
     1612            sage: G= 1
     1613            sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y))
     1614            sage: F= QuasiRationalExpression(G,H)
     1615            sage: c= {z: 4/3, y: 1, x: 1/2}
     1616            sage: F.critical_cone(c)
     1617            ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N)
     1618
     1619        NOTES:
     1620
     1621        The _critical cone_ of a convenient multiple point `c` with
     1622        with `c_k \del_k H_j(c) \neq 0` for all `j=1,\ldots,r` is
     1623        the conical hull of the vectors `\gamma_j(c) =
     1624        \left(\frac{c_1 \del_1 H_j(c)}{c_k \del_k H_j(c)},\ldots,
     1625        \frac{c_d \del_d H_j(c)}{c_k \del_k H_j(c)} \right)`.
     1626        Here `H_1,\ldots,H_r` are the irreducible germs of `self._H` around `c`.
     1627        For more details, see [RaWi2011]_.
     1628
     1629        If this function's optional argument `coordinate` isn't given, then
     1630        this function searches (from `d` down to 1) for the first index `k`
     1631        such that for all `j=1,\ldots,r` we have `c_k \del_k H_j(c) \neq 0`
     1632        and sets `coordinate = k`.
     1633        Almost.
     1634        Since this is Python, all the indices actually start at 0.
     1635
     1636        REFERENCES:
     1637
     1638        .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
     1639        coefficients of multivariate generating functions: improvements
     1640        for smooth points", To appear.
     1641
     1642        AUTHORS:
     1643
     1644        - Alex Raichev (2010-08-25)
     1645        """
     1646        Hs= [SR(h[Integer(0) ]) for h in self._Hfac]  # irreducible factors of H
     1647        X= self._variables
     1648        d= self._d
     1649        # Ensure the variables of `c` lie in SR
     1650        cc= {}
     1651        for x in c.keys():
     1652            cc[SR(x)] = c[x]
     1653        lg= [self._log_grad(h,X,cc) for h in Hs]
     1654        if not coordinate:
     1655            # Search (from d-1 down to 0) for a coordinate k such that
     1656            # for all h in Hs we have cc[k] * diff(h,X[k]) !=0.
     1657            # One is guaranteed to exist in the case of a convenient multiple
     1658            # point.
     1659            for k in reversed(range(d)):
     1660                if Integer(0)  not in [v[k] for v in lg]:
     1661                    coordinate= k
     1662                    break
     1663        gamma= [self._direction(v,coordinate) for v in lg]
     1664        if [[gg in QQ for gg in g] for g in gamma] == \
     1665        [[True for gg in g] for g in gamma]:
     1666            return gamma,Cone(gamma)    # Cone() needs rational vectors
     1667        else:
     1668            return gamma
     1669# D ============================================================================
     1670    def _diff_all(self,f,V,n,ending=[],sub=None,sub_final=None,zero_order=Integer(0) ,rekey=None):
     1671        r"""
     1672        This function returns a dictionary of representative mixed partial
     1673        derivatives of `f` from order 1 up to order `n` with respect to the
     1674        variables in `V`.
     1675        The default is to key the dictionary by all nondecreasing sequences
     1676        in `V` of length 1 up to length `n`.
     1677        For internal use.
     1678
     1679        Does not use `self`.
     1680
     1681        INPUT:
     1682
     1683        - ``f`` - An individual or list of `\mathcal{C}^{n+1}` functions.
     1684        - ``V`` - A list of variables occurring in `f`.
     1685        - ``n`` - A natural number.
     1686        - ``ending`` - A list of variables in `V`.
     1687        - ``sub`` - An individual or list of dictionaries.
     1688        - ``sub_final`` - An individual or list of dictionaries.
     1689        - ``rekey`` - A callable symbolic function in `V` or list thereof.
     1690        - ``zero_order`` - A natural number.
     1691
     1692        OUTPUT:
     1693
     1694        The dictionary `{s_1:deriv_1,...,s_r:deriv_r}`.
     1695        Here `s_1,\ldots,s_r` is a listing of
     1696        all nondecreasing sequences of length 1 up to length `n` over the
     1697        alphabet `V`, where `w > v` in `X`  iff `str(w) > str(v)`, and
     1698        `deriv_j` is the derivative of `f` with respect to the derivative
     1699        sequence `s_j` and simplified with respect to the substitutions in `sub`
     1700        and evaluated at `sub_final`.
     1701        Moreover, all derivatives with respect to sequences of length less than
     1702        `zero_order` (derivatives of order less than `zero_order`) will be made
     1703        zero.
     1704
     1705        If `rekey` is nonempty, then `s_1,\ldots,s_r` will be replaced by the
     1706        symbolic derivatives of the functions in `rekey`.
     1707
     1708        If `ending` is nonempty, then every derivative sequence `s_j` will be
     1709        suffixed by `ending`.
     1710
     1711        EXAMPLES:
     1712
     1713        I'd like to print the entire dictionaries, but that doesn't yield
     1714        consistent output order for doctesting.
     1715        Order of keys changes.
     1716       
     1717        ::
     1718
     1719            sage: R.<x> = PolynomialRing(QQ)
     1720            sage: G,H = 1,1
     1721            sage: F= QuasiRationalExpression(G,H)
     1722            sage: f= function('f',x)
     1723            sage: dd= F._diff_all(f,[x],3)
     1724            sage: dd[(x,x,x)]
     1725            D[0, 0, 0](f)(x)
     1726
     1727            ::
     1728
     1729            sage: d1= {diff(f,x): 4*x^3}
     1730            sage: dd= F._diff_all(f,[x],3,sub=d1)
     1731            sage: dd[(x,x,x)]
     1732            24*x
     1733
     1734            ::
     1735
     1736            sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f)
     1737            sage: dd[diff(f,x,3)]
     1738            24*x
     1739
     1740            ::
     1741
     1742            sage: a= {x:1}
     1743            sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f,sub_final=a)
     1744            sage: dd[diff(f,x,3)]
     1745            24
     1746
     1747            ::
     1748
     1749            sage: X= var('x,y,z')
     1750            sage: f= function('f',*X)
     1751            sage: dd= F._diff_all(f,X,2,ending=[y,y,y])
     1752            sage: dd[(z,y,y,y)]
     1753            D[1, 1, 1, 2](f)(x, y, z)
     1754
     1755            ::
     1756
     1757            sage: g= function('g',*X)
     1758            sage: dd= F._diff_all([f,g],X,2)
     1759            sage: dd[(0,y,z)]
     1760            D[1, 2](f)(x, y, z)
     1761
     1762            ::
     1763
     1764            sage: dd[(1,z,z)]
     1765            D[2, 2](g)(x, y, z)
     1766
     1767            ::
     1768
     1769            sage: f= exp(x*y*z)
     1770            sage: ff= function('ff',*X)
     1771            sage: dd= F._diff_all(f,X,2,rekey=ff)
     1772            sage: dd[diff(ff,x,z)]
     1773            x*y^2*z*e^(x*y*z) + y*e^(x*y*z)
     1774
     1775        AUTHORS:
     1776
     1777        - Alex Raichev (2008-10-01, 2009-04-01, 2010-02-01)
     1778        """
     1779        singleton=False
     1780        if not isinstance(f,list):
     1781            f= [f]
     1782            singleton=True
     1783
     1784        # Build the dictionary of derivatives iteratively from a list of nondecreasing
     1785        # derivative-order sequences.
     1786        derivs= {}
     1787        r= len(f)
     1788        if ending:
     1789            seeds = [ending]
     1790            start = Integer(1)
     1791        else:
     1792            seeds = [[v] for v in V]
     1793            start = Integer(2)
     1794        if singleton:
     1795            for s in seeds:
     1796                derivs[tuple(s)] = self._subs_all(diff(f[Integer(0) ],s),sub)
     1797            for l in (ellipsis_range(start,Ellipsis,n)):
     1798                for t in UnorderedTuples(V,l):
     1799                    s= tuple(t + ending)
     1800                    derivs[s] = self._subs_all(diff(derivs[s[Integer(1) :]],s[Integer(0) ]),sub)
     1801        else:
     1802            # Make the dictionary keys of the form (j,sequence of variables),
     1803            # where j in range(r).
     1804            for s in seeds:
     1805                value= self._subs_all([diff(f[j],s) for j in range(r)],sub)
     1806                derivs.update(dict([(tuple([j]+s),value[j]) for j in range(r)]))
     1807            for l in (ellipsis_range(start,Ellipsis,n)):
     1808                for t in UnorderedTuples(V,l):
     1809                    s= tuple(t + ending)
     1810                    value= self._subs_all(\
     1811                    [diff(derivs[(j,)+s[Integer(1) :]],s[Integer(0) ]) for j in range(r)],sub)
     1812                    derivs.update(dict([((j,)+s,value[j]) for j in range(r)]))
     1813        if zero_order:
     1814            # Zero out all the derivatives of order < zero_order
     1815            if singleton:
     1816                for k in derivs.keys():
     1817                    if len(k) < zero_order:
     1818                        derivs[k]= Integer(0)
     1819            else:
     1820                # Ignore the first of element of k, which is an index.
     1821                for k in derivs.keys():
     1822                    if len(k)-Integer(1)  < zero_order:
     1823                        derivs[k]= Integer(0)
     1824        if sub_final:
     1825            # Substitute sub_final into the values of derivs.
     1826            for k in derivs.keys():
     1827                derivs[k] = self._subs_all(derivs[k],sub_final)
     1828        if rekey:
     1829            # Rekey the derivs dictionary by the value of rekey.
     1830            F= rekey
     1831            if singleton:
     1832                # F must be a singleton.
     1833                derivs= dict( [(diff(F,list(k)), derivs[k]) for k in derivs.keys()] )
     1834            else:
     1835                # F must be a list.
     1836                derivs= dict( [(diff(F[k[Integer(0) ]],list(k)[Integer(1) :]), derivs[k]) for k in derivs.keys()] )
     1837        return derivs
     1838#-------------------------------------------------------------------------------
     1839    def _diff_op(self,A,B,AB_derivs,V,M,r,N):
     1840        r"""
     1841        This function computes the derivatives `DD^(l+k)(A[j] B^l)` evaluated at a
     1842        point `p` for various natural numbers `j,k,l` which depend on `r` and `N`.
     1843        Here `DD` is a specific second-order linear differential operator that depends
     1844        on `M`, `A` is a list of symbolic functions, `B` is symbolic function,
     1845        and `AB_derivs` contains all the derivatives of `A` and `B` evaluated at `p`
     1846        that are necessary for the computation.
     1847        For internal use by the functions _asymptotics_main_multiple() and
     1848        _asymptotics_main_smooth().
     1849
     1850        Does not use `self`.
     1851
     1852        INPUT:
     1853
     1854        - ``A`` - A single or length `r` list of symbolic functions in the
     1855            variables `V`.
     1856        - ``B`` - A symbolic function in the variables `V`.
     1857        - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of
     1858            `A[0],\ldots,A[r-1]` up to order `2N-2` and
     1859            the (symbolic) derivatives of `B` up to order `2N`.
     1860            The values of the dictionary are complex numbers that are
     1861            the keys evaluated at a common point `p`.
     1862        - ``V`` - The variables of the `A[j]` and `B`.
     1863        - ``M`` - A symmetric `l \times l` matrix, where `l` is the length of `V`.
     1864        - ``r,N`` - Natural numbers.
     1865
     1866        OUTPUT:
     1867
     1868        A dictionary whose keys are natural number tuples of the form `(j,k,l)`,
     1869        where `l \le 2k`, `j \le r-1`, and `j+k \le N-1`, and whose values are
     1870        `DD^(l+k)(A[j] B^l)` evaluated at a point `p`, where `DD` is the linear
     1871        second-order differential operator
     1872        `-\sum_{i=0}^{l-1} \sum_{j=0}^{l-1} M[i][j]
     1873        \partial^2 /(\partial V[j] \partial V[i])`.
     1874
     1875        EXAMPLES::
     1876
     1877            sage: R.<x> = PolynomialRing(QQ)
     1878            sage: G,H = 1,1
     1879            sage: F= QuasiRationalExpression(G,H)
     1880            sage: T= var('x,y')
     1881            sage: A= function('A',*tuple(T))
     1882            sage: B= function('B',*tuple(T))
     1883            sage: AB_derivs= {}
     1884            sage: M= matrix([[1,2],[2,1]])
     1885            sage: DD= F._diff_op(A,B,AB_derivs,T,M,1,2)
     1886            sage: DD.keys()
     1887            [(0, 1, 2), (0, 1, 1), (0, 1, 0), (0, 0, 0)]
     1888            sage: len(DD[(0,1,2)])
     1889            246
     1890
     1891        AUTHORS:
     1892
     1893        - Alex Raichev (2008-10-01, 2010-01-12)
     1894        """
     1895        if not isinstance(A,list):
     1896            A= [A]
     1897
     1898        # First, compute the necessary product derivatives of A and B.
     1899        product_derivs= {}
     1900        for (j,k) in mrange([r,N]):
     1901            if j+k <N:
     1902                for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)):
     1903                    for s in UnorderedTuples(V,Integer(2) *(k+l)):
     1904                        product_derivs[tuple([j,k,l]+s)] = \
     1905                        diff(A[j]*B**l,s).subs(AB_derivs)
     1906
     1907        # Second, compute DD^(k+l)(A[j]*B^l)(p) and store values in dictionary.
     1908        DD= {}
     1909        rows= M.nrows()
     1910        for (j,k) in mrange([r,N]):
     1911            if j+k <N:
     1912                for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)):
     1913                    # Take advantage of the symmetry of M by ignoring
     1914                    # the upper-diagonal entries of M and multiplying by
     1915                    # appropriate powers of 2.
     1916                    if k+l == Integer(0) :
     1917                        DD[(j,k,l)] = product_derivs[(j,k,l)]
     1918                        continue
     1919                    S= [(a,b) for (a,b) in mrange([rows,rows]) if b <= a]
     1920                    P=  cartesian_product_iterator([S for i in range(k+l)])
     1921                    diffo= Integer(0)
     1922                    for t in P:
     1923                        if product_derivs[(j,k,l)+self._diff_seq(V,t)] != Integer(0) :
     1924                            MM= Integer(1)
     1925                            for (a,b) in t:
     1926                                MM= MM * M[a][b]
     1927                                if a != b:
     1928                                    MM= Integer(2) *MM
     1929                            diffo= diffo + MM * product_derivs[(j,k,l)+self._diff_seq(V,t)]
     1930                    DD[(j,k,l)] = (-Integer(1) )**(k+l)*diffo
     1931        return DD
     1932#-------------------------------------------------------------------------------
     1933    def _diff_op_simple(self,A,B,AB_derivs,x,v,a,N):
     1934        r"""
     1935        This function computes `DD^(e k + v l)(A B^l)` evaluated at a point `p`
     1936        for various natural numbers `e,k,l` that depend on `v` and `N`.
     1937        Here `DD` is a specific linear differential operator that depends
     1938        on `a` and `v`, `A` and `B` are symbolic functions, and `AB_derivs` contains
     1939        all the derivatives of `A` and `B` evaluated at `p` that are necessary for
     1940        the computation.
     1941        For internal use by the function _asymptotics_main_smooth().
     1942
     1943        Does not use `self`.
     1944
     1945        INPUT:
     1946
     1947        - ``A``,``B`` - Symbolic functions in the variable `x`.
     1948        - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of
     1949            `A` up to order `2N` if `v` is even or `N` if `v` is odd and
     1950            the (symbolic) derivatives of `B` up to order `2N+v` if `v` is even
     1951            or `N+v` if `v` is odd.
     1952            The values of the dictionary are complex numbers that are
     1953            the keys evaluated at a common point `p`.
     1954        - ``x`` - Symbolic variable.
     1955        - ``a`` - A complex number.
     1956        - ``v``,``N`` - Natural numbers.
     1957
     1958        OUTPUT:
     1959
     1960        A dictionary whose keys are natural number pairs of the form `(k,l)`,
     1961        where `k < N` and `l \le 2k` and whose values are
     1962        `DD^(e k + v l)(A B^l)` evaluated at a point `p`.
     1963        Here `e=2` if `v` is even, `e=1` if `v` is odd, and `DD` is a particular
     1964        linear differential operator
     1965        `(a^{-1/v} d/dt)' if `v` is even and `(|a|^{-1/v} i \sgn a d/dt)`
     1966        if `v` is odd.
     1967
     1968        EXAMPLES::
     1969
     1970            sage: R.<x> = PolynomialRing(QQ)
     1971            sage: G,H = 1,1
     1972            sage: F= QuasiRationalExpression(G,H)
     1973            sage: A= function('A',x)
     1974            sage: B= function('B',x)
     1975            sage: AB_derivs= {}
     1976            sage: F._diff_op_simple(A,B,AB_derivs,x,3,2,2)
     1977            {(1, 0): 1/2*I*2^(2/3)*D[0](A)(x), (0, 0): A(x), (1, 1): 1/4*(A(x)*D[0, 0, 0, 0](B)(x) + B(x)*D[0, 0, 0, 0](A)(x) + 4*D[0](A)(x)*D[0, 0, 0](B)(x) + 4*D[0](B)(x)*D[0, 0, 0](A)(x) + 6*D[0, 0](A)(x)*D[0, 0](B)(x))*2^(2/3)}
     1978
     1979        AUTHORS:
     1980
     1981        - Alex Raichev (2010-01-15)
     1982        """
     1983        I= sqrt(-Integer(1) )
     1984        DD= {}
     1985        if v.mod(Integer(2) ) == Integer(0) :
     1986            for k in range(N):
     1987                for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)):
     1988                    DD[(k,l)] = (a**(-Integer(1) /v))**(Integer(2) *k+v*l) \
     1989                    * diff(A*B**l,x,Integer(2) *k+v*l).subs(AB_derivs)
     1990        else:
     1991            for k in range(N):
     1992                for l in (ellipsis_range(Integer(0) ,Ellipsis,k)):
     1993                    DD[(k,l)] = (abs(a)**(-Integer(1) /v)*I*a/abs(a))**(k+v*l) \
     1994                    * diff(A*B**l,x,k+v*l).subs(AB_derivs)
     1995        return DD
     1996#-------------------------------------------------------------------------------
     1997    def _diff_prod(self,f_derivs,u,g,X,interval,end,uderivs,atc):
     1998        r"""
     1999        This function takes various derivatives of the equation `f=ug`,
     2000        evaluates at a point `c`, and solves for the derivatives of `u`.
     2001        For internal use by the function _asymptotics_main_multiple().
     2002
     2003        Does not use `self`.
     2004
     2005        INPUT:
     2006
     2007        - ``f_derivs`` - A dictionary whose keys are tuples \code{s + end} for all
     2008            `X`-variable sequences `s` with length in `interval` and whose
     2009            values are the derivatives of a function `f` evaluated at `c`.
     2010        - ``u`` - A callable symbolic function.
     2011        - ``g`` - An expression or callable symbolic function.
     2012        - ``X`` - A list of symbolic variables.
     2013        - ``interval`` - A list of positive integers.
     2014            Call the first and last values `n` and `nn`, respectively.
     2015        - ``end`` - A possibly empty list of `z`'s where `z` is the last element of
     2016            `X`.
     2017        - ``uderivs`` - A dictionary whose keys are the symbolic
     2018            derivatives of order 0 to order `n-1` of `u` evaluated at `c`
     2019            and whose values are the corresponding derivatives evaluated at `c`.
     2020        - ``atc`` - A dictionary whose keys are the keys of `c` and all the symbolic
     2021            derivatives of order 0 to order `nn` of `g` evaluated `c` and whose
     2022            values are the corresponding derivatives evaluated at `c`.
     2023
     2024        OUTPUT:
     2025
     2026        A dictionary whose keys are the derivatives of `u` up to order
     2027        `nn` and whose values are those derivatives evaluated at `c`.
     2028
     2029        EXAMPLES::
     2030
     2031        I'd like to print out the entire dictionary, but that does not give
     2032        consistent output for doctesting.
     2033        Order of keys changes ::
     2034
     2035            sage: R.<x> = PolynomialRing(QQ)
     2036            sage: G,H = 1,1
     2037            sage: F= QuasiRationalExpression(G,H)
     2038            sage: u= function('u',x)
     2039            sage: g= function('g',x)
     2040            sage: dd= F._diff_prod({(x,):1,(x,x):1},u,g,[x],[1,2],[],{u(x=2):1},{x:2,g(x=2):3,diff(g,x)(x=2):5, diff(g,x,x)(x=2):7})
     2041            sage: dd[diff(u,x,2)(x=2)]
     2042            22/9
     2043
     2044        NOTES:
     2045
     2046        This function works by differentiating the equation `f=ug` with respect
     2047        to the variable sequence \code{s+end}, for all tuples `s` of `X` of
     2048        lengths in `interval`, evaluating at the point `c`,
     2049        and solving for the remaining derivatives of `u`.
     2050        This function assumes that `u` never appears in the differentiations of
     2051        `f=ug` after evaluating at `c`.
     2052
     2053        AUTHORS:
     2054
     2055        - Alex Raichev (2009-05-14, 2010-01-21)
     2056        """
     2057        for l in interval:
     2058            D= {}
     2059            rhs= []
     2060            lhs= []
     2061            for t in UnorderedTuples(X,l):
     2062                s= t+end
     2063                lhs.append(f_derivs[tuple(s)])
     2064                rhs.append(diff(u*g,s).subs(atc).subs(uderivs))
     2065                # Since Sage's solve command can't take derivatives as variable
     2066                # names, i make new variables based on t to stand in for
     2067                # diff(u,t) and store them in D.
     2068                D[diff(u,t).subs(atc)] = self._make_var([var('zing')]+t)
     2069            eqns=[ lhs[i] == rhs[i].subs(uderivs).subs(D) for i in range(len(lhs))]
     2070            vars= D.values()
     2071            sol= solve(eqns,*vars,solution_dict=True)
     2072            uderivs.update(self._subs_all(D,sol[Integer(0) ]))
     2073        return uderivs
     2074#-------------------------------------------------------------------------------
     2075    def _diff_seq(self,V,s):
     2076        r"""
     2077        Given a list `s` of tuples of natural numbers, this function returns the
     2078        list of elements of `V` with indices the elements of the elements of `s`.
     2079        This function is for internal use by the function _diff_op().
     2080
     2081        Does not use `self`.
     2082
     2083        INPUT:
     2084
     2085        - ``V`` - A list.
     2086        - ``s`` - A list of tuples of natural numbers in the interval
     2087            \code{range(len(V))}.
     2088
     2089        OUTPUT:
     2090
     2091        The tuple \code{tuple([V[tt] for tt in sorted(t)])}, where `t` is the
     2092        list of elements of the elements of `s`.
     2093
     2094        EXAMPLES::
     2095
     2096            sage: R.<x> = PolynomialRing(QQ)
     2097            sage: G,H = 1,1
     2098            sage: F= QuasiRationalExpression(G,H)
     2099            sage: V= list(var('x,t,z'))
     2100            sage: F._diff_seq(V,([0,1],[0,2,1],[0,0]))
     2101            (x, x, x, x, t, t, z)
     2102
     2103        AUTHORS:
     2104
     2105        - Alex Raichev (2009.05.19)
     2106        """
     2107        t= []
     2108        for ss in s:
     2109            t.extend(ss)
     2110        return tuple([V[tt] for tt in sorted(t)])
     2111#-------------------------------------------------------------------------------
     2112    def _direction(self,v,coordinate=None):
     2113        r"""
     2114        This function returns \code{[vv/v[coordinate] for vv in v]} where
     2115        coordinate is the last index of v if not specified otherwise.
     2116
     2117        Does not use `self`.
     2118
     2119        INPUT:
     2120
     2121        - ``v`` - A vector.
     2122        - ``coordinate`` - An index for `v` (default: None).
     2123
     2124        OUTPUT:
     2125
     2126        This function returns \code{[vv/v[coordinate] for vv in v]} where
     2127        coordinate is the last index of v if not specified otherwise.
     2128
     2129        EXAMPLES::
     2130
     2131            sage: R.<x> = PolynomialRing(QQ)
     2132            sage: G,H = 1,1
     2133            sage: F= QuasiRationalExpression(G,H)
     2134            sage: F._direction([2,3,5])
     2135            (2/5, 3/5, 1)
     2136            sage: F._direction([2,3,5],0)
     2137            (1, 3/2, 5/2)
     2138
     2139        AUTHORS:
     2140
     2141        - Alex Raichev (2010-08-25)
     2142        """
     2143        if coordinate == None:
     2144            coordinate= len(v)-Integer(1)
     2145        try:
     2146            v[Integer(0) ].variables()
     2147            return tuple([(vv/v[coordinate]).simplify_full() for vv in v])
     2148        except:
     2149            return tuple([vv/v[coordinate] for vv in v])
     2150# E ============================================================================
     2151# F ============================================================================
     2152# G ============================================================================
     2153# H ============================================================================
     2154# I ============================================================================
     2155    def is_cmp(self,points):
     2156        r"""
     2157        Checks if the points in the list `points` are convenient multiple
     2158        points of `V= \{ x\in CC^d : H(x) = 0\}`, where `H=self._H`.
     2159
     2160        INPUT:
     2161
     2162        - ``points`` - An individual or list of dictionaries with keys
     2163            `self._variables` and values in some superfield of
     2164            `self._R.base_ring()`.
     2165
     2166        OUTPUT:
     2167
     2168        A list of tuples `(p,verdict,comment)`, one for each point
     2169        `p` in `points`, where `verdict` is True if `p` is a convenient
     2170        multiple point and False otherwise, and where `comment` is a string
     2171        comment relating to `verdict`, such as 'not a transverse intersection'.
     2172
     2173        EXAMPLES::
     2174
     2175            sage: R.<x,y,z>= PolynomialRing(QQ)
     2176            sage: G= 16
     2177            sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y))
     2178            sage: F= QuasiRationalExpression(G,H)
     2179            sage: points= [{x:1/2,y:1,z:4/3},{x:1/2,y:1,z:2}]
     2180            sage: F.is_cmp(points)
     2181            [({y: 1, z: 4/3, x: 1/2}, True, 'all good'), ({y: 1, z: 2, x: 1/2}, False, 'not a singular point')]
     2182
     2183        NOTES:
     2184
     2185        A point `c` of `V` is a __convenient multiple point__ if `V` is locally
     2186        a union of complex manifolds that intersect transversely at `c`;
     2187        see [RaWi2011]_.
     2188
     2189        REFERENCES:
     2190
     2191        .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
     2192        coefficients of multivariate generating functions: improvements
     2193        for smooth points", To appear.
     2194
     2195        AUTHORS:
     2196
     2197        - Alex Raichev (2011-04-18)
     2198        """
     2199        H= self._H
     2200        d= self._d
     2201        Hs= [SR(h[Integer(0) ]) for h in self._Hfac]  # irreducible factors of H
     2202        X= self._variables
     2203        J= [tuple([diff(h,x) for x in X]) for h in Hs]
     2204        verdicts= []
     2205        if not isinstance(points,list):
     2206            points= [points]
     2207        for p in points:
     2208            # Ensure variables in points lie in SR.
     2209            pp= {}
     2210            for x in p.keys():
     2211                pp[SR(x)]= p[x]
     2212            # Test 1: Is p a zero of all factors of H?
     2213            if [h.subs(pp) for h in Hs] != [Integer(0)  for h in Hs]:
     2214                # Failed test 1.  Move on to next point.
     2215                verdicts.append((p,False,'not a singular point'))
     2216                continue
     2217            # Test 2: Are the factors of H smooth and
     2218            # do theyintersect transversely at p?
     2219            J= [tuple([f.subs(pp) for f in dh]) for dh in J]
     2220            l= len(J)
     2221            if Set(J).cardinality() < l:
     2222                # Fail.  Move on to next point.
     2223                verdicts.append((p,False,'not a transverse intersection'))
     2224                continue
     2225            temp= True
     2226            for S in list(Set(J).subsets())[Integer(1) :]: # Subsets of size >= 1
     2227                k= len(list(S))
     2228                M = matrix(list(S))
     2229                if  M.rank() != min(k,d):
     2230                    # Fail.
     2231                    temp= False
     2232                    verdicts.append((p,False,'not a transvere intersection'))
     2233                    break
     2234            if not temp:
     2235                # Move on to next point
     2236                continue
     2237            # Test 3: Is p convenient?
     2238            Jlog= matrix([self._log_grad(h,X,pp) for h in Hs])
     2239            if [Integer(0)  in f for f in Jlog.columns()] ==\
     2240            [True for f in Jlog.columns()]:
     2241                # Fail.  Move on to next point.
     2242                verdict[p] = (False,'multiple point but not convenient')
     2243                continue
     2244            verdicts.append((p,True,'all good'))
     2245        return verdicts
     2246# J ============================================================================
     2247# K ============================================================================
     2248# L ============================================================================
     2249    def _log_grad(self,f,X,c):
     2250        r"""
     2251        This function returns the logarithmic gradient of `f` with respect to the
     2252        variables of `X` evalutated at `c`.
     2253
     2254        Does not use `self`.
     2255
     2256        INPUT:
     2257
     2258        - ``f`` - An expression in the variables of `X`.
     2259        - ``X`` - A list of variables.
     2260        - ``c`` - A dictionary with keys `X` and values in a field `K`.
     2261
     2262        OUTPUT:
     2263
     2264            \code{[c[x] * diff(f,x).subs(c) for x in X]}.
     2265
     2266        EXAMPLES::
     2267
     2268            sage: R.<x> = PolynomialRing(QQ)
     2269            sage: G,H = 1,1
     2270            sage: F= QuasiRationalExpression(G,H)
     2271            sage: X= var('x,y,z')
     2272            sage: f= x*y*z^2
     2273            sage: c= {x:1,y:2,z:3}
     2274            sage: f.gradient()
     2275            (y*z^2, x*z^2, 2*x*y*z)
     2276            sage: F._log_grad(f,X,c)
     2277            (18, 18, 36)
     2278
     2279            ::
     2280
     2281            sage: R.<x,y,z>= PolynomialRing(QQ)
     2282            sage: f= x*y*z^2
     2283            sage: c= {x:1,y:2,z:3}
     2284            sage: F._log_grad(f,R.gens(),c)
     2285            (18, 18, 36)
     2286
     2287        AUTHORS:
     2288
     2289        - Alex Raichev (2009-03-06)
     2290        """
     2291        return tuple([SR(c[x] * diff(f,x).subs(c)).simplify() for x in X])
     2292# M ============================================================================
     2293    def _make_var(self,L):
     2294        r"""
     2295        This function converts the list `L` to a string (without commas) and returns
     2296        the string as a variable.
     2297        For internal use by the function _diff_op()
     2298
     2299        Does not use `self`.
     2300
     2301        INPUT:
     2302
     2303        - ``L`` - A list.
     2304
     2305        OUTPUT:
     2306
     2307        A variable whose name is the concatenation of the variable names in `L`.
     2308
     2309        EXAMPLES::
     2310
     2311            sage: R.<x> = PolynomialRing(QQ)
     2312            sage: G,H = 1,1
     2313            sage: F= QuasiRationalExpression(G,H)
     2314            sage: L= list(var('x,y,hello'))
     2315            sage: v= F._make_var(L)
     2316            sage: print v, type(v)
     2317            xyhello <type 'sage.symbolic.expression.Expression'>
     2318
     2319        AUTHORS:
     2320
     2321        - Alex Raichev (2010-01-21)
     2322        """
     2323        return var(''.join([str(v) for v in L]))
     2324# N ============================================================================
     2325    def _new_var_name(self,name,V):
     2326        r"""
     2327        This function returns the first string in the sequence `name`, `name+name`,
     2328        `name+name+name`,... that does not appear in the list `V`.
     2329        It is for internal use by the function _asymptotics_main_multiple().
     2330
     2331        Does not use `self`.
     2332
     2333        INPUT:
     2334
     2335        - ``name`` - A string.
     2336        - ``V`` - A list of variables.
     2337
     2338        OUTPUT:
     2339
     2340        The first string in the sequence `name`, `name+name`,
     2341        `name+name+name`,... that does not appear in the list \code{str(V)}.
     2342
     2343        EXAMPLES::
     2344
     2345            sage: R.<x> = PolynomialRing(QQ)
     2346            sage: G,H = 1,1
     2347            sage: F= QuasiRationalExpression(G,H)
     2348            sage: X= var('x,xx,y,z')
     2349            sage: F._new_var_name('x',X)
     2350            'xxx'
     2351
     2352        AUTHORS:
     2353
     2354        - Alex Raichev (2008-10-01)
     2355        """
     2356        newname= name
     2357        while newname in str(V):
     2358            newname= newname +name
     2359        return newname
     2360# O ============================================================================
     2361# P ============================================================================
     2362# Q ============================================================================
     2363# R ============================================================================
     2364    def relative_error(self,approx,alpha,interval,exp_scale=Integer(1) ):
     2365        r"""
     2366        Returns the relative error between the values of the Maclaurin
     2367        coefficients of `self` with multi-indices `m alpha` for `m` in
     2368        `interval` and the values of the functions in `approx`.
     2369
     2370        INPUT:
     2371
     2372        - ``approx`` - An individual or list of symbolic expressions in
     2373            one variable.
     2374        - ``alpha`` - A list of positive integers.
     2375        - ``interval`` - A list of positive integers.
     2376        - ``exp_scale`` - (optional) A number.  Default: 1.
     2377
     2378        OUTPUT:
     2379
     2380        A list whose entries are of the form
     2381        `[m\alpha,a_m,b_m,err_m]` for `m \in interval`.
     2382        Here `m\alpha` is a tuple; `a_m` is the `m alpha` (multi-index)
     2383        coefficient of the Maclaurin series for `F` divided by `exp_scale^m`;
     2384        `b_m` is a list of the values of the functions in `approx` evaluated at
     2385        `m` and divided by `exp_scale^m`; `err_m` is the list of relative errors
     2386        `(a_m-f)/a_m` for `f` in `b_m`.
     2387        All outputs are decimal approximations.
     2388
     2389        EXAMPLES::
     2390
     2391            sage: R.<x,y>= PolynomialRing(QQ)
     2392            sage: G=1
     2393            sage: H= 1-x-y-x*y
     2394            sage: F= QuasiRationalExpression(G,H)
     2395            sage: alpha= [1,1]
     2396            sage: var('n')
     2397            n
     2398            sage: f= (0.573/sqrt(n))*5.83^n
     2399            sage: es= 5.83
     2400            sage: F.relative_error(f,alpha,[1,2,4,8],es)
     2401            Calculating errors table in the form
     2402            exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors...
     2403            [[(1, 1), 0.514579759862779, [0.573000000000000], [-0.113530000000000]], [(2, 2), 0.382477808931739, [0.405172185619892], [-0.0593351461396876]], [(4, 4), 0.277863059517142, [0.286500000000000], [-0.0310834426780842]], [(8, 8), 0.199108827584423, [0.202586092809946], [-0.0174641439443390]]]
     2404            sage: g= (0.573/sqrt(n) - 0.0674/n^(3/2))*5.83^n
     2405            sage: F.relative_error([f,g],alpha,[1,2,4,8],es)
     2406            Calculating errors table in the form
     2407            exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors...
     2408            [[(1, 1), 0.514579759862779, [0.573000000000000, 0.505600000000000], [-0.113530000000000, 0.0174506666666667]], [(2, 2), 0.382477808931739, [0.405172185619892, 0.381342687093905], [-0.0593351461396876, 0.00296781097184384]], [(4, 4), 0.277863059517142, [0.286500000000000, 0.278075000000000], [-0.0310834426780842, -0.000762751562681505]], [(8, 8), 0.199108827584423, [0.202586092809946, 0.199607405494198], [-0.0174641439443390, -0.00250404723800224]]]
     2409
     2410        AUTHORS:
     2411
     2412        - Alex Raichev (2009-05-18, 2011-04-18)
     2413        """
     2414
     2415        if not isinstance(approx,list):
     2416            approx= [approx]
     2417        av= approx[Integer(0) ].variables()[Integer(0) ]
     2418
     2419        print "Calculating errors table in the form"
     2420        print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..."
     2421
     2422        # Get Maclaurin coefficients of self and scale them.
     2423        # Then compute errors.
     2424        n= interval[-Integer(1) ]
     2425        mc= self.maclaurin_coefficients(alpha,n)
     2426        mca={}
     2427        stats=[]
     2428        for m in interval:
     2429            beta= tuple([m*a for a in alpha])
     2430            mc[beta]= mc[beta]/exp_scale**m
     2431            mca[beta]= [f.subs({av:m})/exp_scale**m for f in approx]
     2432            stats_row= [beta, mc[beta].n(), [a.n() for a in mca[beta]]]
     2433            if mc[beta]==Integer(0) :
     2434                 stats_row.extend([None for a in mca[beta]])
     2435            else:
     2436                 stats_row.append([((mc[beta]-a)/mc[beta]).n() for a in mca[beta]])
     2437            stats.append(stats_row)
     2438        return stats
     2439# S ============================================================================
     2440    def singular_points(self):
     2441        r"""
     2442        This function returns a Groebner basis ideal whose variety is the
     2443        set of singular points of the algebraic variety
     2444        `V= \{x\in\CC^d : H(x)=0\}`, where `H=sef._H`.
     2445
     2446        INPUT:
     2447
     2448        OUTPUT:
     2449
     2450        A Groebner basis ideal whose variety is the set of singular points of
     2451        the algebraic variety `V= \{x\in\CC^d : H(x)=0\}`.
     2452
     2453        EXAMPLES::
     2454
     2455            sage: R.<x,y,z>= PolynomialRing(QQ)
     2456            sage: G= 1
     2457            sage: H= (4-2*x-y-z)*(4-x-2*y-z)
     2458            sage: F= QuasiRationalExpression(G,H)
     2459            sage: F.singular_points()
     2460            Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field
     2461
     2462        AUTHORS:
     2463
     2464        - Alex Raichev (2008-10-01, 2008-11-20, 2010-12-03, 2011-04-18)
     2465        """
     2466        H= self._H
     2467        R= self._R
     2468        f= R.ideal(H).radical().gens()[Integer(0) ]     # Compute the reduction of H.
     2469        J= R.ideal([f] + f.gradient())
     2470        return R.ideal(J.groebner_basis())
     2471#-------------------------------------------------------------------------------
     2472    def smooth_critical(self,alpha):
     2473        r"""
     2474        This function returns a Groebner basis ideal whose variety is the set
     2475        of smooth critical points of the algebraic variety
     2476        `V= \{x\in\CC^d : H(x)=0\} for the direction `\alpha` where `H=self._H`.
     2477
     2478        INPUT:
     2479
     2480        - ``alpha`` - A `d`-tuple of positive integers and/or symbolic entries.
     2481
     2482        OUTPUT:
     2483
     2484        A Groebner basis ideal of smooth critical points of `V` for `\alpha`.
     2485
     2486        EXAMPLES::
     2487
     2488            sage: R.<x,y> = PolynomialRing(QQ)
     2489            sage: G=1
     2490            sage: H = (1-x-y-x*y)^2
     2491            sage: F= QuasiRationalExpression(G,H)
     2492            sage: var('a1,a2')
     2493            (a1, a2)
     2494            sage: F.smooth_critical([a1,a2])
     2495            Ideal (y^2 + 2*a1/a2*y - 1, x + (a2/(-a1))*y + (-a2 + a1)/(-a1)) of Multivariate Polynomial Ring in x, y over Fraction Field of Multivariate Polynomial Ring in a2, a1 over Rational Field
     2496
     2497        NOTES:
     2498
     2499        A point `c` of `V` is a __smooth critical point for `alpha`__
     2500        if the gradient of `f` at `c` is not identically zero and `\alpha` is in
     2501        the span of the logarithmic gradient vector
     2502        `(c[0] \partial_1 f(c)),\ldots,c[d-1] \partial_d f(c))`; see [RaWi2008a]_.
     2503
     2504        REFERENCES:
     2505
     2506        .. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
     2507        coefficients of multivariate generating functions: improvements
     2508        for smooth points", Electronic Journal of Combinatorics, Vol. 15
     2509        (2008), R89.
     2510
     2511        AUTHORS:
     2512
     2513        - Alex Raichev (2008-10-01, 2008-11-20, 2009-03-09, 2010-12-02, 2011-04-18)
     2514        """
     2515        H= self._H
     2516        R= H.parent()
     2517        B= R.base_ring()
     2518        d= R.ngens()
     2519        vars= R.gens()
     2520        f= R.ideal(H).radical().gens()[Integer(0) ]   # Compute the reduction of H.
     2521
     2522        # Expand B by the variables of alpha if there are any.
     2523        indets= []
     2524        indets_ind= []
     2525        for a in alpha:
     2526            if not ((a in ZZ) and (a>Integer(0) )):
     2527                try:
     2528                    CC(a)
     2529                except:
     2530                    indets.append(var(a))
     2531                    indets_ind.append(alpha.index(a))
     2532                else:
     2533                    print "The components of", alpha, \
     2534                    "must be positive integers or symbolic variables."
     2535                    return
     2536        indets= list(Set(indets))   # Delete duplicates in indets.
     2537        if indets != []:
     2538            BB= FractionField(PolynomialRing(B,tuple(indets)))
     2539            S= R.change_ring(BB)
     2540            vars= S.gens()
     2541            # Coerce alpha into BB.
     2542            for i in range(len(alpha)):
     2543                alpha[i] = BB(alpha[i])
     2544        else:
     2545            S= R
     2546
     2547        # Find smooth, critical points for alpha.
     2548        f= S(f)
     2549        J= S.ideal([f] +[ alpha[d-Integer(1) ]*vars[i]*diff(f,vars[i]) \
     2550        -alpha[i]*vars[d-Integer(1) ]*diff(f,vars[d-Integer(1) ]) for i in range(d-Integer(1) )])
     2551        return S.ideal(J.groebner_basis())
     2552#-------------------------------------------------------------------------------
     2553    def _subs_all(self,f,sub,simplify=False):
     2554        r"""
     2555        This function returns the items of `f` substituted by the dictionaries of
     2556        `sub` in order of their appearance in `sub`.
     2557
     2558        Does not use `self`.
     2559
     2560        INPUT:
     2561
     2562        - ``f`` - An individual or list of symbolic expressions or dictionaries
     2563        - ``sub`` - An individual or list of dictionaries.
     2564        - ``simplify`` - Boolean (default: False).
     2565
     2566        OUTPUT:
     2567
     2568        The items of `f` substituted by the dictionaries of `sub` in order of
     2569        their appearance in `sub`.  The subs() command is used.
     2570        If simplify is True, then simplify() is used after substitution.
     2571
     2572        EXAMPLES::
     2573
     2574            sage: R.<x> = PolynomialRing(QQ)
     2575            sage: G,H = 1,1
     2576            sage: F= QuasiRationalExpression(G,H)
     2577            sage: var('x,y,z')
     2578            (x, y, z)
     2579            sage: a= {x:1}
     2580            sage: b= {y:2}
     2581            sage: c= {z:3}
     2582            sage: F._subs_all(x+y+z,a)
     2583            y + z + 1
     2584            sage: F._subs_all(x+y+z,[c,a])
     2585            y + 4
     2586            sage: F._subs_all([x+y+z,y^2],b)
     2587            [x + z + 2, 4]
     2588            sage: F._subs_all([x+y+z,y^2],[b,c])
     2589            [x + 5, 4]
     2590
     2591            ::
     2592
     2593            sage: var('x,y')
     2594            (x, y)
     2595            sage: a= {'foo':x^2+y^2, 'bar':x-y}
     2596            sage: b= {x:1,y:2}
     2597            sage: F._subs_all(a,b)
     2598            {'foo': 5, 'bar': -1}
     2599
     2600        AUTHORS:
     2601
     2602        - Alex Raichev (2009-05-05)
     2603        """
     2604        singleton= False
     2605        if not isinstance(f,list):
     2606            f= [f]
     2607            singleton= True
     2608        if not isinstance(sub,list):
     2609            sub= [sub]
     2610        g= []
     2611        for ff in f:
     2612            for D in sub:
     2613                if isinstance(ff,dict):
     2614                    ff= dict( [(k,ff[k].subs(D)) for k in ff.keys()] )
     2615                else:
     2616                    ff= ff.subs(D)
     2617            g.append(ff)
     2618        if singleton and simplify:
     2619            if isinstance(g[Integer(0) ],dict):
     2620                return g[Integer(0) ]
     2621            else:
     2622                return g[Integer(0) ].simplify()
     2623        elif singleton and not simplify:
     2624            return g[Integer(0) ]
     2625        elif not singleton and simplify:
     2626            G= []
     2627            for gg in g:
     2628                if isinstance(gg,dict):
     2629                    G.append(gg)
     2630                else:
     2631                    G.append(gg.simplify())
     2632            return G
     2633        else:
     2634            return g
     2635# T ============================================================================
     2636    def maclaurin_coefficients(self,alpha,n):
     2637        r"""
     2638        Returns the Maclaurin coefficients of self that have multi-indices
     2639        `alpha`, `2*alpha`,...,`n*alpha`.
     2640
     2641        INPUT:
     2642
     2643        - ``n`` - positive integer
     2644        - ``alpha`` - tuple of positive integers representing a multi-index.
     2645
     2646        OUTPUT:
     2647
     2648        A dictionary of the form (beta, beta Maclaurin coefficient of self).
     2649
     2650        AUTHORS:
     2651
     2652        - Alex Raichev (2011-04-08)
     2653        """
     2654        # Do all computations in the Symbolic Ring.
     2655        F= SR(self._G/self._H)
     2656        v= F.variables()
     2657        p= {}
     2658        for x in v:
     2659            p[x]=Integer(0)
     2660        d= len(v)
     2661        coeffs={}
     2662        # Initialize sequence of variables to differentiate with.
     2663        s=[]
     2664        for i in range(d):
     2665            s.extend([v[i] for j in range(alpha[i])])
     2666        F_deriv= diff(F,s)
     2667        coeffs[tuple(alpha)]= F_deriv.subs(p)/ mul([factorial(a) for a in alpha])
     2668        old_beta= alpha
     2669        for k in (ellipsis_range(Integer(2) ,Ellipsis,n)):
     2670            # Update variable sequence to differentiate with.
     2671            beta= [k*a for a in alpha]
     2672            delta= [beta[i]-old_beta[i] for i in range(d)]
     2673            s= []
     2674            for i in range(d):
     2675                s.extend([v[i] for j in range(delta[i])])
     2676            F_deriv= diff(F_deriv,s)
     2677            coeffs[tuple(beta)]= F_deriv.subs(p)/ mul([factorial(b) for b in beta])
     2678            old_beta= copy(beta)
     2679        return coeffs
     2680# U ============================================================================
     2681# V ============================================================================
     2682    def variables(self):
     2683        r"""
     2684        Returns the tuple of variables of `self`.
     2685
     2686        EXAMPLES::
     2687
     2688            sage: R.<x,y>= PolynomialRing(QQ)
     2689            sage: G= exp(x)
     2690            sage: H= 1-y
     2691            sage: F= QuasiRationalExpression(G,H)
     2692            sage: F.variables()
     2693            (x, y)
     2694
     2695            sage: R.<x,y>= PolynomialRing(QQ,order='invlex')
     2696            sage: G= exp(x)
     2697            sage: H= 1-y
     2698            sage: F= QuasiRationalExpression(G,H)
     2699            sage: F.variables()
     2700            (y, x)
     2701
     2702        AUTHORS:
     2703
     2704        - Alex Raichev (2011-04-01)
     2705        """
     2706        return self._variables
     2707# W ============================================================================
     2708# X ============================================================================
     2709# Y ============================================================================
     2710# Z ============================================================================