Ticket #10519: mgf.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 smooth
or multiple points.

For an introduction to this subject, see [PeWi2008]_.
The main algorithms and formulas implemented here come from [RaWi2008]_ 
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.

.. [RaWi2008]  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) : initial version
- Alex Raichev (2010-09-28) : corrected many functions
- Alex Raichev (2010-12-15) : updated documentation

EXAMPLES:

These come from [RaWi2008]_ and [RaWi2011]_.
A smooth point example. ::

    sage: R.<x,y>= PolynomialRing(QQ)
    sage: G= 1
    sage: H= 1-x-y-x*y
    sage: alpha= [3,2]
    sage: I= smooth_crit(H,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(),[SR(g) for g in R.gens()],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: asy= asymptotics(G,H,c,alpha,2)
    Calculating derivatives of implicit functions...
    Calculating derivatives of more implicit functions...
    Calculating actions of the second order differential operator...
    sage: relative_error(SR(G/H),[SR(g) for g in R.gens()],alpha,asy,[1,2,4,8,16])  # optional - slow
    [[1, 25.0000000000000, 24.9440713963148, 0.00223714414740672], [2, 1289.00000000000, 1288.35490191439, 0.000500463991936706], [4, 4.67334500000000e6, 4.67279936981582e6, 0.000116753670910463], [8, 8.52755090867210e13, 8.52731108840740e13, 0.0000281229941947458], [16, 3.97800011417700e28, 3.97797267728392e28, 6.89715743829034e-6]]
     
Another smooth point example.  
Turns out the second coefficient is 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: alpha= list(qq*vector([2,1-q])); alpha
    [4, 1]
    sage: I= smooth_crit(H,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(),[SR(g) for g in R.gens()],solution_dict=true); s
    [{y: 1, x: 1}]
    sage: c= s[0]
    sage: asy= asymptotics(G,H,c,alpha,2); asy
    Calculating derivatives of implicit functions...
    Calculating derivatives of more implicit functions...
    Calculating actions of the second order differential operator...
    1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*n^(1/3))
    sage: relative_error(SR(G/H),[SR(g) for g in R.gens()],alpha,asy,[1,2,4,8,16])  # optional - slow
    [[1, 0.187500000000000, 0.195379467542369, -0.0420238268926343], [2, 0.152343750000000, 0.155072786154872, -0.0179136732217213], [4, 0.122177124023438, 0.123081351936941, -0.00740095922809947], [8, 0.0973967181053013, 0.0976897337711845, -0.00300847576369478], [16, 0.0774425381638039, 0.0775363930774358, -0.00121192972050332]]
    
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: I= singular_points(H); 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: X= [SR(g) for g in R.gens()]
    sage: c= {X[0]: 1/2, X[1]: 1, X[2]: 4/3}
    sage: alpha= [8,3,3]
    sage: Hs= [SR(h[0]) for h in H.factor()];
    sage: Hs[0]= Hs[0]*H.factor().unit(); Hs
    [x*y + x - 1, 2*x^2*y*z + x^2*z - 1]
    sage: cc= critical_cone(Hs,X,c,coordinate=0); cc
    ([(1, 1/2, 0), (1, 1/3, 1/2)], 2-d cone in 3-d lattice N)    
    sage: alpha in cc[1]
    True
    sage: asy= asymptotics(G,H,c,alpha,1); asy
    Calculating derivatives of implicit functions...
    Calculating derivatives of more implicit functions...
    Calculating second-order differential operator actions...
    1/7*sqrt(3)*sqrt(7)*108^n/(sqrt(pi)*sqrt(n))
    sage: relative_error(SR(G/H),[SR(g) for g in R.gens()],alpha,asy,[1,2,4,8],scale=108)  # optional - VERY slow
    
    
Another multiple point example. ::

    sage: R.<x,y,z>= PolynomialRing(QQ)
    sage: G= 16
    sage: H= (4-2*x-y-z)*(4-x-2*y-z)
    sage: I= singular_points(H); 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: X= [SR(g) for g in R.gens()]
    sage: c= {X[0]:1, X[1]:1, X[2]:1}
    sage: alpha= [3,3,2]
    sage: Hs= [SR(h[0]) for h in H.factor()];
    sage: Hs[0]= Hs[0]*H.factor().unit(); Hs
    [x + 2*y + z - 4, 2*x + y + z - 4]
    sage: cc= critical_cone(Hs,X,c); cc
    ([(1, 2, 1), (2, 1, 1)], 2-d cone in 3-d lattice N)    
    sage: alpha in cc[1]
    True
    sage: asy= asymptotics(G,H,c,alpha,2); asy
    Calculating derivatives of implicit functions...
    Calculating derivatives of more implicit functions...
    Calculating second-order differential operator actions...
    4/3*sqrt(3)/(sqrt(pi)*sqrt(n)) - 25/216*sqrt(3)/(sqrt(pi)*n^(3/2))
    sage: relative_error(SR(G/H),[SR(g) for g in R.gens()],alpha,asy,[1,2,4,8,16])  # optional - slow
    [[1, 0.784973144531250, 1.18983759843026, -0.515768541534978], [2, 0.700524947606027, 0.881329983142157, -0.258099352712580], [4, 0.584773265395719, 0.637332211706702, -0.0898791880908176], [8, 0.448554766939249, 0.455660397364105, -0.0158411657808066], [16, 0.323752858725859, 0.323967782414798, -0.000663851092418399]]
      
"""
#*****************************************************************************
#       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/
#*****************************************************************************
#
# Keep functions names in alphabetical order below.
# A ==========================================================================
def algebraic_dep(Bs,leading_term=False):
    r"""
    This function returns an irreducible annihilating polynomial for the
    polynomials in `Bs`, if those polynomials are algebraically dependent.
    Otherwise it returns 0.
    
    INPUT:
 
    - ``Bs`` - A list of the form `[[B_1,e_1],\ldots,[B_r,e_r]]` where each `B_j` 
        is an element of a polynomial ring whose coefficient ring `R` is 
        contaiend in `\CC`.  
        Here `Bs` codes the polynomials `B_1^e_1,\ldots,B_l^e_l`.
    - ``leading_term`` - (default: False) Boolean.
    
    OUTPUT:
  
    If the polynomials coded by `Bs` are algebraically dependent, then the 
    output is an irreducible polynomial `g` in `R[T_1,...,T_r]` such that
    `g(B_1^e_1,...,B_r^e_r) = 0`.
    In this case, if `leading_term` is true, then the output is `[g,glt]`, 
    where `g` is as above and `glt` is the leading term of `g` in the 
    negdeglex order, that is, the monomial of `g` of smallest degree 
    (and largest lex order).
    This optional output is used in the function partial_frac() below.
    If the polynomials in `Bs` are algebraically independent, then the output
    is 0.
    
    EXAMPLES:

        sage: R= QQ['w,z']
        sage: w,z = R.gens()
        sage: Bs= [[w, 2],  [w^2 + z^2 - 1, 2],[w*z - 2, 2]]
        sage: g= algebraic_dep(Bs,leading_term=true)
        sage: g
        [256 - 256*T0 - 256*T2 + 352*T0^2 + 64*T0*T2 + 96*T2^2 - 208*T0^3
        - 32*T0^2*T1 - 48*T0^2*T2 + 16*T0*T2^2 - 16*T2^3 + 145*T0^4 + 16*T0^3*T1
        - 36*T0^3*T2 - 48*T0^2*T1*T2 - 10*T0^2*T2^2 - 4*T0*T2^3 + T2^4 - 52*T0^5
        - 18*T0^4*T1 + 28*T0^4*T2 + 4*T0^3*T1*T2 - 12*T0^3*T2^2 - 2*T0^2*T1*T2^2
        + 4*T0^2*T2^3 + 22*T0^6 + 4*T0^5*T1 - 12*T0^5*T2 + T0^4*T1^2
        - 4*T0^4*T1*T2 + 6*T0^4*T2^2 - 4*T0^7 - 2*T0^6*T1 + 4*T0^6*T2 + T0^8,
        256]
        sage: g[0]([B^e for B,e in Bs])
        0
        
        sage: R= QQ['x,y,z']
        sage: x,y,z = R.gens()
        sage: Bs= [[x*y,1],[y*z,1]]
        sage: g= algebraic_dep(Bs)
        sage: g
        0
    
    AUTHORS:

    - Alex Raichev (2008-10-01)
    """
    fs= [x[0] for x in Bs]
    es= [x[1] for x in Bs]
    l= len(fs)
    RR= fs[0].parent()
    R= RR.base_ring()
    vars= list(RR.gens())
    d= len(vars)
    str_vars= tuple([str(x) for x in vars])
    
    # Extend RR by 2*l new variables.
    Y= 'T'
    while Y in str(vars):
        Y= Y+'T'
    X= Y+'T'
    X_vars= [X + str(j) for j in range(l)]
    Y_vars= [Y + str(j) for j in range(l)]
    RRR= PolynomialRing(R,2*l+d,tuple(Y_vars+X_vars+vars))
    new_vars= list(RRR.gens())
    Y_vars= new_vars[0:l]
    X_vars= new_vars[l:2*l]
    vars=   new_vars[2*l:2*l+d]
    
    # First find an annihilating polynomial g for the fs.
    J= ideal([X_vars[j] -RRR(fs[j]) for j in range(l)])
    JJ= J.elimination_ideal(vars)
    g= JJ.gens()[0]
    
    # Now using g find an annihilating polynomial gg for the fs^es
    J= ideal([g] +[X_vars[j]^es[j] - Y_vars[j] for j in range(l)])
    JJ= J.elimination_ideal(X_vars)
    gg= JJ.gens()[0]
    
    # Shrink the ambient ring of gg to include only Y_vars and order the
    # monomials with negdeglex for the 'leading_term' option.
    RRR= PolynomialRing(R,l,tuple(Y_vars),order='negdeglex')
    gg= RRR(gg)
    if leading_term:
        return [gg,gg.lt()]
    return gg
#-------------------------------------------------------------------------------
def asymptotics(G,H,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`.
    It assumes that `F` is holomorphic in a neighborhood of the origin, 
    that `H` is a polynomial, and that `c` is a strictly minimal smooth 
    or multiple of `F` that is critical for `\alpha`.
    
    INPUT:

    - ``G`` - A function.
    - ``H`` - A `d`-variate polynomial ring element.
    - ``alpha`` - A `d`-tuple of positive natural numbers or, if `c` is a smooth
        point, possibly of symbols.
    - ``c`` - A dictionary with the symbolic variables of `G` and `H` as keys and 
        numbers as values.
    - ``N`` - A positive integer.
    - ``numerical`` - 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`` - Symbolic variable (default: None).
        The variable of the asymptotic expansion.
        If none is given, `var('n')` will be assigned.
        
    OUTPUT:

    The first `N` terms of the asymptotic expansion as of the Maclaurin 
    coefficients `F_{n\alpha}` of the function `F=G/H` as `n\to\infty`. 
        
    EXAMPLES:

    A smooth point example ::
    
        sage: R.<x,y>= PolynomialRing(QQ)
        sage: G= 1
        sage: H= 1-x-y-x*y
        sage: alpha= [3,2]
        sage: c= {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3}
        sage: asymptotics(G,H,c,alpha,2,numerical=3)
        Calculating derivatives of implicit functions...
        Calculating derivatives of more implicit functions...
        Calculating actions of the second order differential operator...
        (0.369/sqrt(n) - 0.0186/n^(3/2))*71.2^n
    
    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: c= {SR(z): 4/3, SR(y): 1, SR(x): 1/2}
        sage: alpha= [8,3,3]
        sage: asymptotics(G,H,c,alpha,1)
        Calculating derivatives of implicit functions...
        Calculating derivatives of more implicit functions...
        Calculating second-order differential operator actions...
        1/7*sqrt(3)*sqrt(7)*108^n/(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:
    
    .. [RaWi2008]  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)
    """
    # The variable for asymptotic expansions.
    if not asy_var:
        asy_var= var('n')
        
    # Create symbolic (non-ring) variables.
    R= H.parent()
    d= R.ngens()
    X= [SR(x) for x in R.gens()]

    # Do steps (1)--(3)
    new_integrands= cohom_integrand(G,H,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])
    asy= []
    for chunk in new_integrands:
        # Convert chunk into Symbolic Ring
        GG= chunk[0]
        HH= [SR(f) for (f,e) in chunk[1:]]
        asy.append(asymptotics_main(GG,HH,X,c,asy_var,alpha,N,numerical))
    # Do step (5).
    return add(asy)
#--------------------------------------------------------------------------------
def asymptotics_main(G,H,X,c,n,alpha,N,numerical):
    r"""
    This function is for internal use by the function asymptotics().
    It finds a variable in `X` to use to calculate asymptotics and decides whether
    to call the function asymptotics_main_smooth() or the function 
    asymptotics_main_multiple().
    
    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().
    
    EXAMPLES:
    
    sage: var('x,y,z,n')
    (x, y, z, n)
    sage: G= 1
    sage: H= [1-x*(1+y),1-z*x^2*(1+2*y)]
    sage: c= {z: 4/3, y: 1, x: 1/2}
    sage: alpha= [8,3,3]
    sage: asymptotics_main(G,H,[x,y,z],c,n,alpha,1,0)
    Calculating derivatives of implicit functions...
    Calculating derivatives of more implicit functions...
    Calculating second-order differential operator actions...
    1/7*sqrt(3)*sqrt(7)*108^n/(sqrt(pi)*sqrt(n))
    
    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 smooth or
    # multiple points.  
    # Search for good x in X from back to front (to be consistent with 
    # [RaWi2008]_ [RaWi2011]_ which use 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.
        return GG.subs(c) / jacobian(H,X).subs(c).determinant().abs()
    elif r==1:
        # So 1 = r < d, and we have a smooth point.
        return 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 asymptotics_main_multiple(G,H,X,c,n,beta,N,numerical)
#--------------------------------------------------------------------------------
def asymptotics_main_multiple(G,H,X,c,n,alpha,N,numerical):
    r"""
    This function is for internal use by the function asymptotics_main().
    It calculates asymptotics in case they depend upon multiple points.
    
    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().
    
    EXAMPLES:
    
    sage: var('x,y,z,n')
    (x, y, z, n)
    sage: G= 1
    sage: H= [1-x*(1+y), 1-z*x^2*(1+2*y)]
    sage: c= {z: 4/3, y: 1, x: 1/2}
    sage: alpha= [3,8,3]
    sage: asymptotics_main_multiple(G,H,[z,x,y],c,n,alpha,1,0)
    Calculating derivatives of implicit functions...
    Calculating derivatives of more implicit functions...
    Calculating second-order differential operator actions...
    1/7*sqrt(3)*sqrt(7)*108^n/(sqrt(pi)*sqrt(n))
    
    NOTES:
    
    The formula used for computing the asymptotic expansion is
    Theorem 3.4 of [RaWi2011]_.

    Currently this function cannot handle symbolic `c`, 
    because aux_multiple() 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(new_var_name('s',X) + str(j)) for j in range(r-1)]
    T= [var(new_var_name('t',X) + str(i)) for i in range(d-1)]
    Sstar= {}
    temp= aux_multiple(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.
    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 = 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= diff_all(Hmul,X,2*N-2+r,ending=end,sub_final=C)
    print "Calculating derivatives of implicit 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= 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= diff_prod(Hmulderivs,U,Hcheck,X,[p+1..2*N-2+m-1],end,Uderivs,atC)
    else:
        Uderivs= 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 implicit functions..."
    AA= [function('A'+str(j),*tuple(T+S)) for j in range(r)]
    At_derivs= diff_all(At,T+S,2*N-2,sub=hderivs1,sub_final=[thetastar,atC],rekey=AA)
    BB= function('BB',*tuple(T+S))
    Phitu_derivs= 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= 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...
        sub_exp = add( [co[0].subs(c).n(digits=numerical)*n^co[1] \
        for co in coeffs] )
        return (1/mul( [(C[X[i]]^alpha[i]).subs(c) for i in range(d)] )) \
        .n(digits=numerical)^n *sub_exp
    else:
        sub_exp = add( [co[0].subs(c)*n^co[1] for co in coeffs] )
        return (1/mul( [(C[X[i]]^alpha[i]).subs(c) for i in range(d)] ))^n \
        *sub_exp    
#--------------------------------------------------------------------------------
def asymptotics_main_smooth(G,H,X,c,n,alpha,N,numerical):
    r"""
    This function is for internal use by the function asymptotics_main().
    It calculates asymptotics in case they depend upon smooth points.
    
    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().
    
    EXAMPLES:
    
    sage: var('x,y,n')
    (x, y, n)
    sage: G= 1
    sage: H= 1-x-y-x*y
    sage: alpha= [3,2]
    sage: c= {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3}
    sage: asymptotics_main_smooth(G,H,[x,y],c,n,alpha,1,0)
    Calculating derivatives of implicit functions...
    Calculating derivatives of more implicit functions...
    Calculating actions of the second order differential operator...
    3*(108/((sqrt(13) - 3)^2*(sqrt(13) - 2)^3))^n/((sqrt(13) - 3)*(sqrt(13) + 1)*sqrt(-(sqrt(13) - 3)^2*(sqrt(13) - 2)^2*(1/(sqrt(13) - 3) + 2/(sqrt(13) - 3)^2)^2/(sqrt(13) + 1)^2 + (sqrt(13) - 3)*(sqrt(13) - 2)^2*(((1/(sqrt(13) - 3) + 2/(sqrt(13) - 3)^2)/(sqrt(13) + 1) + 4*(1/(sqrt(13) - 3) + 2/(sqrt(13) - 3)^2)/((sqrt(13) - 3)*(sqrt(13) + 1)))/(sqrt(13) + 1) - (1/(sqrt(13) - 3) + 2/(sqrt(13) - 3)^2)/(sqrt(13) + 1)^2) + (sqrt(13) - 3)*(sqrt(13) - 2)*(1/(sqrt(13) - 3) + 2/(sqrt(13) - 3)^2)/(sqrt(13) + 1))*sqrt(pi)*sqrt(n))
    
    NOTES:

    The formulas used for computing the asymptotic expansions are
    Theorems 3.2 and 3.3 [RaWi2008]_ with `p=1`.
    Theorem 3.2 is a specialization of Theorem 3.4 of [RaWi2011]_
    with `r=1`.
    
    REFERENCES:
    
    .. [RaWi2008]  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.
    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(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 = 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= diff_all(H,X,2*N,ending=[X[d-1]],sub_final=C)
    print "Calculating derivatives of implicit 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= 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= diff_prod(Hderivs,U,Hcheck,X,[p+1..2*N],end,Uderivs,atC)
    else:
        Uderivs= 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 implicit functions..."
        AA= function('AA',t)
        BB= function('BB',t)
        if v.mod(2)==0:
            At_derivs= diff_all(At,T,2*N-2, \
            sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
            Phitu_derivs= diff_all(Phitu,T,2*N-2+v, \
            sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+1,rekey=BB)
        else:
            At_derivs= diff_all(At,T,N-1,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
            Phitu_derivs= 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= 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 implicit functions..."
        AA= function('AA',*tuple(T))
        At_derivs= diff_all(At,T,2*N-2,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA)
        BB= function('BB',*tuple(T))
        Phitu_derivs= 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= 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...
        sub_exp = add( [co[0].subs(c).n(digits=numerical)*n^co[1] for co in coeffs] )
        return (1/mul( [(C[X[i]]^alpha[i]).subs(c) for i in range(d)] )) \
        .n(digits=numerical)^n *sub_exp
    else:
        sub_exp = add( [co[0].subs(c)*n^co[1] for co in coeffs] )
        return (1/mul( [(C[X[i]]^alpha[i]).subs(c) for i in range(d)] ))^n \
        *sub_exp       
#-----------------------------------------------------------------------------
def aux_multiple(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 the function  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 G = a` for `y`, where `G` 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: 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: aux_multiple(fs,R.gens(),c,alpha)
        (0, 1)
    
    NOTES:

    Use this function only when `G` is well-defined and there is a unique
    solution to the matrix equation `y G = 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)
    Gamma= matrix([direction(log_grad(f,X,c)) for f in fs])
    s= Gamma.solve_left(vector(alpha))
    # The following is a workaround to a Sage 4.6 bug.
    return s/alpha[d-1]
# B ==========================================================================
# C ==========================================================================
def clean_up(As_over_Bs):
    r"""
    This function returns a simplified version of the fractions in `A_over_Bs.`
    
    INPUT:

    - ``As_over_Bs`` -  A list of the form `[chunk_1,\ldots,chunk_r]`, where
      each `chunk_j` is of the form `[A,[B_1,e_1],\ldots,[B_l,e_l]]`, where
      `A` is an expression, `B_1,\ldots,B_l` are elements of some 
      polynomial ring `R`, `e_1,\ldots,e_l` are positive integers, 
      and `A` is some expression in the variables of `R`.  
      Here `[A,[B_1,e_1],\ldots,[B_l,e_l]]` represents 
      `A/(B_1^e_1 *... * B_l^e_l)`.
    
    OUTPUT:

    A list of the same form as `As_over_Bs` but with each fraction written in
    lowest terms and fractions with identical denominators combined.
    
    EXAMPLES:

        sage: R= QQ['w,z']
        sage: w,z= R.gens()
        sage: As_over_Bs= [ [1,[w,1]], [z^2,[w*z-1,1]], [-1,[w,1]], [w,[w,1], \
              [w*z-1,1]] ]
        sage: clean_up(As_over_Bs)
        [[z^2 + 1, [w*z - 1, 1]]]
    
    AUTHORS:

    - Alex Raichev (2008-10-01)
    """
    # Step 1: write each fraction in As_over_Bs in lowest terms and store
    # results in 'chunks' with denominators listed first in preparation for
    # step 2.
    if As_over_Bs == []:
        return As_over_Bs
    chunks= []
    R= As_over_Bs[0][1][0].parent()
    for x in As_over_Bs:
        A=  x[0]
        Bs= x[1:]
        if not A in Frac(R):
            chunks.extend([ Bs + [A] ])
            continue
        # At this point we know we can use quo_rem on numerator and denominator
        new_Bs= []
        if A in R:
            A= R(A)  # Force coercion.
            for (B,e) in Bs:
                quo,rem = A.quo_rem(B)
                while rem == 0 and e > 0:
                    A= quo
                    e= e -1
                    quo,rem = A.quo_rem(B)
                if e != 0:
                    new_Bs.append([B,e])
        else:
            # 'A' is a fraction so cancel from its numerator only.
            A= Frac(R)(A)
            for (B,e) in Bs:
                quo,rem = A.numerator().quo_rem(B)
                while rem == 0 and e > 0:
                    A= quo/A.denominator()
                    e= e -1
                    quo,rem = A.numerator().quo_rem(B)
                if e != 0:
                    new_Bs.append([B,e])
        chunks.extend([ new_Bs +[A] ])
    
    # Step 2: sort fractions in 'chunks' by denominator so that fractions with
    # equal denominators are adjactent.
    # Then sum fractions in 'chunks' with equal denominators and store results in
    # 'new_chunks'.
    chunks.sort()
    new_chunks= []
    j= 0
    L= len(chunks)
    while j < L:
        l= len(chunks[j])
        A= chunks[j][l-1]
        Bs= chunks[j][0:l-1]
        k= j+1
        
        # Look ahead to find fractions with denominators == 'Bs' and add their
        # numerators to 'A'.
        while k < L:
            next_l= len(chunks[k])
            next_A= chunks[k][next_l-1]
            next_Bs= chunks[k][0:next_l-1]
            if Bs == next_Bs:
                A= A +next_A
                k= k+1
            else:
                break
        
        # Update 'new_chunks', which contains fractions with distinct
        # denominators.
        if A != 0:
            new_chunks.extend([ [A] +Bs ])
        j= k
    return new_chunks
#-------------------------------------------------------------------------------
def cohom_equiv(As_over_Bs):
    r"""
    Given each differential form `(A/B) dx_1 \wedge\cdots\wedge dx_d`
    represented in `As_over_Bs`, this function returns a form equivalent in De
    Rham cohomology that has no repeated factors in the denominator.
    
    INPUT:

    - ``As_over_Bs`` - A list of the form `[chunk_1,\ldots,chunk_r]`, where each
        `chunk_j` has the form `[A,[B_1,e_1],\ldots,[B_l,e_l]]`,
        `B_1,\ldots,B_l` are irreducible elements of a predefined polynomial
        ring `R` such that their corresponding algebraic varieties
        `\{x\in\CC^d : B_j(x)=0\}` intersect transversely,
        `e_1,\ldots,e_l` are positive integers, `l \le d`, and
        `A` is a function in some of the indeterminates of `R`.  Here
        `[A,[B_1,e_1],\ldots,[B_l,e_l]]` represents the fraction
        `A/(B_1^e_1 \cdots B_l^e_l)`.
    - ``eBQ`` - An embedding from the base ring of `R` to `\QQbar`.
    
    OUTPUT:

    A list of the form `[chunk_1,\ldots,chunk_s]`, where each `chunk_j` has
    the form `[A',[B'_1,e'_1],\ldots,[B'_k,e'_k]]` as above, where each
    `B'_j` is one of `B_1,\ldots,B_l`.  Moreover the sum of the fractions
    represented by `chunk_1,\ldots,chunk_r` is equivalent in De Rham
    cohomology to the sum of the fractions represented by `As_over_Bs`.

    EXAMPLES:

        sage: R.<w,z>= PolynomialRing(QQ)
        sage: cohom_equiv([[ 1, [w,3] ]])
        []
        sage: cohom_equiv([[ 1, [w,3], [w*z-1,2] ]])
        [[-3*z^2, [w, 1], [w*z - 1, 1]], [-5*z^3, [w*z - 1, 1]]]
    
    NOTES:
    
    This is a recursive function and stops calling itself when all the
    `e_j` equal 1.  
    It does most of the work for the function cohom_equiv().  
    The algorithm used here is that of Theorem 17.4 of
    [AiYu1983]_.  
    The varieties corresponding to `B_1,\ldots,B_l` 
    __intersect transversely__ iff for each point `c` of their intersection
    and for all `k \le l`, the Jacobian matrix of any `k` polynomials from
    `\{B_1,\ldots,B_l\}` has rank equal to `\min\{k,d\}` when evaluated at
    `c`.
    
    The embedding `eBQ` is used to coerce the coefficients of the numerators
    of the output of cohom_equiv() into `\QQbar`.
    
    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)
    """
    import copy     # Will need this to make copies of a nested list.
    if As_over_Bs == []:
        return As_over_Bs
    R= As_over_Bs[0][1][0].parent()
    V= list(R.gens())
    d= R.ngens()
    all_reduced= true
    chunks= []
    for x in As_over_Bs:
        A=  x[0]
        Bs= x[1:]
        l= len(Bs)
        if A==0:
            # Then do not include this element in final list.
            continue
        if add([e for B,e in Bs]) == l:
            # Then nothing to be done.
            chunks.append([A] +Bs)
            continue
        # At this point we know some reducing needs to be done.
        all_reduced= false
        dets= []
        vars_subsets= Set(V).subsets(l)
        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([B for B,e in Bs],v)
            dets.append(jac.determinant())
        # Get a Nullstellensatz certificate for J.
        J= [B for B,e in Bs] + dets
        L= R(1).lift(ideal(J))
        for j in range(l):
           if L[j] != 0:
                # Make a copy of (and not a reference to) the nested list Bs.
                new_Bs = copy.deepcopy(Bs)
                if new_Bs[j][1] > 1:
                    new_Bs[j][1]= new_Bs[j][1] -1
                else:
                    del new_Bs[j]
                chunks.append( [A*L[j]] +new_Bs )
        for k in range(vars_subsets.cardinality()):
           if L[l+k] != 0:
                new_Bs = copy.deepcopy(Bs)
                for j in range(l):
                    if new_Bs[j][1] > 1:
                        new_Bs[j][1]= new_Bs[j][1] -1
                        v=  sorted(vars_subsets[k],reverse=true)
                        jac= jacobian([SR(A*L[l+k])] +[ SR(Bs[jj][0]) for \
                        jj in range(l) if jj !=j], [SR(vv) for vv in v])
                        det= jac.determinant()
                        if det != 0:
                            chunks.append([perm_sign(v,V) \
                              *(-1)^j/new_Bs[j][1] *det] +new_Bs)
                        break
    if all_reduced:
        return clean_up(chunks)
    return cohom_equiv(chunks)    
#-------------------------------------------------------------------------------
def cohom_integrand(G,H,alpha,asy_var=None):
    r"""
    This function takes an integral and breaks it up into a list of 
    nicer integrals for the purposes of computing asymptotics of the original 
    integral.
    The sum of the nicer integrals is de Rham cohomologous to the original
    integral.
    
    INPUT:
    
    - ``G`` - A symbolic expression.
    - ``H`` - An element of a polynomial ring.
    - ``alpha`` - A list of positive integers or symbolic variables.
    - ``asy_var`` - A symbolic variable (default: None).
        Eventually set to `var('n')` is None is given.
    
    OUTPUT:

    A list of the form `[chunk_1,\ldots,chunk_r]`, where each
    `chunk_j` has the form `[A,[B_1,1],\ldots,[B_l,1]]`,
    `B_1,\ldots,B_l` are irreducible elements of a predefined polynomial
    ring `R` such that their corresponding algebraic varieties
    `\{x\in\CC^d : B_j(x)=0\}` intersect transversely, `l \le d`, and
    `A` is a function in some of the indeterminates of `R` and `asy_var`. 
    Here `[A,[B_1,1],\ldots,[B_l,1]]` represents the fraction
    `A/(B_1 \cdots B_l)`.
    
    EXAMPLES:

    sage: R.<x,y,z>= PolynomialRing(QQ)
    sage: G= 16
    sage: H= (4-2*x-y-z)^2*(4-x-2*y-z)
    sage: c= {SR(x):1, SR(y):1, SR(z):1}
    sage: alpha= [3,3,2]
    sage: cohom_integrand(G,H,alpha)
    [[16*(4*x^(-3*n - 1)*y^(-3*n - 1)*z^(-2*n - 2) - 3*x^(-3*n - 1)*y^(-3*n - 2)*z^(-2*n - 1))*n*x^(3*n + 1)*y^(3*n + 1)*z^(2*n + 1) + 16*(2*x^(-3*n - 1)*y^(-3*n - 1)*z^(-2*n - 2) - x^(-3*n - 1)*y^(-3*n - 2)*z^(-2*n - 1))*x^(3*n + 1)*y^(3*n + 1)*z^(2*n + 1), [x + 2*y + z - 4, 1], [2*x + y + z - 4, 1]]]
    
    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) 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')
    
    # Put G,H into useable format.
    G_over_Hs= format(G,H)
    
    # Create symbolic (non-ring) variables.
    R= H.parent()
    d= R.ngens()
    X= [SR(x) for x in R.gens()]
    n_part=  1/mul([X[j]^(alpha[j]*asy_var+1) for j in range(d)])
    
    # Do steps (2) and (1).
    pf= partial_frac(G_over_Hs)  
    integrand= []
    for chunk in pf:
        integrand.append( [chunk[0]*n_part] + chunk[1:] )
       
    # Do step (3).
    ce= cohom_equiv(integrand)
    ce_new= []
    for chunk in ce:
        ce_new.append( [(chunk[0]/n_part).simplify().collect(n)] + chunk[1:] )
    return ce_new
#-------------------------------------------------------------------------------
def critical_cone(fs,X,c,coordinate=None):
    r"""
    This function returns the convex hull of vectors derived from `fs` with
    respect to the variables `X` and evaluated at `c`.

    INPUT:

    - ``fs`` - A list of symbolic expressions in the variables `X`.
    - ``X`` - A list of variables.
    - ``c`` - A dictionary with keys `X` and values in a field.

    OUTPUT:

    The directions of the logarithmic gradients of the 
    functions in `fs` evaluated at `c` and their convex hull. 
    See the function direction() above.

    EXAMPLES:

        sage: var('x,y,z')
        (x, y, z)
        sage: fs= [x + 2*y + z - 4, 2*x + y + z - 4]
        sage: a= {z: 1, y: 1, x: 1}
        sage: print critical_cone(fs,[x,y,z],a)
        ([(1, 2, 1), (2, 1, 1)], 2-d cone in 3-d lattice N)

    AUTHORS:

    - Alex Raichev (2010-08-25)
    """
    vecs= [direction(log_grad(f,X,c),coordinate) for f in fs]
    return vecs,Cone(vecs)
# D ============================================================================
def diff_all(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 `X`. 
    The default is to key the dictionary by all nondecreasing sequences
    in `V` of length 1 up to length `n`.
    For internal use.
    
    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: f= function('f',x)
        sage: dd= 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= diff_all(f,[x],3,sub=d1)
        sage: dd[(x,x,x)]
        24*x
        sage: dd= diff_all(f,[x],3,sub=d1,rekey=f)
        sage: dd[diff(f,x,3)]
        24*x
        sage: a= {x:1}
        sage: dd= 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= 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= 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: F= function('F',*X)
        sage: dd= diff_all(f,X,2,rekey=F)
        sage: dd[diff(F,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)] = subs_all(diff(f[0],s),sub)
        for l in [start..n]:
            for t in UnorderedTuples(V,l):
                s= tuple(t + ending)
                derivs[s] = 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= 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= 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] = 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_implicit(F,y,f_name,vars,v):
    r"""
    This function assumes `F=0` defines `y` as an implicit function `f` of the
    variables `vars`.  
    It returns `f` and its derivative with respect to `v`.
    
    INPUT:

    - ``F`` - A function.
    - ``y`` - A symbolic variable occurring in `F`.
    - ``f_name`` - A string.
    - ``vars`` - A list of some symbolic variables different from `y` that occur -
        in `F`.
    - ``v`` - A symbolic variable occurring in `vars`.

    OUTPUT:

    The pair `(f,\partial f/\partial t)`, where `f` is a callable
    symbolic function with name `f_name` and arguments `vars`.
    Here `f(vars) = y` is defined implicitly by `F=0`.

    EXAMPLES:
    
        sage: var('a,b,c')
        (a, b, c)
        sage: F= a^2+b^2+c^2-1
        sage: v= [a,b]
        sage: diff_implicit(F,c,'f',v,a)
        (f(a, b), -a/f(a, b))
    
    WARNING:
    
    Choosing `\code{f_name = str(y)}` could lead to confusion.
    
    AUTHORS:

    - Alex Raichev (2008-10-01)
    """
    f = function(f_name,*tuple(vars))
    # Replace y by f in F.
    # Use ** trickery since we don't know the name of the symbolic variable
    # that y references.
    FF= F.substitute(**{str(y):f})
    s= solve(diff(FF,v) == 0, diff(f, v, 1),solution_dict=True)
    return f, s[0][diff(f, v, 1)]
#-------------------------------------------------------------------------------
def diff_op(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().
    
    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: 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= 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)+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)+diff_seq(V,t)]
                DD[(j,k,l)] = (-1)^(k+l)*diffo
    return DD
#-------------------------------------------------------------------------------
def diff_op_simple(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().
    
    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: A= function('A',x)
        sage: B= function('B',x)
        sage: AB_derivs= {}
        sage: 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(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().
    
    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: u= function('u',x)
        sage: g= function('g',x)
        sage: dd= 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)] = 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(subs_all(D,sol[0]))
    return uderivs
#-------------------------------------------------------------------------------
def diff_seq(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().
    
    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: V= list(var('x,t,z'))
        sage: 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(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.

    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: direction([2,3,5])
        (2/5, 3/5, 1)

        sage: direction([2,3,5],0)
        (1, 3/2, 5/2)

    AUTHORS:

    - Alex Raichev (2010-08-25)
    """
    if coordinate == None:
        coordinate= len(v)-1
    return tuple([SR(vv/v[coordinate]).simplify() for vv in v])
# E ============================================================================
# F ============================================================================
def format(A,B):
    r"""
    This functions puts `A/B` into input format used by the function 
    cohom_equiv() and related functions.
    
    INPUT:

    - ``A`` - An expression.
    - ``B`` - An element of a polynomial ring UFD `R`.
    
    OUTPUT:

    A list `[ [A/u, [B_1,e_1],\ldots,[B_r,e_r] ]`, where `u` is a unit and
    `B= u B_1^{e_1} \cdots B_r^{e_r}` is the unique factorization of 
    `B` in `R`.
    
    EXAMPLES:

        sage: R.<x,y,z>= PolynomialRing(QQ)
        sage: G= x^2 +y^3
        sage: H= y^3*(1-y*z)*(y^2+z^2-1)^2
        sage: format(G,H)
        [[-y^3 - x^2, [y, 3], [y*z - 1, 1], [y^2 + z^2 - 1, 2]]]
        
    AUTHORS:

    - Alex Raichev (2008-10-01)
    """
    fac= B.factor()
    return [ [A/fac.unit()] +[[f,e] for f,e in fac] ]
# G ============================================================================
# H ============================================================================
# I ============================================================================
# J ============================================================================
# K ============================================================================
# L ============================================================================
def log_grad(f,X,c):
    r"""
    This function returns the logarithmic gradient of `f` with respect to the
    variables of `X` evalutated at `c`.
    
    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: 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: 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: 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(L):
    r"""
    This function converts the list `L` to a string (without commas) and returns
    the string as a variable.
    It is for internal use by the function diff_op()
    
    INPUT:

    - ``L`` - A list.
        
    OUTPUT:

    A variable whose name is the concatenation of the variable names in `L`.
    
    EXAMPLES:

        sage: L= list(var('x,y,hello'))
        sage: v= 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(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().
    
    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: X= var('x,xx,y,z')
        sage: 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 ============================================================================
def partial_frac(A_over_Bs):
    r"""
    This function returns the partial fraction decomposition of the fractional
    expression represented by `A_over_Bs`.
    
    INPUT:

    - ``As_over_Bs`` - List of the form `[chunk_1,\ldots,chunk_r]`, where each
        `chunk_j` has the form `[A,[B_1,e_1],\ldots,[B_l,e_l]]`, where
        `A, B_1,\ldots,B_l` are elements of the predefined polynomial ring `R`
        and `e_1,\ldots,e_l` are positive integers.  Here
        `[A,[B_1,e_1],\ldots,[B_l,e_l]]` represents the fraction
        `A/(B_1^e_1 \cdots B_l^e_l)`.
    - ``eBQ`` - An embedding from the base ring of `R` into `\QQbar`.
    
    OUTPUT:

    A list of the form `[chunk_1,\ldots,chunk_r]`, where each `chunk_j` has
    the form `[A',[B'_1,e'_1],\ldots,[B'_k,e'_k]]` as above, each `B'_j` is
    one of `B_1,\ldots,B_l`, and `k \le d`.  
    Moreover the sum of the fractions represented  by `chunk_1,\ldots,chunk_r` 
    equals the fraction represented by A_over_Bs.
    
    EXAMPLES:

        sage: R= QQ['w,z']
        sage: w,z = R.gens()
        sage: A_over_Bs= [[1,[w,2], [w*z-1,1], [w^2+z^2-1,2],[w*z-2,2]]]
        sage: pf= partial_frac_main(A_over_Bs,R,len(R.gens()))
        sage: if unformat(A_over_Bs) != unformat(pf): print "fail" \n

                  
        sage: R= QQ['w,z']
        sage: w,z = R.gens()
        sage: Bs= [[w,2], [w*z-1,1], [w^2+z^2-1,2],[w*z-2,2]]
        sage: A_over_Bs= [[exp(w)/z] +Bs]
        sage: pf= partial_frac(A_over_Bs)
        sage: if unformat(A_over_Bs) != unformat(pf): print "fail" \n

            
    NOTES:

    This function calls partial_frac_main() to do most of the work.
    
    The embedding `eBQ` is used to coerce the coefficients of the numerator
    of the output of partial_frac() into `\QQbar`.
    
    AUTHORS:

    - Alex Raichev (2008-10-01)
    """
    
    if A_over_Bs == []:
        return A_over_Bs
    A= A_over_Bs[0][0]
    Bs= A_over_Bs[0][1:]
    R= Bs[0][0].parent()
    d= len(R.gens())
    
    # It is possible that A is not in Frac(R).
    # Since partial_frac_main() can not deal with this scenario, decompose
    # [[1,Bs]] into partial fractions and then multiply result by A.
    pf= partial_frac_main([ [R(1)] +Bs],R,d)
    for x in pf:
        x[0]= A*x[0]
    return pf
#-------------------------------------------------------------------------------
def partial_frac_main(As_over_Bs,R,d):
    r"""
    This function returns the partial fraction decomposition of the fractional
    expression coded by `As_over_Bs`.
    
    INPUT:

    - ``As_over_Bs`` - List of the form `[chunk_1,\ldots,chunk_r]`, where each
      `chunk_j` has the form `[A,[B_1,e_1],\ldots,[B_l,e_l]]`, where
      `A, B_1,\ldots,B_l` are elements of the predefined polynomial ring `R`
      and `e_1,\ldots,e_l` are positive integers.  Here
      `[A,[B_1,e_1],\ldots,[B_l,e_l]]` represents the fraction
      `A/(B_1^e_1 \cdots B_l^e_l)`.
    - ``R`` - A polynomial ring.
    - ``d`` - The number of indeterminates of `R`.

    OUTPUT:

    A list of the form `[chunk_1,\ldots,chunk_r]`, where each `chunk_j` has
    the form `[A',[B'_1,e'_1],\ldots,[B'_k,e'_k]]` as above, each `B'_j` is
    one of `B_1,\ldots,B_l`.  
    Moreover the sum of the fractions represented
    by `chunk_1,\ldots,chunk_r` equals the fraction represented by
    `A_over_Bs`.  Stops recursively calling itself when `k \le d`.
    
    EXAMPLES:

        sage: R= QQ['w,z']
        sage: w,z = R.gens()
        sage: A_over_Bs= [[1,[w,2], [w*z-1,1], [w^2+z^2-1,2],[w*z-2,2]]]
        sage: pf= partial_frac_main(A_over_Bs,R,len(R.gens()))
        sage: unformat(A_over_Bs)==unformat(pf)
        True
                  
    NOTES:
    
    This is a recursive function and does most of the work for the function
    partial_frac().  
    The algorithm used here is that of Theorem 1 [Lein1978]_.
    
    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 (2008-10-01)
    """
    # First reduce all fractions.
    As_over_Bs= clean_up(As_over_Bs)
    if As_over_Bs == []:
        return As_over_Bs
    
    # Proceed with main loop.
    chunks= []
    all_reduced= true
    for x in As_over_Bs:
        A=  x[0]
        Bs= x[1:]
        l= len(Bs)
        # If the number of terms in the denominator is less than or equal to the
        # dimension then no partial fractionizing needed.
        if l <= d:
            chunks.append([A]+Bs)
            continue    # Go to next iteration of loop.
        all_reduced= false
        g,glt= algebraic_dep(Bs,leading_term=true)
        new_vars= g.parent().gens()
        gg= (glt -g)/(glt.lc())
        numer= map(mul,zip(gg.coefficients(),gg.monomials()))
        subby= [new_vars[j]^Bs[j][1] for j in range(l)]
        numer= [x(subby) for x in numer]
        e= list(glt.exponents())[0:l]
        denom= [[new_vars[j],(Bs[j][1]*(e[0][j]+1))] for j in range(l)]
        temp_chunks= []
        for y in numer:
            temp_chunks.append([y]+denom)
        temp_chunks= clean_up(temp_chunks)
        L= [B for B,e in Bs]
        for y in temp_chunks:
            y[0]= A*R(y[0](L))
            for yy in y[1:]:
                yy[0]= R(yy[0](L))
        chunks= chunks +temp_chunks
    if all_reduced:
        return chunks
    return partial_frac_main(chunks,R,d)
#-------------------------------------------------------------------------------
def perm_sign(v,vars):
    r"""
    This function returns the sign of the permutation on `1,\ldots,len(vars)`
    that is induced by the subtlist `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: perm_sign(v,vars)
        -1
        sage: v= ['d','b']
        sage: perm_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()
# Q ============================================================================
# R ============================================================================
def relative_error(F,V,alpha,approx,interval,scale=1):
    r"""
    This function returns the relative error of the approximate values in
    `approx` compared to the true values of the Maclaurin coefficients of `F`
    in the direction `m alpha` for `m` in `interval`.
    
    INPUT:

    - ``F`` - A symbolic expression.
    - ``V`` - A list of variables occurring in F.
    - ``alpha`` - A list of positive integers of length the length of `V`.
    - ``approx`` - An individual or list of symbolic expressions in one variable.
    - ``interval`` - A list of positive integers.
    - ``scale`` - An integer (default: 1).
    
    OUTPUT:

    A list whose entries are of the form
    `[m,a,b_0,\dots,b_k,err_0,\ldots,err_k]` for `m \in interval`.
    Here `a` is the `m alpha` coefficient (in multi-index notation) of the
    Maclaurin series for `F` divided by `scale^m`;
    `b_i` is the value of `approx[i]` evaluated at `m` (for `i=0,\ldots,k`)
    divided by `scale^m`;
    `err_i` is (a-b_i)/a (for `i=0,\ldots,k`).
    All outputs are decimal approximations.
    
    EXAMPLES:

        sage: var('x,y,n')
        (x, y, n)
        sage: relative_error(1/(1-x-y),[x,y],[1,1],4^n/sqrt(pi*n),[1,2,4,8,16])     
        [[1, 2.00000000000000, 2.25675833419103, -0.128379167095513], [2, 6.00000000000000, 6.38307648642292, -0.0638460810704871], [4, 70.0000000000000, 72.2162666941128, -0.0316609527730400], [8, 12870.0000000000, 13072.5406441941, -0.0157374237913090], [16, 6.01080390000000e8, 6.05793952520368e8, -0.00784181716586718]]
    
    AUTHORS:

    - Alex Raichev (2009.05.18)
    """
    if not isinstance(approx,list):
        approx= [approx]
    av= approx[0].variables()[0]
    # Scale approx according to the value of `scale`.
    approx= [app/scale^av for app in approx]
    interval.sort()
    d= len(V)  # = len(alpha)
    p= dict([(x,0) for x in V])
    stats= []
    s= []
    for i in range(d):
        s.extend([V[i] for j in range(alpha[i])])
    # Store last F derivative taken to calculate next F derivative.
    F_deriv= diff(F,s)
    old_beta= alpha
    for m in interval:
        beta= [m*a for a in alpha]
        delta= [beta[i]-old_beta[i] for i in range(len(beta))]
        s= []
        for i in range(d):
            s.extend([V[i] for j in range(delta[i])])
        F_deriv= diff(F_deriv,s)
        a= F_deriv.subs(p)/ (mul([factorial(b) for b in beta]) * scale^m)
        b= [f.subs({av:m}) for f in approx]
        stats_row= [m, a.n()]
        stats_row.extend([bb.n() for bb in b])
        stats_row.extend([((a-bb)/a).n() for bb in b])
        stats.append(stats_row)
        old_beta= beta
    return stats
# S ============================================================================
def singular_points(f):
    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 : f(x)=0\}`.
    
    INPUT:

    - ``f`` - An element of a polynomial ring over an algebraic number field.
    
    OUTPUT:

    A Groebner basis ideal whose variety is the set of singular points of 
    the algebraic variety `V= \{x\in\CC^d : f(x)=0\}`.
    
    EXAMPLES:

        sage: R.<x,y,z>= PolynomialRing(QQ)
        sage: H= (4-2*x-y-z)*(4-x-2*y-z)
        sage: singular_points(H)
        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)
    """
    R= f.parent()
    fred= R.ideal(f).radical().gens()[0]     # Compute the reduction of f.
    J= R.ideal([fred] + fred.gradient())
    return R.ideal(J.groebner_basis())
#-------------------------------------------------------------------------------
def smooth_crit(f,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 : f(x)=0\} for the direction `\alpha`.
    
    INPUT:

    - ``f`` - An element of a `d`-variate polynomial ring over a number field.
    - ``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: f = (1-x-y-x*y)^2
        sage: var('a1,a2')
        (a1, a2)
        sage: I = smooth_crit(f,[a1,a2]); I
        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
        sage: solve(I.gens(),[var(g) for g in R.gens()],solution_dict=True)
        [{y: -(a1 + sqrt(a1^2 + a2^2))/a2, x: -(a2 + sqrt(a1^2 + a2^2))/a1}, {y: -(a1 - sqrt(a1^2 + a2^2))/a2, x: -(a2 - sqrt(a1^2 + a2^2))/a1}]
        
    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 [RaWi2008]_.
            
    REFERENCES:
    
    .. [RaWi2008]  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)
    """
    R= f.parent()
    B= R.base_ring()
    d= R.ngens()
    vars= R.gens()
    f= R.ideal(f).radical().gens()[0]   # Compute the reduction of f.
    
    # 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(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`.
    
    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: var('x,y,z')
        (x, y, z)
        sage: a= {x:1}
        sage: b= {y:2}
        sage: c= {z:3}
        sage: subs_all(x+y+z,a)
        y + z + 1
        sage: subs_all(x+y+z,[c,a])
        y + 4
        sage: subs_all([x+y+z,y^2],b)
        [x + z + 2, 4]
        sage: 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: 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 ============================================================================
# U ============================================================================
def unformat(As_over_Bs):
    r"""
    Converts the output format of the function cohom_equiv() into the fractions
    that it represents.
    
    INPUT:

    - ``As_over_Bs`` - A list of the form `[chunk_1,\ldots,chunk_r]` where each
      `chunk_j` has the form `[A,[B_1,e_1],\ldots,[B_l,e_l]]`.  Here
      `[A,[B_1,e_1],\ldots,[B_l,e_l]]` represents the fraction
      `A/(B_1^e_1 \cdots B_l^e_l)`.

    OUTPUT:

    The sum of the fractions coded by `chunk_1,\dots,chunk_r`.
    
    EXAMPLES:

        sage: R= QQ['w,z']
        sage: w,z = R.gens()
        sage: unformat([[1, [w, 1], [w*z - 1, 1]],[x, [w, 3]]])
        1/(w^2*z - w) + x/w^3
    
    AUTHORS:

    - Alex Raichev (2008-10-01)
    """
    sum= 0
    for x in As_over_Bs:
        sum= sum+ x[0]/prod([B^e for B,e in x[1:]])
    return sum
# V ============================================================================
# W ============================================================================
# X ============================================================================
# Y ============================================================================
# Z ============================================================================