Ticket #10519: trac_10519-fixed.patch

sage/combinat/all.py

@@ -116 +116 @@
 from cyclic_sieving_phenomenon import CyclicSievingPolynomial, CyclicSievingCheck
 
 from sidon_sets import sidon_sets
+
+from amgf import QuasiRationalExpression

new file sage/combinat/amgf.py

@@ +1 @@
+r"""
+Asymptotics of Maclaurin coefficients of generating functions
+
+This code relates to analytic combinatorics. More specifically, it is a
+collection of functions designed to compute asymptotics of Maclaurin
+coefficients of certain classes of multivariate generating functions.
+
+The main function asymptotics() returns the first `N` terms of
+the asymptotic expansion of the Maclaurin coefficients `F_{n\alpha}`
+of the multivariate meromorphic function `F=G/H` as `n\to\infty`.
+It assumes that `F` is holomorphic in a neighborhood of the origin,
+that `H` is a polynomial, and that asymptotics in the direction of
+`\alpha` (a tuple of positive integers) are controlled by convenient
+smooth or multiple points.
+
+For an introduction to this subject, see [PeWi2008]_.
+The main algorithms and formulas implemented here come from [RaWi2008a]_
+and [RaWi2011]_.
+
+REFERENCES:
+
+.. [AiYu1983]  I. A. A\u\izenberg and A. P. Yuzhakov, "Integral
+representations and residues in multidimensional complex analysis",
+Translations of Mathematical Monographs, 58. American Mathematical
+Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X.
+
+.. [DeLo2006]  Wolfram Decker and Christoph Lossen, "Computing in
+algebraic geometry", Chapter 7.1, Springer-Verlag, 2006.
+
+.. [DiEm2005]  Alicia Dickenstein and Ioannis Z. Emiris (editors),
+"Solving polynomial equations", Chapter 9.0, Springer-Verlag, 2005.
+
+.. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
+several variables into partial fractions", Soviet Math. (Iz. VUZ)
+22 (1978), no. 10, 35--38.
+
+.. [PeWi2008]  Robin Pemantle and Mark C. Wilson, "Twenty combinatorial
+examples of asymptotics derived from multivariate generating
+functions", SIAM Rev. (2008) vol. 50 (2) pp. 199-272.
+
+.. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
+coefficients of multivariate generating functions: improvements
+for smooth points", Electronic Journal of Combinatorics, Vol. 15
+(2008), R89.
+
+.. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
+coefficients of multivariate generating functions: improvements
+for smooth points", To appear.
+
+AUTHORS:
+
+- Alex Raichev (2008-10-01) : Initial version
+- Alex Raichev (2010-09-28) : Corrected many functions
+- Alex Raichev (2010-12-15) : Updated documentation
+- Alex Raichev (2011-03-09) : Fixed a division by zero bug in relative_error()
+- Alex Raichev (2011-04-26) : Rewrote in object-oreinted style
+- Alex Raichev (2011-05-06) : Fixed bug in cohomologous_integrand() and
+    fixed _crit_cone_combo() to work in SR
+
+EXAMPLES:
+
+These come from [RaWi2008a]_ and [RaWi2011]_.
+A smooth point example::
+
+    sage: R.<x,y>= PolynomialRing(QQ)
+    sage: G= 1
+    sage: H= 1-x-y-x*y
+    sage: F= QuasiRationalExpression(G,H)
+    sage: alpha= [3,2]
+    sage: F.maclaurin_coefficients(alpha,8)
+    {(21, 14): 1279919346549, (6, 4): 1289, (3, 2): 25, (18, 12): 19403906105, (24, 16): 85275509086721, (15, 10): 298199265, (9, 6): 75517, (12, 8): 4673345}
+    sage: I= F.smooth_critical(alpha); I
+    Ideal (y^2 + 3*y - 1, x - 2/3*y - 1/3) of Multivariate Polynomial Ring in x, y over Rational Field
+    sage: s= solve(I.gens(),F.variables(),solution_dict=true); s
+    [{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}]
+    sage: c= s[1]
+    sage: asys= [F.asymptotics(c,alpha,n,numerical=3) for n in [1..2]]; asys
+    Initializing auxiliary functions...
+    Calculating derivatives of auxiallary functions...
+    Calculating derivatives of more auxiliary functions...
+    Calculating actions of the second order differential operator...
+    Initializing auxiliary functions...
+    Calculating derivatives of auxiallary functions...
+    Calculating derivatives of more auxiliary functions...
+    Calculating actions of the second order differential operator...
+    [(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))]
+    sage: F.relative_error([f[0] for f in asys],alpha,[1,2,4,8,16],asys[0][1])
+    Calculating errors table in the form
+    exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors...
+    [[(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]]]
+
+Another smooth point example. Here it turns out that the terms 2--4 of the
+asymptotic expansion are 0.
+
+::
+
+    sage: R.<x,y> = PolynomialRing(QQ)
+    sage: q= 1/2
+    sage: qq= q.denominator()
+    sage: G= 1-q*x
+    sage: H= 1-q*x +q*x*y -x^2*y
+    sage: F= QuasiRationalExpression(G,H)
+    sage: alpha= list(qq*vector([2,1-q])); alpha
+    [4, 1]
+    sage: I= F.smooth_critical(alpha); I
+    Ideal (y^2 - 2*y + 1, x + 1/4*y - 5/4) of Multivariate Polynomial Ring in x, y over Rational Field
+    sage: s= solve(I.gens(),F.variables(),solution_dict=true); s
+    [{y: 1, x: 1}]
+    sage: c= s[0]
+    sage: asy= F.asymptotics(c,alpha,5); asy
+    Initializing auxiliary functions...
+    Calculating derivatives of auxiallary functions...
+    Calculating derivatives of more auxiliary functions...
+    Calculating actions of the second order differential operator...
+    (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)))
+    sage: F.relative_error(asy[0],alpha,[1,2],asy[1])
+    Calculating errors table in the form
+    exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors...
+    [[(4, 1), 0.187500000000000, [0.185581424449526], [0.0102324029358597]], [(8, 2), 0.152343750000000, [0.151986595969759], [0.00234439568568301]]]
+
+A multiple point example.
+
+::
+
+    sage: R.<x,y,z>= PolynomialRing(QQ)
+    sage: G= 1
+    sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y))
+    sage: F= QuasiRationalExpression(G,H)
+    sage: I= F.singular_points(); I
+    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
+    sage: c= {x: 1/2, y: 1, z: 4/3}
+    sage: F.is_cmp(c)
+    [({y: 1, z: 4/3, x: 1/2}, True, 'all good')]
+    sage: cc= F.critical_cone(c); cc
+    ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N)
+    sage: alpha= [8,3,3]
+    sage: alpha in cc[1]
+    True
+    sage: asy= F.asymptotics(c,alpha,2); asy
+    Initializing auxiliary functions...
+    Calculating derivatives of auxiliary functions...
+    Calculating derivatives of more auxiliary functions...
+    Calculating second-order differential operator actions...
+    (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)))
+
+Another multiple point example. 
+
+::
+
+    sage: R.<x,y,z>= PolynomialRing(QQ)
+    sage: G= 16
+    sage: H= (4-2*x-y-z)^2*(4-x-2*y-z)
+    sage: F= QuasiRationalExpression(G,H)
+    sage: var('a1,a2,a3')
+    (a1, a2, a3)
+    sage: F.cohomologous_integrand([a1,a2,a3])
+    [[-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]]]
+    sage: I= F.singular_points(); I
+    Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field
+    sage: c= {x:1, y:1, z:1}
+    sage: F.is_cmp(c)
+    [({y: 1, z: 1, x: 1}, True, 'all good')]
+    sage: alpha= [3,3,2]
+    sage: cc= F.critical_cone(c); cc
+    ([(1, 2, 1), (2, 1, 1)], 2-d cone in 3-d lattice N)
+    sage: alpha in cc[1]
+    True
+    sage: asy= F.asymptotics(c,alpha,2); asy
+    Initializing auxiliary functions...
+    Calculating derivatives of auxiliary functions...
+    Calculating derivatives of more auxiliary functions...
+    Calculating second-order differential operator actions...
+    (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)))
+
+"""
+#*****************************************************************************
+#       Copyright (C) 2008 Alexander Raichev <tortoise.said@gmail.com>
+#
+#  Distributed under the terms of the GNU General Public License (GPL)
+#                  http://www.gnu.org/licenses/
+#*****************************************************************************
+
+from copy import copy, deepcopy
+
+# Sage libraries
+from sage.calculus.functional                           import diff
+from sage.calculus.functions                            import jacobian
+from sage.calculus.var                                  import function, var
+from sage.combinat.cartesian_product                    import CartesianProduct
+from sage.combinat.combinat                             import stirling_number1
+from sage.combinat.permutation                          import Permutation
+from sage.combinat.tuple                                import UnorderedTuples
+from sage.functions.log                                 import exp, log
+from sage.functions.other                               import factorial, gamma, sqrt
+from sage.geometry.cone                                 import Cone
+from sage.matrix.constructor                            import matrix
+from sage.misc.misc                                     import add, ellipsis_range, mul
+from sage.misc.misc_c                                   import prod
+from sage.misc.mrange                                   import cartesian_product_iterator, mrange
+from sage.modules.free_module_element                   import vector
+from sage.rings.arith                                   import binomial
+from sage.rings.all                                     import CC
+from sage.rings.fraction_field                          import FractionField
+from sage.rings.integer                                 import Integer
+from sage.rings.integer_ring                            import ZZ
+from sage.rings.polynomial.polynomial_ring_constructor  import PolynomialRing
+from sage.rings.rational_field                          import QQ
+from sage.sets.set                                      import Set
+from sage.structure.sage_object                         import SageObject
+from sage.symbolic.constants                            import pi
+from sage.symbolic.relation                             import solve
+from sage.symbolic.ring                                 import SR
+
+
+# Functions to incorporate later into existing Sage classes.
+def algebraic_dependence(fs):
+    r"""
+    This function returns an irreducible annihilating polynomial for the
+    polynomials in `fs`, if those polynomials are algebraically dependent.
+    Otherwise it returns 0.
+
+    INPUT:
+
+    - ``fs`` - A list of polynomials `f_1,\ldots,f_r' from a common polynomial
+        ring.
+
+    OUTPUT:
+
+    If the polynomials in `fs` are algebraically dependent, then the
+    output is an irreducible polynomial `g` in `K[T_1,\ldots,T_r]` such that
+    `g(f_1,\ldots,f_r) = 0`.
+    Here `K` is the coefficient ring of self and `T_1,\ldots,T_r` are
+    indeterminates different from those of self.
+    If the polynomials in `fs` are algebraically independent, then the output
+    is the zero polynomial.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.amgf import *
+        sage: R.<x>= PolynomialRing(QQ)
+        sage: fs= [x-3, x^2+1]
+        sage: g= algebraic_dependence(fs); g
+        10 + 6*T0 - T1 + T0^2
+        sage: g(fs)
+        0
+
+    ::
+
+        sage: R.<w,z>= PolynomialRing(QQ)
+        sage: fs= [w,  (w^2 + z^2 - 1)^2, w*z - 2]
+        sage: g= algebraic_dependence(fs); g
+        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
+        sage: g(fs)
+        0
+        sage: g.parent()
+        Multivariate Polynomial Ring in T0, T1, T2 over Rational Field
+
+    ::
+
+        sage: R.<x,y,z>= PolynomialRing(GF(7))
+        sage: fs= [x*y,y*z]
+        sage: g= algebraic_dependence(fs); g
+        0
+
+    AUTHORS:
+
+    - Alex Raichev (2011-01-06)
+    """
+    r= len(fs)
+    R= fs[Integer(0) ].parent()
+    B= R.base_ring()
+    Xs= list(R.gens())
+    d= len(Xs)
+
+    # Expand R by r new variables.
+    T= 'T'
+    while T in [str(x) for x in Xs]:
+        T= T+'T'
+    Ts= [T + str(j) for j in range(r)]
+    RR= PolynomialRing(B,d+r,tuple(Xs+Ts))
+    Vs= list(RR.gens())
+    Xs= Vs[Integer(0) :d]
+    Ts= Vs[d:]
+
+    # Find an irreducible annihilating polynomial g for the fs if there is one.
+    J= RR.ideal([ Ts[j] -RR(fs[j]) for j in range(r)])
+    JJ= J.elimination_ideal(Xs)
+    g= JJ.gens()[Integer(0) ]
+
+    # Shrink the ambient ring of g to include only Ts.
+    # I chose the negdeglex order simply because i find it useful in my work.
+    RRR= PolynomialRing(B,r,tuple(Ts),order='negdeglex')
+    return RRR(g)
+#-------------------------------------------------------------------------------
+def partial_fraction_decomposition(f):
+    r"""
+    Return a partial fraction decomposition of the rational expression `f`.
+    Works for univariate and mulitivariate `f`.
+
+    INPUT:
+
+    - ``f`` - An element of the field of fractions `F` of a polynomial ring
+        `R` whose coefficients ring is a field.
+        In the univariate case, the coefficient ring doesn't have to be a field.
+
+    OUTPUT:
+
+    A tuple `whole,parts` where `whole \in R` and `parts` is the list of
+    terms (in `F`) of a partial fraction decomposition of `f - whole`.
+    See the notes below for more details.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.amgf import *
+        sage: S.<t> = QQ[]
+        sage: q = 1/(t+1) + 2/(t+2) + 3/(t-3)
+        sage: whole, parts = partial_fraction_decomposition(q); parts
+        [3/(t - 3), 1/(t + 1), 2/(t + 2)]
+        sage: whole +sum(parts) == q
+        True
+
+    Notice that there is a whole part below despite the appearance of q ::
+
+        sage: R.<x,y>= PolynomialRing(QQ)
+        sage: q= 1/( x *y *(x*y-1) )
+        sage: whole,parts= partial_fraction_decomposition(q)
+        sage: whole, parts
+        (-1, [x*y/(x*y - 1), (-1)/(x*y)])
+        sage: q == whole +sum(parts)
+        True
+
+    ::
+
+        sage: R.<x,y>= PolynomialRing(QQ)
+        sage: q= x +1/x +1/(x*y-2) +1/(x^2+y^2-1)
+        sage: whole,parts= partial_fraction_decomposition(q)
+        sage: whole, parts
+        (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)])
+        sage: whole +sum(parts)==q
+        True
+
+    ::
+
+        sage: R.<x,y,z>= PolynomialRing(QQ)
+        sage: q= 1/x +1/(x*y-z-2)^2 +1/(x^2+y^2 +z^2-1)^3 +1/(x*y-3)
+        sage: whole,parts= partial_fraction_decomposition(q)
+        sage: whole +sum(parts)==q
+        True
+
+    NOTES:
+
+    In the case of univariate `f` this function calls the old univariate
+    partial fraction decomposition function.
+    In the multivariate case, it uses a different notion of and algorithm for
+    partial fraction decompositions.
+
+    Let `f= P/Q` where `P,Q \in R`, let `Q_1^{e_1} \cdots Q_m^{e_m}` be the
+    unique factorization of `Q` in `R` into irreducible factors, and let `d` be
+    the number of indeterminates of `R`.
+    Then `f` can be written as a sum `\sum P_A/\prod_{j\in A} Q_j^{b_j}`,
+    where the `b_j \le e_j` are positive integers, the `P_A` are in `R`, and
+    the sum is taken over all size `\le m` subsets `A \subseteq \{ 1, \ldots, d \}`
+    such that `S:= \{ Q_j : j\in A \}` is algebraically independent and the
+    ideal generated by `S` is not all of `R` (that is, by the Nullstellensatz,
+    the polynomials of `S` have no common roots in the algebraic closure of the
+    coefficient field of `R`).
+    Following [Lein1978]_, i call any such decomposition of `f` a
+    `\emph{partial fraction decomposition}`.
+
+    The algorithm used below comes from Theorem 1, Lemma 2, and Lemma 3 of
+    [Lein1978]_.
+    By the way, that paper has a typo in equation (c) on the
+    second page and equation (b) on the third page.
+    The right sides of (c) and (b) should be multiplied by `P`.
+
+    REFERENCES:
+
+    .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
+    several variables into partial fractions", Soviet Math. (Iz. VUZ)
+    22 (1978), no. 10, 35--38.
+
+    AUTHORS:
+
+    - Alex Raichev (2011-01-10)
+    """
+    R= f.denominator().parent()
+    d= len(R.gens())
+    if d==Integer(1) :
+        return f.partial_fraction_decomposition()
+    Q= f.denominator()
+    whole,P= f.numerator().quo_rem(Q)
+    parts= [format_quotient(Integer(1) /Q)]
+    # Decompose via nullstellensatz trick
+    # (which is faster than the algebraic dependence trick)
+    parts= decompose_via_nullstellensatz(parts)
+    # Decompose via algebraic dependence trick
+    parts= decompose_via_algebraic_dependence(parts)
+    # Rewrite parts back in terms of rational expressions
+    new_parts=[]
+    for p in parts:
+        f= unformat_quotient(p)
+        if f.denominator() == Integer(1) :
+            whole= whole +f
+        else:
+            new_parts.append(f)
+    # Put P back in
+    new_parts= [P*f for f in new_parts]
+    return whole, new_parts
+#-------------------------------------------------------------------------------
+def combine_like_terms(parts,rationomial=True):
+    r"""
+    Combines like terms of the fractions represented by parts.
+    For use by partial_fraction_decomposition() above.
+
+    EXAMPLES:
+
+    sage: from sage.combinat.amgf import *
+    sage: R.<x,y>= PolynomialRing(QQ)
+    sage: parts =[[1, [x*y, 1]], [x, [x*y, 1]]]
+    sage: combine_like_terms(parts)
+    [[x + 1, [y, 1], [x, 1]]]
+
+    ::
+
+    sage: R.<x>= PolynomialRing(QQ)
+    sage: parts =[[1, [x, 1]], [x-1, [x, 1]]]
+    sage: combine_like_terms(parts)
+    [[1]]
+
+    AUTHORS:
+
+    - Alex Raichev (2011-01-10)
+    """
+
+    if parts == []: return parts
+    # Sort parts by denominators
+    new_parts= sorted(parts, key= lambda p:p[Integer(1) :])
+    # Combine items of parts with same denominators.
+    newnew_parts=[]
+    left,right=Integer(0) ,Integer(1)
+    glom= new_parts[left][Integer(0) ]
+    while right <= len(new_parts)-Integer(1) :
+        if new_parts[left][Integer(1) :] == new_parts[right][Integer(1) :]:
+            glom= glom +new_parts[right][Integer(0) ]
+        else:
+            newnew_parts.append([glom]+new_parts[left][Integer(1) :])
+            left= right
+            glom= new_parts[left][Integer(0) ]
+        right= right +Integer(1)
+    if glom != Integer(0) :
+        newnew_parts.append([glom]+new_parts[left][Integer(1) :])
+    if rationomial:
+        # Reduce fractions in newnew_parts.
+        # Todo: speed up below by working directly with newnew_parts and
+        # thereby make fewer calls to format_quotient() which in turn
+        # calls factor().
+        newnew_parts= [format_quotient(unformat_quotient(part)) for part in newnew_parts]
+    return newnew_parts
+#-------------------------------------------------------------------------------
+def decompose_via_algebraic_dependence(parts):
+    r"""
+    Returns a decomposition of parts.
+    Used by partial_fraction_decomposition() above.
+    Implements Lemma 2 of [Lein1978]_.
+    Recursive.
+
+    REFERENCES:
+
+    .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
+    several variables into partial fractions", Soviet Math. (Iz. VUZ)
+    22 (1978), no. 10, 35--38.
+
+    AUTHORS:
+
+    - Alex Raichev (2011-01-10)
+    """
+    decomposing_done= True
+    new_parts= []
+    for p in parts:
+        p_parts= [p]
+        P= p[Integer(0) ]
+        Qs= p[Integer(1) :]
+        m= len(Qs)
+        G= algebraic_dependence([q for q,e in Qs])
+        if G:
+            # Then the denominator factors are algebraically dependent
+            # and so we can decompose p.
+            decomposing_done= False
+            P= p[Integer(0) ]
+            Qs= p[Integer(1) :]
+
+            # Todo: speed up step below by using
+            # G to calculate F.  See [Lein1978]_ Lemma 1.
+            F= algebraic_dependence([q**e for q,e in Qs])
+            new_vars= F.parent().gens()
+
+            # Note that new_vars[j] corresponds to Qs[j] so that
+            # F([q^e for q,e in Qs]) = 0.
+            # Assuming below that F.parent() has negdeglex term order
+            # so that F.lt() is indeed the monomial we want.
+            FF= (F.lt() -F)/(F.lc())
+            numers= map(mul,zip(FF.coefficients(),FF.monomials()))
+            e= list(F.lt().exponents())[Integer(0) :m]
+            denom= [[new_vars[j], e[Integer(0) ][j]+Integer(1) ] for j in range(m)]
+
+            # Before making things messy by substituting in Qs,
+            # reduce terms and combine like terms.
+            p_parts_temp= [format_quotient(unformat_quotient([a]+denom)) for a in numers]
+            p_parts_temp= combine_like_terms(p_parts_temp)
+
+            # Substitute Qs into new_p.
+            Qpowsub= dict([(new_vars[j],Qs[j][Integer(0) ]**Qs[j][Integer(1) ]) for j in range(m)])
+            p_parts=[]
+            for x in p_parts_temp:
+                y= P*F.parent()(x[Integer(0) ]).subs(Qpowsub)
+                yy=[]
+                for xx in x[Integer(1) :]:
+                    if xx[Integer(0) ] in new_vars:
+                        j= new_vars.index(xx[Integer(0) ])
+                        yy.append([Qs[j][Integer(0) ],Qs[j][Integer(1) ]*xx[Integer(1) ]])
+                    else:
+                        # Occasionally xx[0] is an integer.
+                        yy.append(xx)
+                p_parts.append([y]+yy)
+        # Done one step of decomposing p.  Add it to new_parts.
+        new_parts.extend(p_parts)
+    if decomposing_done:
+        return new_parts
+    else:
+        return decompose_via_algebraic_dependence(new_parts)
+#-------------------------------------------------------------------------------
+def decompose_via_cohomology(parts):
+    r"""
+    Given each nice (described below) differential form
+    `(P/Q) dx_1 \wedge\cdots\wedge dx_d` in `parts`,
+    this function returns a differential form equivalent in De Rham
+    cohomology that has no repeated factors in the denominator.
+
+    INPUT:
+
+    - ``parts`` - A list of the form `[chunk_1,\ldots,chunk_r]`, where each
+        `chunk_j` has the form `[P,[Q_1,e_1],\ldots,[Q_m,e_m]]`,
+        `Q_1,\ldots,Q_m` are irreducible elements of a common polynomial
+        ring `R` such that their corresponding algebraic varieties
+        `\{x\in F^d : B_j(x)=0\}` intersect transversely (where `F` is the
+        algebraic closure of the field of coefficients of `R`),
+        `e_1,\ldots,e_m` are positive integers, `m \le d`, and
+        `P` is a symbolic expression in some of the indeterminates of `R`.
+        Here `[P,[Q_1,e_1],\ldots,[Q_m,e_m]]` represents the fraction
+        `P/(Q_1^e_1 \cdots Q_m^e_m)`.
+
+    OUTPUT:
+
+    A list of the form `[chunky_1,\ldots,chunky_s]`, where each
+    `chunky_j` has the form `[P,[Q_1,1],\ldots,[Q_m,1]]`.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.amgf import *
+        sage: R.<x,y>= PolynomialRing(QQ)
+        sage: decompose_via_cohomology([[ 1, [x,3] ]])
+        []
+
+    ::
+
+        sage: R.<x,y>= PolynomialRing(QQ)
+        sage: decompose_via_cohomology([[ 1, [x,3], [x*y-1,2] ]])
+        [[-3*y^2, [x, 1], [x*y - 1, 1]], [-5*y^3, [x*y - 1, 1]]]
+
+    NOTES:
+
+    This is a recursive function thats stops calling itself when all the
+    `e_j` equal 1 or `parts == []`.
+    The algorithm used here is that of Theorem 17.4 of
+    [AiYu1983]_.
+    The algebraic varieties `\{x\in F^d : Q_j(x)=0\}`
+    (where `F` is the algebraic closure of the field of coefficients of `R`)
+    corresponding to the `Q_j` __intersect transversely__ iff for each
+    point `c` of their intersection and for all `k \le m`,
+    the Jacobian matrix of any `k` polynomials from
+    `\{Q_1,\ldots,Q_m\}` has rank equal to `\min\{k,d\}` when evaluated at
+    `c`.
+
+    REFERENCES:
+
+    .. [AiYu1983]  I. A. A\u\izenberg and A. P. Yuzhakov, "Integral
+    representations and residues in multidimensional complex analysis",
+    Translations of Mathematical Monographs, 58. American Mathematical
+    Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X.
+
+    AUTHORS:
+
+    - Alex Raichev (2008-10-01, 2011-01-15)
+    """
+    if parts == []: return parts
+    decomposing_done= True
+    new_parts= []
+    R= parts[Integer(0) ][Integer(1) ][Integer(0) ].parent()
+    V= list(R.gens())
+    for p in parts:
+        p_parts= [p]
+        P= p[Integer(0) ]
+        Qs= p[Integer(1) :]
+        m= len(Qs)
+        if sum([e for q,e in Qs]) > m:
+            # Then we can decompose p
+            p_parts= []
+            decomposing_done=False
+            dets= []
+            vars_subsets= Set(V).subsets(m)
+            for v in vars_subsets:
+                # Sort variables so that first polynomial ring indeterminate
+                # declared is first in vars list.
+                v= sorted(v,reverse=True)
+                jac= jacobian([q for q,e in Qs],v)
+                dets.append(jac.determinant())
+            # Get a Nullstellensatz certificate.
+            L= R(Integer(1) ).lift(R.ideal([q for q,e in Qs] +dets))
+            for j in range(m):
+                if L[j] != Integer(0) :
+                    # Make a copy of (and not a reference to) the nested list Qs.
+                    new_Qs = deepcopy(Qs)
+                    if new_Qs[j][Integer(1) ] > Integer(1) :
+                        new_Qs[j][Integer(1) ]= new_Qs[j][Integer(1) ] -Integer(1)
+                    else:
+                        del new_Qs[j]
+                    p_parts.append( [P*L[j]] +new_Qs )
+            for k in range(vars_subsets.cardinality()):
+                if L[m+k] != Integer(0) :
+                    new_Qs = deepcopy(Qs)
+                    for j in range(m):
+                        if new_Qs[j][Integer(1) ] > Integer(1) :
+                            new_Qs[j][Integer(1) ]= new_Qs[j][Integer(1) ] -Integer(1)
+                            v=sorted(vars_subsets[k],reverse=True)
+                            jac= jacobian([SR(P*L[m+k])] +[ SR(Qs[jj][Integer(0) ]) for \
+                            jj in range(m) if jj !=j], [SR(vv) for vv in v])
+                            det= jac.determinant()
+                            if det != Integer(0) :
+                                p_parts.append([permutation_sign(v,V) \
+                                  *(-Integer(1) )**j/new_Qs[j][Integer(1) ] *det] +new_Qs)
+                            break
+        # Done one step of decomposing p.  Add it to new_parts.
+        new_parts.extend(p_parts)
+        new_parts= combine_like_terms(new_parts,rationomial=False)
+    if decomposing_done:
+        return new_parts
+    else:
+        return decompose_via_cohomology(new_parts)
+#-------------------------------------------------------------------------------
+def decompose_via_nullstellensatz(parts):
+    r"""
+    Returns a decomposition of parts.
+    Used by partial_fraction_decomposition() above.
+    Implements Lemma 3 of [Lein1978]_.
+    Recursive.
+
+    REFERENCES:
+
+    .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of
+    several variables into partial fractions", Soviet Math. (Iz. VUZ)
+    22 (1978), no. 10, 35--38.
+
+    AUTHORS:
+
+    - Alex Raichev (2011-01-10)
+    """
+    decomposing_done= True
+    new_parts= []
+    R= parts[Integer(0) ][Integer(0) ].parent()
+    for p in parts:
+        p_parts= [p]
+        P= p[Integer(0) ]
+        Qs= p[Integer(1) :]
+        m= len(Qs)
+        if R(Integer(1) ) in R.ideal([q for q,e in Qs]):
+            # Then we can decompose p.
+            decomposing_done= False
+            L= R(Integer(1) ).lift(R.ideal([q**e for q,e in Qs]))
+            p_parts= [ [P*L[i]] + \
+            [[Qs[j][Integer(0) ],Qs[j][Integer(1) ]] for j in range(m) if j != i] \
+            for i in range(m) if L[i]!=Integer(0) ]
+            # The procedure above yields no like terms to combine.
+        # Done one step of decomposing p.  Add it to new_parts.
+        new_parts.extend(p_parts)
+    if decomposing_done:
+        return new_parts
+    else:
+        return decompose_via_nullstellensatz(new_parts)
+#-------------------------------------------------------------------------------
+def format_quotient(f):
+    r"""
+    Formats `f` for use by partial_fraction_decomposition() above.
+
+    AUTHORS:
+
+    - Alex Raichev (2011-01-10)
+    """
+    P= f.numerator()
+    Q= f.denominator()
+    Qs= Q.factor()
+    if Qs.unit() != Integer(1) :
+        P= P/Qs.unit()
+    Qs= sorted([[q,e] for q,e in Qs])  # sorting for future bookkeeping
+    return [P]+Qs
+#-------------------------------------------------------------------------------
+def permutation_sign(v,vars):
+    r"""
+    This function returns the sign of the permutation on `1,\ldots,len(vars)`
+    that is induced by the sublist `v` of `vars`.
+    For internal use by _cohom_equiv_main().
+
+    INPUT:
+
+    - ``v`` - A sublist of `vars`.
+    - ``vars`` - A list of predefined variables or numbers.
+
+    OUTPUT:
+
+    The sign of the permutation obtained by taking indices (and adding 1)
+    within `vars` of the list `v,w`, where `w` is the list `vars` with the
+    elements of `v` removed.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.amgf import *
+        sage: vars= ['a','b','c','d','e']
+        sage: v= ['b','d']
+        sage: permutation_sign(v,vars)
+        -1
+        sage: v= ['d','b']
+        sage: permutation_sign(v,vars)
+        1
+
+    AUTHORS:
+
+    - Alex Raichev (2008-10-01)
+    """
+    # Convert variable lists to lists of numbers in {1,...,len(vars)}
+    A= [x+Integer(1)  for x in range(len(vars))]
+    B= [vars.index(x)+Integer(1)  for x in v]
+    C= list(Set(A).difference(Set(B)))
+    C.sort()
+    P= Permutation(B+C)
+    return P.signature()
+#-------------------------------------------------------------------------------
+def unformat_quotient(part):
+    r"""
+    Unformats `f` for use by partial_fraction_decomposition() above.
+    Inverse of format_quotient() above.
+
+    AUTHORS:
+
+    - Alex Raichev (2011-01-10)
+    """
+    P= part[Integer(0) ]
+    Qs= part[Integer(1) :]
+    Q= prod([q**e for q,e in Qs])
+    return P/Q
+#===============================================================================
+# Class for calculation of asymptotics of multivariate generating functions.
+class QuasiRationalExpression(SageObject):
+    "Store an expression G/H, where H comes from a polynomial ring R and \
+    G comes from R or the Symbolic Ring."
+    def __init__(self, G, H):
+        # Store important information about object as attributes of self.
+        # G, H, H's ring, ring dimension, H's factorization.
+        self._G = G
+        self._H = H
+        R= H.parent()
+        self._R = R
+        self._d = len(R.gens())
+        self._Hfac= list(H.factor())
+
+        # Variables of self as elements of the SR.
+        # Remember that G might be in SR and not in R.
+        try:
+            # This fails if G is a constant, for example.
+            Gv= Set([R(x) for x in G.variables()])
+        except:
+            Gv= Set([])
+        try:
+            # This fails if H is a constant, for example.
+            Hv= Set(H.variables())
+        except:
+            Hv= Set([])
+        # Preserve the ring ordering of the variables which some methods below
+        # depends upon.
+        V= sorted(list(Gv.union(Hv)),reverse=True)
+        self._variables = tuple([SR(x) for x in V])
+#-------------------------------------------------------------------------------
+# Keeping methods in alphabetical order (ignoring initial single underscores)
+    def asymptotics(self,c,alpha,N,numerical=Integer(0) ,asy_var=None):
+        r"""
+        This function returns the first `N` terms of the asymptotic expansion
+        of the Maclaurin coefficients `F_{n\alpha}` of the
+        multivariate meromorphic function `F=G/H` as `n\to\infty`,
+        where `F = self`.
+        It assumes that `F` is holomorphic in a neighborhood of the origin,
+        that `H` is a polynomial, and that `c` is a convenient strictly minimal
+        smooth or multiple of `F` that is critical for `\alpha`.
+
+        INPUT:
+
+        - ``alpha`` - A `d`-tuple of positive integers or, if `c` is a smooth
+            point, possibly of symbolic variables.
+        - ``c`` - A dictionary with keys `self._variables` and values from a
+            superfield of the field of `self._R.base_ring()`.
+        - ``N`` - A positive integer.
+        - ``numerical`` - A natural number (default: 0).
+            If k=numerical > 0, then a numerical approximation of the coefficients
+            of `F_{n\alpha}` with k digits of precision will be returned.
+            Otherwise exact values will be returned.
+        - ``asy_var`` - A  symbolic variable (default: None).
+            The variable of the asymptotic expansion.
+            If none is given, `var('n')` will be assigned.
+
+        OUTPUT:
+
+        The tuple `(asy,exp_part,subexp_part)`, where `asy` is first `N` terms
+        of the asymptotic expansion of the Maclaurin coefficients `F_{n\alpha}`
+        of the function `F=self` as `n\to\infty`, `exp_part^n` is the exponential
+        factor of `asy`, and `subexp_part` is the subexponential factor of
+        `asy`.
+
+        EXAMPLES::
+
+        A smooth point example ::
+
+            sage: R.<x,y>= PolynomialRing(QQ)
+            sage: G= 1
+            sage: H= 1-x-y-x*y
+            sage: F= QuasiRationalExpression(G,H)
+            sage: alpha= [3,2]
+            sage: c= {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3}
+            sage: F.asymptotics(c,alpha,2,numerical=3)
+            Initializing auxiliary functions...
+            Calculating derivatives of auxiallary functions...
+            Calculating derivatives of more auxiliary functions...
+            Calculating actions of the second order differential operator...
+            ((0.369/sqrt(n) - 0.0186/n^(3/2))*71.2^n, 71.2, 0.369/sqrt(n) - 0.0186/n^(3/2))
+
+        A multiple point example ::
+
+            sage: R.<x,y,z>= PolynomialRing(QQ)
+            sage: G= 1
+            sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y))
+            sage: F= QuasiRationalExpression(G,H)
+            sage: c= {z: 4/3, y: 1, x: 1/2}
+            sage: alpha= [8,3,3]
+            sage: F.asymptotics(c,alpha,1)
+            Initializing auxiliary functions...
+            Calculating derivatives of auxiliary functions...
+            Calculating derivatives of more auxiliary functions...
+            Calculating second-order differential operator actions...
+            (1/7*sqrt(3)*sqrt(7)*108^n/(sqrt(pi)*sqrt(n)), 108, 1/7*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n)))
+
+        NOTES:
+
+        A zero `c` of `H` is __strictly minimal__ if there is no zero `x` of `H`
+        such that `x_j < c_j` for all `0 \le j < d`.
+        For definitions of the terms "smooth critical point for `\alpha`" and
+        "multiple critical point for `\alpha`",
+        see the documentation for _asymptotics_main_smooth() and
+        _asymptotics_main_multiple(), which are the functions that do most of the
+        work.
+
+        ALGORITHM:
+
+        The algorithm used here comes from [RaWi2011]_.
+
+        (1) Use Cauchy's multivariate integral formula to write `F_{n\alpha}` as
+        an integral around a polycirle centered at the origin of the
+        differential form `\frac{G(x) dx[0] \wedge\cdots\wedge
+        dx[d-1]}{H(x)x^\alpha}`.
+
+        (2) Decompose `G/H` into a sum of partial fractions `P[0] +\cdots+ P[r]`
+        so that each term of the sum has at most `d` irreducible factors of `H`
+        in the denominator.
+
+        (3) For each differential form `P[j] dx[0] \wedge\cdots\wedge dx[d-1]`,
+        find an equivalent form `\omega[j]` in de Rham cohomology with no
+        repeated irreducible factors of `H` in its denominator.
+
+        (4) Compute an asymptotic expansion for each integral `\omega[j]`.
+
+        (5) Add the expansions.
+
+        REFERENCES:
+
+        .. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
+        coefficients of multivariate generating functions: improvements
+        for smooth points", Electronic Journal of Combinatorics, Vol. 15
+        (2008), R89.
+
+        .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
+        coefficients of multivariate generating functions: improvements
+        for smooth points", To appear.
+
+        AUTHORS:
+
+        - Alex Raichev (2008-10-01, 2010-09-28, 2011-04-27)
+        """
+        # The variable for asymptotic expansions.
+        if not asy_var:
+            asy_var= var('n')
+
+        # Create symbolic (non-ring) variables.
+        R= self._R
+        d= self._d
+        X= list(self._variables)
+
+        # Do steps (1)--(3)
+        new_integrands= self.cohomologous_integrand(alpha,asy_var)
+
+        # Coerce everything into the Symbolic Ring, as polynomial ring
+        # calculations are no longer needed.
+        # Calculate asymptotics.
+        cc={}
+        for k in c.keys():
+            cc[SR(k)] = SR(c[k])
+        c= cc
+        for i in range(len(alpha)):
+            alpha[i] = SR(alpha[i])
+        subexp_parts= []
+        for chunk in new_integrands:
+            # Convert chunk into Symbolic Ring
+            GG= SR(chunk[Integer(0) ])
+            HH= [SR(f) for (f,e) in chunk[Integer(1) :]]
+            asy= self._asymptotics_main(GG,HH,X,c,asy_var,alpha,N,numerical)
+            subexp_parts.append(asy[Integer(2) ])
+        exp_scale= asy[Integer(1) ]  # same for all chunk in new_integrands
+
+        # Do step (5).
+        subexp_part= add(subexp_parts)
+        return exp_scale**asy_var *subexp_part, exp_scale, subexp_part
+#--------------------------------------------------------------------------------
+    def _asymptotics_main(self,G,H,X,c,n,alpha,N,numerical):
+        r"""
+        This function is for internal use by asymptotics().
+        It finds a variable in `X` to use to calculate asymptotics and decides
+        whether to call _asymptotics_main_smooth() or
+        _asymptotics_main_multiple().
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``G`` - A symbolic expression.
+        - ``H`` - A list of symbolic expressions.
+        - ``X`` - The list of variables occurring in `G` and `H`.
+            Call its length `d`.
+        - ``c`` - A dictionary with `X` as keys and numbers as values.
+        - ``n`` - The variable of the asymptotic expansion.
+        - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic
+            variables if `c` is a smooth point.
+        - ``N`` - A positive integer.
+        - ``numerical`` - Natural number.
+            If k=numerical > 0, then a numerical approximation of the coefficients
+            of `F_{n\alpha}` with k digits of precision will be returned.
+            Otherwise exact values will be returned.
+
+        OUTPUT:
+
+        The same as the function asymptotics().
+
+        AUTHORS:
+
+        - Alex Raichev (2008-10-01, 2010-09-28)
+        """
+        d= len(X)
+        r= len(H)  # We know 1 <= r <= d.
+
+        # Find a good variable x in X to do asymptotics calculations with, that is,
+        # a variable x in X such that for all h in H, diff(h,x).subs(c) != 0.
+        # A good variable is guaranteed to exist since we are dealing with
+        # convenient smooth or multiple points.
+        # Search for good x in X from back to front (to be consistent with
+        # [RaWi2008a]_ [RaWi2011]_ which uses X[d-1] as a good variable).
+        # Put each not good x found at the beginning of X and reshuffle alpha
+        # in the same way.
+        x= X[d-Integer(1) ]
+        beta= copy(alpha)
+        while Integer(0)  in [diff(h,x).subs(c) for h in H]:
+            X.pop()
+            X.insert(Integer(0) ,x)
+            x= X[d-Integer(1) ]
+            a= beta.pop()
+            beta.insert(Integer(0) ,a)
+        if d==r:
+            # This is the case of a 'simple' multiple point.
+            A= G.subs(c) / jacobian(H,X).subs(c).determinant().abs()
+            return A,Integer(1) ,A
+        elif r==Integer(1) :
+            # So 1 = r < d, and we have a smooth point.
+            return self._asymptotics_main_smooth(G,H[Integer(0) ],X,c,n,beta,N,numerical)
+        else:
+            # So 1 < r < d, and we have a non-smooth multiple point.
+            return self._asymptotics_main_multiple(G,H,X,c,n,beta,N,numerical)
+#--------------------------------------------------------------------------------
+    def _asymptotics_main_multiple(self,G,H,X,c,n,alpha,N,numerical):
+        r"""
+        This function is for internal use by _asymptotics_main().
+        It calculates asymptotics in case they depend upon multiple points.
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``G`` - A symbolic expression.
+        - ``H`` - A list of symbolic expressions.
+        - ``X`` - The list of variables occurring in `G` and `H`.
+            Call its length `d`.
+        - ``c`` - A dictionary with `X` as keys and numbers as values.
+        - ``n`` - The variable of the asymptotic expansion.
+        - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic
+            variables if `c` is a smooth point.
+        - ``N`` - A positive integer.
+        - ``numerical`` - Natural number.
+            If k=numerical > 0, then a numerical approximation of the coefficients
+            of `F_{n\alpha}` with k digits of precision will be returned.
+            Otherwise exact values will be returned.
+
+        OUTPUT:
+
+        The same as the function asymptotics().
+
+        NOTES:
+
+        The formula used for computing the asymptotic expansion is
+        Theorem 3.4 of [RaWi2011]_.
+
+        Currently this function cannot handle `c` with symbolic variable keys,
+        because _crit_cone_combo() crashes.
+
+        REFERENCES:
+
+        .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
+        coefficients of multivariate generating functions: improvements
+        for smooth points", To appear.
+
+        AUTHORS:
+
+        - Alex Raichev (2008-10-01, 2010-09-28)
+        """
+        I= sqrt(-Integer(1) )
+        d= len(X)   # > r > 1
+        r= len(H)   # > 1
+        C= copy(c)
+
+        S= [var(self._new_var_name('s',X) + str(j)) for j in range(r-Integer(1) )]
+        T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-Integer(1) )]
+        Sstar= {}
+        temp= self._crit_cone_combo(H,X,c,alpha)
+        for j in range(r-Integer(1) ):
+            Sstar[S[j]]= temp[j]
+        thetastar= dict([(t,Integer(0) ) for t in T])
+        thetastar.update(Sstar)
+
+        # Setup.
+        print "Initializing auxiliary functions..."
+        Hmul= mul(H)
+        h= [function('h'+str(j),*tuple(X[:d-Integer(1) ])) for j in range(r)]    # Implicit functions
+        U = function('U',*tuple(X))
+        # All other functions are defined in terms of h, U, and explicit functions.
+        Hcheck=  mul([ X[d-Integer(1) ] -Integer(1) /h[j] for j in range(r)] )
+        Gcheck= -G/U *mul( [-h[j]/X[d-Integer(1) ] for j in range(r)] )
+        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)]
+        e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-Integer(1) )])
+        ht= [hh.subs(e) for hh in h]
+        Ht= [H[j].subs(e).subs({X[d-Integer(1) ]:Integer(1) /ht[j]}) for j in range(r)]
+        hsumt= add([S[j]*ht[j] for j in range(r-Integer(1) )]) +(Integer(1) -add(S))*ht[r-Integer(1) ]
+        At= [AA.subs(e).subs({X[d-Integer(1) ]:hsumt}) for AA in A]
+        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) )])
+        # atC Stores h and U and all their derivatives evaluated at C.
+        atC = copy(C)
+        atC.update(dict( [(hh.subs(C),Integer(1) /C[X[d-Integer(1) ]]) for hh in h ]))
+
+        # Compute the derivatives of h up to order 2*N and evaluate at C.
+        hderivs1= {}   # First derivatives of h.
+        for (i,j) in mrange([d-Integer(1) ,r]):
+            s= solve( diff(H[j].subs({X[d-Integer(1) ]:Integer(1) /h[j]}),X[i]), diff(h[j],X[i]) )[Integer(0) ].rhs()\
+            .simplify()
+            hderivs1.update({diff(h[j],X[i]):s})
+            atC.update({diff(h[j],X[i]).subs(C):s.subs(C).subs(atC)})
+        hderivs = self._diff_all(h,X[Integer(0) :d-Integer(1) ],Integer(2) *N,sub=hderivs1,rekey=h)
+        for k in hderivs.keys():
+            atC.update({k.subs(C):hderivs[k].subs(atC)})
+
+        # Compute the derivatives of U up to order 2*N-2+ min{r,N}-1 and evaluate at C.
+        # To do this, differentiate H = U*Hcheck over and over, evaluate at C,
+        # and solve for the derivatives of U at C.
+        # Need the derivatives of H with short keys to pass on to diff_prod later.
+        m= min(r,N)
+        end= [X[d-Integer(1) ] for j in range(r)]
+        Hmulderivs= self._diff_all(Hmul,X,Integer(2) *N-Integer(2) +r,ending=end,sub_final=C)
+        print "Calculating derivatives of auxiliary functions..."
+        atC.update({U.subs(C):diff(Hmul,X[d-Integer(1) ],r).subs(C)/factorial(r)})
+        Uderivs={}
+        p= Hmul.polynomial(CC).degree()-r
+        if p == Integer(0) :
+            # Then, using a proposition at the end of [RaWi2011], we can
+            # conclude that all higher derivatives of U are zero.
+            for l in (ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )):
+                for s in UnorderedTuples(X,l):
+                    Uderivs[diff(U,s).subs(C)] = Integer(0)
+        elif p > Integer(0)  and p < Integer(2) *N-Integer(2) +m-Integer(1) :
+            all_zero= True
+            Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,p)),end,Uderivs,atC)
+            # Check for a nonzero U derivative.
+            if Uderivs.values() != [Integer(0)  for i in range(len(Uderivs))]:
+                all_zero= False
+            if all_zero:
+                # Then, using a proposition at the end of [RaWi2011], we can
+                # conclude that all higher derivatives of U are zero.
+                for l in (ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )):
+                    for s in UnorderedTuples(X,l):
+                        Uderivs.update({diff(U,s).subs(C):Integer(0) })
+            else:
+                # Have to compute the rest of the derivatives.
+                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)
+        else:
+            Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )),end,Uderivs,atC)
+        atC.update(Uderivs)
+        Phit1= jacobian(Phit,T+S).subs(hderivs1)
+        a= jacobian(Phit1,T+S).subs(hderivs1).subs(thetastar).subs(atC)
+        a_inv= a.inverse()
+        Phitu= Phit -(Integer(1) /Integer(2) ) *matrix([T+S]) *a *matrix([T+S]).transpose()
+        Phitu= Phitu[Integer(0) ][Integer(0) ]
+
+        # Compute all partial derivatives of At and Phitu up to orders 2*N-2
+        # and 2*N, respectively.  Take advantage of the fact that At and Phitu
+        # are sufficiently differentiable functions so that mixed partials
+        # are equal.  Thus only need to compute representative partials.
+        # Choose nondecreasing sequences as representative differentiation-
+        # order sequences.
+        # To speed up later computations, create symbolic functions AA and BB
+        # to stand in for the expressions At and Phitu respectively.
+        print "Calculating derivatives of more auxiliary functions..."
+        AA= [function('A'+str(j),*tuple(T+S)) for j in range(r)]
+        At_derivs= self._diff_all(At,T+S,Integer(2) *N-Integer(2) ,sub=hderivs1,sub_final=[thetastar,atC],rekey=AA)
+        BB= function('BB',*tuple(T+S))
+        Phitu_derivs= self._diff_all(Phitu,T+S,Integer(2) *N,sub=hderivs1,sub_final=[thetastar,atC],rekey=BB,zero_order=Integer(3) )
+        AABB_derivs= At_derivs
+        AABB_derivs.update(Phitu_derivs)
+        for j in range(r):
+            AABB_derivs[AA[j]] = At[j].subs(thetastar).subs(atC)
+        AABB_derivs[BB] = Phitu.subs(thetastar).subs(atC)
+        print "Calculating second-order differential operator actions..."
+        DD= self._diff_op(AA,BB,AABB_derivs,T+S,a_inv,r,N)
+
+        L={}
+        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) ))):
+            if j+k <= N-Integer(1) :
+                L[(j,k)] = add([ \
+                DD[(j,k,l)] /( (-Integer(1) )**k *Integer(2) **(k+l) *factorial(l) *factorial(k+l) ) \
+                for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k))] )
+        # The next line's QQ coercion is a workaround for the Sage 4.6 bug reported
+        # on http://trac.sagemath.org/sage_trac/ticket/8659.
+        # Once the bug is fixed, the QQ can be removed.
+        det= QQ(a.determinant())**(-Integer(1) /Integer(2) ) *(Integer(2) *pi)**((r-d)/Integer(2) )
+        chunk= det *add([ (alpha[d-Integer(1) ]*n)**((r-d)/Integer(2) -q) *add([ \
+        L[(j,k)] *binomial(r-Integer(1) ,j) *stirling_number1(r-j,r+k-q) *(-Integer(1) )**(q-j-k) \
+        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 ]) \
+        for q in range(N)])
+
+        chunk= chunk.subs(C).simplify()
+        coeffs= chunk.coefficients(n)
+        coeffs.reverse()
+        coeffs= coeffs[:N]
+        if numerical: # If a numerical approximation is desired...
+            subexp_part = add( [co[Integer(0) ].subs(c).n(digits=numerical)*n**co[Integer(1) ] \
+            for co in coeffs] )
+            exp_scale= (Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )) \
+            .n(digits=numerical)
+        else:
+            subexp_part = add( [co[Integer(0) ].subs(c)*n**co[Integer(1) ] for co in coeffs] )
+            exp_scale= Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )
+        return exp_scale**n*subexp_part, exp_scale, subexp_part
+#--------------------------------------------------------------------------------
+    def _asymptotics_main_smooth(self,G,H,X,c,n,alpha,N,numerical):
+        r"""
+        This function is for internal use by _asymptotics_main().
+        It calculates asymptotics in case they depend upon smooth points.
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``G`` - A symbolic expression.
+        - ``H`` - A list of symbolic expressions.
+        - ``X`` - The list of variables occurring in `G` and `H`.
+            Call its length `d`.
+        - ``c`` - A dictionary with `X` as keys and numbers as values.
+        - ``n`` - The variable of the asymptotic expansion.
+        - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic
+            variables if `c` is a smooth point.
+        - ``N`` - A positive integer.
+        - ``numerical`` - Natural number.
+            If k=numerical > 0, then a numerical approximation of the coefficients
+            of `F_{n\alpha}` with k digits of precision will be returned.
+            Otherwise exact values will be returned.
+
+        OUTPUT:
+
+        The same as the function asymptotics().
+
+        NOTES:
+
+        The formulas used for computing the asymptotic expansions are
+        Theorems 3.2 and 3.3 [RaWi2008a]_ with `p=1`.
+        Theorem 3.2 is a specialization of Theorem 3.4 of [RaWi2011]_
+        with `r=1`.
+
+        REFERENCES:
+
+        .. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
+        coefficients of multivariate generating functions: improvements
+        for smooth points", Electronic Journal of Combinatorics, Vol. 15
+        (2008), R89.
+
+        .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
+        coefficients of multivariate generating functions: improvements
+        for smooth points", To appear.
+
+        AUTHORS:
+
+        - Alex Raichev (2008-10-01, 2010-09-28)
+        """
+        I= sqrt(-Integer(1) )
+        d= len(X) # > 1
+
+        # If c is a tuple of rationals, then compute with it directly.
+        # Otherwise, compute symbolically and plug in c at the end.
+        if vector(c.values()) in QQ**d:
+            C= c
+        else:
+            Cs= [var('cs' +str(j)) for j in range(d)]
+            C= dict( [(X[j],Cs[j]) for j in range(d)] )
+            c= dict( [(Cs[j],c[X[j]]) for j in range(d)] )
+
+        # Setup.
+        print "Initializing auxiliary functions..."
+        h= function('h',*tuple(X[:d-Integer(1) ]))    # Implicit functions
+        U = function('U',*tuple(X))         #
+        # All other functions are defined in terms of h, U, and explicit functions.
+        Gcheck = -G/U *(h/X[d-Integer(1) ])
+        A= Gcheck.subs({X[d-Integer(1) ]:Integer(1) /h})/h
+        T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-Integer(1) )]
+        e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-Integer(1) )])
+        ht= h.subs(e)
+        Ht= H.subs(e).subs({X[d-Integer(1) ]:Integer(1) /ht})
+        At= A.subs(e)
+        Phit = -log(C[X[d-Integer(1) ]] *ht)\
+        + I* add([alpha[i]/alpha[d-Integer(1) ] *T[i] for i in range(d-Integer(1) )])
+        Tstar= dict([(t,Integer(0) ) for t in T])
+        atC = copy(C)
+        atC.update({h.subs(C):Integer(1) /C[X[d-Integer(1) ]]}) # Stores h and U and all their derivatives
+                                            # evaluated at C.
+
+        # Compute the derivatives of h up to order 2*N and evaluate at C and store
+        # in atC.  Keep a copy of unevaluated h derivatives for use in the case
+        # d=2 and v > 2 below.
+        hderivs1= {}   # First derivatives of h.
+        for i in range(d-Integer(1) ):
+            s= solve( diff(H.subs({X[d-Integer(1) ]:Integer(1) /h}),X[i]), diff(h,X[i]) )[Integer(0) ].rhs()\
+            .simplify()
+            hderivs1.update({diff(h,X[i]):s})
+            atC.update({diff(h,X[i]).subs(C):s.subs(C).subs(atC)})
+        hderivs = self._diff_all(h,X[Integer(0) :d-Integer(1) ],Integer(2) *N,sub=hderivs1,rekey=h)
+        for k in hderivs.keys():
+            atC.update({k.subs(C):hderivs[k].subs(atC)})
+
+        # Compute the derivatives of U up to order 2*N and evaluate at C.
+        # To do this, differentiate H = U*Hcheck over and over, evaluate at C,
+        # and solve for the derivatives of U at C.
+        # Need the derivatives of H with short keys to pass on to diff_prod later.
+        Hderivs= self._diff_all(H,X,Integer(2) *N,ending=[X[d-Integer(1) ]],sub_final=C)
+        print "Calculating derivatives of auxiallary functions..."
+        # For convenience in checking if all the nontrivial derivatives of U at c
+        # are zero a few line below, store the value of U(c) in atC instead of in
+        # Uderivs.
+        Uderivs={}
+        atC.update({U.subs(C):diff(H,X[d-Integer(1) ]).subs(C)})
+        end= [X[d-Integer(1) ]]
+        Hcheck= X[d-Integer(1) ] - Integer(1) /h
+        p= H.polynomial(CC).degree()-Integer(1)
+        if p == Integer(0) :
+            # Then, using a proposition at the end of [RaWi2011], we can
+            # conclude that all higher derivatives of U are zero.
+            for l in (ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N)):
+                for s in UnorderedTuples(X,l):
+                    Uderivs[diff(U,s).subs(C)] = Integer(0)
+        elif p > Integer(0)  and p < Integer(2) *N:
+            all_zero= True
+            Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,p)),end,Uderivs,atC)
+            # Check for a nonzero U derivative.
+            if Uderivs.values() != [Integer(0)  for i in range(len(Uderivs))]:
+                all_zero= False
+            if all_zero:
+                # Then, using a proposition at the end of [RaWi2011], we can
+                # conclude that all higher derivatives of U are zero.
+                for l in (ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N)):
+                    for s in UnorderedTuples(X,l):
+                        Uderivs.update({diff(U,s).subs(C):Integer(0) })
+            else:
+                # Have to compute the rest of the derivatives.
+                Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N)),end,Uderivs,atC)
+        else:
+            Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N)),end,Uderivs,atC)
+        atC.update(Uderivs)
+
+        # In general, this algorithm is not designed to handle the case of a
+        # singular Phit''(Tstar).  However, when d=2 the algorithm can cope.
+        if d==Integer(2) :
+            # Compute v, the order of vanishing at Tstar of Phit. It is at least 2.
+            v=Integer(2)
+            Phitderiv= diff(Phit,T[Integer(0) ],Integer(2) )
+            splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify()
+            while splat==Integer(0) :
+                v= v+Integer(1)
+                if v > Integer(2) *N:  # Then need to compute more derivatives of h for atC.
+                    hderivs.update({diff(h,X[Integer(0) ],v) \
+                    :diff(hderivs[diff(h,X[Integer(0) ],v-Integer(1) )],X[Integer(0) ]).subs(hderivs1)})
+                    atC.update({diff(h,X[Integer(0) ],v).subs(C) \
+                    :hderivs[diff(h,X[Integer(0) ],v)].subs(atC)})
+                Phitderiv= diff(Phitderiv,T[Integer(0) ])
+                splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify()
+        if d==Integer(2)  and v>Integer(2) :
+            t= T[Integer(0) ]  # Simplify variable names.
+            a= splat/factorial(v)
+            Phitu= Phit -a*t**v
+
+            # Compute all partial derivatives of At and Phitu up to orders 2*(N-1)
+            # and 2*(N-1)+v, respectively, in case v is even.
+            # Otherwise, compute up to orders N-1 and N-1+v, respectively.
+            # To speed up later computations, create symbolic functions AA and BB
+            # to stand in for the expressions At and Phitu, respectively.
+            print "Calculating derivatives of more auxiliary functions..."
+            AA= function('AA',t)
+            BB= function('BB',t)
+            if v.mod(Integer(2) )==Integer(0) :
+                At_derivs= self._diff_all(At,T,Integer(2) *N-Integer(2) , \
+                sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
+                Phitu_derivs= self._diff_all(Phitu,T,Integer(2) *N-Integer(2) +v, \
+                sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+Integer(1) ,rekey=BB)
+            else:
+                At_derivs= self._diff_all(At,T,N-Integer(1) ,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
+                Phitu_derivs= self._diff_all(Phitu,T,N-Integer(1) +v,sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+Integer(1) ,rekey=BB)
+            AABB_derivs= At_derivs
+            AABB_derivs.update(Phitu_derivs)
+            AABB_derivs[AA] = At.subs(Tstar).subs(atC)
+            AABB_derivs[BB] = Phitu.subs(Tstar).subs(atC)
+            print "Calculating actions of the second order differential operator..."
+            DD= self._diff_op_simple(AA,BB,AABB_derivs,t,v,a,N)
+            # Plug above into asymptotic formula.
+            L = []
+            if v.mod(Integer(2) ) == Integer(0) :
+                for k in range(N):
+                    L.append( add([ \
+                    (-Integer(1) )**l *gamma((Integer(2) *k+v*l+Integer(1) )/v) \
+                    / (factorial(l) *factorial(Integer(2) *k+v*l)) \
+                    * DD[(k,l)] for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)) ]) )
+                chunk= a**(-Integer(1) /v) /(pi*v) *add([ alpha[d-Integer(1) ]**(-(Integer(2) *k+Integer(1) )/v) \
+                * L[k] *n**(-(Integer(2) *k+Integer(1) )/v) for k in range(N) ])
+            else:
+                zeta= exp(I*pi/(Integer(2) *v))
+                for k in range(N):
+                    L.append( add([ \
+                    (-Integer(1) )**l *gamma((k+v*l+Integer(1) )/v) \
+                    / (factorial(l) *factorial(k+v*l)) \
+                    * (zeta**(k+v*l+Integer(1) ) +(-Integer(1) )**(k+v*l)*zeta**(-(k+v*l+Integer(1) ))) \
+                    * DD[(k,l)] for l in (ellipsis_range(Integer(0) ,Ellipsis,k)) ]) )
+                chunk= abs(a)**(-Integer(1) /v) /(Integer(2) *pi*v) *add([ alpha[d-Integer(1) ]**(-(k+Integer(1) )/v) \
+                * L[k] *n**(-(k+Integer(1) )/v) for k in range(N) ])
+        # Asymptotics for d>=2 case.  A singular Phit''(Tstar) will cause a crash
+        # in this case.
+        else:
+            Phit1= jacobian(Phit,T).subs(hderivs1)
+            a= jacobian(Phit1,T).subs(hderivs1).subs(Tstar).subs(atC)
+            a_inv= a.inverse()
+            Phitu= Phit -(Integer(1) /Integer(2) ) *matrix([T]) *a *matrix([T]).transpose()
+            Phitu= Phitu[Integer(0) ][Integer(0) ]
+            # Compute all partial derivatives of At and Phitu up to orders 2*N-2
+            # and 2*N, respectively.  Take advantage of the fact that At and Phitu
+            # are sufficiently differentiable functions so that mixed partials
+            # are equal.  Thus only need to compute representative partials.
+            # Choose nondecreasing sequences as representative differentiation-
+            # order sequences.
+            # To speed up later computations, create symbolic functions AA and BB
+            # to stand in for the expressions At and Phitu respectively.
+            print "Calculating derivatives of more auxiliary functions..."
+            AA= function('AA',*tuple(T))
+            At_derivs= self._diff_all(At,T,Integer(2) *N-Integer(2) ,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
+            BB= function('BB',*tuple(T))
+            Phitu_derivs= self._diff_all(Phitu,T,Integer(2) *N,sub=hderivs1,sub_final=[Tstar,atC],rekey=BB,zero_order=Integer(3) )
+            AABB_derivs= At_derivs
+            AABB_derivs.update(Phitu_derivs)
+            AABB_derivs[AA] = At.subs(Tstar).subs(atC)
+            AABB_derivs[BB] = Phitu.subs(Tstar).subs(atC)
+            print "Calculating actions of the second order differential operator..."
+            DD= self._diff_op(AA,BB,AABB_derivs,T,a_inv,Integer(1) ,N)
+
+            # Plug above into asymptotic formula.
+            L=[]
+            for k in range(N):
+                L.append( add([ \
+                  DD[(Integer(0) ,k,l)] / ( (-Integer(1) )**k *Integer(2) **(l+k) *factorial(l) *factorial(l+k) ) \
+                  for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k))]) )
+            chunk= add([ (Integer(2) *pi)**((Integer(1) -d)/Integer(2) ) *a.determinant()**(-Integer(1) /Integer(2) ) \
+                  *alpha[d-Integer(1) ]**((Integer(1) -d)/Integer(2)  -k) *L[k] \
+                  *n**((Integer(1) -d)/Integer(2) -k) for k in range(N) ])
+
+        chunk= chunk.subs(c).simplify()
+        coeffs= chunk.coefficients(n)
+        coeffs.reverse()
+        coeffs= coeffs[:N]
+        if numerical: # If a numerical approximation is desired...
+            subexp_part = add( [co[Integer(0) ].subs(c).n(digits=numerical)*n**co[Integer(1) ] for co in coeffs] )
+            exp_scale=(Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )) \
+            .n(digits=numerical)
+        else:
+            subexp_part = add( [co[Integer(0) ].subs(c)*n**co[Integer(1) ] for co in coeffs] )
+            exp_scale= Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )
+        return exp_scale**n*subexp_part, exp_scale, subexp_part
+#-----------------------------------------------------------------------------
+    def _crit_cone_combo(self,fs,X,c,alpha):
+        r"""
+        This function returns an auxiliary point associated to the multiple
+        point `c` of the factors `fs`.
+        It is for internal use by _asymptotics_main_multiple().
+
+        INPUT:
+
+        - ``fs`` - A list of expressions in the variables of `X`.
+        - ``X`` - A list of variables.
+        - ``c`` - A dictionary with keys `X` and values in some field.
+        - ``alpha`` - A list of rationals.
+
+        OUTPUT:
+
+        A solution of the matrix equation `y Gamma = a` for `y`,
+        where `Gamma` is the matrix whose `j`th row is
+        _direction(_log_grad(fj,X,c)) where `fj`
+        is the `j`th item in `fs` and where `a` is _direction(alpha).
+
+        EXAMPLES::
+
+            sage: R.<x,y,z>= PolynomialRing(QQ)
+            sage: G,H= 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: fs= [x + 2*y + z - 4, 2*x + y + z - 4]
+            sage: c= {x:1,y:1,z:1}
+            sage: alpha= [2,1,1]
+            sage: F._crit_cone_combo(fs,R.gens(),c,alpha)
+            [0, 1]
+
+        NOTES:
+
+        Use this function only when `Gamma` is well-defined and
+        there is a unique solution to the matrix equation `y Gamma = a`.
+        Fails otherwise.
+
+        AUTHORS:
+
+        - Alex Raichev (2008-10-01, 2008-11-25, 2009-03-04, 2010-09-08,
+            2010-12-02)
+        """
+        # Assuming here that each _log_grad(f) has nonzero final component.
+        # Then 'direction' will not throw a division by zero error.
+        d= len(X)
+        r= len(fs)
+        Gamma= matrix([self._direction(self._log_grad(f,X,c)) for f in fs])
+        # solve_left() fails when working in SR :-(. So use solve() instead.
+        #s= Gamma.solve_left(vector(alpha)/alpha[d-1])
+        V= [var('sss'+str(i)) for i in range(r)]
+        M= matrix(V)*Gamma
+        eqns= [M[Integer(0) ][i]== alpha[i]/alpha[d-Integer(1) ] for i in range(d)]
+        s= solve(eqns,V,solution_dict=True)[Integer(0) ]  # Assuming a unique solution.
+        return [s[v] for v in V]
+# B ==========================================================================
+# C ==========================================================================
+    def cohomologous_integrand(self,alpha,asy_var=None):
+        r"""
+        This function takes a multivariate Cauchy type integral
+        `\int F / x^{asy_var\alpha+1} dx`, where `F=self`, and breaks it up
+        into a list of nicer Cauchy type integrals for the purposes of
+        computing asymptotics of the original integral as `asy_var\to\infty`.
+        The sum of the nicer integrals is de Rham cohomologous to the original
+        integral.
+        It assumes that algebraic varieties corresponding to the irreducible
+        factors of `self._H` intersect transversely (see notes below).
+
+        INPUT:
+
+        - ``alpha`` - A list of positive integers or symbolic variables.
+        - ``asy_var`` - A symbolic variable (default: None).
+            Eventually set to `var('n')` if None is given.
+
+        OUTPUT:
+
+        A list of the form `[chunk_1,\ldots,chunk_r]`, where each
+        `chunk_j` has the form `[P,[B_1,1],\ldots,[B_l,1]]`.
+        Here `l \le d`, `P` is a symbolic expression in the indeterminates of
+        `R` and `asy_var`, `\{B_1,\ldots,B_l\} \subseteq \{Q_1,\ldots,Q_m\}`,
+        and `[P,[B_1,1],\ldots,[B_l,1]]` represents the integral
+        `\int P/(B_1 \cdots B_l) dx`.
+
+        EXAMPLES::
+
+        sage: R.<x,y>= PolynomialRing(QQ)
+        sage: G= 9*exp(x+y)
+        sage: H= (3-2*x-y)*(3-x-2*y)
+        sage: F= QuasiRationalExpression(G,H)
+        sage: alpha= [4,3]
+        sage: F.cohomologous_integrand(alpha)
+        [[9*e^(x + y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]]
+
+        sage: R.<x,y>= PolynomialRing(QQ)
+        sage: G= 9*exp(x+y)
+        sage: H= (3-2*x-y)^2*(3-x-2*y)
+        sage: F= QuasiRationalExpression(G,H)
+        sage: alpha= [4,3]
+        sage: F.cohomologous_integrand(alpha)
+        [[-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]]]
+
+        sage: R.<x,y,z>= PolynomialRing(QQ)
+        sage: G= 16
+        sage: H= (4-2*x-y-z)^2*(4-x-2*y-z)
+        sage: F= QuasiRationalExpression(G,H)
+        sage: alpha= [3,3,2]
+        sage: F.cohomologous_integrand(alpha)
+        [[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]]]
+
+        NOTES:
+
+        The varieties  corresponding to `Q_1,\ldots,Q_m`
+        __intersect transversely__ iff for each point `c` of their intersection
+        and for all `k \le l`, the Jacobian matrix of any `k` polynomials from
+        `\{Q_1,\ldots,Q_m\}` has rank equal to `\min\{k,d\}` when evaluated at
+        `c`.
+
+        ALGORITHM:
+
+        Let `asy_var= n` and consider the integral around a polycirle centered
+        at the origin of the `d`-variate differential form
+        `\frac{G(x) dx_1 \wedge\cdots\wedge dx_d}{H(x) x^{n\alpha+1}}`, where
+        `G=self._G` and `H=self._H`.
+
+        (1) Decompose `G/H` into a sum of partial fractions `P_1 +\cdots+ P_r`
+        so that each term of the sum has at most `d` irreducible factors of `H`
+        in the denominator.
+
+        (2) For each differential form `P_j dx_1 \wedge\cdots\wedge dx_d`,
+        find an equivalent form `\omega_j` in de Rham cohomology with no
+        repeated irreducible factors of `H` in its denominator.
+
+        AUTHOR:
+
+        - Alex Raichev (2010-09-22)
+        """
+        if not asy_var:
+            asy_var = var('n')
+
+        # Create symbolic (non-ring) variables.
+        G= self._G
+        H= self._H
+        R= self._R
+        d= self._d
+        X= self._variables
+        # Prepare input for partial_fraction_decomposition() which only works
+        # for functions in the field of fractions of R.
+        if G in R:
+            numer=Integer(1)
+            F= G/H
+        else:
+            numer= G
+            F= R(Integer(1) )/H
+        nstuff= Integer(1) /mul([X[j]**(alpha[j]*asy_var+Integer(1) ) for j in range(d)])
+
+        # Do steps (2) and (1).
+        integrands= []
+        whole,parts= partial_fraction_decomposition(F)
+        for f in parts:
+            a= format_quotient(f)
+            integrands.append( [a[Integer(0) ]*numer*nstuff] + a[Integer(1) :] )
+
+        # Do step (3).
+        ce= decompose_via_cohomology(integrands)
+        ce_new= []
+        for a in ce:
+            ce_new.append( [(a[Integer(0) ]/nstuff).simplify_full().collect(n)] + a[Integer(1) :] )
+        return ce_new
+#-------------------------------------------------------------------------------
+    def critical_cone(self,c,coordinate=None):
+        r"""
+        Returns the critical cone of a convenient multiple point `c`.
+
+        INPUT:
+
+        - ``c`` - A dictionary with keys `self.variables()` and values
+            in a field.
+        - ``coordinate`` - (optional) A natural number.
+
+        OUTPUT:
+
+        A list of vectors that generate the critical cone of `c` and
+        the cone itself if the values of `c` lie in QQ.
+
+        EXAMPLES::
+
+            sage: R.<x,y,z>= PolynomialRing(QQ)
+            sage: G= 1
+            sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y))
+            sage: F= QuasiRationalExpression(G,H)
+            sage: c= {z: 4/3, y: 1, x: 1/2}
+            sage: F.critical_cone(c)
+            ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N)
+
+        NOTES:
+
+        The _critical cone_ of a convenient multiple point `c` with
+        with `c_k \del_k H_j(c) \neq 0` for all `j=1,\ldots,r` is
+        the conical hull of the vectors `\gamma_j(c) =
+        \left(\frac{c_1 \del_1 H_j(c)}{c_k \del_k H_j(c)},\ldots,
+        \frac{c_d \del_d H_j(c)}{c_k \del_k H_j(c)} \right)`.
+        Here `H_1,\ldots,H_r` are the irreducible germs of `self._H` around `c`.
+        For more details, see [RaWi2011]_.
+
+        If this function's optional argument `coordinate` isn't given, then
+        this function searches (from `d` down to 1) for the first index `k`
+        such that for all `j=1,\ldots,r` we have `c_k \del_k H_j(c) \neq 0`
+        and sets `coordinate = k`.
+        Almost.
+        Since this is Python, all the indices actually start at 0.
+
+        REFERENCES:
+
+        .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
+        coefficients of multivariate generating functions: improvements
+        for smooth points", To appear.
+
+        AUTHORS:
+
+        - Alex Raichev (2010-08-25)
+        """
+        Hs= [SR(h[Integer(0) ]) for h in self._Hfac]  # irreducible factors of H
+        X= self._variables
+        d= self._d
+        # Ensure the variables of `c` lie in SR
+        cc= {}
+        for x in c.keys():
+            cc[SR(x)] = c[x]
+        lg= [self._log_grad(h,X,cc) for h in Hs]
+        if not coordinate:
+            # Search (from d-1 down to 0) for a coordinate k such that
+            # for all h in Hs we have cc[k] * diff(h,X[k]) !=0.
+            # One is guaranteed to exist in the case of a convenient multiple
+            # point.
+            for k in reversed(range(d)):
+                if Integer(0)  not in [v[k] for v in lg]:
+                    coordinate= k
+                    break
+        gamma= [self._direction(v,coordinate) for v in lg]
+        if [[gg in QQ for gg in g] for g in gamma] == \
+        [[True for gg in g] for g in gamma]:
+            return gamma,Cone(gamma)    # Cone() needs rational vectors
+        else:
+            return gamma
+# D ============================================================================
+    def _diff_all(self,f,V,n,ending=[],sub=None,sub_final=None,zero_order=Integer(0) ,rekey=None):
+        r"""
+        This function returns a dictionary of representative mixed partial
+        derivatives of `f` from order 1 up to order `n` with respect to the
+        variables in `V`.
+        The default is to key the dictionary by all nondecreasing sequences
+        in `V` of length 1 up to length `n`.
+        For internal use.
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``f`` - An individual or list of `\mathcal{C}^{n+1}` functions.
+        - ``V`` - A list of variables occurring in `f`.
+        - ``n`` - A natural number.
+        - ``ending`` - A list of variables in `V`.
+        - ``sub`` - An individual or list of dictionaries.
+        - ``sub_final`` - An individual or list of dictionaries.
+        - ``rekey`` - A callable symbolic function in `V` or list thereof.
+        - ``zero_order`` - A natural number.
+
+        OUTPUT:
+
+        The dictionary `{s_1:deriv_1,...,s_r:deriv_r}`.
+        Here `s_1,\ldots,s_r` is a listing of
+        all nondecreasing sequences of length 1 up to length `n` over the
+        alphabet `V`, where `w > v` in `X`  iff `str(w) > str(v)`, and
+        `deriv_j` is the derivative of `f` with respect to the derivative
+        sequence `s_j` and simplified with respect to the substitutions in `sub`
+        and evaluated at `sub_final`.
+        Moreover, all derivatives with respect to sequences of length less than
+        `zero_order` (derivatives of order less than `zero_order`) will be made
+        zero.
+
+        If `rekey` is nonempty, then `s_1,\ldots,s_r` will be replaced by the
+        symbolic derivatives of the functions in `rekey`.
+
+        If `ending` is nonempty, then every derivative sequence `s_j` will be
+        suffixed by `ending`.
+
+        EXAMPLES:
+
+        I'd like to print the entire dictionaries, but that doesn't yield
+        consistent output order for doctesting.
+        Order of keys changes. 
+        
+        ::
+
+            sage: R.<x> = PolynomialRing(QQ)
+            sage: G,H = 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: f= function('f',x)
+            sage: dd= F._diff_all(f,[x],3)
+            sage: dd[(x,x,x)]
+            D[0, 0, 0](f)(x)
+
+            ::
+
+            sage: d1= {diff(f,x): 4*x^3}
+            sage: dd= F._diff_all(f,[x],3,sub=d1)
+            sage: dd[(x,x,x)]
+            24*x
+
+            ::
+
+            sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f)
+            sage: dd[diff(f,x,3)]
+            24*x
+
+            ::
+
+            sage: a= {x:1}
+            sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f,sub_final=a)
+            sage: dd[diff(f,x,3)]
+            24
+
+            ::
+
+            sage: X= var('x,y,z')
+            sage: f= function('f',*X)
+            sage: dd= F._diff_all(f,X,2,ending=[y,y,y])
+            sage: dd[(z,y,y,y)]
+            D[1, 1, 1, 2](f)(x, y, z)
+
+            ::
+
+            sage: g= function('g',*X)
+            sage: dd= F._diff_all([f,g],X,2)
+            sage: dd[(0,y,z)]
+            D[1, 2](f)(x, y, z)
+
+            ::
+
+            sage: dd[(1,z,z)]
+            D[2, 2](g)(x, y, z)
+
+            ::
+
+            sage: f= exp(x*y*z)
+            sage: ff= function('ff',*X)
+            sage: dd= F._diff_all(f,X,2,rekey=ff)
+            sage: dd[diff(ff,x,z)]
+            x*y^2*z*e^(x*y*z) + y*e^(x*y*z)
+
+        AUTHORS:
+
+        - Alex Raichev (2008-10-01, 2009-04-01, 2010-02-01)
+        """
+        singleton=False
+        if not isinstance(f,list):
+            f= [f]
+            singleton=True
+
+        # Build the dictionary of derivatives iteratively from a list of nondecreasing
+        # derivative-order sequences.
+        derivs= {}
+        r= len(f)
+        if ending:
+            seeds = [ending]
+            start = Integer(1)
+        else:
+            seeds = [[v] for v in V]
+            start = Integer(2)
+        if singleton:
+            for s in seeds:
+                derivs[tuple(s)] = self._subs_all(diff(f[Integer(0) ],s),sub)
+            for l in (ellipsis_range(start,Ellipsis,n)):
+                for t in UnorderedTuples(V,l):
+                    s= tuple(t + ending)
+                    derivs[s] = self._subs_all(diff(derivs[s[Integer(1) :]],s[Integer(0) ]),sub)
+        else:
+            # Make the dictionary keys of the form (j,sequence of variables),
+            # where j in range(r).
+            for s in seeds:
+                value= self._subs_all([diff(f[j],s) for j in range(r)],sub)
+                derivs.update(dict([(tuple([j]+s),value[j]) for j in range(r)]))
+            for l in (ellipsis_range(start,Ellipsis,n)):
+                for t in UnorderedTuples(V,l):
+                    s= tuple(t + ending)
+                    value= self._subs_all(\
+                    [diff(derivs[(j,)+s[Integer(1) :]],s[Integer(0) ]) for j in range(r)],sub)
+                    derivs.update(dict([((j,)+s,value[j]) for j in range(r)]))
+        if zero_order:
+            # Zero out all the derivatives of order < zero_order
+            if singleton:
+                for k in derivs.keys():
+                    if len(k) < zero_order:
+                        derivs[k]= Integer(0)
+            else:
+                # Ignore the first of element of k, which is an index.
+                for k in derivs.keys():
+                    if len(k)-Integer(1)  < zero_order:
+                        derivs[k]= Integer(0)
+        if sub_final:
+            # Substitute sub_final into the values of derivs.
+            for k in derivs.keys():
+                derivs[k] = self._subs_all(derivs[k],sub_final)
+        if rekey:
+            # Rekey the derivs dictionary by the value of rekey.
+            F= rekey
+            if singleton:
+                # F must be a singleton.
+                derivs= dict( [(diff(F,list(k)), derivs[k]) for k in derivs.keys()] )
+            else:
+                # F must be a list.
+                derivs= dict( [(diff(F[k[Integer(0) ]],list(k)[Integer(1) :]), derivs[k]) for k in derivs.keys()] )
+        return derivs
+#-------------------------------------------------------------------------------
+    def _diff_op(self,A,B,AB_derivs,V,M,r,N):
+        r"""
+        This function computes the derivatives `DD^(l+k)(A[j] B^l)` evaluated at a
+        point `p` for various natural numbers `j,k,l` which depend on `r` and `N`.
+        Here `DD` is a specific second-order linear differential operator that depends
+        on `M`, `A` is a list of symbolic functions, `B` is symbolic function,
+        and `AB_derivs` contains all the derivatives of `A` and `B` evaluated at `p`
+        that are necessary for the computation.
+        For internal use by the functions _asymptotics_main_multiple() and
+        _asymptotics_main_smooth().
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``A`` - A single or length `r` list of symbolic functions in the
+            variables `V`.
+        - ``B`` - A symbolic function in the variables `V`.
+        - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of
+            `A[0],\ldots,A[r-1]` up to order `2N-2` and
+            the (symbolic) derivatives of `B` up to order `2N`.
+            The values of the dictionary are complex numbers that are
+            the keys evaluated at a common point `p`.
+        - ``V`` - The variables of the `A[j]` and `B`.
+        - ``M`` - A symmetric `l \times l` matrix, where `l` is the length of `V`.
+        - ``r,N`` - Natural numbers.
+
+        OUTPUT:
+
+        A dictionary whose keys are natural number tuples of the form `(j,k,l)`,
+        where `l \le 2k`, `j \le r-1`, and `j+k \le N-1`, and whose values are
+        `DD^(l+k)(A[j] B^l)` evaluated at a point `p`, where `DD` is the linear
+        second-order differential operator
+        `-\sum_{i=0}^{l-1} \sum_{j=0}^{l-1} M[i][j]
+        \partial^2 /(\partial V[j] \partial V[i])`.
+
+        EXAMPLES::
+
+            sage: R.<x> = PolynomialRing(QQ)
+            sage: G,H = 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: T= var('x,y')
+            sage: A= function('A',*tuple(T))
+            sage: B= function('B',*tuple(T))
+            sage: AB_derivs= {}
+            sage: M= matrix([[1,2],[2,1]])
+            sage: DD= F._diff_op(A,B,AB_derivs,T,M,1,2)
+            sage: DD.keys()
+            [(0, 1, 2), (0, 1, 1), (0, 1, 0), (0, 0, 0)]
+            sage: len(DD[(0,1,2)])
+            246
+
+        AUTHORS:
+
+        - Alex Raichev (2008-10-01, 2010-01-12)
+        """
+        if not isinstance(A,list):
+            A= [A]
+
+        # First, compute the necessary product derivatives of A and B.
+        product_derivs= {}
+        for (j,k) in mrange([r,N]):
+            if j+k <N:
+                for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)):
+                    for s in UnorderedTuples(V,Integer(2) *(k+l)):
+                        product_derivs[tuple([j,k,l]+s)] = \
+                        diff(A[j]*B**l,s).subs(AB_derivs)
+
+        # Second, compute DD^(k+l)(A[j]*B^l)(p) and store values in dictionary.
+        DD= {}
+        rows= M.nrows()
+        for (j,k) in mrange([r,N]):
+            if j+k <N:
+                for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)):
+                    # Take advantage of the symmetry of M by ignoring
+                    # the upper-diagonal entries of M and multiplying by
+                    # appropriate powers of 2.
+                    if k+l == Integer(0) :
+                        DD[(j,k,l)] = product_derivs[(j,k,l)]
+                        continue
+                    S= [(a,b) for (a,b) in mrange([rows,rows]) if b <= a]
+                    P=  cartesian_product_iterator([S for i in range(k+l)])
+                    diffo= Integer(0)
+                    for t in P:
+                        if product_derivs[(j,k,l)+self._diff_seq(V,t)] != Integer(0) :
+                            MM= Integer(1)
+                            for (a,b) in t:
+                                MM= MM * M[a][b]
+                                if a != b:
+                                    MM= Integer(2) *MM
+                            diffo= diffo + MM * product_derivs[(j,k,l)+self._diff_seq(V,t)]
+                    DD[(j,k,l)] = (-Integer(1) )**(k+l)*diffo
+        return DD
+#-------------------------------------------------------------------------------
+    def _diff_op_simple(self,A,B,AB_derivs,x,v,a,N):
+        r"""
+        This function computes `DD^(e k + v l)(A B^l)` evaluated at a point `p`
+        for various natural numbers `e,k,l` that depend on `v` and `N`.
+        Here `DD` is a specific linear differential operator that depends
+        on `a` and `v`, `A` and `B` are symbolic functions, and `AB_derivs` contains
+        all the derivatives of `A` and `B` evaluated at `p` that are necessary for
+        the computation.
+        For internal use by the function _asymptotics_main_smooth().
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``A``,``B`` - Symbolic functions in the variable `x`.
+        - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of
+            `A` up to order `2N` if `v` is even or `N` if `v` is odd and
+            the (symbolic) derivatives of `B` up to order `2N+v` if `v` is even
+            or `N+v` if `v` is odd.
+            The values of the dictionary are complex numbers that are
+            the keys evaluated at a common point `p`.
+        - ``x`` - Symbolic variable.
+        - ``a`` - A complex number.
+        - ``v``,``N`` - Natural numbers.
+
+        OUTPUT:
+
+        A dictionary whose keys are natural number pairs of the form `(k,l)`,
+        where `k < N` and `l \le 2k` and whose values are
+        `DD^(e k + v l)(A B^l)` evaluated at a point `p`.
+        Here `e=2` if `v` is even, `e=1` if `v` is odd, and `DD` is a particular
+        linear differential operator
+        `(a^{-1/v} d/dt)' if `v` is even and `(|a|^{-1/v} i \sgn a d/dt)`
+        if `v` is odd.
+
+        EXAMPLES::
+
+            sage: R.<x> = PolynomialRing(QQ)
+            sage: G,H = 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: A= function('A',x)
+            sage: B= function('B',x)
+            sage: AB_derivs= {}
+            sage: F._diff_op_simple(A,B,AB_derivs,x,3,2,2)
+            {(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)}
+
+        AUTHORS:
+
+        - Alex Raichev (2010-01-15)
+        """
+        I= sqrt(-Integer(1) )
+        DD= {}
+        if v.mod(Integer(2) ) == Integer(0) :
+            for k in range(N):
+                for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)):
+                    DD[(k,l)] = (a**(-Integer(1) /v))**(Integer(2) *k+v*l) \
+                    * diff(A*B**l,x,Integer(2) *k+v*l).subs(AB_derivs)
+        else:
+            for k in range(N):
+                for l in (ellipsis_range(Integer(0) ,Ellipsis,k)):
+                    DD[(k,l)] = (abs(a)**(-Integer(1) /v)*I*a/abs(a))**(k+v*l) \
+                    * diff(A*B**l,x,k+v*l).subs(AB_derivs)
+        return DD
+#-------------------------------------------------------------------------------
+    def _diff_prod(self,f_derivs,u,g,X,interval,end,uderivs,atc):
+        r"""
+        This function takes various derivatives of the equation `f=ug`,
+        evaluates at a point `c`, and solves for the derivatives of `u`.
+        For internal use by the function _asymptotics_main_multiple().
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``f_derivs`` - A dictionary whose keys are tuples \code{s + end} for all
+            `X`-variable sequences `s` with length in `interval` and whose
+            values are the derivatives of a function `f` evaluated at `c`.
+        - ``u`` - A callable symbolic function.
+        - ``g`` - An expression or callable symbolic function.
+        - ``X`` - A list of symbolic variables.
+        - ``interval`` - A list of positive integers.
+            Call the first and last values `n` and `nn`, respectively.
+        - ``end`` - A possibly empty list of `z`'s where `z` is the last element of
+            `X`.
+        - ``uderivs`` - A dictionary whose keys are the symbolic
+            derivatives of order 0 to order `n-1` of `u` evaluated at `c`
+            and whose values are the corresponding derivatives evaluated at `c`.
+        - ``atc`` - A dictionary whose keys are the keys of `c` and all the symbolic
+            derivatives of order 0 to order `nn` of `g` evaluated `c` and whose
+            values are the corresponding derivatives evaluated at `c`.
+
+        OUTPUT:
+
+        A dictionary whose keys are the derivatives of `u` up to order
+        `nn` and whose values are those derivatives evaluated at `c`.
+
+        EXAMPLES::
+
+        I'd like to print out the entire dictionary, but that does not give
+        consistent output for doctesting.
+        Order of keys changes ::
+
+            sage: R.<x> = PolynomialRing(QQ)
+            sage: G,H = 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: u= function('u',x)
+            sage: g= function('g',x)
+            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})
+            sage: dd[diff(u,x,2)(x=2)]
+            22/9
+
+        NOTES:
+
+        This function works by differentiating the equation `f=ug` with respect
+        to the variable sequence \code{s+end}, for all tuples `s` of `X` of
+        lengths in `interval`, evaluating at the point `c`,
+        and solving for the remaining derivatives of `u`.
+        This function assumes that `u` never appears in the differentiations of
+        `f=ug` after evaluating at `c`.
+
+        AUTHORS:
+
+        - Alex Raichev (2009-05-14, 2010-01-21)
+        """
+        for l in interval:
+            D= {}
+            rhs= []
+            lhs= []
+            for t in UnorderedTuples(X,l):
+                s= t+end
+                lhs.append(f_derivs[tuple(s)])
+                rhs.append(diff(u*g,s).subs(atc).subs(uderivs))
+                # Since Sage's solve command can't take derivatives as variable
+                # names, i make new variables based on t to stand in for
+                # diff(u,t) and store them in D.
+                D[diff(u,t).subs(atc)] = self._make_var([var('zing')]+t)
+            eqns=[ lhs[i] == rhs[i].subs(uderivs).subs(D) for i in range(len(lhs))]
+            vars= D.values()
+            sol= solve(eqns,*vars,solution_dict=True)
+            uderivs.update(self._subs_all(D,sol[Integer(0) ]))
+        return uderivs
+#-------------------------------------------------------------------------------
+    def _diff_seq(self,V,s):
+        r"""
+        Given a list `s` of tuples of natural numbers, this function returns the
+        list of elements of `V` with indices the elements of the elements of `s`.
+        This function is for internal use by the function _diff_op().
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``V`` - A list.
+        - ``s`` - A list of tuples of natural numbers in the interval
+            \code{range(len(V))}.
+
+        OUTPUT:
+
+        The tuple \code{tuple([V[tt] for tt in sorted(t)])}, where `t` is the
+        list of elements of the elements of `s`.
+
+        EXAMPLES::
+
+            sage: R.<x> = PolynomialRing(QQ)
+            sage: G,H = 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: V= list(var('x,t,z'))
+            sage: F._diff_seq(V,([0,1],[0,2,1],[0,0]))
+            (x, x, x, x, t, t, z)
+
+        AUTHORS:
+
+        - Alex Raichev (2009.05.19)
+        """
+        t= []
+        for ss in s:
+            t.extend(ss)
+        return tuple([V[tt] for tt in sorted(t)])
+#-------------------------------------------------------------------------------
+    def _direction(self,v,coordinate=None):
+        r"""
+        This function returns \code{[vv/v[coordinate] for vv in v]} where
+        coordinate is the last index of v if not specified otherwise.
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``v`` - A vector.
+        - ``coordinate`` - An index for `v` (default: None).
+
+        OUTPUT:
+
+        This function returns \code{[vv/v[coordinate] for vv in v]} where
+        coordinate is the last index of v if not specified otherwise.
+
+        EXAMPLES::
+
+            sage: R.<x> = PolynomialRing(QQ)
+            sage: G,H = 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: F._direction([2,3,5])
+            (2/5, 3/5, 1)
+            sage: F._direction([2,3,5],0)
+            (1, 3/2, 5/2)
+
+        AUTHORS:
+
+        - Alex Raichev (2010-08-25)
+        """
+        if coordinate == None:
+            coordinate= len(v)-Integer(1)
+        try:
+            v[Integer(0) ].variables()
+            return tuple([(vv/v[coordinate]).simplify_full() for vv in v])
+        except:
+            return tuple([vv/v[coordinate] for vv in v])
+# E ============================================================================
+# F ============================================================================
+# G ============================================================================
+# H ============================================================================
+# I ============================================================================
+    def is_cmp(self,points):
+        r"""
+        Checks if the points in the list `points` are convenient multiple
+        points of `V= \{ x\in CC^d : H(x) = 0\}`, where `H=self._H`.
+
+        INPUT:
+
+        - ``points`` - An individual or list of dictionaries with keys
+            `self._variables` and values in some superfield of
+            `self._R.base_ring()`.
+
+        OUTPUT:
+
+        A list of tuples `(p,verdict,comment)`, one for each point
+        `p` in `points`, where `verdict` is True if `p` is a convenient
+        multiple point and False otherwise, and where `comment` is a string
+        comment relating to `verdict`, such as 'not a transverse intersection'.
+
+        EXAMPLES::
+
+            sage: R.<x,y,z>= PolynomialRing(QQ)
+            sage: G= 16
+            sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y))
+            sage: F= QuasiRationalExpression(G,H)
+            sage: points= [{x:1/2,y:1,z:4/3},{x:1/2,y:1,z:2}]
+            sage: F.is_cmp(points)
+            [({y: 1, z: 4/3, x: 1/2}, True, 'all good'), ({y: 1, z: 2, x: 1/2}, False, 'not a singular point')]
+
+        NOTES:
+
+        A point `c` of `V` is a __convenient multiple point__ if `V` is locally
+        a union of complex manifolds that intersect transversely at `c`;
+        see [RaWi2011]_.
+
+        REFERENCES:
+
+        .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
+        coefficients of multivariate generating functions: improvements
+        for smooth points", To appear.
+
+        AUTHORS:
+
+        - Alex Raichev (2011-04-18)
+        """
+        H= self._H
+        d= self._d
+        Hs= [SR(h[Integer(0) ]) for h in self._Hfac]  # irreducible factors of H
+        X= self._variables
+        J= [tuple([diff(h,x) for x in X]) for h in Hs]
+        verdicts= []
+        if not isinstance(points,list):
+            points= [points]
+        for p in points:
+            # Ensure variables in points lie in SR.
+            pp= {}
+            for x in p.keys():
+                pp[SR(x)]= p[x]
+            # Test 1: Is p a zero of all factors of H?
+            if [h.subs(pp) for h in Hs] != [Integer(0)  for h in Hs]:
+                # Failed test 1.  Move on to next point.
+                verdicts.append((p,False,'not a singular point'))
+                continue
+            # Test 2: Are the factors of H smooth and
+            # do theyintersect transversely at p?
+            J= [tuple([f.subs(pp) for f in dh]) for dh in J]
+            l= len(J)
+            if Set(J).cardinality() < l:
+                # Fail.  Move on to next point.
+                verdicts.append((p,False,'not a transverse intersection'))
+                continue
+            temp= True
+            for S in list(Set(J).subsets())[Integer(1) :]: # Subsets of size >= 1
+                k= len(list(S))
+                M = matrix(list(S))
+                if  M.rank() != min(k,d):
+                    # Fail.
+                    temp= False
+                    verdicts.append((p,False,'not a transvere intersection'))
+                    break
+            if not temp:
+                # Move on to next point
+                continue
+            # Test 3: Is p convenient?
+            Jlog= matrix([self._log_grad(h,X,pp) for h in Hs])
+            if [Integer(0)  in f for f in Jlog.columns()] ==\
+            [True for f in Jlog.columns()]:
+                # Fail.  Move on to next point.
+                verdict[p] = (False,'multiple point but not convenient')
+                continue
+            verdicts.append((p,True,'all good'))
+        return verdicts
+# J ============================================================================
+# K ============================================================================
+# L ============================================================================
+    def _log_grad(self,f,X,c):
+        r"""
+        This function returns the logarithmic gradient of `f` with respect to the
+        variables of `X` evalutated at `c`.
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``f`` - An expression in the variables of `X`.
+        - ``X`` - A list of variables.
+        - ``c`` - A dictionary with keys `X` and values in a field `K`.
+
+        OUTPUT:
+
+            \code{[c[x] * diff(f,x).subs(c) for x in X]}.
+
+        EXAMPLES::
+
+            sage: R.<x> = PolynomialRing(QQ)
+            sage: G,H = 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: X= var('x,y,z')
+            sage: f= x*y*z^2
+            sage: c= {x:1,y:2,z:3}
+            sage: f.gradient()
+            (y*z^2, x*z^2, 2*x*y*z)
+            sage: F._log_grad(f,X,c)
+            (18, 18, 36)
+
+            ::
+
+            sage: R.<x,y,z>= PolynomialRing(QQ)
+            sage: f= x*y*z^2
+            sage: c= {x:1,y:2,z:3}
+            sage: F._log_grad(f,R.gens(),c)
+            (18, 18, 36)
+
+        AUTHORS:
+
+        - Alex Raichev (2009-03-06)
+        """
+        return tuple([SR(c[x] * diff(f,x).subs(c)).simplify() for x in X])
+# M ============================================================================
+    def _make_var(self,L):
+        r"""
+        This function converts the list `L` to a string (without commas) and returns
+        the string as a variable.
+        For internal use by the function _diff_op()
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``L`` - A list.
+
+        OUTPUT:
+
+        A variable whose name is the concatenation of the variable names in `L`.
+
+        EXAMPLES::
+
+            sage: R.<x> = PolynomialRing(QQ)
+            sage: G,H = 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: L= list(var('x,y,hello'))
+            sage: v= F._make_var(L)
+            sage: print v, type(v)
+            xyhello <type 'sage.symbolic.expression.Expression'>
+
+        AUTHORS:
+
+        - Alex Raichev (2010-01-21)
+        """
+        return var(''.join([str(v) for v in L]))
+# N ============================================================================
+    def _new_var_name(self,name,V):
+        r"""
+        This function returns the first string in the sequence `name`, `name+name`,
+        `name+name+name`,... that does not appear in the list `V`.
+        It is for internal use by the function _asymptotics_main_multiple().
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``name`` - A string.
+        - ``V`` - A list of variables.
+
+        OUTPUT:
+
+        The first string in the sequence `name`, `name+name`,
+        `name+name+name`,... that does not appear in the list \code{str(V)}.
+
+        EXAMPLES::
+
+            sage: R.<x> = PolynomialRing(QQ)
+            sage: G,H = 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: X= var('x,xx,y,z')
+            sage: F._new_var_name('x',X)
+            'xxx'
+
+        AUTHORS:
+
+        - Alex Raichev (2008-10-01)
+        """
+        newname= name
+        while newname in str(V):
+            newname= newname +name
+        return newname
+# O ============================================================================
+# P ============================================================================
+# Q ============================================================================
+# R ============================================================================
+    def relative_error(self,approx,alpha,interval,exp_scale=Integer(1) ):
+        r"""
+        Returns the relative error between the values of the Maclaurin
+        coefficients of `self` with multi-indices `m alpha` for `m` in
+        `interval` and the values of the functions in `approx`.
+
+        INPUT:
+
+        - ``approx`` - An individual or list of symbolic expressions in
+            one variable.
+        - ``alpha`` - A list of positive integers.
+        - ``interval`` - A list of positive integers.
+        - ``exp_scale`` - (optional) A number.  Default: 1.
+
+        OUTPUT:
+
+        A list whose entries are of the form
+        `[m\alpha,a_m,b_m,err_m]` for `m \in interval`.
+        Here `m\alpha` is a tuple; `a_m` is the `m alpha` (multi-index)
+        coefficient of the Maclaurin series for `F` divided by `exp_scale^m`;
+        `b_m` is a list of the values of the functions in `approx` evaluated at
+        `m` and divided by `exp_scale^m`; `err_m` is the list of relative errors
+        `(a_m-f)/a_m` for `f` in `b_m`.
+        All outputs are decimal approximations.
+
+        EXAMPLES::
+
+            sage: R.<x,y>= PolynomialRing(QQ)
+            sage: G=1
+            sage: H= 1-x-y-x*y
+            sage: F= QuasiRationalExpression(G,H)
+            sage: alpha= [1,1]
+            sage: var('n')
+            n
+            sage: f= (0.573/sqrt(n))*5.83^n
+            sage: es= 5.83
+            sage: F.relative_error(f,alpha,[1,2,4,8],es)
+            Calculating errors table in the form
+            exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors...
+            [[(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]]]
+            sage: g= (0.573/sqrt(n) - 0.0674/n^(3/2))*5.83^n
+            sage: F.relative_error([f,g],alpha,[1,2,4,8],es)
+            Calculating errors table in the form
+            exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors...
+            [[(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]]]
+
+        AUTHORS:
+
+        - Alex Raichev (2009-05-18, 2011-04-18)
+        """
+
+        if not isinstance(approx,list):
+            approx= [approx]
+        av= approx[Integer(0) ].variables()[Integer(0) ]
+
+        print "Calculating errors table in the form"
+        print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..."
+
+        # Get Maclaurin coefficients of self and scale them.
+        # Then compute errors.
+        n= interval[-Integer(1) ]
+        mc= self.maclaurin_coefficients(alpha,n)
+        mca={}
+        stats=[]
+        for m in interval:
+            beta= tuple([m*a for a in alpha])
+            mc[beta]= mc[beta]/exp_scale**m
+            mca[beta]= [f.subs({av:m})/exp_scale**m for f in approx]
+            stats_row= [beta, mc[beta].n(), [a.n() for a in mca[beta]]]
+            if mc[beta]==Integer(0) :
+                 stats_row.extend([None for a in mca[beta]])
+            else:
+                 stats_row.append([((mc[beta]-a)/mc[beta]).n() for a in mca[beta]])
+            stats.append(stats_row)
+        return stats
+# S ============================================================================
+    def singular_points(self):
+        r"""
+        This function returns a Groebner basis ideal whose variety is the
+        set of singular points of the algebraic variety
+        `V= \{x\in\CC^d : H(x)=0\}`, where `H=sef._H`.
+
+        INPUT:
+
+        OUTPUT:
+
+        A Groebner basis ideal whose variety is the set of singular points of
+        the algebraic variety `V= \{x\in\CC^d : H(x)=0\}`.
+
+        EXAMPLES::
+
+            sage: R.<x,y,z>= PolynomialRing(QQ)
+            sage: G= 1
+            sage: H= (4-2*x-y-z)*(4-x-2*y-z)
+            sage: F= QuasiRationalExpression(G,H)
+            sage: F.singular_points()
+            Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field
+
+        AUTHORS:
+
+        - Alex Raichev (2008-10-01, 2008-11-20, 2010-12-03, 2011-04-18)
+        """
+        H= self._H
+        R= self._R
+        f= R.ideal(H).radical().gens()[Integer(0) ]     # Compute the reduction of H.
+        J= R.ideal([f] + f.gradient())
+        return R.ideal(J.groebner_basis())
+#-------------------------------------------------------------------------------
+    def smooth_critical(self,alpha):
+        r"""
+        This function returns a Groebner basis ideal whose variety is the set
+        of smooth critical points of the algebraic variety
+        `V= \{x\in\CC^d : H(x)=0\} for the direction `\alpha` where `H=self._H`.
+
+        INPUT:
+
+        - ``alpha`` - A `d`-tuple of positive integers and/or symbolic entries.
+
+        OUTPUT:
+
+        A Groebner basis ideal of smooth critical points of `V` for `\alpha`.
+
+        EXAMPLES::
+
+            sage: R.<x,y> = PolynomialRing(QQ)
+            sage: G=1
+            sage: H = (1-x-y-x*y)^2
+            sage: F= QuasiRationalExpression(G,H)
+            sage: var('a1,a2')
+            (a1, a2)
+            sage: F.smooth_critical([a1,a2])
+            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
+
+        NOTES:
+
+        A point `c` of `V` is a __smooth critical point for `alpha`__
+        if the gradient of `f` at `c` is not identically zero and `\alpha` is in
+        the span of the logarithmic gradient vector
+        `(c[0] \partial_1 f(c)),\ldots,c[d-1] \partial_d f(c))`; see [RaWi2008a]_.
+
+        REFERENCES:
+
+        .. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of
+        coefficients of multivariate generating functions: improvements
+        for smooth points", Electronic Journal of Combinatorics, Vol. 15
+        (2008), R89.
+
+        AUTHORS:
+
+        - Alex Raichev (2008-10-01, 2008-11-20, 2009-03-09, 2010-12-02, 2011-04-18)
+        """
+        H= self._H
+        R= H.parent()
+        B= R.base_ring()
+        d= R.ngens()
+        vars= R.gens()
+        f= R.ideal(H).radical().gens()[Integer(0) ]   # Compute the reduction of H.
+
+        # Expand B by the variables of alpha if there are any.
+        indets= []
+        indets_ind= []
+        for a in alpha:
+            if not ((a in ZZ) and (a>Integer(0) )):
+                try:
+                    CC(a)
+                except:
+                    indets.append(var(a))
+                    indets_ind.append(alpha.index(a))
+                else:
+                    print "The components of", alpha, \
+                    "must be positive integers or symbolic variables."
+                    return
+        indets= list(Set(indets))   # Delete duplicates in indets.
+        if indets != []:
+            BB= FractionField(PolynomialRing(B,tuple(indets)))
+            S= R.change_ring(BB)
+            vars= S.gens()
+            # Coerce alpha into BB.
+            for i in range(len(alpha)):
+                alpha[i] = BB(alpha[i])
+        else:
+            S= R
+
+        # Find smooth, critical points for alpha.
+        f= S(f)
+        J= S.ideal([f] +[ alpha[d-Integer(1) ]*vars[i]*diff(f,vars[i]) \
+        -alpha[i]*vars[d-Integer(1) ]*diff(f,vars[d-Integer(1) ]) for i in range(d-Integer(1) )])
+        return S.ideal(J.groebner_basis())
+#-------------------------------------------------------------------------------
+    def _subs_all(self,f,sub,simplify=False):
+        r"""
+        This function returns the items of `f` substituted by the dictionaries of
+        `sub` in order of their appearance in `sub`.
+
+        Does not use `self`.
+
+        INPUT:
+
+        - ``f`` - An individual or list of symbolic expressions or dictionaries
+        - ``sub`` - An individual or list of dictionaries.
+        - ``simplify`` - Boolean (default: False).
+
+        OUTPUT:
+
+        The items of `f` substituted by the dictionaries of `sub` in order of
+        their appearance in `sub`.  The subs() command is used.
+        If simplify is True, then simplify() is used after substitution.
+
+        EXAMPLES::
+
+            sage: R.<x> = PolynomialRing(QQ)
+            sage: G,H = 1,1
+            sage: F= QuasiRationalExpression(G,H)
+            sage: var('x,y,z')
+            (x, y, z)
+            sage: a= {x:1}
+            sage: b= {y:2}
+            sage: c= {z:3}
+            sage: F._subs_all(x+y+z,a)
+            y + z + 1
+            sage: F._subs_all(x+y+z,[c,a])
+            y + 4
+            sage: F._subs_all([x+y+z,y^2],b)
+            [x + z + 2, 4]
+            sage: F._subs_all([x+y+z,y^2],[b,c])
+            [x + 5, 4]
+
+            ::
+
+            sage: var('x,y')
+            (x, y)
+            sage: a= {'foo':x^2+y^2, 'bar':x-y}
+            sage: b= {x:1,y:2}
+            sage: F._subs_all(a,b)
+            {'foo': 5, 'bar': -1}
+
+        AUTHORS:
+
+        - Alex Raichev (2009-05-05)
+        """
+        singleton= False
+        if not isinstance(f,list):
+            f= [f]
+            singleton= True
+        if not isinstance(sub,list):
+            sub= [sub]
+        g= []
+        for ff in f:
+            for D in sub:
+                if isinstance(ff,dict):
+                    ff= dict( [(k,ff[k].subs(D)) for k in ff.keys()] )
+                else:
+                    ff= ff.subs(D)
+            g.append(ff)
+        if singleton and simplify:
+            if isinstance(g[Integer(0) ],dict):
+                return g[Integer(0) ]
+            else:
+                return g[Integer(0) ].simplify()
+        elif singleton and not simplify:
+            return g[Integer(0) ]
+        elif not singleton and simplify:
+            G= []
+            for gg in g:
+                if isinstance(gg,dict):
+                    G.append(gg)
+                else:
+                    G.append(gg.simplify())
+            return G
+        else:
+            return g
+# T ============================================================================
+    def maclaurin_coefficients(self,alpha,n):
+        r"""
+        Returns the Maclaurin coefficients of self that have multi-indices
+        `alpha`, `2*alpha`,...,`n*alpha`.
+
+        INPUT:
+
+        - ``n`` - positive integer
+        - ``alpha`` - tuple of positive integers representing a multi-index.
+
+        OUTPUT:
+
+        A dictionary of the form (beta, beta Maclaurin coefficient of self).
+
+        AUTHORS:
+
+        - Alex Raichev (2011-04-08)
+        """
+        # Do all computations in the Symbolic Ring.
+        F= SR(self._G/self._H)
+        v= F.variables()
+        p= {}
+        for x in v:
+            p[x]=Integer(0)
+        d= len(v)
+        coeffs={}
+        # Initialize sequence of variables to differentiate with.
+        s=[]
+        for i in range(d):
+            s.extend([v[i] for j in range(alpha[i])])
+        F_deriv= diff(F,s)
+        coeffs[tuple(alpha)]= F_deriv.subs(p)/ mul([factorial(a) for a in alpha])
+        old_beta= alpha
+        for k in (ellipsis_range(Integer(2) ,Ellipsis,n)):
+            # Update variable sequence to differentiate with.
+            beta= [k*a for a in alpha]
+            delta= [beta[i]-old_beta[i] for i in range(d)]
+            s= []
+            for i in range(d):
+                s.extend([v[i] for j in range(delta[i])])
+            F_deriv= diff(F_deriv,s)
+            coeffs[tuple(beta)]= F_deriv.subs(p)/ mul([factorial(b) for b in beta])
+            old_beta= copy(beta)
+        return coeffs
+# U ============================================================================
+# V ============================================================================
+    def variables(self):
+        r"""
+        Returns the tuple of variables of `self`.
+
+        EXAMPLES::
+
+            sage: R.<x,y>= PolynomialRing(QQ)
+            sage: G= exp(x)
+            sage: H= 1-y
+            sage: F= QuasiRationalExpression(G,H)
+            sage: F.variables()
+            (x, y)
+
+            sage: R.<x,y>= PolynomialRing(QQ,order='invlex')
+            sage: G= exp(x)
+            sage: H= 1-y
+            sage: F= QuasiRationalExpression(G,H)
+            sage: F.variables()
+            (y, x)
+
+        AUTHORS:
+
+        - Alex Raichev (2011-04-01)
+        """
+        return self._variables
+# W ============================================================================
+# X ============================================================================
+# Y ============================================================================
+# Z ============================================================================