Ticket #10519: mgf.sage

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