Ticket #10519: amgf.sage

File amgf.sage, 105.5 KB (added by araichev, 6 years ago)

Fixed a bug in cohomologous_integrand()

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