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