# Ticket #10519: trac_10519-fixed.patch

File trac_10519-fixed.patch, 114.0 KB (added by davidloeffler, 7 years ago)

Patch against 5.0

• ## sage/combinat/all.py

# HG changeset patch
# User Alex Raichev <tortoise.said@gmail.com>
# Date 1331399002 0
# Node ID b7819b98b3c6ed2078a83c4e08e0c59fe4b91357
# Parent  68e89260148df131fec1708338ceba3ea964b2bb
#10519: initial version

diff --git a/sage/combinat/all.py b/sage/combinat/all.py
 a from cyclic_sieving_phenomenon import CyclicSievingPolynomial, CyclicSievingCheck from sidon_sets import sidon_sets from amgf import QuasiRationalExpression
• ## new file sage/combinat/amgf.py

diff --git a/sage/combinat/amgf.py b/sage/combinat/amgf.py
new file mode 100644
 - r""" Asymptotics of Maclaurin coefficients of generating functions This code relates to analytic combinatorics. More specifically, it is a collection of functions designed to compute asymptotics of Maclaurin coefficients of certain classes of multivariate generating functions. The main function asymptotics() returns the first N terms of the asymptotic expansion of the Maclaurin coefficients F_{n\alpha} of the multivariate meromorphic function F=G/H as n\to\infty. It assumes that F is holomorphic in a neighborhood of the origin, that H is a polynomial, and that asymptotics in the direction of \alpha (a tuple of positive integers) are controlled by convenient smooth or multiple points. For an introduction to this subject, see [PeWi2008]_. The main algorithms and formulas implemented here come from [RaWi2008a]_ and [RaWi2011]_. REFERENCES: .. [AiYu1983]  I. A. A\u\izenberg and A. P. Yuzhakov, "Integral representations and residues in multidimensional complex analysis", Translations of Mathematical Monographs, 58. American Mathematical Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X. .. [DeLo2006]  Wolfram Decker and Christoph Lossen, "Computing in algebraic geometry", Chapter 7.1, Springer-Verlag, 2006. .. [DiEm2005]  Alicia Dickenstein and Ioannis Z. Emiris (editors), "Solving polynomial equations", Chapter 9.0, Springer-Verlag, 2005. .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of several variables into partial fractions", Soviet Math. (Iz. VUZ) 22 (1978), no. 10, 35--38. .. [PeWi2008]  Robin Pemantle and Mark C. Wilson, "Twenty combinatorial examples of asymptotics derived from multivariate generating functions", SIAM Rev. (2008) vol. 50 (2) pp. 199-272. .. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", Electronic Journal of Combinatorics, Vol. 15 (2008), R89. .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", To appear. AUTHORS: - Alex Raichev (2008-10-01) : Initial version - Alex Raichev (2010-09-28) : Corrected many functions - Alex Raichev (2010-12-15) : Updated documentation - Alex Raichev (2011-03-09) : Fixed a division by zero bug in relative_error() - Alex Raichev (2011-04-26) : Rewrote in object-oreinted style - Alex Raichev (2011-05-06) : Fixed bug in cohomologous_integrand() and fixed _crit_cone_combo() to work in SR EXAMPLES: These come from [RaWi2008a]_ and [RaWi2011]_. A smooth point example:: sage: R.= PolynomialRing(QQ) sage: G= 1 sage: H= 1-x-y-x*y sage: F= QuasiRationalExpression(G,H) sage: alpha= [3,2] sage: F.maclaurin_coefficients(alpha,8) {(21, 14): 1279919346549, (6, 4): 1289, (3, 2): 25, (18, 12): 19403906105, (24, 16): 85275509086721, (15, 10): 298199265, (9, 6): 75517, (12, 8): 4673345} sage: I= F.smooth_critical(alpha); I Ideal (y^2 + 3*y - 1, x - 2/3*y - 1/3) of Multivariate Polynomial Ring in x, y over Rational Field sage: s= solve(I.gens(),F.variables(),solution_dict=true); s [{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}] sage: c= s[1] sage: asys= [F.asymptotics(c,alpha,n,numerical=3) for n in [1..2]]; asys Initializing auxiliary functions... Calculating derivatives of auxiallary functions... Calculating derivatives of more auxiliary functions... Calculating actions of the second order differential operator... Initializing auxiliary functions... Calculating derivatives of auxiallary functions... Calculating derivatives of more auxiliary functions... Calculating actions of the second order differential operator... [(0.369*71.2^n/sqrt(n), 71.2, 0.369/sqrt(n)), ((0.369/sqrt(n) - 0.0186/n^(3/2))*71.2^n, 71.2, 0.369/sqrt(n) - 0.0186/n^(3/2))] sage: F.relative_error([f[0] for f in asys],alpha,[1,2,4,8,16],asys[0][1]) Calculating errors table in the form exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... [[(3, 2), 0.351196289062500, [0.369140625000000, 0.351000000000000], [-0.0510940551757812, 0.00174000000000000]], [(6, 4), 0.254364013671875, [0.261021839148940, 0.254558441227157], [-0.0262197902254473, 0.000151011402221735]], [(12, 8), 0.181976318359375, [0.184570312500000, 0.182000000000000], [-0.0142545700073242, -0.00151000000000000]], [(24, 16), 0.129287719726562, [0.130510919574470, 0.129683383669613], [-0.00947434154988702, -0.00267741572252445]], [(48, 32), 0.0914535522460938, [0.0922851562500000, 0.0920000000000000], [-0.00909328460693359, -0.00592000000000000]]] Another smooth point example. Here it turns out that the terms 2--4 of the asymptotic expansion are 0. :: sage: R. = PolynomialRing(QQ) sage: q= 1/2 sage: qq= q.denominator() sage: G= 1-q*x sage: H= 1-q*x +q*x*y -x^2*y sage: F= QuasiRationalExpression(G,H) sage: alpha= list(qq*vector([2,1-q])); alpha [4, 1] sage: I= F.smooth_critical(alpha); I Ideal (y^2 - 2*y + 1, x + 1/4*y - 5/4) of Multivariate Polynomial Ring in x, y over Rational Field sage: s= solve(I.gens(),F.variables(),solution_dict=true); s [{y: 1, x: 1}] sage: c= s[0] sage: asy= F.asymptotics(c,alpha,5); asy Initializing auxiliary functions... Calculating derivatives of auxiallary functions... Calculating derivatives of more auxiliary functions... Calculating actions of the second order differential operator... (1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*n^(1/3)) - 1/96*2^(1/3)*sqrt(3)*gamma(2/3)/(pi*n^(5/3)), 1, 1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*n^(1/3)) - 1/96*2^(1/3)*sqrt(3)*gamma(2/3)/(pi*n^(5/3))) sage: F.relative_error(asy[0],alpha,[1,2],asy[1]) Calculating errors table in the form exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... [[(4, 1), 0.187500000000000, [0.185581424449526], [0.0102324029358597]], [(8, 2), 0.152343750000000, [0.151986595969759], [0.00234439568568301]]] A multiple point example. :: sage: R.= PolynomialRing(QQ) sage: G= 1 sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) sage: F= QuasiRationalExpression(G,H) sage: I= F.singular_points(); I 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 sage: c= {x: 1/2, y: 1, z: 4/3} sage: F.is_cmp(c) [({y: 1, z: 4/3, x: 1/2}, True, 'all good')] sage: cc= F.critical_cone(c); cc ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N) sage: alpha= [8,3,3] sage: alpha in cc[1] True sage: asy= F.asymptotics(c,alpha,2); asy Initializing auxiliary functions... Calculating derivatives of auxiliary functions... Calculating derivatives of more auxiliary functions... Calculating second-order differential operator actions... (1/172872*(24696*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n)) - 1231*sqrt(3)*sqrt(7)/(sqrt(pi)*n^(3/2)))*108^n, 108, 1/7*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n)) - 1231/172872*sqrt(3)*sqrt(7)/(sqrt(pi)*n^(3/2))) Another multiple point example. :: sage: R.= PolynomialRing(QQ) sage: G= 16 sage: H= (4-2*x-y-z)^2*(4-x-2*y-z) sage: F= QuasiRationalExpression(G,H) sage: var('a1,a2,a3') (a1, a2, a3) sage: F.cohomologous_integrand([a1,a2,a3]) [[-16*(a2*z - 2*a3*y)*n/(y*z) + 16*(2*y - z)/(y*z), [x + 2*y + z - 4, 1], [2*x + y + z - 4, 1]]] sage: I= F.singular_points(); I Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field sage: c= {x:1, y:1, z:1} sage: F.is_cmp(c) [({y: 1, z: 1, x: 1}, True, 'all good')] sage: alpha= [3,3,2] sage: cc= F.critical_cone(c); cc ([(1, 2, 1), (2, 1, 1)], 2-d cone in 3-d lattice N) sage: alpha in cc[1] True sage: asy= F.asymptotics(c,alpha,2); asy Initializing auxiliary functions... Calculating derivatives of auxiliary functions... Calculating derivatives of more auxiliary functions... Calculating second-order differential operator actions... (4/3*sqrt(3)*sqrt(n)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(n)), 1, 4/3*sqrt(3)*sqrt(n)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(n))) """ #***************************************************************************** #       Copyright (C) 2008 Alexander Raichev # #  Distributed under the terms of the GNU General Public License (GPL) #                  http://www.gnu.org/licenses/ #***************************************************************************** from copy import copy, deepcopy # Sage libraries from sage.calculus.functional                           import diff from sage.calculus.functions                            import jacobian from sage.calculus.var                                  import function, var from sage.combinat.cartesian_product                    import CartesianProduct from sage.combinat.combinat                             import stirling_number1 from sage.combinat.permutation                          import Permutation from sage.combinat.tuple                                import UnorderedTuples from sage.functions.log                                 import exp, log from sage.functions.other                               import factorial, gamma, sqrt from sage.geometry.cone                                 import Cone from sage.matrix.constructor                            import matrix from sage.misc.misc                                     import add, ellipsis_range, mul from sage.misc.misc_c                                   import prod from sage.misc.mrange                                   import cartesian_product_iterator, mrange from sage.modules.free_module_element                   import vector from sage.rings.arith                                   import binomial from sage.rings.all                                     import CC from sage.rings.fraction_field                          import FractionField from sage.rings.integer                                 import Integer from sage.rings.integer_ring                            import ZZ from sage.rings.polynomial.polynomial_ring_constructor  import PolynomialRing from sage.rings.rational_field                          import QQ from sage.sets.set                                      import Set from sage.structure.sage_object                         import SageObject from sage.symbolic.constants                            import pi from sage.symbolic.relation                             import solve from sage.symbolic.ring                                 import SR # Functions to incorporate later into existing Sage classes. def algebraic_dependence(fs): r""" This function returns an irreducible annihilating polynomial for the polynomials in fs, if those polynomials are algebraically dependent. Otherwise it returns 0. INPUT: - fs - A list of polynomials f_1,\ldots,f_r' from a common polynomial ring. OUTPUT: If the polynomials in fs are algebraically dependent, then the output is an irreducible polynomial g in K[T_1,\ldots,T_r] such that g(f_1,\ldots,f_r) = 0. Here K is the coefficient ring of self and T_1,\ldots,T_r are indeterminates different from those of self. If the polynomials in fs are algebraically independent, then the output is the zero polynomial. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: fs= [x-3, x^2+1] sage: g= algebraic_dependence(fs); g 10 + 6*T0 - T1 + T0^2 sage: g(fs) 0 :: sage: R.= PolynomialRing(QQ) sage: fs= [w,  (w^2 + z^2 - 1)^2, w*z - 2] sage: g= algebraic_dependence(fs); g 16 + 32*T2 - 8*T0^2 + 24*T2^2 - 8*T0^2*T2 + 8*T2^3 + 9*T0^4 - 2*T0^2*T2^2 + T2^4 - T0^4*T1 + 8*T0^4*T2 - 2*T0^6 + 2*T0^4*T2^2 + T0^8 sage: g(fs) 0 sage: g.parent() Multivariate Polynomial Ring in T0, T1, T2 over Rational Field :: sage: R.= PolynomialRing(GF(7)) sage: fs= [x*y,y*z] sage: g= algebraic_dependence(fs); g 0 AUTHORS: - Alex Raichev (2011-01-06) """ r= len(fs) R= fs[Integer(0) ].parent() B= R.base_ring() Xs= list(R.gens()) d= len(Xs) # Expand R by r new variables. T= 'T' while T in [str(x) for x in Xs]: T= T+'T' Ts= [T + str(j) for j in range(r)] RR= PolynomialRing(B,d+r,tuple(Xs+Ts)) Vs= list(RR.gens()) Xs= Vs[Integer(0) :d] Ts= Vs[d:] # Find an irreducible annihilating polynomial g for the fs if there is one. J= RR.ideal([ Ts[j] -RR(fs[j]) for j in range(r)]) JJ= J.elimination_ideal(Xs) g= JJ.gens()[Integer(0) ] # Shrink the ambient ring of g to include only Ts. # I chose the negdeglex order simply because i find it useful in my work. RRR= PolynomialRing(B,r,tuple(Ts),order='negdeglex') return RRR(g) #------------------------------------------------------------------------------- def partial_fraction_decomposition(f): r""" Return a partial fraction decomposition of the rational expression f. Works for univariate and mulitivariate f. INPUT: - f - An element of the field of fractions F of a polynomial ring R whose coefficients ring is a field. In the univariate case, the coefficient ring doesn't have to be a field. OUTPUT: A tuple whole,parts where whole \in R and parts is the list of terms (in F) of a partial fraction decomposition of f - whole. See the notes below for more details. EXAMPLES:: sage: from sage.combinat.amgf import * sage: S. = QQ[] sage: q = 1/(t+1) + 2/(t+2) + 3/(t-3) sage: whole, parts = partial_fraction_decomposition(q); parts [3/(t - 3), 1/(t + 1), 2/(t + 2)] sage: whole +sum(parts) == q True Notice that there is a whole part below despite the appearance of q :: sage: R.= PolynomialRing(QQ) sage: q= 1/( x *y *(x*y-1) ) sage: whole,parts= partial_fraction_decomposition(q) sage: whole, parts (-1, [x*y/(x*y - 1), (-1)/(x*y)]) sage: q == whole +sum(parts) True :: sage: R.= PolynomialRing(QQ) sage: q= x +1/x +1/(x*y-2) +1/(x^2+y^2-1) sage: whole,parts= partial_fraction_decomposition(q) sage: whole, parts (x, [(1/2*x^3*y^2 + 1/2*x*y^4 + 1/2*x^3*y + 1/2*x^2*y^2 + 1/2*x*y^3 - x^2*y - 1/2*x*y^2 - y^3 - 3/2*x*y + y)/(x^3*y + x*y^3 - 2*x^2 - x*y - 2*y^2 + 2), (-1/2*x^3*y - 1/2*x*y^3 - 1/2*x^3 - 1/2*x^2*y - 1/2*x*y^2 + x^2 + 1/2*x*y + y^2 + 3/2*x - 1)/(x^3 + x*y^2 - x)]) sage: whole +sum(parts)==q True :: sage: R.= PolynomialRing(QQ) sage: q= 1/x +1/(x*y-z-2)^2 +1/(x^2+y^2 +z^2-1)^3 +1/(x*y-3) sage: whole,parts= partial_fraction_decomposition(q) sage: whole +sum(parts)==q True NOTES: In the case of univariate f this function calls the old univariate partial fraction decomposition function. In the multivariate case, it uses a different notion of and algorithm for partial fraction decompositions. Let f= P/Q where P,Q \in R, let Q_1^{e_1} \cdots Q_m^{e_m} be the unique factorization of Q in R into irreducible factors, and let d be the number of indeterminates of R. Then f can be written as a sum \sum P_A/\prod_{j\in A} Q_j^{b_j}, where the b_j \le e_j are positive integers, the P_A are in R, and the sum is taken over all size \le m subsets A \subseteq \{ 1, \ldots, d \} such that S:= \{ Q_j : j\in A \} is algebraically independent and the ideal generated by S is not all of R (that is, by the Nullstellensatz, the polynomials of S have no common roots in the algebraic closure of the coefficient field of R). Following [Lein1978]_, i call any such decomposition of f a \emph{partial fraction decomposition}. The algorithm used below comes from Theorem 1, Lemma 2, and Lemma 3 of [Lein1978]_. By the way, that paper has a typo in equation (c) on the second page and equation (b) on the third page. The right sides of (c) and (b) should be multiplied by P. REFERENCES: .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of several variables into partial fractions", Soviet Math. (Iz. VUZ) 22 (1978), no. 10, 35--38. AUTHORS: - Alex Raichev (2011-01-10) """ R= f.denominator().parent() d= len(R.gens()) if d==Integer(1) : return f.partial_fraction_decomposition() Q= f.denominator() whole,P= f.numerator().quo_rem(Q) parts= [format_quotient(Integer(1) /Q)] # Decompose via nullstellensatz trick # (which is faster than the algebraic dependence trick) parts= decompose_via_nullstellensatz(parts) # Decompose via algebraic dependence trick parts= decompose_via_algebraic_dependence(parts) # Rewrite parts back in terms of rational expressions new_parts=[] for p in parts: f= unformat_quotient(p) if f.denominator() == Integer(1) : whole= whole +f else: new_parts.append(f) # Put P back in new_parts= [P*f for f in new_parts] return whole, new_parts #------------------------------------------------------------------------------- def combine_like_terms(parts,rationomial=True): r""" Combines like terms of the fractions represented by parts. For use by partial_fraction_decomposition() above. EXAMPLES: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: parts =[[1, [x*y, 1]], [x, [x*y, 1]]] sage: combine_like_terms(parts) [[x + 1, [y, 1], [x, 1]]] :: sage: R.= PolynomialRing(QQ) sage: parts =[[1, [x, 1]], [x-1, [x, 1]]] sage: combine_like_terms(parts) [[1]] AUTHORS: - Alex Raichev (2011-01-10) """ if parts == []: return parts # Sort parts by denominators new_parts= sorted(parts, key= lambda p:p[Integer(1) :]) # Combine items of parts with same denominators. newnew_parts=[] left,right=Integer(0) ,Integer(1) glom= new_parts[left][Integer(0) ] while right <= len(new_parts)-Integer(1) : if new_parts[left][Integer(1) :] == new_parts[right][Integer(1) :]: glom= glom +new_parts[right][Integer(0) ] else: newnew_parts.append([glom]+new_parts[left][Integer(1) :]) left= right glom= new_parts[left][Integer(0) ] right= right +Integer(1) if glom != Integer(0) : newnew_parts.append([glom]+new_parts[left][Integer(1) :]) if rationomial: # Reduce fractions in newnew_parts. # Todo: speed up below by working directly with newnew_parts and # thereby make fewer calls to format_quotient() which in turn # calls factor(). newnew_parts= [format_quotient(unformat_quotient(part)) for part in newnew_parts] return newnew_parts #------------------------------------------------------------------------------- def decompose_via_algebraic_dependence(parts): r""" Returns a decomposition of parts. Used by partial_fraction_decomposition() above. Implements Lemma 2 of [Lein1978]_. Recursive. REFERENCES: .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of several variables into partial fractions", Soviet Math. (Iz. VUZ) 22 (1978), no. 10, 35--38. AUTHORS: - Alex Raichev (2011-01-10) """ decomposing_done= True new_parts= [] for p in parts: p_parts= [p] P= p[Integer(0) ] Qs= p[Integer(1) :] m= len(Qs) G= algebraic_dependence([q for q,e in Qs]) if G: # Then the denominator factors are algebraically dependent # and so we can decompose p. decomposing_done= False P= p[Integer(0) ] Qs= p[Integer(1) :] # Todo: speed up step below by using # G to calculate F.  See [Lein1978]_ Lemma 1. F= algebraic_dependence([q**e for q,e in Qs]) new_vars= F.parent().gens() # Note that new_vars[j] corresponds to Qs[j] so that # F([q^e for q,e in Qs]) = 0. # Assuming below that F.parent() has negdeglex term order # so that F.lt() is indeed the monomial we want. FF= (F.lt() -F)/(F.lc()) numers= map(mul,zip(FF.coefficients(),FF.monomials())) e= list(F.lt().exponents())[Integer(0) :m] denom= [[new_vars[j], e[Integer(0) ][j]+Integer(1) ] for j in range(m)] # Before making things messy by substituting in Qs, # reduce terms and combine like terms. p_parts_temp= [format_quotient(unformat_quotient([a]+denom)) for a in numers] p_parts_temp= combine_like_terms(p_parts_temp) # Substitute Qs into new_p. Qpowsub= dict([(new_vars[j],Qs[j][Integer(0) ]**Qs[j][Integer(1) ]) for j in range(m)]) p_parts=[] for x in p_parts_temp: y= P*F.parent()(x[Integer(0) ]).subs(Qpowsub) yy=[] for xx in x[Integer(1) :]: if xx[Integer(0) ] in new_vars: j= new_vars.index(xx[Integer(0) ]) yy.append([Qs[j][Integer(0) ],Qs[j][Integer(1) ]*xx[Integer(1) ]]) else: # Occasionally xx[0] is an integer. yy.append(xx) p_parts.append([y]+yy) # Done one step of decomposing p.  Add it to new_parts. new_parts.extend(p_parts) if decomposing_done: return new_parts else: return decompose_via_algebraic_dependence(new_parts) #------------------------------------------------------------------------------- def decompose_via_cohomology(parts): r""" Given each nice (described below) differential form (P/Q) dx_1 \wedge\cdots\wedge dx_d in parts, this function returns a differential form equivalent in De Rham cohomology that has no repeated factors in the denominator. INPUT: - parts - A list of the form [chunk_1,\ldots,chunk_r], where each chunk_j has the form [P,[Q_1,e_1],\ldots,[Q_m,e_m]], Q_1,\ldots,Q_m are irreducible elements of a common polynomial ring R such that their corresponding algebraic varieties \{x\in F^d : B_j(x)=0\} intersect transversely (where F is the algebraic closure of the field of coefficients of R), e_1,\ldots,e_m are positive integers, m \le d, and P is a symbolic expression in some of the indeterminates of R. Here [P,[Q_1,e_1],\ldots,[Q_m,e_m]] represents the fraction P/(Q_1^e_1 \cdots Q_m^e_m). OUTPUT: A list of the form [chunky_1,\ldots,chunky_s], where each chunky_j has the form [P,[Q_1,1],\ldots,[Q_m,1]]. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: decompose_via_cohomology([[ 1, [x,3] ]]) [] :: sage: R.= PolynomialRing(QQ) sage: decompose_via_cohomology([[ 1, [x,3], [x*y-1,2] ]]) [[-3*y^2, [x, 1], [x*y - 1, 1]], [-5*y^3, [x*y - 1, 1]]] NOTES: This is a recursive function thats stops calling itself when all the e_j equal 1 or parts == []. The algorithm used here is that of Theorem 17.4 of [AiYu1983]_. The algebraic varieties \{x\in F^d : Q_j(x)=0\} (where F is the algebraic closure of the field of coefficients of R) corresponding to the Q_j __intersect transversely__ iff for each point c of their intersection and for all k \le m, the Jacobian matrix of any k polynomials from \{Q_1,\ldots,Q_m\} has rank equal to \min\{k,d\} when evaluated at c. REFERENCES: .. [AiYu1983]  I. A. A\u\izenberg and A. P. Yuzhakov, "Integral representations and residues in multidimensional complex analysis", Translations of Mathematical Monographs, 58. American Mathematical Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X. AUTHORS: - Alex Raichev (2008-10-01, 2011-01-15) """ if parts == []: return parts decomposing_done= True new_parts= [] R= parts[Integer(0) ][Integer(1) ][Integer(0) ].parent() V= list(R.gens()) for p in parts: p_parts= [p] P= p[Integer(0) ] Qs= p[Integer(1) :] m= len(Qs) if sum([e for q,e in Qs]) > m: # Then we can decompose p p_parts= [] decomposing_done=False dets= [] vars_subsets= Set(V).subsets(m) for v in vars_subsets: # Sort variables so that first polynomial ring indeterminate # declared is first in vars list. v= sorted(v,reverse=True) jac= jacobian([q for q,e in Qs],v) dets.append(jac.determinant()) # Get a Nullstellensatz certificate. L= R(Integer(1) ).lift(R.ideal([q for q,e in Qs] +dets)) for j in range(m): if L[j] != Integer(0) : # Make a copy of (and not a reference to) the nested list Qs. new_Qs = deepcopy(Qs) if new_Qs[j][Integer(1) ] > Integer(1) : new_Qs[j][Integer(1) ]= new_Qs[j][Integer(1) ] -Integer(1) else: del new_Qs[j] p_parts.append( [P*L[j]] +new_Qs ) for k in range(vars_subsets.cardinality()): if L[m+k] != Integer(0) : new_Qs = deepcopy(Qs) for j in range(m): if new_Qs[j][Integer(1) ] > Integer(1) : new_Qs[j][Integer(1) ]= new_Qs[j][Integer(1) ] -Integer(1) v=sorted(vars_subsets[k],reverse=True) jac= jacobian([SR(P*L[m+k])] +[ SR(Qs[jj][Integer(0) ]) for \ jj in range(m) if jj !=j], [SR(vv) for vv in v]) det= jac.determinant() if det != Integer(0) : p_parts.append([permutation_sign(v,V) \ *(-Integer(1) )**j/new_Qs[j][Integer(1) ] *det] +new_Qs) break # Done one step of decomposing p.  Add it to new_parts. new_parts.extend(p_parts) new_parts= combine_like_terms(new_parts,rationomial=False) if decomposing_done: return new_parts else: return decompose_via_cohomology(new_parts) #------------------------------------------------------------------------------- def decompose_via_nullstellensatz(parts): r""" Returns a decomposition of parts. Used by partial_fraction_decomposition() above. Implements Lemma 3 of [Lein1978]_. Recursive. REFERENCES: .. [Lein1978]  E. K. Leinartas, "On expansion of rational functions of several variables into partial fractions", Soviet Math. (Iz. VUZ) 22 (1978), no. 10, 35--38. AUTHORS: - Alex Raichev (2011-01-10) """ decomposing_done= True new_parts= [] R= parts[Integer(0) ][Integer(0) ].parent() for p in parts: p_parts= [p] P= p[Integer(0) ] Qs= p[Integer(1) :] m= len(Qs) if R(Integer(1) ) in R.ideal([q for q,e in Qs]): # Then we can decompose p. decomposing_done= False L= R(Integer(1) ).lift(R.ideal([q**e for q,e in Qs])) p_parts= [ [P*L[i]] + \ [[Qs[j][Integer(0) ],Qs[j][Integer(1) ]] for j in range(m) if j != i] \ for i in range(m) if L[i]!=Integer(0) ] # The procedure above yields no like terms to combine. # Done one step of decomposing p.  Add it to new_parts. new_parts.extend(p_parts) if decomposing_done: return new_parts else: return decompose_via_nullstellensatz(new_parts) #------------------------------------------------------------------------------- def format_quotient(f): r""" Formats f for use by partial_fraction_decomposition() above. AUTHORS: - Alex Raichev (2011-01-10) """ P= f.numerator() Q= f.denominator() Qs= Q.factor() if Qs.unit() != Integer(1) : P= P/Qs.unit() Qs= sorted([[q,e] for q,e in Qs])  # sorting for future bookkeeping return [P]+Qs #------------------------------------------------------------------------------- def permutation_sign(v,vars): r""" This function returns the sign of the permutation on 1,\ldots,len(vars) that is induced by the sublist v of vars. For internal use by _cohom_equiv_main(). INPUT: - v - A sublist of vars. - vars - A list of predefined variables or numbers. OUTPUT: The sign of the permutation obtained by taking indices (and adding 1) within vars of the list v,w, where w is the list vars with the elements of v removed. EXAMPLES:: sage: from sage.combinat.amgf import * sage: vars= ['a','b','c','d','e'] sage: v= ['b','d'] sage: permutation_sign(v,vars) -1 sage: v= ['d','b'] sage: permutation_sign(v,vars) 1 AUTHORS: - Alex Raichev (2008-10-01) """ # Convert variable lists to lists of numbers in {1,...,len(vars)} A= [x+Integer(1)  for x in range(len(vars))] B= [vars.index(x)+Integer(1)  for x in v] C= list(Set(A).difference(Set(B))) C.sort() P= Permutation(B+C) return P.signature() #------------------------------------------------------------------------------- def unformat_quotient(part): r""" Unformats f for use by partial_fraction_decomposition() above. Inverse of format_quotient() above. AUTHORS: - Alex Raichev (2011-01-10) """ P= part[Integer(0) ] Qs= part[Integer(1) :] Q= prod([q**e for q,e in Qs]) return P/Q #=============================================================================== # Class for calculation of asymptotics of multivariate generating functions. class QuasiRationalExpression(SageObject): "Store an expression G/H, where H comes from a polynomial ring R and \ G comes from R or the Symbolic Ring." def __init__(self, G, H): # Store important information about object as attributes of self. # G, H, H's ring, ring dimension, H's factorization. self._G = G self._H = H R= H.parent() self._R = R self._d = len(R.gens()) self._Hfac= list(H.factor()) # Variables of self as elements of the SR. # Remember that G might be in SR and not in R. try: # This fails if G is a constant, for example. Gv= Set([R(x) for x in G.variables()]) except: Gv= Set([]) try: # This fails if H is a constant, for example. Hv= Set(H.variables()) except: Hv= Set([]) # Preserve the ring ordering of the variables which some methods below # depends upon. V= sorted(list(Gv.union(Hv)),reverse=True) self._variables = tuple([SR(x) for x in V]) #------------------------------------------------------------------------------- # Keeping methods in alphabetical order (ignoring initial single underscores) def asymptotics(self,c,alpha,N,numerical=Integer(0) ,asy_var=None): r""" This function returns the first N terms of the asymptotic expansion of the Maclaurin coefficients F_{n\alpha} of the multivariate meromorphic function F=G/H as n\to\infty, where F = self. It assumes that F is holomorphic in a neighborhood of the origin, that H is a polynomial, and that c is a convenient strictly minimal smooth or multiple of F that is critical for \alpha. INPUT: - alpha - A d-tuple of positive integers or, if c is a smooth point, possibly of symbolic variables. - c - A dictionary with keys self._variables and values from a superfield of the field of self._R.base_ring(). - N - A positive integer. - numerical - A natural number (default: 0). If k=numerical > 0, then a numerical approximation of the coefficients of F_{n\alpha} with k digits of precision will be returned. Otherwise exact values will be returned. - asy_var - A  symbolic variable (default: None). The variable of the asymptotic expansion. If none is given, var('n') will be assigned. OUTPUT: The tuple (asy,exp_part,subexp_part), where asy is first N terms of the asymptotic expansion of the Maclaurin coefficients F_{n\alpha} of the function F=self as n\to\infty, exp_part^n is the exponential factor of asy, and subexp_part is the subexponential factor of asy. EXAMPLES:: A smooth point example :: sage: R.= PolynomialRing(QQ) sage: G= 1 sage: H= 1-x-y-x*y sage: F= QuasiRationalExpression(G,H) sage: alpha= [3,2] sage: c= {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3} sage: F.asymptotics(c,alpha,2,numerical=3) Initializing auxiliary functions... Calculating derivatives of auxiallary functions... Calculating derivatives of more auxiliary functions... Calculating actions of the second order differential operator... ((0.369/sqrt(n) - 0.0186/n^(3/2))*71.2^n, 71.2, 0.369/sqrt(n) - 0.0186/n^(3/2)) A multiple point example :: sage: R.= PolynomialRing(QQ) sage: G= 1 sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) sage: F= QuasiRationalExpression(G,H) sage: c= {z: 4/3, y: 1, x: 1/2} sage: alpha= [8,3,3] sage: F.asymptotics(c,alpha,1) Initializing auxiliary functions... Calculating derivatives of auxiliary functions... Calculating derivatives of more auxiliary functions... Calculating second-order differential operator actions... (1/7*sqrt(3)*sqrt(7)*108^n/(sqrt(pi)*sqrt(n)), 108, 1/7*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n))) NOTES: A zero c of H is __strictly minimal__ if there is no zero x of H such that x_j < c_j for all 0 \le j < d. For definitions of the terms "smooth critical point for \alpha" and "multiple critical point for \alpha", see the documentation for _asymptotics_main_smooth() and _asymptotics_main_multiple(), which are the functions that do most of the work. ALGORITHM: The algorithm used here comes from [RaWi2011]_. (1) Use Cauchy's multivariate integral formula to write F_{n\alpha} as an integral around a polycirle centered at the origin of the differential form \frac{G(x) dx[0] \wedge\cdots\wedge dx[d-1]}{H(x)x^\alpha}. (2) Decompose G/H into a sum of partial fractions P[0] +\cdots+ P[r] so that each term of the sum has at most d irreducible factors of H in the denominator. (3) For each differential form P[j] dx[0] \wedge\cdots\wedge dx[d-1], find an equivalent form \omega[j] in de Rham cohomology with no repeated irreducible factors of H in its denominator. (4) Compute an asymptotic expansion for each integral \omega[j]. (5) Add the expansions. REFERENCES: .. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", Electronic Journal of Combinatorics, Vol. 15 (2008), R89. .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", To appear. AUTHORS: - Alex Raichev (2008-10-01, 2010-09-28, 2011-04-27) """ # The variable for asymptotic expansions. if not asy_var: asy_var= var('n') # Create symbolic (non-ring) variables. R= self._R d= self._d X= list(self._variables) # Do steps (1)--(3) new_integrands= self.cohomologous_integrand(alpha,asy_var) # Coerce everything into the Symbolic Ring, as polynomial ring # calculations are no longer needed. # Calculate asymptotics. cc={} for k in c.keys(): cc[SR(k)] = SR(c[k]) c= cc for i in range(len(alpha)): alpha[i] = SR(alpha[i]) subexp_parts= [] for chunk in new_integrands: # Convert chunk into Symbolic Ring GG= SR(chunk[Integer(0) ]) HH= [SR(f) for (f,e) in chunk[Integer(1) :]] asy= self._asymptotics_main(GG,HH,X,c,asy_var,alpha,N,numerical) subexp_parts.append(asy[Integer(2) ]) exp_scale= asy[Integer(1) ]  # same for all chunk in new_integrands # Do step (5). subexp_part= add(subexp_parts) return exp_scale**asy_var *subexp_part, exp_scale, subexp_part #-------------------------------------------------------------------------------- def _asymptotics_main(self,G,H,X,c,n,alpha,N,numerical): r""" This function is for internal use by asymptotics(). It finds a variable in X to use to calculate asymptotics and decides whether to call _asymptotics_main_smooth() or _asymptotics_main_multiple(). Does not use self. INPUT: - G - A symbolic expression. - H - A list of symbolic expressions. - X - The list of variables occurring in G and H. Call its length d. - c - A dictionary with X as keys and numbers as values. - n - The variable of the asymptotic expansion. - alpha - A d-tuple of positive natural numbers or possibly of symbolic variables if c is a smooth point. - N - A positive integer. - numerical - Natural number. If k=numerical > 0, then a numerical approximation of the coefficients of F_{n\alpha} with k digits of precision will be returned. Otherwise exact values will be returned. OUTPUT: The same as the function asymptotics(). AUTHORS: - Alex Raichev (2008-10-01, 2010-09-28) """ d= len(X) r= len(H)  # We know 1 <= r <= d. # Find a good variable x in X to do asymptotics calculations with, that is, # a variable x in X such that for all h in H, diff(h,x).subs(c) != 0. # A good variable is guaranteed to exist since we are dealing with # convenient smooth or multiple points. # Search for good x in X from back to front (to be consistent with # [RaWi2008a]_ [RaWi2011]_ which uses X[d-1] as a good variable). # Put each not good x found at the beginning of X and reshuffle alpha # in the same way. x= X[d-Integer(1) ] beta= copy(alpha) while Integer(0)  in [diff(h,x).subs(c) for h in H]: X.pop() X.insert(Integer(0) ,x) x= X[d-Integer(1) ] a= beta.pop() beta.insert(Integer(0) ,a) if d==r: # This is the case of a 'simple' multiple point. A= G.subs(c) / jacobian(H,X).subs(c).determinant().abs() return A,Integer(1) ,A elif r==Integer(1) : # So 1 = r < d, and we have a smooth point. return self._asymptotics_main_smooth(G,H[Integer(0) ],X,c,n,beta,N,numerical) else: # So 1 < r < d, and we have a non-smooth multiple point. return self._asymptotics_main_multiple(G,H,X,c,n,beta,N,numerical) #-------------------------------------------------------------------------------- def _asymptotics_main_multiple(self,G,H,X,c,n,alpha,N,numerical): r""" This function is for internal use by _asymptotics_main(). It calculates asymptotics in case they depend upon multiple points. Does not use self. INPUT: - G - A symbolic expression. - H - A list of symbolic expressions. - X - The list of variables occurring in G and H. Call its length d. - c - A dictionary with X as keys and numbers as values. - n - The variable of the asymptotic expansion. - alpha - A d-tuple of positive natural numbers or possibly of symbolic variables if c is a smooth point. - N - A positive integer. - numerical - Natural number. If k=numerical > 0, then a numerical approximation of the coefficients of F_{n\alpha} with k digits of precision will be returned. Otherwise exact values will be returned. OUTPUT: The same as the function asymptotics(). NOTES: The formula used for computing the asymptotic expansion is Theorem 3.4 of [RaWi2011]_. Currently this function cannot handle c with symbolic variable keys, because _crit_cone_combo() crashes. REFERENCES: .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", To appear. AUTHORS: - Alex Raichev (2008-10-01, 2010-09-28) """ I= sqrt(-Integer(1) ) d= len(X)   # > r > 1 r= len(H)   # > 1 C= copy(c) S= [var(self._new_var_name('s',X) + str(j)) for j in range(r-Integer(1) )] T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-Integer(1) )] Sstar= {} temp= self._crit_cone_combo(H,X,c,alpha) for j in range(r-Integer(1) ): Sstar[S[j]]= temp[j] thetastar= dict([(t,Integer(0) ) for t in T]) thetastar.update(Sstar) # Setup. print "Initializing auxiliary functions..." Hmul= mul(H) h= [function('h'+str(j),*tuple(X[:d-Integer(1) ])) for j in range(r)]    # Implicit functions U = function('U',*tuple(X)) # All other functions are defined in terms of h, U, and explicit functions. Hcheck=  mul([ X[d-Integer(1) ] -Integer(1) /h[j] for j in range(r)] ) Gcheck= -G/U *mul( [-h[j]/X[d-Integer(1) ] for j in range(r)] ) A= [(-Integer(1) )**(r-Integer(1) ) *X[d-Integer(1) ]**(-r+j)*diff(Gcheck.subs({X[d-Integer(1) ]:Integer(1) /X[d-Integer(1) ]}),X[d-Integer(1) ],j) for j in range(r)] e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-Integer(1) )]) ht= [hh.subs(e) for hh in h] Ht= [H[j].subs(e).subs({X[d-Integer(1) ]:Integer(1) /ht[j]}) for j in range(r)] hsumt= add([S[j]*ht[j] for j in range(r-Integer(1) )]) +(Integer(1) -add(S))*ht[r-Integer(1) ] At= [AA.subs(e).subs({X[d-Integer(1) ]:hsumt}) for AA in A] Phit = -log(C[X[d-Integer(1) ]] *hsumt)+ I*add([alpha[i]/alpha[d-Integer(1) ] *T[i] for i in range(d-Integer(1) )]) # atC Stores h and U and all their derivatives evaluated at C. atC = copy(C) atC.update(dict( [(hh.subs(C),Integer(1) /C[X[d-Integer(1) ]]) for hh in h ])) # Compute the derivatives of h up to order 2*N and evaluate at C. hderivs1= {}   # First derivatives of h. for (i,j) in mrange([d-Integer(1) ,r]): s= solve( diff(H[j].subs({X[d-Integer(1) ]:Integer(1) /h[j]}),X[i]), diff(h[j],X[i]) )[Integer(0) ].rhs()\ .simplify() hderivs1.update({diff(h[j],X[i]):s}) atC.update({diff(h[j],X[i]).subs(C):s.subs(C).subs(atC)}) hderivs = self._diff_all(h,X[Integer(0) :d-Integer(1) ],Integer(2) *N,sub=hderivs1,rekey=h) for k in hderivs.keys(): atC.update({k.subs(C):hderivs[k].subs(atC)}) # Compute the derivatives of U up to order 2*N-2+ min{r,N}-1 and evaluate at C. # To do this, differentiate H = U*Hcheck over and over, evaluate at C, # and solve for the derivatives of U at C. # Need the derivatives of H with short keys to pass on to diff_prod later. m= min(r,N) end= [X[d-Integer(1) ] for j in range(r)] Hmulderivs= self._diff_all(Hmul,X,Integer(2) *N-Integer(2) +r,ending=end,sub_final=C) print "Calculating derivatives of auxiliary functions..." atC.update({U.subs(C):diff(Hmul,X[d-Integer(1) ],r).subs(C)/factorial(r)}) Uderivs={} p= Hmul.polynomial(CC).degree()-r if p == Integer(0) : # Then, using a proposition at the end of [RaWi2011], we can # conclude that all higher derivatives of U are zero. for l in (ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )): for s in UnorderedTuples(X,l): Uderivs[diff(U,s).subs(C)] = Integer(0) elif p > Integer(0)  and p < Integer(2) *N-Integer(2) +m-Integer(1) : all_zero= True Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,p)),end,Uderivs,atC) # Check for a nonzero U derivative. if Uderivs.values() != [Integer(0)  for i in range(len(Uderivs))]: all_zero= False if all_zero: # Then, using a proposition at the end of [RaWi2011], we can # conclude that all higher derivatives of U are zero. for l in (ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )): for s in UnorderedTuples(X,l): Uderivs.update({diff(U,s).subs(C):Integer(0) }) else: # Have to compute the rest of the derivatives. Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )),end,Uderivs,atC) else: Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )),end,Uderivs,atC) atC.update(Uderivs) Phit1= jacobian(Phit,T+S).subs(hderivs1) a= jacobian(Phit1,T+S).subs(hderivs1).subs(thetastar).subs(atC) a_inv= a.inverse() Phitu= Phit -(Integer(1) /Integer(2) ) *matrix([T+S]) *a *matrix([T+S]).transpose() Phitu= Phitu[Integer(0) ][Integer(0) ] # Compute all partial derivatives of At and Phitu up to orders 2*N-2 # and 2*N, respectively.  Take advantage of the fact that At and Phitu # are sufficiently differentiable functions so that mixed partials # are equal.  Thus only need to compute representative partials. # Choose nondecreasing sequences as representative differentiation- # order sequences. # To speed up later computations, create symbolic functions AA and BB # to stand in for the expressions At and Phitu respectively. print "Calculating derivatives of more auxiliary functions..." AA= [function('A'+str(j),*tuple(T+S)) for j in range(r)] At_derivs= self._diff_all(At,T+S,Integer(2) *N-Integer(2) ,sub=hderivs1,sub_final=[thetastar,atC],rekey=AA) BB= function('BB',*tuple(T+S)) Phitu_derivs= self._diff_all(Phitu,T+S,Integer(2) *N,sub=hderivs1,sub_final=[thetastar,atC],rekey=BB,zero_order=Integer(3) ) AABB_derivs= At_derivs AABB_derivs.update(Phitu_derivs) for j in range(r): AABB_derivs[AA[j]] = At[j].subs(thetastar).subs(atC) AABB_derivs[BB] = Phitu.subs(thetastar).subs(atC) print "Calculating second-order differential operator actions..." DD= self._diff_op(AA,BB,AABB_derivs,T+S,a_inv,r,N) L={} for (j,k) in  CartesianProduct((ellipsis_range(Integer(0) ,Ellipsis,min(r-Integer(1) ,N-Integer(1) ))), (ellipsis_range(max(Integer(0) ,N-Integer(1) -r),Ellipsis,N-Integer(1) ))): if j+k <= N-Integer(1) : L[(j,k)] = add([ \ DD[(j,k,l)] /( (-Integer(1) )**k *Integer(2) **(k+l) *factorial(l) *factorial(k+l) ) \ for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k))] ) # The next line's QQ coercion is a workaround for the Sage 4.6 bug reported # on http://trac.sagemath.org/sage_trac/ticket/8659. # Once the bug is fixed, the QQ can be removed. det= QQ(a.determinant())**(-Integer(1) /Integer(2) ) *(Integer(2) *pi)**((r-d)/Integer(2) ) chunk= det *add([ (alpha[d-Integer(1) ]*n)**((r-d)/Integer(2) -q) *add([ \ L[(j,k)] *binomial(r-Integer(1) ,j) *stirling_number1(r-j,r+k-q) *(-Integer(1) )**(q-j-k) \ for (j,k) in CartesianProduct((ellipsis_range(Integer(0) ,Ellipsis,min(r-Integer(1) ,q))), (ellipsis_range(max(Integer(0) ,q-r),Ellipsis,q))) if j+k <= q ]) \ for q in range(N)]) chunk= chunk.subs(C).simplify() coeffs= chunk.coefficients(n) coeffs.reverse() coeffs= coeffs[:N] if numerical: # If a numerical approximation is desired... subexp_part = add( [co[Integer(0) ].subs(c).n(digits=numerical)*n**co[Integer(1) ] \ for co in coeffs] ) exp_scale= (Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )) \ .n(digits=numerical) else: subexp_part = add( [co[Integer(0) ].subs(c)*n**co[Integer(1) ] for co in coeffs] ) exp_scale= Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] ) return exp_scale**n*subexp_part, exp_scale, subexp_part #-------------------------------------------------------------------------------- def _asymptotics_main_smooth(self,G,H,X,c,n,alpha,N,numerical): r""" This function is for internal use by _asymptotics_main(). It calculates asymptotics in case they depend upon smooth points. Does not use self. INPUT: - G - A symbolic expression. - H - A list of symbolic expressions. - X - The list of variables occurring in G and H. Call its length d. - c - A dictionary with X as keys and numbers as values. - n - The variable of the asymptotic expansion. - alpha - A d-tuple of positive natural numbers or possibly of symbolic variables if c is a smooth point. - N - A positive integer. - numerical - Natural number. If k=numerical > 0, then a numerical approximation of the coefficients of F_{n\alpha} with k digits of precision will be returned. Otherwise exact values will be returned. OUTPUT: The same as the function asymptotics(). NOTES: The formulas used for computing the asymptotic expansions are Theorems 3.2 and 3.3 [RaWi2008a]_ with p=1. Theorem 3.2 is a specialization of Theorem 3.4 of [RaWi2011]_ with r=1. REFERENCES: .. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", Electronic Journal of Combinatorics, Vol. 15 (2008), R89. .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", To appear. AUTHORS: - Alex Raichev (2008-10-01, 2010-09-28) """ I= sqrt(-Integer(1) ) d= len(X) # > 1 # If c is a tuple of rationals, then compute with it directly. # Otherwise, compute symbolically and plug in c at the end. if vector(c.values()) in QQ**d: C= c else: Cs= [var('cs' +str(j)) for j in range(d)] C= dict( [(X[j],Cs[j]) for j in range(d)] ) c= dict( [(Cs[j],c[X[j]]) for j in range(d)] ) # Setup. print "Initializing auxiliary functions..." h= function('h',*tuple(X[:d-Integer(1) ]))    # Implicit functions U = function('U',*tuple(X))         # # All other functions are defined in terms of h, U, and explicit functions. Gcheck = -G/U *(h/X[d-Integer(1) ]) A= Gcheck.subs({X[d-Integer(1) ]:Integer(1) /h})/h T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-Integer(1) )] e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-Integer(1) )]) ht= h.subs(e) Ht= H.subs(e).subs({X[d-Integer(1) ]:Integer(1) /ht}) At= A.subs(e) Phit = -log(C[X[d-Integer(1) ]] *ht)\ + I* add([alpha[i]/alpha[d-Integer(1) ] *T[i] for i in range(d-Integer(1) )]) Tstar= dict([(t,Integer(0) ) for t in T]) atC = copy(C) atC.update({h.subs(C):Integer(1) /C[X[d-Integer(1) ]]}) # Stores h and U and all their derivatives # evaluated at C. # Compute the derivatives of h up to order 2*N and evaluate at C and store # in atC.  Keep a copy of unevaluated h derivatives for use in the case # d=2 and v > 2 below. hderivs1= {}   # First derivatives of h. for i in range(d-Integer(1) ): s= solve( diff(H.subs({X[d-Integer(1) ]:Integer(1) /h}),X[i]), diff(h,X[i]) )[Integer(0) ].rhs()\ .simplify() hderivs1.update({diff(h,X[i]):s}) atC.update({diff(h,X[i]).subs(C):s.subs(C).subs(atC)}) hderivs = self._diff_all(h,X[Integer(0) :d-Integer(1) ],Integer(2) *N,sub=hderivs1,rekey=h) for k in hderivs.keys(): atC.update({k.subs(C):hderivs[k].subs(atC)}) # Compute the derivatives of U up to order 2*N and evaluate at C. # To do this, differentiate H = U*Hcheck over and over, evaluate at C, # and solve for the derivatives of U at C. # Need the derivatives of H with short keys to pass on to diff_prod later. Hderivs= self._diff_all(H,X,Integer(2) *N,ending=[X[d-Integer(1) ]],sub_final=C) print "Calculating derivatives of auxiallary functions..." # For convenience in checking if all the nontrivial derivatives of U at c # are zero a few line below, store the value of U(c) in atC instead of in # Uderivs. Uderivs={} atC.update({U.subs(C):diff(H,X[d-Integer(1) ]).subs(C)}) end= [X[d-Integer(1) ]] Hcheck= X[d-Integer(1) ] - Integer(1) /h p= H.polynomial(CC).degree()-Integer(1) if p == Integer(0) : # Then, using a proposition at the end of [RaWi2011], we can # conclude that all higher derivatives of U are zero. for l in (ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N)): for s in UnorderedTuples(X,l): Uderivs[diff(U,s).subs(C)] = Integer(0) elif p > Integer(0)  and p < Integer(2) *N: all_zero= True Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,p)),end,Uderivs,atC) # Check for a nonzero U derivative. if Uderivs.values() != [Integer(0)  for i in range(len(Uderivs))]: all_zero= False if all_zero: # Then, using a proposition at the end of [RaWi2011], we can # conclude that all higher derivatives of U are zero. for l in (ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N)): for s in UnorderedTuples(X,l): Uderivs.update({diff(U,s).subs(C):Integer(0) }) else: # Have to compute the rest of the derivatives. Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N)),end,Uderivs,atC) else: Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N)),end,Uderivs,atC) atC.update(Uderivs) # In general, this algorithm is not designed to handle the case of a # singular Phit''(Tstar).  However, when d=2 the algorithm can cope. if d==Integer(2) : # Compute v, the order of vanishing at Tstar of Phit. It is at least 2. v=Integer(2) Phitderiv= diff(Phit,T[Integer(0) ],Integer(2) ) splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify() while splat==Integer(0) : v= v+Integer(1) if v > Integer(2) *N:  # Then need to compute more derivatives of h for atC. hderivs.update({diff(h,X[Integer(0) ],v) \ :diff(hderivs[diff(h,X[Integer(0) ],v-Integer(1) )],X[Integer(0) ]).subs(hderivs1)}) atC.update({diff(h,X[Integer(0) ],v).subs(C) \ :hderivs[diff(h,X[Integer(0) ],v)].subs(atC)}) Phitderiv= diff(Phitderiv,T[Integer(0) ]) splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify() if d==Integer(2)  and v>Integer(2) : t= T[Integer(0) ]  # Simplify variable names. a= splat/factorial(v) Phitu= Phit -a*t**v # Compute all partial derivatives of At and Phitu up to orders 2*(N-1) # and 2*(N-1)+v, respectively, in case v is even. # Otherwise, compute up to orders N-1 and N-1+v, respectively. # To speed up later computations, create symbolic functions AA and BB # to stand in for the expressions At and Phitu, respectively. print "Calculating derivatives of more auxiliary functions..." AA= function('AA',t) BB= function('BB',t) if v.mod(Integer(2) )==Integer(0) : At_derivs= self._diff_all(At,T,Integer(2) *N-Integer(2) , \ sub=hderivs1,sub_final=[Tstar,atC],rekey=AA) Phitu_derivs= self._diff_all(Phitu,T,Integer(2) *N-Integer(2) +v, \ sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+Integer(1) ,rekey=BB) else: At_derivs= self._diff_all(At,T,N-Integer(1) ,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA) Phitu_derivs= self._diff_all(Phitu,T,N-Integer(1) +v,sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+Integer(1) ,rekey=BB) AABB_derivs= At_derivs AABB_derivs.update(Phitu_derivs) AABB_derivs[AA] = At.subs(Tstar).subs(atC) AABB_derivs[BB] = Phitu.subs(Tstar).subs(atC) print "Calculating actions of the second order differential operator..." DD= self._diff_op_simple(AA,BB,AABB_derivs,t,v,a,N) # Plug above into asymptotic formula. L = [] if v.mod(Integer(2) ) == Integer(0) : for k in range(N): L.append( add([ \ (-Integer(1) )**l *gamma((Integer(2) *k+v*l+Integer(1) )/v) \ / (factorial(l) *factorial(Integer(2) *k+v*l)) \ * DD[(k,l)] for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)) ]) ) chunk= a**(-Integer(1) /v) /(pi*v) *add([ alpha[d-Integer(1) ]**(-(Integer(2) *k+Integer(1) )/v) \ * L[k] *n**(-(Integer(2) *k+Integer(1) )/v) for k in range(N) ]) else: zeta= exp(I*pi/(Integer(2) *v)) for k in range(N): L.append( add([ \ (-Integer(1) )**l *gamma((k+v*l+Integer(1) )/v) \ / (factorial(l) *factorial(k+v*l)) \ * (zeta**(k+v*l+Integer(1) ) +(-Integer(1) )**(k+v*l)*zeta**(-(k+v*l+Integer(1) ))) \ * DD[(k,l)] for l in (ellipsis_range(Integer(0) ,Ellipsis,k)) ]) ) chunk= abs(a)**(-Integer(1) /v) /(Integer(2) *pi*v) *add([ alpha[d-Integer(1) ]**(-(k+Integer(1) )/v) \ * L[k] *n**(-(k+Integer(1) )/v) for k in range(N) ]) # Asymptotics for d>=2 case.  A singular Phit''(Tstar) will cause a crash # in this case. else: Phit1= jacobian(Phit,T).subs(hderivs1) a= jacobian(Phit1,T).subs(hderivs1).subs(Tstar).subs(atC) a_inv= a.inverse() Phitu= Phit -(Integer(1) /Integer(2) ) *matrix([T]) *a *matrix([T]).transpose() Phitu= Phitu[Integer(0) ][Integer(0) ] # Compute all partial derivatives of At and Phitu up to orders 2*N-2 # and 2*N, respectively.  Take advantage of the fact that At and Phitu # are sufficiently differentiable functions so that mixed partials # are equal.  Thus only need to compute representative partials. # Choose nondecreasing sequences as representative differentiation- # order sequences. # To speed up later computations, create symbolic functions AA and BB # to stand in for the expressions At and Phitu respectively. print "Calculating derivatives of more auxiliary functions..." AA= function('AA',*tuple(T)) At_derivs= self._diff_all(At,T,Integer(2) *N-Integer(2) ,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA) BB= function('BB',*tuple(T)) Phitu_derivs= self._diff_all(Phitu,T,Integer(2) *N,sub=hderivs1,sub_final=[Tstar,atC],rekey=BB,zero_order=Integer(3) ) AABB_derivs= At_derivs AABB_derivs.update(Phitu_derivs) AABB_derivs[AA] = At.subs(Tstar).subs(atC) AABB_derivs[BB] = Phitu.subs(Tstar).subs(atC) print "Calculating actions of the second order differential operator..." DD= self._diff_op(AA,BB,AABB_derivs,T,a_inv,Integer(1) ,N) # Plug above into asymptotic formula. L=[] for k in range(N): L.append( add([ \ DD[(Integer(0) ,k,l)] / ( (-Integer(1) )**k *Integer(2) **(l+k) *factorial(l) *factorial(l+k) ) \ for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k))]) ) chunk= add([ (Integer(2) *pi)**((Integer(1) -d)/Integer(2) ) *a.determinant()**(-Integer(1) /Integer(2) ) \ *alpha[d-Integer(1) ]**((Integer(1) -d)/Integer(2)  -k) *L[k] \ *n**((Integer(1) -d)/Integer(2) -k) for k in range(N) ]) chunk= chunk.subs(c).simplify() coeffs= chunk.coefficients(n) coeffs.reverse() coeffs= coeffs[:N] if numerical: # If a numerical approximation is desired... subexp_part = add( [co[Integer(0) ].subs(c).n(digits=numerical)*n**co[Integer(1) ] for co in coeffs] ) exp_scale=(Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )) \ .n(digits=numerical) else: subexp_part = add( [co[Integer(0) ].subs(c)*n**co[Integer(1) ] for co in coeffs] ) exp_scale= Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] ) return exp_scale**n*subexp_part, exp_scale, subexp_part #----------------------------------------------------------------------------- def _crit_cone_combo(self,fs,X,c,alpha): r""" This function returns an auxiliary point associated to the multiple point c of the factors fs. It is for internal use by _asymptotics_main_multiple(). INPUT: - fs - A list of expressions in the variables of X. - X - A list of variables. - c - A dictionary with keys X and values in some field. - alpha - A list of rationals. OUTPUT: A solution of the matrix equation y Gamma = a for y, where Gamma is the matrix whose jth row is _direction(_log_grad(fj,X,c)) where fj is the jth item in fs and where a is _direction(alpha). EXAMPLES:: sage: R.= PolynomialRing(QQ) sage: G,H= 1,1 sage: F= QuasiRationalExpression(G,H) sage: fs= [x + 2*y + z - 4, 2*x + y + z - 4] sage: c= {x:1,y:1,z:1} sage: alpha= [2,1,1] sage: F._crit_cone_combo(fs,R.gens(),c,alpha) [0, 1] NOTES: Use this function only when Gamma is well-defined and there is a unique solution to the matrix equation y Gamma = a. Fails otherwise. AUTHORS: - Alex Raichev (2008-10-01, 2008-11-25, 2009-03-04, 2010-09-08, 2010-12-02) """ # Assuming here that each _log_grad(f) has nonzero final component. # Then 'direction' will not throw a division by zero error. d= len(X) r= len(fs) Gamma= matrix([self._direction(self._log_grad(f,X,c)) for f in fs]) # solve_left() fails when working in SR :-(. So use solve() instead. #s= Gamma.solve_left(vector(alpha)/alpha[d-1]) V= [var('sss'+str(i)) for i in range(r)] M= matrix(V)*Gamma eqns= [M[Integer(0) ][i]== alpha[i]/alpha[d-Integer(1) ] for i in range(d)] s= solve(eqns,V,solution_dict=True)[Integer(0) ]  # Assuming a unique solution. return [s[v] for v in V] # B ========================================================================== # C ========================================================================== def cohomologous_integrand(self,alpha,asy_var=None): r""" This function takes a multivariate Cauchy type integral \int F / x^{asy_var\alpha+1} dx, where F=self, and breaks it up into a list of nicer Cauchy type integrals for the purposes of computing asymptotics of the original integral as asy_var\to\infty. The sum of the nicer integrals is de Rham cohomologous to the original integral. It assumes that algebraic varieties corresponding to the irreducible factors of self._H intersect transversely (see notes below). INPUT: - alpha - A list of positive integers or symbolic variables. - asy_var - A symbolic variable (default: None). Eventually set to var('n') if None is given. OUTPUT: A list of the form [chunk_1,\ldots,chunk_r], where each chunk_j has the form [P,[B_1,1],\ldots,[B_l,1]]. Here l \le d, P is a symbolic expression in the indeterminates of R and asy_var, \{B_1,\ldots,B_l\} \subseteq \{Q_1,\ldots,Q_m\}, and [P,[B_1,1],\ldots,[B_l,1]] represents the integral \int P/(B_1 \cdots B_l) dx. EXAMPLES:: sage: R.= PolynomialRing(QQ) sage: G= 9*exp(x+y) sage: H= (3-2*x-y)*(3-x-2*y) sage: F= QuasiRationalExpression(G,H) sage: alpha= [4,3] sage: F.cohomologous_integrand(alpha) [[9*e^(x + y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]] sage: R.= PolynomialRing(QQ) sage: G= 9*exp(x+y) sage: H= (3-2*x-y)^2*(3-x-2*y) sage: F= QuasiRationalExpression(G,H) sage: alpha= [4,3] sage: F.cohomologous_integrand(alpha) [[-3*(3*x*e^x - 8*y*e^x)*n*e^y/(x*y) - 3*((x - 2)*y*e^x + x*e^x)*e^y/(x*y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]] sage: R.= PolynomialRing(QQ) sage: G= 16 sage: H= (4-2*x-y-z)^2*(4-x-2*y-z) sage: F= QuasiRationalExpression(G,H) sage: alpha= [3,3,2] sage: F.cohomologous_integrand(alpha) [[16*(4*y - 3*z)*n/(y*z) + 16*(2*y - z)/(y*z), [x + 2*y + z - 4, 1], [2*x + y + z - 4, 1]]] NOTES: The varieties  corresponding to Q_1,\ldots,Q_m __intersect transversely__ iff for each point c of their intersection and for all k \le l, the Jacobian matrix of any k polynomials from \{Q_1,\ldots,Q_m\} has rank equal to \min\{k,d\} when evaluated at c. ALGORITHM: Let asy_var= n and consider the integral around a polycirle centered at the origin of the d-variate differential form \frac{G(x) dx_1 \wedge\cdots\wedge dx_d}{H(x) x^{n\alpha+1}}, where G=self._G and H=self._H. (1) Decompose G/H into a sum of partial fractions P_1 +\cdots+ P_r so that each term of the sum has at most d irreducible factors of H in the denominator. (2) For each differential form P_j dx_1 \wedge\cdots\wedge dx_d, find an equivalent form \omega_j in de Rham cohomology with no repeated irreducible factors of H in its denominator. AUTHOR: - Alex Raichev (2010-09-22) """ if not asy_var: asy_var = var('n') # Create symbolic (non-ring) variables. G= self._G H= self._H R= self._R d= self._d X= self._variables # Prepare input for partial_fraction_decomposition() which only works # for functions in the field of fractions of R. if G in R: numer=Integer(1) F= G/H else: numer= G F= R(Integer(1) )/H nstuff= Integer(1) /mul([X[j]**(alpha[j]*asy_var+Integer(1) ) for j in range(d)]) # Do steps (2) and (1). integrands= [] whole,parts= partial_fraction_decomposition(F) for f in parts: a= format_quotient(f) integrands.append( [a[Integer(0) ]*numer*nstuff] + a[Integer(1) :] ) # Do step (3). ce= decompose_via_cohomology(integrands) ce_new= [] for a in ce: ce_new.append( [(a[Integer(0) ]/nstuff).simplify_full().collect(n)] + a[Integer(1) :] ) return ce_new #------------------------------------------------------------------------------- def critical_cone(self,c,coordinate=None): r""" Returns the critical cone of a convenient multiple point c. INPUT: - c - A dictionary with keys self.variables() and values in a field. - coordinate - (optional) A natural number. OUTPUT: A list of vectors that generate the critical cone of c and the cone itself if the values of c lie in QQ. EXAMPLES:: sage: R.= PolynomialRing(QQ) sage: G= 1 sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) sage: F= QuasiRationalExpression(G,H) sage: c= {z: 4/3, y: 1, x: 1/2} sage: F.critical_cone(c) ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N) NOTES: The _critical cone_ of a convenient multiple point c with with c_k \del_k H_j(c) \neq 0 for all j=1,\ldots,r is the conical hull of the vectors \gamma_j(c) = \left(\frac{c_1 \del_1 H_j(c)}{c_k \del_k H_j(c)},\ldots, \frac{c_d \del_d H_j(c)}{c_k \del_k H_j(c)} \right). Here H_1,\ldots,H_r are the irreducible germs of self._H around c. For more details, see [RaWi2011]_. If this function's optional argument coordinate isn't given, then this function searches (from d down to 1) for the first index k such that for all j=1,\ldots,r we have c_k \del_k H_j(c) \neq 0 and sets coordinate = k. Almost. Since this is Python, all the indices actually start at 0. REFERENCES: .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", To appear. AUTHORS: - Alex Raichev (2010-08-25) """ Hs= [SR(h[Integer(0) ]) for h in self._Hfac]  # irreducible factors of H X= self._variables d= self._d # Ensure the variables of c lie in SR cc= {} for x in c.keys(): cc[SR(x)] = c[x] lg= [self._log_grad(h,X,cc) for h in Hs] if not coordinate: # Search (from d-1 down to 0) for a coordinate k such that # for all h in Hs we have cc[k] * diff(h,X[k]) !=0. # One is guaranteed to exist in the case of a convenient multiple # point. for k in reversed(range(d)): if Integer(0)  not in [v[k] for v in lg]: coordinate= k break gamma= [self._direction(v,coordinate) for v in lg] if [[gg in QQ for gg in g] for g in gamma] == \ [[True for gg in g] for g in gamma]: return gamma,Cone(gamma)    # Cone() needs rational vectors else: return gamma # D ============================================================================ def _diff_all(self,f,V,n,ending=[],sub=None,sub_final=None,zero_order=Integer(0) ,rekey=None): r""" This function returns a dictionary of representative mixed partial derivatives of f from order 1 up to order n with respect to the variables in V. The default is to key the dictionary by all nondecreasing sequences in V of length 1 up to length n. For internal use. Does not use self. INPUT: - f - An individual or list of \mathcal{C}^{n+1} functions. - V - A list of variables occurring in f. - n - A natural number. - ending - A list of variables in V. - sub - An individual or list of dictionaries. - sub_final - An individual or list of dictionaries. - rekey - A callable symbolic function in V or list thereof. - zero_order - A natural number. OUTPUT: The dictionary {s_1:deriv_1,...,s_r:deriv_r}. Here s_1,\ldots,s_r is a listing of all nondecreasing sequences of length 1 up to length n over the alphabet V, where w > v in X  iff str(w) > str(v), and deriv_j is the derivative of f with respect to the derivative sequence s_j and simplified with respect to the substitutions in sub and evaluated at sub_final. Moreover, all derivatives with respect to sequences of length less than zero_order (derivatives of order less than zero_order) will be made zero. If rekey is nonempty, then s_1,\ldots,s_r will be replaced by the symbolic derivatives of the functions in rekey. If ending is nonempty, then every derivative sequence s_j will be suffixed by ending. EXAMPLES: I'd like to print the entire dictionaries, but that doesn't yield consistent output order for doctesting. Order of keys changes. :: sage: R. = PolynomialRing(QQ) sage: G,H = 1,1 sage: F= QuasiRationalExpression(G,H) sage: f= function('f',x) sage: dd= F._diff_all(f,[x],3) sage: dd[(x,x,x)] D[0, 0, 0](f)(x) :: sage: d1= {diff(f,x): 4*x^3} sage: dd= F._diff_all(f,[x],3,sub=d1) sage: dd[(x,x,x)] 24*x :: sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f) sage: dd[diff(f,x,3)] 24*x :: sage: a= {x:1} sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f,sub_final=a) sage: dd[diff(f,x,3)] 24 :: sage: X= var('x,y,z') sage: f= function('f',*X) sage: dd= F._diff_all(f,X,2,ending=[y,y,y]) sage: dd[(z,y,y,y)] D[1, 1, 1, 2](f)(x, y, z) :: sage: g= function('g',*X) sage: dd= F._diff_all([f,g],X,2) sage: dd[(0,y,z)] D[1, 2](f)(x, y, z) :: sage: dd[(1,z,z)] D[2, 2](g)(x, y, z) :: sage: f= exp(x*y*z) sage: ff= function('ff',*X) sage: dd= F._diff_all(f,X,2,rekey=ff) sage: dd[diff(ff,x,z)] x*y^2*z*e^(x*y*z) + y*e^(x*y*z) AUTHORS: - Alex Raichev (2008-10-01, 2009-04-01, 2010-02-01) """ singleton=False if not isinstance(f,list): f= [f] singleton=True # Build the dictionary of derivatives iteratively from a list of nondecreasing # derivative-order sequences. derivs= {} r= len(f) if ending: seeds = [ending] start = Integer(1) else: seeds = [[v] for v in V] start = Integer(2) if singleton: for s in seeds: derivs[tuple(s)] = self._subs_all(diff(f[Integer(0) ],s),sub) for l in (ellipsis_range(start,Ellipsis,n)): for t in UnorderedTuples(V,l): s= tuple(t + ending) derivs[s] = self._subs_all(diff(derivs[s[Integer(1) :]],s[Integer(0) ]),sub) else: # Make the dictionary keys of the form (j,sequence of variables), # where j in range(r). for s in seeds: value= self._subs_all([diff(f[j],s) for j in range(r)],sub) derivs.update(dict([(tuple([j]+s),value[j]) for j in range(r)])) for l in (ellipsis_range(start,Ellipsis,n)): for t in UnorderedTuples(V,l): s= tuple(t + ending) value= self._subs_all(\ [diff(derivs[(j,)+s[Integer(1) :]],s[Integer(0) ]) for j in range(r)],sub) derivs.update(dict([((j,)+s,value[j]) for j in range(r)])) if zero_order: # Zero out all the derivatives of order < zero_order if singleton: for k in derivs.keys(): if len(k) < zero_order: derivs[k]= Integer(0) else: # Ignore the first of element of k, which is an index. for k in derivs.keys(): if len(k)-Integer(1)  < zero_order: derivs[k]= Integer(0) if sub_final: # Substitute sub_final into the values of derivs. for k in derivs.keys(): derivs[k] = self._subs_all(derivs[k],sub_final) if rekey: # Rekey the derivs dictionary by the value of rekey. F= rekey if singleton: # F must be a singleton. derivs= dict( [(diff(F,list(k)), derivs[k]) for k in derivs.keys()] ) else: # F must be a list. derivs= dict( [(diff(F[k[Integer(0) ]],list(k)[Integer(1) :]), derivs[k]) for k in derivs.keys()] ) return derivs #------------------------------------------------------------------------------- def _diff_op(self,A,B,AB_derivs,V,M,r,N): r""" This function computes the derivatives DD^(l+k)(A[j] B^l) evaluated at a point p for various natural numbers j,k,l which depend on r and N. Here DD is a specific second-order linear differential operator that depends on M, A is a list of symbolic functions, B is symbolic function, and AB_derivs contains all the derivatives of A and B evaluated at p that are necessary for the computation. For internal use by the functions _asymptotics_main_multiple() and _asymptotics_main_smooth(). Does not use self. INPUT: - A - A single or length r list of symbolic functions in the variables V. - B - A symbolic function in the variables V. - AB_derivs - A dictionary whose keys are the (symbolic) derivatives of A[0],\ldots,A[r-1] up to order 2N-2 and the (symbolic) derivatives of B up to order 2N. The values of the dictionary are complex numbers that are the keys evaluated at a common point p. - V - The variables of the A[j] and B. - M - A symmetric l \times l matrix, where l is the length of V. - r,N - Natural numbers. OUTPUT: A dictionary whose keys are natural number tuples of the form (j,k,l), where l \le 2k, j \le r-1, and j+k \le N-1, and whose values are DD^(l+k)(A[j] B^l) evaluated at a point p, where DD is the linear second-order differential operator -\sum_{i=0}^{l-1} \sum_{j=0}^{l-1} M[i][j] \partial^2 /(\partial V[j] \partial V[i]). EXAMPLES:: sage: R. = PolynomialRing(QQ) sage: G,H = 1,1 sage: F= QuasiRationalExpression(G,H) sage: T= var('x,y') sage: A= function('A',*tuple(T)) sage: B= function('B',*tuple(T)) sage: AB_derivs= {} sage: M= matrix([[1,2],[2,1]]) sage: DD= F._diff_op(A,B,AB_derivs,T,M,1,2) sage: DD.keys() [(0, 1, 2), (0, 1, 1), (0, 1, 0), (0, 0, 0)] sage: len(DD[(0,1,2)]) 246 AUTHORS: - Alex Raichev (2008-10-01, 2010-01-12) """ if not isinstance(A,list): A= [A] # First, compute the necessary product derivatives of A and B. product_derivs= {} for (j,k) in mrange([r,N]): if j+k = PolynomialRing(QQ) sage: G,H = 1,1 sage: F= QuasiRationalExpression(G,H) sage: A= function('A',x) sage: B= function('B',x) sage: AB_derivs= {} sage: F._diff_op_simple(A,B,AB_derivs,x,3,2,2) {(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)} AUTHORS: - Alex Raichev (2010-01-15) """ I= sqrt(-Integer(1) ) DD= {} if v.mod(Integer(2) ) == Integer(0) : for k in range(N): for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)): DD[(k,l)] = (a**(-Integer(1) /v))**(Integer(2) *k+v*l) \ * diff(A*B**l,x,Integer(2) *k+v*l).subs(AB_derivs) else: for k in range(N): for l in (ellipsis_range(Integer(0) ,Ellipsis,k)): DD[(k,l)] = (abs(a)**(-Integer(1) /v)*I*a/abs(a))**(k+v*l) \ * diff(A*B**l,x,k+v*l).subs(AB_derivs) return DD #------------------------------------------------------------------------------- def _diff_prod(self,f_derivs,u,g,X,interval,end,uderivs,atc): r""" This function takes various derivatives of the equation f=ug, evaluates at a point c, and solves for the derivatives of u. For internal use by the function _asymptotics_main_multiple(). Does not use self. INPUT: - f_derivs - A dictionary whose keys are tuples \code{s + end} for all X-variable sequences s with length in interval and whose values are the derivatives of a function f evaluated at c. - u - A callable symbolic function. - g - An expression or callable symbolic function. - X - A list of symbolic variables. - interval - A list of positive integers. Call the first and last values n and nn, respectively. - end - A possibly empty list of z's where z is the last element of X. - uderivs - A dictionary whose keys are the symbolic derivatives of order 0 to order n-1 of u evaluated at c and whose values are the corresponding derivatives evaluated at c. - atc - A dictionary whose keys are the keys of c and all the symbolic derivatives of order 0 to order nn of g evaluated c and whose values are the corresponding derivatives evaluated at c. OUTPUT: A dictionary whose keys are the derivatives of u up to order nn and whose values are those derivatives evaluated at c. EXAMPLES:: I'd like to print out the entire dictionary, but that does not give consistent output for doctesting. Order of keys changes :: sage: R. = PolynomialRing(QQ) sage: G,H = 1,1 sage: F= QuasiRationalExpression(G,H) sage: u= function('u',x) sage: g= function('g',x) sage: dd= F._diff_prod({(x,):1,(x,x):1},u,g,[x],[1,2],[],{u(x=2):1},{x:2,g(x=2):3,diff(g,x)(x=2):5, diff(g,x,x)(x=2):7}) sage: dd[diff(u,x,2)(x=2)] 22/9 NOTES: This function works by differentiating the equation f=ug with respect to the variable sequence \code{s+end}, for all tuples s of X of lengths in interval, evaluating at the point c, and solving for the remaining derivatives of u. This function assumes that u never appears in the differentiations of f=ug after evaluating at c. AUTHORS: - Alex Raichev (2009-05-14, 2010-01-21) """ for l in interval: D= {} rhs= [] lhs= [] for t in UnorderedTuples(X,l): s= t+end lhs.append(f_derivs[tuple(s)]) rhs.append(diff(u*g,s).subs(atc).subs(uderivs)) # Since Sage's solve command can't take derivatives as variable # names, i make new variables based on t to stand in for # diff(u,t) and store them in D. D[diff(u,t).subs(atc)] = self._make_var([var('zing')]+t) eqns=[ lhs[i] == rhs[i].subs(uderivs).subs(D) for i in range(len(lhs))] vars= D.values() sol= solve(eqns,*vars,solution_dict=True) uderivs.update(self._subs_all(D,sol[Integer(0) ])) return uderivs #------------------------------------------------------------------------------- def _diff_seq(self,V,s): r""" Given a list s of tuples of natural numbers, this function returns the list of elements of V with indices the elements of the elements of s. This function is for internal use by the function _diff_op(). Does not use self. INPUT: - V - A list. - s - A list of tuples of natural numbers in the interval \code{range(len(V))}. OUTPUT: The tuple \code{tuple([V[tt] for tt in sorted(t)])}, where t is the list of elements of the elements of s. EXAMPLES:: sage: R. = PolynomialRing(QQ) sage: G,H = 1,1 sage: F= QuasiRationalExpression(G,H) sage: V= list(var('x,t,z')) sage: F._diff_seq(V,([0,1],[0,2,1],[0,0])) (x, x, x, x, t, t, z) AUTHORS: - Alex Raichev (2009.05.19) """ t= [] for ss in s: t.extend(ss) return tuple([V[tt] for tt in sorted(t)]) #------------------------------------------------------------------------------- def _direction(self,v,coordinate=None): r""" This function returns \code{[vv/v[coordinate] for vv in v]} where coordinate is the last index of v if not specified otherwise. Does not use self. INPUT: - v - A vector. - coordinate - An index for v (default: None). OUTPUT: This function returns \code{[vv/v[coordinate] for vv in v]} where coordinate is the last index of v if not specified otherwise. EXAMPLES:: sage: R. = PolynomialRing(QQ) sage: G,H = 1,1 sage: F= QuasiRationalExpression(G,H) sage: F._direction([2,3,5]) (2/5, 3/5, 1) sage: F._direction([2,3,5],0) (1, 3/2, 5/2) AUTHORS: - Alex Raichev (2010-08-25) """ if coordinate == None: coordinate= len(v)-Integer(1) try: v[Integer(0) ].variables() return tuple([(vv/v[coordinate]).simplify_full() for vv in v]) except: return tuple([vv/v[coordinate] for vv in v]) # E ============================================================================ # F ============================================================================ # G ============================================================================ # H ============================================================================ # I ============================================================================ def is_cmp(self,points): r""" Checks if the points in the list points are convenient multiple points of V= \{ x\in CC^d : H(x) = 0\}, where H=self._H. INPUT: - points - An individual or list of dictionaries with keys self._variables and values in some superfield of self._R.base_ring(). OUTPUT: A list of tuples (p,verdict,comment), one for each point p in points, where verdict is True if p is a convenient multiple point and False otherwise, and where comment is a string comment relating to verdict, such as 'not a transverse intersection'. EXAMPLES:: sage: R.= PolynomialRing(QQ) sage: G= 16 sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) sage: F= QuasiRationalExpression(G,H) sage: points= [{x:1/2,y:1,z:4/3},{x:1/2,y:1,z:2}] sage: F.is_cmp(points) [({y: 1, z: 4/3, x: 1/2}, True, 'all good'), ({y: 1, z: 2, x: 1/2}, False, 'not a singular point')] NOTES: A point c of V is a __convenient multiple point__ if V is locally a union of complex manifolds that intersect transversely at c; see [RaWi2011]_. REFERENCES: .. [RaWi2011]  Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", To appear. AUTHORS: - Alex Raichev (2011-04-18) """ H= self._H d= self._d Hs= [SR(h[Integer(0) ]) for h in self._Hfac]  # irreducible factors of H X= self._variables J= [tuple([diff(h,x) for x in X]) for h in Hs] verdicts= [] if not isinstance(points,list): points= [points] for p in points: # Ensure variables in points lie in SR. pp= {} for x in p.keys(): pp[SR(x)]= p[x] # Test 1: Is p a zero of all factors of H? if [h.subs(pp) for h in Hs] != [Integer(0)  for h in Hs]: # Failed test 1.  Move on to next point. verdicts.append((p,False,'not a singular point')) continue # Test 2: Are the factors of H smooth and # do theyintersect transversely at p? J= [tuple([f.subs(pp) for f in dh]) for dh in J] l= len(J) if Set(J).cardinality() < l: # Fail.  Move on to next point. verdicts.append((p,False,'not a transverse intersection')) continue temp= True for S in list(Set(J).subsets())[Integer(1) :]: # Subsets of size >= 1 k= len(list(S)) M = matrix(list(S)) if  M.rank() != min(k,d): # Fail. temp= False verdicts.append((p,False,'not a transvere intersection')) break if not temp: # Move on to next point continue # Test 3: Is p convenient? Jlog= matrix([self._log_grad(h,X,pp) for h in Hs]) if [Integer(0)  in f for f in Jlog.columns()] ==\ [True for f in Jlog.columns()]: # Fail.  Move on to next point. verdict[p] = (False,'multiple point but not convenient') continue verdicts.append((p,True,'all good')) return verdicts # J ============================================================================ # K ============================================================================ # L ============================================================================ def _log_grad(self,f,X,c): r""" This function returns the logarithmic gradient of f with respect to the variables of X evalutated at c. Does not use self. INPUT: - f - An expression in the variables of X. - X - A list of variables. - c - A dictionary with keys X and values in a field K. OUTPUT: \code{[c[x] * diff(f,x).subs(c) for x in X]}. EXAMPLES:: sage: R. = PolynomialRing(QQ) sage: G,H = 1,1 sage: F= QuasiRationalExpression(G,H) sage: X= var('x,y,z') sage: f= x*y*z^2 sage: c= {x:1,y:2,z:3} sage: f.gradient() (y*z^2, x*z^2, 2*x*y*z) sage: F._log_grad(f,X,c) (18, 18, 36) :: sage: R.= PolynomialRing(QQ) sage: f= x*y*z^2 sage: c= {x:1,y:2,z:3} sage: F._log_grad(f,R.gens(),c) (18, 18, 36) AUTHORS: - Alex Raichev (2009-03-06) """ return tuple([SR(c[x] * diff(f,x).subs(c)).simplify() for x in X]) # M ============================================================================ def _make_var(self,L): r""" This function converts the list L to a string (without commas) and returns the string as a variable. For internal use by the function _diff_op() Does not use self. INPUT: - L - A list. OUTPUT: A variable whose name is the concatenation of the variable names in L. EXAMPLES:: sage: R. = PolynomialRing(QQ) sage: G,H = 1,1 sage: F= QuasiRationalExpression(G,H) sage: L= list(var('x,y,hello')) sage: v= F._make_var(L) sage: print v, type(v) xyhello AUTHORS: - Alex Raichev (2010-01-21) """ return var(''.join([str(v) for v in L])) # N ============================================================================ def _new_var_name(self,name,V): r""" This function returns the first string in the sequence name, name+name, name+name+name,... that does not appear in the list V. It is for internal use by the function _asymptotics_main_multiple(). Does not use self. INPUT: - name - A string. - V - A list of variables. OUTPUT: The first string in the sequence name, name+name, name+name+name,... that does not appear in the list \code{str(V)}. EXAMPLES:: sage: R. = PolynomialRing(QQ) sage: G,H = 1,1 sage: F= QuasiRationalExpression(G,H) sage: X= var('x,xx,y,z') sage: F._new_var_name('x',X) 'xxx' AUTHORS: - Alex Raichev (2008-10-01) """ newname= name while newname in str(V): newname= newname +name return newname # O ============================================================================ # P ============================================================================ # Q ============================================================================ # R ============================================================================ def relative_error(self,approx,alpha,interval,exp_scale=Integer(1) ): r""" Returns the relative error between the values of the Maclaurin coefficients of self with multi-indices m alpha for m in interval and the values of the functions in approx. INPUT: - approx - An individual or list of symbolic expressions in one variable. - alpha - A list of positive integers. - interval - A list of positive integers. - exp_scale - (optional) A number.  Default: 1. OUTPUT: A list whose entries are of the form [m\alpha,a_m,b_m,err_m] for m \in interval. Here m\alpha is a tuple; a_m is the m alpha (multi-index) coefficient of the Maclaurin series for F divided by exp_scale^m; b_m is a list of the values of the functions in approx evaluated at m and divided by exp_scale^m; err_m is the list of relative errors (a_m-f)/a_m for f in b_m. All outputs are decimal approximations. EXAMPLES:: sage: R.= PolynomialRing(QQ) sage: G=1 sage: H= 1-x-y-x*y sage: F= QuasiRationalExpression(G,H) sage: alpha= [1,1] sage: var('n') n sage: f= (0.573/sqrt(n))*5.83^n sage: es= 5.83 sage: F.relative_error(f,alpha,[1,2,4,8],es) Calculating errors table in the form exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... [[(1, 1), 0.514579759862779, [0.573000000000000], [-0.113530000000000]], [(2, 2), 0.382477808931739, [0.405172185619892], [-0.0593351461396876]], [(4, 4), 0.277863059517142, [0.286500000000000], [-0.0310834426780842]], [(8, 8), 0.199108827584423, [0.202586092809946], [-0.0174641439443390]]] sage: g= (0.573/sqrt(n) - 0.0674/n^(3/2))*5.83^n sage: F.relative_error([f,g],alpha,[1,2,4,8],es) Calculating errors table in the form exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... [[(1, 1), 0.514579759862779, [0.573000000000000, 0.505600000000000], [-0.113530000000000, 0.0174506666666667]], [(2, 2), 0.382477808931739, [0.405172185619892, 0.381342687093905], [-0.0593351461396876, 0.00296781097184384]], [(4, 4), 0.277863059517142, [0.286500000000000, 0.278075000000000], [-0.0310834426780842, -0.000762751562681505]], [(8, 8), 0.199108827584423, [0.202586092809946, 0.199607405494198], [-0.0174641439443390, -0.00250404723800224]]] AUTHORS: - Alex Raichev (2009-05-18, 2011-04-18) """ if not isinstance(approx,list): approx= [approx] av= approx[Integer(0) ].variables()[Integer(0) ] print "Calculating errors table in the form" print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..." # Get Maclaurin coefficients of self and scale them. # Then compute errors. n= interval[-Integer(1) ] mc= self.maclaurin_coefficients(alpha,n) mca={} stats=[] for m in interval: beta= tuple([m*a for a in alpha]) mc[beta]= mc[beta]/exp_scale**m mca[beta]= [f.subs({av:m})/exp_scale**m for f in approx] stats_row= [beta, mc[beta].n(), [a.n() for a in mca[beta]]] if mc[beta]==Integer(0) : stats_row.extend([None for a in mca[beta]]) else: stats_row.append([((mc[beta]-a)/mc[beta]).n() for a in mca[beta]]) stats.append(stats_row) return stats # S ============================================================================ def singular_points(self): r""" This function returns a Groebner basis ideal whose variety is the set of singular points of the algebraic variety V= \{x\in\CC^d : H(x)=0\}, where H=sef._H. INPUT: OUTPUT: A Groebner basis ideal whose variety is the set of singular points of the algebraic variety V= \{x\in\CC^d : H(x)=0\}. EXAMPLES:: sage: R.= PolynomialRing(QQ) sage: G= 1 sage: H= (4-2*x-y-z)*(4-x-2*y-z) sage: F= QuasiRationalExpression(G,H) sage: F.singular_points() Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field AUTHORS: - Alex Raichev (2008-10-01, 2008-11-20, 2010-12-03, 2011-04-18) """ H= self._H R= self._R f= R.ideal(H).radical().gens()[Integer(0) ]     # Compute the reduction of H. J= R.ideal([f] + f.gradient()) return R.ideal(J.groebner_basis()) #------------------------------------------------------------------------------- def smooth_critical(self,alpha): r""" This function returns a Groebner basis ideal whose variety is the set of smooth critical points of the algebraic variety V= \{x\in\CC^d : H(x)=0\} for the direction \alpha where H=self._H. INPUT: - alpha - A d-tuple of positive integers and/or symbolic entries. OUTPUT: A Groebner basis ideal of smooth critical points of V for \alpha. EXAMPLES:: sage: R. = PolynomialRing(QQ) sage: G=1 sage: H = (1-x-y-x*y)^2 sage: F= QuasiRationalExpression(G,H) sage: var('a1,a2') (a1, a2) sage: F.smooth_critical([a1,a2]) 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 NOTES: A point c of V is a __smooth critical point for alpha__ if the gradient of f at c is not identically zero and \alpha is in the span of the logarithmic gradient vector (c[0] \partial_1 f(c)),\ldots,c[d-1] \partial_d f(c)); see [RaWi2008a]_. REFERENCES: .. [RaWi2008a]  Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", Electronic Journal of Combinatorics, Vol. 15 (2008), R89. AUTHORS: - Alex Raichev (2008-10-01, 2008-11-20, 2009-03-09, 2010-12-02, 2011-04-18) """ H= self._H R= H.parent() B= R.base_ring() d= R.ngens() vars= R.gens() f= R.ideal(H).radical().gens()[Integer(0) ]   # Compute the reduction of H. # Expand B by the variables of alpha if there are any. indets= [] indets_ind= [] for a in alpha: if not ((a in ZZ) and (a>Integer(0) )): try: CC(a) except: indets.append(var(a)) indets_ind.append(alpha.index(a)) else: print "The components of", alpha, \ "must be positive integers or symbolic variables." return indets= list(Set(indets))   # Delete duplicates in indets. if indets != []: BB= FractionField(PolynomialRing(B,tuple(indets))) S= R.change_ring(BB) vars= S.gens() # Coerce alpha into BB. for i in range(len(alpha)): alpha[i] = BB(alpha[i]) else: S= R # Find smooth, critical points for alpha. f= S(f) J= S.ideal([f] +[ alpha[d-Integer(1) ]*vars[i]*diff(f,vars[i]) \ -alpha[i]*vars[d-Integer(1) ]*diff(f,vars[d-Integer(1) ]) for i in range(d-Integer(1) )]) return S.ideal(J.groebner_basis()) #------------------------------------------------------------------------------- def _subs_all(self,f,sub,simplify=False): r""" This function returns the items of f substituted by the dictionaries of sub in order of their appearance in sub. Does not use self. INPUT: - f - An individual or list of symbolic expressions or dictionaries - sub - An individual or list of dictionaries. - simplify - Boolean (default: False). OUTPUT: The items of f substituted by the dictionaries of sub in order of their appearance in sub.  The subs() command is used. If simplify is True, then simplify() is used after substitution. EXAMPLES:: sage: R. = PolynomialRing(QQ) sage: G,H = 1,1 sage: F= QuasiRationalExpression(G,H) sage: var('x,y,z') (x, y, z) sage: a= {x:1} sage: b= {y:2} sage: c= {z:3} sage: F._subs_all(x+y+z,a) y + z + 1 sage: F._subs_all(x+y+z,[c,a]) y + 4 sage: F._subs_all([x+y+z,y^2],b) [x + z + 2, 4] sage: F._subs_all([x+y+z,y^2],[b,c]) [x + 5, 4] :: sage: var('x,y') (x, y) sage: a= {'foo':x^2+y^2, 'bar':x-y} sage: b= {x:1,y:2} sage: F._subs_all(a,b) {'foo': 5, 'bar': -1} AUTHORS: - Alex Raichev (2009-05-05) """ singleton= False if not isinstance(f,list): f= [f] singleton= True if not isinstance(sub,list): sub= [sub] g= [] for ff in f: for D in sub: if isinstance(ff,dict): ff= dict( [(k,ff[k].subs(D)) for k in ff.keys()] ) else: ff= ff.subs(D) g.append(ff) if singleton and simplify: if isinstance(g[Integer(0) ],dict): return g[Integer(0) ] else: return g[Integer(0) ].simplify() elif singleton and not simplify: return g[Integer(0) ] elif not singleton and simplify: G= [] for gg in g: if isinstance(gg,dict): G.append(gg) else: G.append(gg.simplify()) return G else: return g # T ============================================================================ def maclaurin_coefficients(self,alpha,n): r""" Returns the Maclaurin coefficients of self that have multi-indices alpha, 2*alpha,...,n*alpha. INPUT: - n - positive integer - alpha - tuple of positive integers representing a multi-index. OUTPUT: A dictionary of the form (beta, beta Maclaurin coefficient of self). AUTHORS: - Alex Raichev (2011-04-08) """ # Do all computations in the Symbolic Ring. F= SR(self._G/self._H) v= F.variables() p= {} for x in v: p[x]=Integer(0) d= len(v) coeffs={} # Initialize sequence of variables to differentiate with. s=[] for i in range(d): s.extend([v[i] for j in range(alpha[i])]) F_deriv= diff(F,s) coeffs[tuple(alpha)]= F_deriv.subs(p)/ mul([factorial(a) for a in alpha]) old_beta= alpha for k in (ellipsis_range(Integer(2) ,Ellipsis,n)): # Update variable sequence to differentiate with. beta= [k*a for a in alpha] delta= [beta[i]-old_beta[i] for i in range(d)] s= [] for i in range(d): s.extend([v[i] for j in range(delta[i])]) F_deriv= diff(F_deriv,s) coeffs[tuple(beta)]= F_deriv.subs(p)/ mul([factorial(b) for b in beta]) old_beta= copy(beta) return coeffs # U ============================================================================ # V ============================================================================ def variables(self): r""" Returns the tuple of variables of self. EXAMPLES:: sage: R.= PolynomialRing(QQ) sage: G= exp(x) sage: H= 1-y sage: F= QuasiRationalExpression(G,H) sage: F.variables() (x, y) sage: R.= PolynomialRing(QQ,order='invlex') sage: G= exp(x) sage: H= 1-y sage: F= QuasiRationalExpression(G,H) sage: F.variables() (y, x) AUTHORS: - Alex Raichev (2011-04-01) """ return self._variables # W ============================================================================ # X ============================================================================ # Y ============================================================================ # Z ============================================================================