Ticket #10519: trac_10519.patch

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