Ticket #10519: amgf.py

File amgf.py, 111.3 KB (added by araichev, 6 years ago)

Pythonified amgf.sage

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