Ticket #10519: amgf.sage

r"""
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.  
Turns out 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/
#*****************************************************************************
#
# 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: 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[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[0:d]
    Ts= Vs[d:]
    
    # Find an irreducible annihilating polynomial g for the fs if there is one.
    J= ideal([ Ts[j] -RR(fs[j]) for j in range(r)])
    JJ= J.elimination_ideal(Xs)
    g= JJ.gens()[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: 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==1:
        return f.partial_fraction_decomposition()
    Q= f.denominator()
    whole,P= f.numerator().quo_rem(Q)
    parts= [format_quotient(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() == 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: 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[1:])
    # Combine items of parts with same denominators.
    newnew_parts=[]
    left,right=0,1
    glom= new_parts[left][0]
    while right <= len(new_parts)-1:
        if new_parts[left][1:] == new_parts[right][1:]:
            glom= glom +new_parts[right][0]
        else:
            newnew_parts.append([glom]+new_parts[left][1:])
            left= right
            glom= new_parts[left][0]
        right= right +1  
    if glom != 0:
        newnew_parts.append([glom]+new_parts[left][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[0]
        Qs= p[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[0]
            Qs= p[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())[0:m]
            denom= [[new_vars[j], e[0][j]+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][0]^Qs[j][1]) for j in range(m)])
            p_parts=[]
            for x in p_parts_temp:
                y= P*F.parent()(x[0]).subs(Qpowsub)
                yy=[]
                for xx in x[1:]:
                    if xx[0] in new_vars:
                        j= new_vars.index(xx[0])
                        yy.append([Qs[j][0],Qs[j][1]*xx[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: 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
    import copy     # Will need this to make copies of a nested list.
    decomposing_done= True
    new_parts= []
    R= parts[0][1][0].parent()
    V= list(R.gens())
    for p in parts:
        p_parts= [p]
        P= p[0]
        Qs= p[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(1).lift(R.ideal([q for q,e in Qs] +dets))
            for j in range(m):
                if L[j] != 0:
                    # Make a copy of (and not a reference to) the nested list Qs.
                    new_Qs = copy.deepcopy(Qs)
                    if new_Qs[j][1] > 1:
                        new_Qs[j][1]= new_Qs[j][1] -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] != 0:
                    new_Qs = copy.deepcopy(Qs)
                    for j in range(m):
                        if new_Qs[j][1] > 1:
                            new_Qs[j][1]= new_Qs[j][1] -1
                            v=  sorted(vars_subsets[k],reverse=true)
                            jac= jacobian([SR(P*L[m+k])] +[ SR(Qs[jj][0]) for \
                            jj in range(m) if jj !=j], [SR(vv) for vv in v])
                            det= jac.determinant()
                            if det != 0:
                                p_parts.append([permutation_sign(v,V) \
                                  *(-1)^j/new_Qs[j][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[0][0].parent()
    for p in parts:
        p_parts= [p]
        P= p[0]
        Qs= p[1:]
        m= len(Qs)
        if R(1) in ideal([q for q,e in Qs]):
            # Then we can decompose p.
            decomposing_done= False
            L= R(1).lift(R.ideal([q^e for q,e in Qs]))
            p_parts= [ [P*L[i]] + \
            [[Qs[j][0],Qs[j][1]] for j in range(m) if j != i] \
            for i in range(m) if L[i]!=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() != 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: 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+1 for x in range(len(vars))]
    B= [vars.index(x)+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[0]
    Qs= part[1:]
    Q= prod([q^e for q,e in Qs])
    return P/Q
#===============================================================================
# Class for calculation of asymptotics of multivariate generating functions.
class QuasiRationalExpression(object):
    "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=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", submitted.

        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[0])
            HH= [SR(f) for (f,e) in chunk[1:]]
            asy= self._asymptotics_main(GG,HH,X,c,asy_var,alpha,N,numerical)
            subexp_parts.append(asy[2])
        exp_scale= asy[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-1]
        beta= copy(alpha)
        while 0 in [diff(h,x).subs(c) for h in H]:
            X.pop()
            X.insert(0,x)
            x= X[d-1]
            a= beta.pop()
            beta.insert(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,1,A
        elif r==1:
            # So 1 = r < d, and we have a smooth point.
            return self._asymptotics_main_smooth(G,H[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", submitted.

        AUTHORS:

        - Alex Raichev (2008-10-01, 2010-09-28)        
        """
        I= sqrt(-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-1)]
        T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-1)]
        Sstar= {}
        temp= self._crit_cone_combo(H,X,c,alpha) 
        for j in range(r-1):
            Sstar[S[j]]= temp[j]
        thetastar= dict([(t,0) for t in T])
        thetastar.update(Sstar)

        # Setup.
        print "Initializing auxiliary functions..."
        Hmul= mul(H)
        h= [function('h'+str(j),*tuple(X[:d-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-1] -1/h[j] for j in range(r)] )
        Gcheck= -G/U *mul( [-h[j]/X[d-1] for j in range(r)] )
        A= [(-1)^(r-1) *X[d-1]^(-r+j)*diff(Gcheck.subs({X[d-1]:1/X[d-1]}),X[d-1],j) for j in range(r)]
        e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-1)])
        ht= [hh.subs(e) for hh in h]
        Ht= [H[j].subs(e).subs({X[d-1]:1/ht[j]}) for j in range(r)]
        hsumt= add([S[j]*ht[j] for j in range(r-1)]) +(1-add(S))*ht[r-1]
        At= [AA.subs(e).subs({X[d-1]:hsumt}) for AA in A]
        Phit = -log(C[X[d-1]] *hsumt)+ I*add([alpha[i]/alpha[d-1] *T[i] for i in range(d-1)]) 
        # atC Stores h and U and all their derivatives evaluated at C.
        atC = copy(C)
        atC.update(dict( [(hh.subs(C),1/C[X[d-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-1,r]):
            s= solve( diff(H[j].subs({X[d-1]:1/h[j]}),X[i]), diff(h[j],X[i]) )[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[0:d-1],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-1] for j in range(r)]
        Hmulderivs= self._diff_all(Hmul,X,2*N-2+r,ending=end,sub_final=C)
        print "Calculating derivatives of auxiliary functions..."
        atC.update({U.subs(C):diff(Hmul,X[d-1],r).subs(C)/factorial(r)})
        Uderivs={}
        p= Hmul.polynomial(CC).degree()-r
        if p == 0:
            # Then, using a proposition at the end of [RaWi2011], we can
            # conclude that all higher derivatives of U are zero.
            for l in [1..2*N-2+m-1]:
                for s in UnorderedTuples(X,l):
                    Uderivs[diff(U,s).subs(C)] = 0
        elif p > 0 and p < 2*N-2+m-1:
            all_zero= True
            Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,[1..p],end,Uderivs,atC)
            # Check for a nonzero U derivative.
            if Uderivs.values() != [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 [p+1..2*N-2+m-1]:
                    for s in UnorderedTuples(X,l):
                        Uderivs.update({diff(U,s).subs(C):0})
            else:
                # Have to compute the rest of the derivatives.
                Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,[p+1..2*N-2+m-1],end,Uderivs,atC)
        else:
            Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,[1..2*N-2+m-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 -(1/2) *matrix([T+S]) *a *transpose(matrix([T+S]))
        Phitu= Phitu[0][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,2*N-2,sub=hderivs1,sub_final=[thetastar,atC],rekey=AA)
        BB= function('BB',*tuple(T+S))
        Phitu_derivs= self._diff_all(Phitu,T+S,2*N,sub=hderivs1,sub_final=[thetastar,atC],rekey=BB,zero_order=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([0..min(r-1,N-1)], [max(0,N-1-r)..N-1]):
            if j+k <= N-1: 
                L[(j,k)] = add([ \
                DD[(j,k,l)] /( (-1)^k *2^(k+l) *factorial(l) *factorial(k+l) ) \
                for l in [0..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())^(-1/2) *(2*pi)^((r-d)/2)   
        chunk= det *add([ (alpha[d-1]*n)^((r-d)/2-q) *add([ \
        L[(j,k)] *binomial(r-1,j) *stirling_number1(r-j,r+k-q) *(-1)^(q-j-k) \
        for (j,k) in CartesianProduct([0..min(r-1,q)], [max(0,q-r)..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[0].subs(c).n(digits=numerical)*n^co[1] \
            for co in coeffs] )
            exp_scale= (1/mul( [(C[X[i]]^alpha[i]).subs(c) for i in range(d)] )) \
            .n(digits=numerical)
        else:
            subexp_part = add( [co[0].subs(c)*n^co[1] for co in coeffs] )
            exp_scale= 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", submitted.

        AUTHORS:

        - Alex Raichev (2008-10-01, 2010-09-28)
        """
        I= sqrt(-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-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-1])
        A= Gcheck.subs({X[d-1]:1/h})/h
        T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-1)]
        e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-1)])
        ht= h.subs(e)
        Ht= H.subs(e).subs({X[d-1]:1/ht})
        At= A.subs(e)  
        Phit = -log(C[X[d-1]] *ht)\
        + I* add([alpha[i]/alpha[d-1] *T[i] for i in range(d-1)]) 
        Tstar= dict([(t,0) for t in T]) 
        atC = copy(C)
        atC.update({h.subs(C):1/C[X[d-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-1):
            s= solve( diff(H.subs({X[d-1]:1/h}),X[i]), diff(h,X[i]) )[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[0:d-1],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,2*N,ending=[X[d-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-1]).subs(C)})
        end= [X[d-1]]
        Hcheck= X[d-1] - 1/h
        p= H.polynomial(CC).degree()-1
        if p == 0:
            # Then, using a proposition at the end of [RaWi2011], we can
            # conclude that all higher derivatives of U are zero.
            for l in [1..2*N]:
                for s in UnorderedTuples(X,l):
                    Uderivs[diff(U,s).subs(C)] = 0
        elif p > 0 and p < 2*N:
            all_zero= True
            Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,[1..p],end,Uderivs,atC)
            # Check for a nonzero U derivative.
            if Uderivs.values() != [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 [p+1..2*N]:
                    for s in UnorderedTuples(X,l):
                        Uderivs.update({diff(U,s).subs(C):0})
            else:
                # Have to compute the rest of the derivatives.
                Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,[p+1..2*N],end,Uderivs,atC)
        else:
            Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,[1..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==2:
            # Compute v, the order of vanishing at Tstar of Phit. It is at least 2.
            v=2
            Phitderiv= diff(Phit,T[0],2) 
            splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify()
            while splat==0:
                v= v+1
                if v > 2*N:  # Then need to compute more derivatives of h for atC.
                    hderivs.update({diff(h,X[0],v) \
                    :diff(hderivs[diff(h,X[0],v-1)],X[0]).subs(hderivs1)})
                    atC.update({diff(h,X[0],v).subs(C) \
                    :hderivs[diff(h,X[0],v)].subs(atC)})
                Phitderiv= diff(Phitderiv,T[0])
                splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify()    
        if d==2 and v>2:     
            t= T[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(2)==0:
                At_derivs= self._diff_all(At,T,2*N-2, \
                sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
                Phitu_derivs= self._diff_all(Phitu,T,2*N-2+v, \
                sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+1,rekey=BB)
            else:
                At_derivs= self._diff_all(At,T,N-1,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
                Phitu_derivs= self._diff_all(Phitu,T,N-1+v,sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+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(2) == 0:
                for k in range(N):
                    L.append( add([ \
                    (-1)^l *gamma((2*k+v*l+1)/v) \
                    / (factorial(l) *factorial(2*k+v*l)) \
                    * DD[(k,l)] for l in [0..2*k] ]) )
                chunk= a^(-1/v) /(pi*v) *add([ alpha[d-1]^(-(2*k+1)/v) \
                * L[k] *n^(-(2*k+1)/v) for k in range(N) ])
            else:
                zeta= exp(I*pi/(2*v))
                for k in range(N):
                    L.append( add([ \
                    (-1)^l *gamma((k+v*l+1)/v) \
                    / (factorial(l) *factorial(k+v*l)) \
                    * (zeta^(k+v*l+1) +(-1)^(k+v*l)*zeta^(-(k+v*l+1))) \
                    * DD[(k,l)] for l in [0..k] ]) ) 
                chunk= abs(a)^(-1/v) /(2*pi*v) *add([ alpha[d-1]^(-(k+1)/v) \
                * L[k] *n^(-(k+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 -(1/2) *matrix([T]) *a *transpose(matrix([T]))
            Phitu= Phitu[0][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,2*N-2,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
            BB= function('BB',*tuple(T))
            Phitu_derivs= self._diff_all(Phitu,T,2*N,sub=hderivs1,sub_final=[Tstar,atC],rekey=BB,zero_order=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,1,N)
    
            # Plug above into asymptotic formula.
            L=[]
            for k in range(N):
                L.append( add([ \
                  DD[(0,k,l)] / ( (-1)^k *2^(l+k) *factorial(l) *factorial(l+k) ) \
                  for l in [0..2*k]]) )
            chunk= add([ (2*pi)^((1-d)/2) *a.determinant()^(-1/2) \
                  *alpha[d-1]^((1-d)/2 -k) *L[k] \
                  *n^((1-d)/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[0].subs(c).n(digits=numerical)*n^co[1] for co in coeffs] )
            exp_scale=(1/mul( [(C[X[i]]^alpha[i]).subs(c) for i in range(d)] )) \
            .n(digits=numerical) 
        else:
            subexp_part = add( [co[0].subs(c)*n^co[1] for co in coeffs] )
            exp_scale= 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[0][i]== alpha[i]/alpha[d-1] for i in range(d)]
        s= solve(eqns,V,solution_dict=True)[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=1
            F= G/H
        else:
            numer= G
            F= R(1)/H
        nstuff= 1/mul([X[j]^(alpha[j]*asy_var+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[0]*numer*nstuff] + a[1:] )

        # Do step (3).
        ce= decompose_via_cohomology(integrands)
        ce_new= []
        for a in ce:
            ce_new.append( [(a[0]/nstuff).simplify_full().collect(n)] + a[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", submitted.

        AUTHORS:

        - Alex Raichev (2010-08-25)
        """
        Hs= [SR(h[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 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=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 = 1
        else:
            seeds = [[v] for v in V]
            start = 2
        if singleton:
            for s in seeds:
                derivs[tuple(s)] = self._subs_all(diff(f[0],s),sub)
            for l in [start..n]:
                for t in UnorderedTuples(V,l):
                    s= tuple(t + ending)
                    derivs[s] = self._subs_all(diff(derivs[s[1:]],s[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 [start..n]:
                for t in UnorderedTuples(V,l):
                    s= tuple(t + ending)
                    value= self._subs_all(\
                    [diff(derivs[(j,)+s[1:]],s[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]= 0
            else:
                # Ignore the first of element of k, which is an index.
                for k in derivs.keys():
                    if len(k)-1 < zero_order:
                        derivs[k]= 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[0]],list(k)[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 [0..2*k]:
                    for s in UnorderedTuples(V,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 [0..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 == 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= 0
                    for t in P:
                        if product_derivs[(j,k,l)+self._diff_seq(V,t)] != 0:
                            MM= 1
                            for (a,b) in t:
                                MM= MM * M[a][b]
                                if a != b:
                                    MM= 2*MM
                            diffo= diffo + MM * product_derivs[(j,k,l)+self._diff_seq(V,t)]
                    DD[(j,k,l)] = (-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(-1)
        DD= {}
        if v.mod(2) == 0:
            for k in range(N):
                for l in [0..2*k]:
                    DD[(k,l)] = (a^(-1/v))^(2*k+v*l) \
                    * diff(A*B^l,x,2*k+v*l).subs(AB_derivs)
        else:
            for k in range(N):
                for l in [0..k]:
                    DD[(k,l)] = (abs(a)^(-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[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)-1    
        try:
            v[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", submitted.

        AUTHORS:

        - Alex Raichev (2011-04-18)
        """
        H= self._H
        d= self._d
        Hs= [SR(h[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] != [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())[1:]: # Subsets of size >= 1 
                k= len(list(S))
                M = matrix(list(S))
                if  rank(M) != 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 [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=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[0].variables()[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[-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]==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()[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()[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>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-1]*vars[i]*diff(f,vars[i]) \
        -alpha[i]*vars[d-1]*diff(f,vars[d-1]) for i in range(d-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[0],dict):
                return g[0]
            else:
                return g[0].simplify()
        elif singleton and not simplify:
            return g[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]=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 [2..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 ============================================================================