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