Ticket #10519: amgf_alpha_0_5.sage

File amgf_alpha_0_5.sage, 104.5 KB (added by araichev, 6 years ago)

New object oriented version

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