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