# Ticket #10519: mgf.sage

File mgf.sage, 89.3 KB (added by , 12 years ago) |
---|

Line | |
---|---|

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 |