# Ticket #10519: amgf.py

File amgf.py, 111.3 KB (added by , 11 years ago) |
---|

Line | |
---|---|

1 | from sage.all_cmdline import * # import sage library |

2 | r""" |

3 | This code relates to analytic combinatorics. |

4 | More specifically, it is a collection of functions designed |

5 | to compute asymptotics of Maclaurin coefficients of certain classes of |

6 | multivariate generating functions. |

7 | |

8 | The main function asymptotics() returns the first `N` terms of |

9 | the asymptotic expansion of the Maclaurin coefficients `F_{n\alpha}` |

10 | of the multivariate meromorphic function `F=G/H` as `n\to\infty`. |

11 | It assumes that `F` is holomorphic in a neighborhood of the origin, |

12 | that `H` is a polynomial, and that asymptotics in the direction of |

13 | `\alpha` (a tuple of positive integers) are controlled by convenient |

14 | smooth or multiple points. |

15 | |

16 | For an introduction to this subject, see [PeWi2008]_. |

17 | The main algorithms and formulas implemented here come from [RaWi2008a]_ |

18 | and [RaWi2011]_. |

19 | |

20 | REFERENCES: |

21 | |

22 | .. [AiYu1983] I. A. A\u\izenberg and A. P. Yuzhakov, "Integral |

23 | representations and residues in multidimensional complex analysis", |

24 | Translations of Mathematical Monographs, 58. American Mathematical |

25 | Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X. |

26 | |

27 | .. [DeLo2006] Wolfram Decker and Christoph Lossen, "Computing in |

28 | algebraic geometry", Chapter 7.1, Springer-Verlag, 2006. |

29 | |

30 | .. [DiEm2005] Alicia Dickenstein and Ioannis Z. Emiris (editors), |

31 | "Solving polynomial equations", Chapter 9.0, Springer-Verlag, 2005. |

32 | |

33 | .. [Lein1978] E. K. Leinartas, "On expansion of rational functions of |

34 | several variables into partial fractions", Soviet Math. (Iz. VUZ) |

35 | 22 (1978), no. 10, 35--38. |

36 | |

37 | .. [PeWi2008] Robin Pemantle and Mark C. Wilson, "Twenty combinatorial |

38 | examples of asymptotics derived from multivariate generating |

39 | functions", SIAM Rev. (2008) vol. 50 (2) pp. 199-272. |

40 | |

41 | .. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

42 | coefficients of multivariate generating functions: improvements |

43 | for smooth points", Electronic Journal of Combinatorics, Vol. 15 |

44 | (2008), R89. |

45 | |

46 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

47 | coefficients of multivariate generating functions: improvements |

48 | for smooth points", To appear. |

49 | |

50 | AUTHORS: |

51 | |

52 | - Alex Raichev (2008-10-01) : Initial version |

53 | - Alex Raichev (2010-09-28) : Corrected many functions |

54 | - Alex Raichev (2010-12-15) : Updated documentation |

55 | - Alex Raichev (2011-03-09) : Fixed a division by zero bug in relative_error() |

56 | - Alex Raichev (2011-04-26) : Rewrote in object-oreinted style |

57 | - Alex Raichev (2011-05-06) : Fixed bug in cohomologous_integrand() and |

58 | fixed _crit_cone_combo() to work in SR |

59 | |

60 | EXAMPLES:: |

61 | |

62 | These come from [RaWi2008a]_ and [RaWi2011]_. |

63 | A smooth point example. :: |

64 | |

65 | sage: from sage.combinat.amgf import * |

66 | sage: R.<x,y>= PolynomialRing(QQ) |

67 | sage: G= 1 |

68 | sage: H= 1-x-y-x*y |

69 | sage: F= QuasiRationalExpression(G,H) |

70 | sage: alpha= [3,2] |

71 | sage: F.maclaurin_coefficients(alpha,8) |

72 | {(21, 14): 1279919346549, (6, 4): 1289, (3, 2): 25, (18, 12): 19403906105, (24, 16): 85275509086721, (15, 10): 298199265, (9, 6): 75517, (12, 8): 4673345} |

73 | sage: I= F.smooth_critical(alpha); I |

74 | Ideal (y^2 + 3*y - 1, x - 2/3*y - 1/3) of Multivariate Polynomial Ring in x, y over Rational Field |

75 | sage: s= solve(I.gens(),F.variables(),solution_dict=true); s |

76 | [{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}] |

77 | sage: c= s[1] |

78 | sage: asys= [F.asymptotics(c,alpha,n,numerical=3) for n in [1..2]]; asys |

79 | Initializing auxiliary functions... |

80 | Calculating derivatives of auxiallary functions... |

81 | Calculating derivatives of more auxiliary functions... |

82 | Calculating actions of the second order differential operator... |

83 | Initializing auxiliary functions... |

84 | Calculating derivatives of auxiallary functions... |

85 | Calculating derivatives of more auxiliary functions... |

86 | Calculating actions of the second order differential operator... |

87 | [(0.369*71.2^n/sqrt(n), 71.2, 0.369/sqrt(n)), ((0.369/sqrt(n) - 0.0186/n^(3/2))*71.2^n, 71.2, 0.369/sqrt(n) - 0.0186/n^(3/2))] |

88 | sage: F.relative_error([f[0] for f in asys],alpha,[1,2,4,8,16],asys[0][1]) |

89 | Calculating errors table in the form |

90 | exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... |

91 | [[(3, 2), 0.351196289062500, [0.369140625000000, 0.351000000000000], [-0.0510940551757812, 0.00174000000000000]], [(6, 4), 0.254364013671875, [0.261021839148940, 0.254558441227157], [-0.0262197902254473, 0.000151011402221735]], [(12, 8), 0.181976318359375, [0.184570312500000, 0.182000000000000], [-0.0142545700073242, -0.00151000000000000]], [(24, 16), 0.129287719726562, [0.130510919574470, 0.129683383669613], [-0.00947434154988702, -0.00267741572252445]], [(48, 32), 0.0914535522460938, [0.0922851562500000, 0.0920000000000000], [-0.00909328460693359, -0.00592000000000000]]] |

92 | |

93 | Another smooth point example. |

94 | Turns out the terms 2--4 of the asymptotic expansion are 0. :: |

95 | |

96 | sage: R.<x,y> = PolynomialRing(QQ) |

97 | sage: q= 1/2 |

98 | sage: qq= q.denominator() |

99 | sage: G= 1-q*x |

100 | sage: H= 1-q*x +q*x*y -x^2*y |

101 | sage: F= QuasiRationalExpression(G,H) |

102 | sage: alpha= list(qq*vector([2,1-q])); alpha |

103 | [4, 1] |

104 | sage: I= F.smooth_critical(alpha); I |

105 | Ideal (y^2 - 2*y + 1, x + 1/4*y - 5/4) of Multivariate Polynomial Ring in x, y over Rational Field |

106 | sage: s= solve(I.gens(),F.variables(),solution_dict=true); s |

107 | [{y: 1, x: 1}] |

108 | sage: c= s[0] |

109 | sage: asy= F.asymptotics(c,alpha,5); asy |

110 | Initializing auxiliary functions... |

111 | Calculating derivatives of auxiallary functions... |

112 | Calculating derivatives of more auxiliary functions... |

113 | Calculating actions of the second order differential operator... |

114 | (1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*n^(1/3)) - 1/96*2^(1/3)*sqrt(3)*gamma(2/3)/(pi*n^(5/3)), 1, 1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*n^(1/3)) - 1/96*2^(1/3)*sqrt(3)*gamma(2/3)/(pi*n^(5/3))) |

115 | sage: F.relative_error(asy[0],alpha,[1,2],asy[1]) |

116 | Calculating errors table in the form |

117 | exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... |

118 | [[(4, 1), 0.187500000000000, [0.185581424449526], [0.0102324029358597]], [(8, 2), 0.152343750000000, [0.151986595969759], [0.00234439568568301]]] |

119 | |

120 | A multiple point example. :: |

121 | |

122 | sage: R.<x,y,z>= PolynomialRing(QQ) |

123 | sage: G= 1 |

124 | sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) |

125 | sage: F= QuasiRationalExpression(G,H) |

126 | sage: I= F.singular_points(); I |

127 | 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 |

128 | sage: c= {x: 1/2, y: 1, z: 4/3} |

129 | sage: F.is_cmp(c) |

130 | [({y: 1, z: 4/3, x: 1/2}, True, 'all good')] |

131 | sage: cc= F.critical_cone(c); cc |

132 | ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N) |

133 | sage: alpha= [8,3,3] |

134 | sage: alpha in cc[1] |

135 | True |

136 | sage: asy= F.asymptotics(c,alpha,2); asy |

137 | Initializing auxiliary functions... |

138 | Calculating derivatives of auxiliary functions... |

139 | Calculating derivatives of more auxiliary functions... |

140 | Calculating second-order differential operator actions... |

141 | (1/172872*(24696*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n)) - 1231*sqrt(3)*sqrt(7)/(sqrt(pi)*n^(3/2)))*108^n, 108, 1/7*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n)) - 1231/172872*sqrt(3)*sqrt(7)/(sqrt(pi)*n^(3/2))) |

142 | |

143 | Another multiple point example. :: |

144 | |

145 | sage: R.<x,y,z>= PolynomialRing(QQ) |

146 | sage: G= 16 |

147 | sage: H= (4-2*x-y-z)^2*(4-x-2*y-z) |

148 | sage: F= QuasiRationalExpression(G,H) |

149 | sage: var('a1,a2,a3') |

150 | (a1, a2, a3) |

151 | sage: F.cohomologous_integrand([a1,a2,a3]) |

152 | [[-16*(a2*z - 2*a3*y)*n/(y*z) + 16*(2*y - z)/(y*z), [x + 2*y + z - 4, 1], [2*x + y + z - 4, 1]]] |

153 | sage: I= F.singular_points(); I |

154 | Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field |

155 | sage: c= {x:1, y:1, z:1} |

156 | sage: F.is_cmp(c) |

157 | [({y: 1, z: 1, x: 1}, True, 'all good')] |

158 | sage: alpha= [3,3,2] |

159 | sage: cc= F.critical_cone(c); cc |

160 | ([(1, 2, 1), (2, 1, 1)], 2-d cone in 3-d lattice N) |

161 | sage: alpha in cc[1] |

162 | True |

163 | sage: asy= F.asymptotics(c,alpha,2); asy |

164 | Initializing auxiliary functions... |

165 | Calculating derivatives of auxiliary functions... |

166 | Calculating derivatives of more auxiliary functions... |

167 | Calculating second-order differential operator actions... |

168 | (4/3*sqrt(3)*sqrt(n)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(n)), 1, 4/3*sqrt(3)*sqrt(n)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(n))) |

169 | |

170 | """ |

171 | #***************************************************************************** |

172 | # Copyright (C) 2008 Alexander Raichev <tortoise.said@gmail.com> |

173 | # |

174 | # Distributed under the terms of the GNU General Public License (GPL) |

175 | # http://www.gnu.org/licenses/ |

176 | #***************************************************************************** |

177 | from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing |

178 | |

179 | # Functions to incorporate later into existing Sage classes. |

180 | def algebraic_dependence(fs): |

181 | r""" |

182 | This function returns an irreducible annihilating polynomial for the |

183 | polynomials in `fs`, if those polynomials are algebraically dependent. |

184 | Otherwise it returns 0. |

185 | |

186 | INPUT: |

187 | |

188 | - ``fs`` - A list of polynomials `f_1,\ldots,f_r' from a common polynomial |

189 | ring. |

190 | |

191 | OUTPUT: |

192 | |

193 | If the polynomials in `fs` are algebraically dependent, then the |

194 | output is an irreducible polynomial `g` in `K[T_1,\ldots,T_r]` such that |

195 | `g(f_1,\ldots,f_r) = 0`. |

196 | Here `K` is the coefficient ring of self and `T_1,\ldots,T_r` are |

197 | indeterminates different from those of self. |

198 | If the polynomials in `fs` are algebraically independent, then the output |

199 | is the zero polynomial. |

200 | |

201 | EXAMPLES:: |

202 | |

203 | sage: from sage.combinat.amgf import * |

204 | sage: R.<x>= PolynomialRing(QQ) |

205 | sage: fs= [x-3, x^2+1] |

206 | sage: g= algebraic_dependence(fs); g |

207 | 10 + 6*T0 - T1 + T0^2 |

208 | sage: g(fs) |

209 | 0 |

210 | |

211 | :: |

212 | |

213 | sage: R.<w,z>= PolynomialRing(QQ) |

214 | sage: fs= [w, (w^2 + z^2 - 1)^2, w*z - 2] |

215 | sage: g= algebraic_dependence(fs); g |

216 | 16 + 32*T2 - 8*T0^2 + 24*T2^2 - 8*T0^2*T2 + 8*T2^3 + 9*T0^4 - 2*T0^2*T2^2 + T2^4 - T0^4*T1 + 8*T0^4*T2 - 2*T0^6 + 2*T0^4*T2^2 + T0^8 |

217 | sage: g(fs) |

218 | 0 |

219 | sage: g.parent() |

220 | Multivariate Polynomial Ring in T0, T1, T2 over Rational Field |

221 | |

222 | :: |

223 | |

224 | sage: R.<x,y,z>= PolynomialRing(GF(7)) |

225 | sage: fs= [x*y,y*z] |

226 | sage: g= algebraic_dependence(fs); g |

227 | 0 |

228 | |

229 | AUTHORS: |

230 | |

231 | - Alex Raichev (2011-01-06) |

232 | """ |

233 | r= len(fs) |

234 | R= fs[Integer(0) ].parent() |

235 | B= R.base_ring() |

236 | Xs= list(R.gens()) |

237 | d= len(Xs) |

238 | |

239 | # Expand R by r new variables. |

240 | T= 'T' |

241 | while T in [str(x) for x in Xs]: |

242 | T= T+'T' |

243 | Ts= [T + str(j) for j in range(r)] |

244 | RR= PolynomialRing(B,d+r,tuple(Xs+Ts)) |

245 | Vs= list(RR.gens()) |

246 | Xs= Vs[Integer(0) :d] |

247 | Ts= Vs[d:] |

248 | |

249 | # Find an irreducible annihilating polynomial g for the fs if there is one. |

250 | J= ideal([ Ts[j] -RR(fs[j]) for j in range(r)]) |

251 | JJ= J.elimination_ideal(Xs) |

252 | g= JJ.gens()[Integer(0) ] |

253 | |

254 | # Shrink the ambient ring of g to include only Ts. |

255 | # I chose the negdeglex order simply because i find it useful in my work. |

256 | RRR= PolynomialRing(B,r,tuple(Ts),order='negdeglex') |

257 | return RRR(g) |

258 | #------------------------------------------------------------------------------- |

259 | def partial_fraction_decomposition(f): |

260 | r""" |

261 | Return a partial fraction decomposition of the rational expression `f`. |

262 | Works for univariate and mulitivariate `f`. |

263 | |

264 | INPUT: |

265 | |

266 | - ``f`` - An element of the field of fractions `F` of a polynomial ring |

267 | `R` whose coefficients ring is a field. |

268 | In the univariate case, the coefficient ring doesn't have to be a field. |

269 | |

270 | OUTPUT: |

271 | |

272 | A tuple `whole,parts` where `whole \in R` and `parts` is the list of |

273 | terms (in `F`) of a partial fraction decomposition of `f - whole`. |

274 | See the notes below for more details. |

275 | |

276 | EXAMPLES:: |

277 | |

278 | sage: from sage.combinat.amgf import * |

279 | sage: S.<t> = QQ[] |

280 | sage: q = 1/(t+1) + 2/(t+2) + 3/(t-3) |

281 | sage: whole, parts = partial_fraction_decomposition(q); parts |

282 | [3/(t - 3), 1/(t + 1), 2/(t + 2)] |

283 | sage: whole +sum(parts) == q |

284 | True |

285 | |

286 | Notice that there is a whole part below despite the appearance of q :: |

287 | |

288 | sage: R.<x,y>= PolynomialRing(QQ) |

289 | sage: q= 1/( x *y *(x*y-1) ) |

290 | sage: whole,parts= partial_fraction_decomposition(q) |

291 | sage: whole, parts |

292 | (-1, [x*y/(x*y - 1), (-1)/(x*y)]) |

293 | sage: q == whole +sum(parts) |

294 | True |

295 | |

296 | :: |

297 | |

298 | sage: R.<x,y>= PolynomialRing(QQ) |

299 | sage: q= x +1/x +1/(x*y-2) +1/(x^2+y^2-1) |

300 | sage: whole,parts= partial_fraction_decomposition(q) |

301 | sage: whole, parts |

302 | (x, [(1/2*x^3*y^2 + 1/2*x*y^4 + 1/2*x^3*y + 1/2*x^2*y^2 + 1/2*x*y^3 - x^2*y - 1/2*x*y^2 - y^3 - 3/2*x*y + y)/(x^3*y + x*y^3 - 2*x^2 - x*y - 2*y^2 + 2), (-1/2*x^3*y - 1/2*x*y^3 - 1/2*x^3 - 1/2*x^2*y - 1/2*x*y^2 + x^2 + 1/2*x*y + y^2 + 3/2*x - 1)/(x^3 + x*y^2 - x)]) |

303 | sage: whole +sum(parts)==q |

304 | True |

305 | |

306 | :: |

307 | |

308 | sage: R.<x,y,z>= PolynomialRing(QQ) |

309 | sage: q= 1/x +1/(x*y-z-2)^2 +1/(x^2+y^2 +z^2-1)^3 +1/(x*y-3) |

310 | sage: whole,parts= partial_fraction_decomposition(q) |

311 | sage: whole +sum(parts)==q |

312 | True |

313 | |

314 | NOTES: |

315 | |

316 | In the case of univariate `f` this function calls the old univariate |

317 | partial fraction decomposition function. |

318 | In the multivariate case, it uses a different notion of and algorithm for |

319 | partial fraction decompositions. |

320 | |

321 | Let `f= P/Q` where `P,Q \in R`, let `Q_1^{e_1} \cdots Q_m^{e_m}` be the |

322 | unique factorization of `Q` in `R` into irreducible factors, and let `d` be |

323 | the number of indeterminates of `R`. |

324 | Then `f` can be written as a sum `\sum P_A/\prod_{j\in A} Q_j^{b_j}`, |

325 | where the `b_j \le e_j` are positive integers, the `P_A` are in `R`, and |

326 | the sum is taken over all size `\le m` subsets `A \subseteq \{ 1, \ldots, d \}` |

327 | such that `S:= \{ Q_j : j\in A \}` is algebraically independent and the |

328 | ideal generated by `S` is not all of `R` (that is, by the Nullstellensatz, |

329 | the polynomials of `S` have no common roots in the algebraic closure of the |

330 | coefficient field of `R`). |

331 | Following [Lein1978]_, i call any such decomposition of `f` a |

332 | `\emph{partial fraction decomposition}`. |

333 | |

334 | The algorithm used below comes from Theorem 1, Lemma 2, and Lemma 3 of |

335 | [Lein1978]_. |

336 | By the way, that paper has a typo in equation (c) on the |

337 | second page and equation (b) on the third page. |

338 | The right sides of (c) and (b) should be multiplied by `P`. |

339 | |

340 | REFERENCES: |

341 | |

342 | .. [Lein1978] E. K. Leinartas, "On expansion of rational functions of |

343 | several variables into partial fractions", Soviet Math. (Iz. VUZ) |

344 | 22 (1978), no. 10, 35--38. |

345 | |

346 | AUTHORS: |

347 | |

348 | - Alex Raichev (2011-01-10) |

349 | """ |

350 | R= f.denominator().parent() |

351 | d= len(R.gens()) |

352 | if d==Integer(1) : |

353 | return f.partial_fraction_decomposition() |

354 | Q= f.denominator() |

355 | whole,P= f.numerator().quo_rem(Q) |

356 | parts= [format_quotient(Integer(1) /Q)] |

357 | # Decompose via nullstellensatz trick |

358 | # (which is faster than the algebraic dependence trick) |

359 | parts= decompose_via_nullstellensatz(parts) |

360 | # Decompose via algebraic dependence trick |

361 | parts= decompose_via_algebraic_dependence(parts) |

362 | # Rewrite parts back in terms of rational expressions |

363 | new_parts=[] |

364 | for p in parts: |

365 | f= unformat_quotient(p) |

366 | if f.denominator() == Integer(1) : |

367 | whole= whole +f |

368 | else: |

369 | new_parts.append(f) |

370 | # Put P back in |

371 | new_parts= [P*f for f in new_parts] |

372 | return whole, new_parts |

373 | #------------------------------------------------------------------------------- |

374 | def combine_like_terms(parts,rationomial=True): |

375 | r""" |

376 | Combines like terms of the fractions represented by parts. |

377 | For use by partial_fraction_decomposition() above. |

378 | |

379 | EXAMPLES: |

380 | |

381 | sage: from sage.combinat.amgf import * |

382 | sage: R.<x,y>= PolynomialRing(QQ) |

383 | sage: parts =[[1, [x*y, 1]], [x, [x*y, 1]]] |

384 | sage: combine_like_terms(parts) |

385 | [[x + 1, [y, 1], [x, 1]]] |

386 | |

387 | :: |

388 | |

389 | sage: R.<x>= PolynomialRing(QQ) |

390 | sage: parts =[[1, [x, 1]], [x-1, [x, 1]]] |

391 | sage: combine_like_terms(parts) |

392 | [[1]] |

393 | |

394 | AUTHORS: |

395 | |

396 | - Alex Raichev (2011-01-10) |

397 | """ |

398 | |

399 | if parts == []: return parts |

400 | # Sort parts by denominators |

401 | new_parts= sorted(parts, key= lambda p:p[Integer(1) :]) |

402 | # Combine items of parts with same denominators. |

403 | newnew_parts=[] |

404 | left,right=Integer(0) ,Integer(1) |

405 | glom= new_parts[left][Integer(0) ] |

406 | while right <= len(new_parts)-Integer(1) : |

407 | if new_parts[left][Integer(1) :] == new_parts[right][Integer(1) :]: |

408 | glom= glom +new_parts[right][Integer(0) ] |

409 | else: |

410 | newnew_parts.append([glom]+new_parts[left][Integer(1) :]) |

411 | left= right |

412 | glom= new_parts[left][Integer(0) ] |

413 | right= right +Integer(1) |

414 | if glom != Integer(0) : |

415 | newnew_parts.append([glom]+new_parts[left][Integer(1) :]) |

416 | if rationomial: |

417 | # Reduce fractions in newnew_parts. |

418 | # Todo: speed up below by working directly with newnew_parts and |

419 | # thereby make fewer calls to format_quotient() which in turn |

420 | # calls factor(). |

421 | newnew_parts= [format_quotient(unformat_quotient(part)) for part in newnew_parts] |

422 | return newnew_parts |

423 | #------------------------------------------------------------------------------- |

424 | def decompose_via_algebraic_dependence(parts): |

425 | r""" |

426 | Returns a decomposition of parts. |

427 | Used by partial_fraction_decomposition() above. |

428 | Implements Lemma 2 of [Lein1978]_. |

429 | Recursive. |

430 | |

431 | REFERENCES: |

432 | |

433 | .. [Lein1978] E. K. Leinartas, "On expansion of rational functions of |

434 | several variables into partial fractions", Soviet Math. (Iz. VUZ) |

435 | 22 (1978), no. 10, 35--38. |

436 | |

437 | AUTHORS: |

438 | |

439 | - Alex Raichev (2011-01-10) |

440 | """ |

441 | decomposing_done= True |

442 | new_parts= [] |

443 | for p in parts: |

444 | p_parts= [p] |

445 | P= p[Integer(0) ] |

446 | Qs= p[Integer(1) :] |

447 | m= len(Qs) |

448 | G= algebraic_dependence([q for q,e in Qs]) |

449 | if G: |

450 | # Then the denominator factors are algebraically dependent |

451 | # and so we can decompose p. |

452 | decomposing_done= False |

453 | P= p[Integer(0) ] |

454 | Qs= p[Integer(1) :] |

455 | |

456 | # Todo: speed up step below by using |

457 | # G to calculate F. See [Lein1978]_ Lemma 1. |

458 | F= algebraic_dependence([q**e for q,e in Qs]) |

459 | new_vars= F.parent().gens() |

460 | |

461 | # Note that new_vars[j] corresponds to Qs[j] so that |

462 | # F([q^e for q,e in Qs]) = 0. |

463 | # Assuming below that F.parent() has negdeglex term order |

464 | # so that F.lt() is indeed the monomial we want. |

465 | FF= (F.lt() -F)/(F.lc()) |

466 | numers= map(mul,zip(FF.coefficients(),FF.monomials())) |

467 | e= list(F.lt().exponents())[Integer(0) :m] |

468 | denom= [[new_vars[j], e[Integer(0) ][j]+Integer(1) ] for j in range(m)] |

469 | |

470 | # Before making things messy by substituting in Qs, |

471 | # reduce terms and combine like terms. |

472 | p_parts_temp= [format_quotient(unformat_quotient([a]+denom)) for a in numers] |

473 | p_parts_temp= combine_like_terms(p_parts_temp) |

474 | |

475 | # Substitute Qs into new_p. |

476 | Qpowsub= dict([(new_vars[j],Qs[j][Integer(0) ]**Qs[j][Integer(1) ]) for j in range(m)]) |

477 | p_parts=[] |

478 | for x in p_parts_temp: |

479 | y= P*F.parent()(x[Integer(0) ]).subs(Qpowsub) |

480 | yy=[] |

481 | for xx in x[Integer(1) :]: |

482 | if xx[Integer(0) ] in new_vars: |

483 | j= new_vars.index(xx[Integer(0) ]) |

484 | yy.append([Qs[j][Integer(0) ],Qs[j][Integer(1) ]*xx[Integer(1) ]]) |

485 | else: |

486 | # Occasionally xx[0] is an integer. |

487 | yy.append(xx) |

488 | p_parts.append([y]+yy) |

489 | # Done one step of decomposing p. Add it to new_parts. |

490 | new_parts.extend(p_parts) |

491 | if decomposing_done: |

492 | return new_parts |

493 | else: |

494 | return decompose_via_algebraic_dependence(new_parts) |

495 | #------------------------------------------------------------------------------- |

496 | def decompose_via_cohomology(parts): |

497 | r""" |

498 | Given each nice (described below) differential form |

499 | `(P/Q) dx_1 \wedge\cdots\wedge dx_d` in `parts`, |

500 | this function returns a differential form equivalent in De Rham |

501 | cohomology that has no repeated factors in the denominator. |

502 | |

503 | INPUT: |

504 | |

505 | - ``parts`` - A list of the form `[chunk_1,\ldots,chunk_r]`, where each |

506 | `chunk_j` has the form `[P,[Q_1,e_1],\ldots,[Q_m,e_m]]`, |

507 | `Q_1,\ldots,Q_m` are irreducible elements of a common polynomial |

508 | ring `R` such that their corresponding algebraic varieties |

509 | `\{x\in F^d : B_j(x)=0\}` intersect transversely (where `F` is the |

510 | algebraic closure of the field of coefficients of `R`), |

511 | `e_1,\ldots,e_m` are positive integers, `m \le d`, and |

512 | `P` is a symbolic expression in some of the indeterminates of `R`. |

513 | Here `[P,[Q_1,e_1],\ldots,[Q_m,e_m]]` represents the fraction |

514 | `P/(Q_1^e_1 \cdots Q_m^e_m)`. |

515 | |

516 | OUTPUT: |

517 | |

518 | A list of the form `[chunky_1,\ldots,chunky_s]`, where each |

519 | `chunky_j` has the form `[P,[Q_1,1],\ldots,[Q_m,1]]`. |

520 | |

521 | EXAMPLES:: |

522 | |

523 | sage: from sage.combinat.amgf import * |

524 | sage: R.<x,y>= PolynomialRing(QQ) |

525 | sage: decompose_via_cohomology([[ 1, [x,3] ]]) |

526 | [] |

527 | |

528 | :: |

529 | |

530 | sage: R.<x,y>= PolynomialRing(QQ) |

531 | sage: decompose_via_cohomology([[ 1, [x,3], [x*y-1,2] ]]) |

532 | [[-3*y^2, [x, 1], [x*y - 1, 1]], [-5*y^3, [x*y - 1, 1]]] |

533 | |

534 | NOTES: |

535 | |

536 | This is a recursive function thats stops calling itself when all the |

537 | `e_j` equal 1 or `parts == []`. |

538 | The algorithm used here is that of Theorem 17.4 of |

539 | [AiYu1983]_. |

540 | The algebraic varieties `\{x\in F^d : Q_j(x)=0\}` |

541 | (where `F` is the algebraic closure of the field of coefficients of `R`) |

542 | corresponding to the `Q_j` __intersect transversely__ iff for each |

543 | point `c` of their intersection and for all `k \le m`, |

544 | the Jacobian matrix of any `k` polynomials from |

545 | `\{Q_1,\ldots,Q_m\}` has rank equal to `\min\{k,d\}` when evaluated at |

546 | `c`. |

547 | |

548 | REFERENCES: |

549 | |

550 | .. [AiYu1983] I. A. A\u\izenberg and A. P. Yuzhakov, "Integral |

551 | representations and residues in multidimensional complex analysis", |

552 | Translations of Mathematical Monographs, 58. American Mathematical |

553 | Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X. |

554 | |

555 | AUTHORS: |

556 | |

557 | - Alex Raichev (2008-10-01, 2011-01-15) |

558 | """ |

559 | if parts == []: return parts |

560 | import copy # Will need this to make copies of a nested list. |

561 | decomposing_done= True |

562 | new_parts= [] |

563 | R= parts[Integer(0) ][Integer(1) ][Integer(0) ].parent() |

564 | V= list(R.gens()) |

565 | for p in parts: |

566 | p_parts= [p] |

567 | P= p[Integer(0) ] |

568 | Qs= p[Integer(1) :] |

569 | m= len(Qs) |

570 | if sum([e for q,e in Qs]) > m: |

571 | # Then we can decompose p |

572 | p_parts= [] |

573 | decomposing_done= False |

574 | dets= [] |

575 | vars_subsets= Set(V).subsets(m) |

576 | for v in vars_subsets: |

577 | # Sort variables so that first polynomial ring indeterminate |

578 | # declared is first in vars list. |

579 | v= sorted(v,reverse=true) |

580 | jac= jacobian([q for q,e in Qs],v) |

581 | dets.append(jac.determinant()) |

582 | # Get a Nullstellensatz certificate. |

583 | L= R(Integer(1) ).lift(R.ideal([q for q,e in Qs] +dets)) |

584 | for j in range(m): |

585 | if L[j] != Integer(0) : |

586 | # Make a copy of (and not a reference to) the nested list Qs. |

587 | new_Qs = copy.deepcopy(Qs) |

588 | if new_Qs[j][Integer(1) ] > Integer(1) : |

589 | new_Qs[j][Integer(1) ]= new_Qs[j][Integer(1) ] -Integer(1) |

590 | else: |

591 | del new_Qs[j] |

592 | p_parts.append( [P*L[j]] +new_Qs ) |

593 | for k in range(vars_subsets.cardinality()): |

594 | if L[m+k] != Integer(0) : |

595 | new_Qs = copy.deepcopy(Qs) |

596 | for j in range(m): |

597 | if new_Qs[j][Integer(1) ] > Integer(1) : |

598 | new_Qs[j][Integer(1) ]= new_Qs[j][Integer(1) ] -Integer(1) |

599 | v= sorted(vars_subsets[k],reverse=true) |

600 | jac= jacobian([SR(P*L[m+k])] +[ SR(Qs[jj][Integer(0) ]) for \ |

601 | jj in range(m) if jj !=j], [SR(vv) for vv in v]) |

602 | det= jac.determinant() |

603 | if det != Integer(0) : |

604 | p_parts.append([permutation_sign(v,V) \ |

605 | *(-Integer(1) )**j/new_Qs[j][Integer(1) ] *det] +new_Qs) |

606 | break |

607 | # Done one step of decomposing p. Add it to new_parts. |

608 | new_parts.extend(p_parts) |

609 | new_parts= combine_like_terms(new_parts,rationomial=False) |

610 | if decomposing_done: |

611 | return new_parts |

612 | else: |

613 | return decompose_via_cohomology(new_parts) |

614 | #------------------------------------------------------------------------------- |

615 | def decompose_via_nullstellensatz(parts): |

616 | r""" |

617 | Returns a decomposition of parts. |

618 | Used by partial_fraction_decomposition() above. |

619 | Implements Lemma 3 of [Lein1978]_. |

620 | Recursive. |

621 | |

622 | REFERENCES: |

623 | |

624 | .. [Lein1978] E. K. Leinartas, "On expansion of rational functions of |

625 | several variables into partial fractions", Soviet Math. (Iz. VUZ) |

626 | 22 (1978), no. 10, 35--38. |

627 | |

628 | AUTHORS: |

629 | |

630 | - Alex Raichev (2011-01-10) |

631 | """ |

632 | decomposing_done= True |

633 | new_parts= [] |

634 | R= parts[Integer(0) ][Integer(0) ].parent() |

635 | for p in parts: |

636 | p_parts= [p] |

637 | P= p[Integer(0) ] |

638 | Qs= p[Integer(1) :] |

639 | m= len(Qs) |

640 | if R(Integer(1) ) in ideal([q for q,e in Qs]): |

641 | # Then we can decompose p. |

642 | decomposing_done= False |

643 | L= R(Integer(1) ).lift(R.ideal([q**e for q,e in Qs])) |

644 | p_parts= [ [P*L[i]] + \ |

645 | [[Qs[j][Integer(0) ],Qs[j][Integer(1) ]] for j in range(m) if j != i] \ |

646 | for i in range(m) if L[i]!=Integer(0) ] |

647 | # The procedure above yields no like terms to combine. |

648 | # Done one step of decomposing p. Add it to new_parts. |

649 | new_parts.extend(p_parts) |

650 | if decomposing_done: |

651 | return new_parts |

652 | else: |

653 | return decompose_via_nullstellensatz(new_parts) |

654 | #------------------------------------------------------------------------------- |

655 | def format_quotient(f): |

656 | r""" |

657 | Formats `f` for use by partial_fraction_decomposition() above. |

658 | |

659 | AUTHORS: |

660 | |

661 | - Alex Raichev (2011-01-10) |

662 | """ |

663 | P= f.numerator() |

664 | Q= f.denominator() |

665 | Qs= Q.factor() |

666 | if Qs.unit() != Integer(1) : |

667 | P= P/Qs.unit() |

668 | Qs= sorted([[q,e] for q,e in Qs]) # sorting for future bookkeeping |

669 | return [P]+Qs |

670 | #------------------------------------------------------------------------------- |

671 | def permutation_sign(v,vars): |

672 | r""" |

673 | This function returns the sign of the permutation on `1,\ldots,len(vars)` |

674 | that is induced by the sublist `v` of `vars`. |

675 | For internal use by _cohom_equiv_main(). |

676 | |

677 | INPUT: |

678 | |

679 | - ``v`` - A sublist of `vars`. |

680 | - ``vars`` - A list of predefined variables or numbers. |

681 | |

682 | OUTPUT: |

683 | |

684 | The sign of the permutation obtained by taking indices (and adding 1) |

685 | within `vars` of the list `v,w`, where `w` is the list `vars` with the |

686 | elements of `v` removed. |

687 | |

688 | EXAMPLES:: |

689 | |

690 | sage: from sage.combinat.amgf import * |

691 | sage: vars= ['a','b','c','d','e'] |

692 | sage: v= ['b','d'] |

693 | sage: permutation_sign(v,vars) |

694 | -1 |

695 | sage: v= ['d','b'] |

696 | sage: permutation_sign(v,vars) |

697 | 1 |

698 | |

699 | AUTHORS: |

700 | |

701 | - Alex Raichev (2008-10-01) |

702 | """ |

703 | # Convert variable lists to lists of numbers in {1,...,len(vars)} |

704 | A= [x+Integer(1) for x in range(len(vars))] |

705 | B= [vars.index(x)+Integer(1) for x in v] |

706 | C= list(Set(A).difference(Set(B))) |

707 | C.sort() |

708 | P= Permutation(B+C) |

709 | return P.signature() |

710 | #------------------------------------------------------------------------------- |

711 | def unformat_quotient(part): |

712 | r""" |

713 | Unformats `f` for use by partial_fraction_decomposition() above. |

714 | Inverse of format_quotient() above. |

715 | |

716 | AUTHORS: |

717 | |

718 | - Alex Raichev (2011-01-10) |

719 | """ |

720 | P= part[Integer(0) ] |

721 | Qs= part[Integer(1) :] |

722 | Q= prod([q**e for q,e in Qs]) |

723 | return P/Q |

724 | #=============================================================================== |

725 | # Class for calculation of asymptotics of multivariate generating functions. |

726 | class QuasiRationalExpression(object): |

727 | "Store an expression G/H, where H comes from a polynomial ring R and \ |

728 | G comes from R or the Symbolic Ring." |

729 | def __init__(self, G, H): |

730 | # Store important information about object as attributes of self. |

731 | # G, H, H's ring, ring dimension, H's factorization. |

732 | self._G = G |

733 | self._H = H |

734 | R= H.parent() |

735 | self._R = R |

736 | self._d = len(R.gens()) |

737 | self._Hfac= list(H.factor()) |

738 | |

739 | # Variables of self as elements of the SR. |

740 | # Remember that G might be in SR and not in R. |

741 | try: |

742 | # This fails if G is a constant, for example. |

743 | Gv= Set([R(x) for x in G.variables()]) |

744 | except: |

745 | Gv= Set([]) |

746 | try: |

747 | # This fails if H is a constant, for example. |

748 | Hv= Set(H.variables()) |

749 | except: |

750 | Hv= Set([]) |

751 | # Preserve the ring ordering of the variables which some methods below |

752 | # depends upon. |

753 | V= sorted(list(Gv.union(Hv)),reverse=True) |

754 | self._variables = tuple([SR(x) for x in V]) |

755 | #------------------------------------------------------------------------------- |

756 | # Keeping methods in alphabetical order (ignoring initial single underscores) |

757 | def asymptotics(self,c,alpha,N,numerical=Integer(0) ,asy_var=None): |

758 | r""" |

759 | This function returns the first `N` terms of the asymptotic expansion |

760 | of the Maclaurin coefficients `F_{n\alpha}` of the |

761 | multivariate meromorphic function `F=G/H` as `n\to\infty`, |

762 | where `F = self`. |

763 | It assumes that `F` is holomorphic in a neighborhood of the origin, |

764 | that `H` is a polynomial, and that `c` is a convenient strictly minimal |

765 | smooth or multiple of `F` that is critical for `\alpha`. |

766 | |

767 | INPUT: |

768 | |

769 | - ``alpha`` - A `d`-tuple of positive integers or, if `c` is a smooth |

770 | point, possibly of symbolic variables. |

771 | - ``c`` - A dictionary with keys `self._variables` and values from a |

772 | superfield of the field of `self._R.base_ring()`. |

773 | - ``N`` - A positive integer. |

774 | - ``numerical`` - A natural number (default: 0). |

775 | If k=numerical > 0, then a numerical approximation of the coefficients |

776 | of `F_{n\alpha}` with k digits of precision will be returned. |

777 | Otherwise exact values will be returned. |

778 | - ``asy_var`` - A symbolic variable (default: None). |

779 | The variable of the asymptotic expansion. |

780 | If none is given, `var('n')` will be assigned. |

781 | |

782 | OUTPUT: |

783 | |

784 | The tuple `(asy,exp_part,subexp_part)`, where `asy` is first `N` terms |

785 | of the asymptotic expansion of the Maclaurin coefficients `F_{n\alpha}` |

786 | of the function `F=self` as `n\to\infty`, `exp_part^n` is the exponential |

787 | factor of `asy`, and `subexp_part` is the subexponential factor of |

788 | `asy`. |

789 | |

790 | EXAMPLES:: |

791 | |

792 | A smooth point example :: |

793 | |

794 | sage: from sage.combinat.amgf import * |

795 | sage: R.<x,y>= PolynomialRing(QQ) |

796 | sage: G= 1 |

797 | sage: H= 1-x-y-x*y |

798 | sage: F= QuasiRationalExpression(G,H) |

799 | sage: alpha= [3,2] |

800 | sage: c= {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3} |

801 | sage: F.asymptotics(c,alpha,2,numerical=3) |

802 | Initializing auxiliary functions... |

803 | Calculating derivatives of auxiallary functions... |

804 | Calculating derivatives of more auxiliary functions... |

805 | Calculating actions of the second order differential operator... |

806 | ((0.369/sqrt(n) - 0.0186/n^(3/2))*71.2^n, 71.2, 0.369/sqrt(n) - 0.0186/n^(3/2)) |

807 | |

808 | A multiple point example :: |

809 | |

810 | sage: R.<x,y,z>= PolynomialRing(QQ) |

811 | sage: G= 1 |

812 | sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) |

813 | sage: F= QuasiRationalExpression(G,H) |

814 | sage: c= {z: 4/3, y: 1, x: 1/2} |

815 | sage: alpha= [8,3,3] |

816 | sage: F.asymptotics(c,alpha,1) |

817 | Initializing auxiliary functions... |

818 | Calculating derivatives of auxiliary functions... |

819 | Calculating derivatives of more auxiliary functions... |

820 | Calculating second-order differential operator actions... |

821 | (1/7*sqrt(3)*sqrt(7)*108^n/(sqrt(pi)*sqrt(n)), 108, 1/7*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n))) |

822 | |

823 | NOTES: |

824 | |

825 | A zero `c` of `H` is __strictly minimal__ if there is no zero `x` of `H` |

826 | such that `x_j < c_j` for all `0 \le j < d`. |

827 | For definitions of the terms "smooth critical point for `\alpha`" and |

828 | "multiple critical point for `\alpha`", |

829 | see the documentation for _asymptotics_main_smooth() and |

830 | _asymptotics_main_multiple(), which are the functions that do most of the |

831 | work. |

832 | |

833 | ALGORITHM: |

834 | |

835 | The algorithm used here comes from [RaWi2011]_. |

836 | |

837 | (1) Use Cauchy's multivariate integral formula to write `F_{n\alpha}` as |

838 | an integral around a polycirle centered at the origin of the |

839 | differential form `\frac{G(x) dx[0] \wedge\cdots\wedge |

840 | dx[d-1]}{H(x)x^\alpha}`. |

841 | |

842 | (2) Decompose `G/H` into a sum of partial fractions `P[0] +\cdots+ P[r]` |

843 | so that each term of the sum has at most `d` irreducible factors of `H` |

844 | in the denominator. |

845 | |

846 | (3) For each differential form `P[j] dx[0] \wedge\cdots\wedge dx[d-1]`, |

847 | find an equivalent form `\omega[j]` in de Rham cohomology with no |

848 | repeated irreducible factors of `H` in its denominator. |

849 | |

850 | (4) Compute an asymptotic expansion for each integral `\omega[j]`. |

851 | |

852 | (5) Add the expansions. |

853 | |

854 | REFERENCES: |

855 | |

856 | .. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

857 | coefficients of multivariate generating functions: improvements |

858 | for smooth points", Electronic Journal of Combinatorics, Vol. 15 |

859 | (2008), R89. |

860 | |

861 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

862 | coefficients of multivariate generating functions: improvements |

863 | for smooth points", To appear. |

864 | |

865 | AUTHORS: |

866 | |

867 | - Alex Raichev (2008-10-01, 2010-09-28, 2011-04-27) |

868 | """ |

869 | # The variable for asymptotic expansions. |

870 | if not asy_var: |

871 | asy_var= var('n') |

872 | |

873 | # Create symbolic (non-ring) variables. |

874 | R= self._R |

875 | d= self._d |

876 | X= list(self._variables) |

877 | |

878 | # Do steps (1)--(3) |

879 | new_integrands= self.cohomologous_integrand(alpha,asy_var) |

880 | |

881 | # Coerce everything into the Symbolic Ring, as polynomial ring |

882 | # calculations are no longer needed. |

883 | # Calculate asymptotics. |

884 | cc={} |

885 | for k in c.keys(): |

886 | cc[SR(k)] = SR(c[k]) |

887 | c= cc |

888 | for i in range(len(alpha)): |

889 | alpha[i] = SR(alpha[i]) |

890 | subexp_parts= [] |

891 | for chunk in new_integrands: |

892 | # Convert chunk into Symbolic Ring |

893 | GG= SR(chunk[Integer(0) ]) |

894 | HH= [SR(f) for (f,e) in chunk[Integer(1) :]] |

895 | asy= self._asymptotics_main(GG,HH,X,c,asy_var,alpha,N,numerical) |

896 | subexp_parts.append(asy[Integer(2) ]) |

897 | exp_scale= asy[Integer(1) ] # same for all chunk in new_integrands |

898 | |

899 | # Do step (5). |

900 | subexp_part= add(subexp_parts) |

901 | return exp_scale**asy_var *subexp_part, exp_scale, subexp_part |

902 | #-------------------------------------------------------------------------------- |

903 | def _asymptotics_main(self,G,H,X,c,n,alpha,N,numerical): |

904 | r""" |

905 | This function is for internal use by asymptotics(). |

906 | It finds a variable in `X` to use to calculate asymptotics and decides |

907 | whether to call _asymptotics_main_smooth() or |

908 | _asymptotics_main_multiple(). |

909 | |

910 | Does not use `self`. |

911 | |

912 | INPUT: |

913 | |

914 | - ``G`` - A symbolic expression. |

915 | - ``H`` - A list of symbolic expressions. |

916 | - ``X`` - The list of variables occurring in `G` and `H`. |

917 | Call its length `d`. |

918 | - ``c`` - A dictionary with `X` as keys and numbers as values. |

919 | - ``n`` - The variable of the asymptotic expansion. |

920 | - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic |

921 | variables if `c` is a smooth point. |

922 | - ``N`` - A positive integer. |

923 | - ``numerical`` - Natural number. |

924 | If k=numerical > 0, then a numerical approximation of the coefficients |

925 | of `F_{n\alpha}` with k digits of precision will be returned. |

926 | Otherwise exact values will be returned. |

927 | |

928 | OUTPUT: |

929 | |

930 | The same as the function asymptotics(). |

931 | |

932 | AUTHORS: |

933 | |

934 | - Alex Raichev (2008-10-01, 2010-09-28) |

935 | """ |

936 | d= len(X) |

937 | r= len(H) # We know 1 <= r <= d. |

938 | |

939 | # Find a good variable x in X to do asymptotics calculations with, that is, |

940 | # a variable x in X such that for all h in H, diff(h,x).subs(c) != 0. |

941 | # A good variable is guaranteed to exist since we are dealing with |

942 | # convenient smooth or multiple points. |

943 | # Search for good x in X from back to front (to be consistent with |

944 | # [RaWi2008a]_ [RaWi2011]_ which uses X[d-1] as a good variable). |

945 | # Put each not good x found at the beginning of X and reshuffle alpha |

946 | # in the same way. |

947 | x= X[d-Integer(1) ] |

948 | beta= copy(alpha) |

949 | while Integer(0) in [diff(h,x).subs(c) for h in H]: |

950 | X.pop() |

951 | X.insert(Integer(0) ,x) |

952 | x= X[d-Integer(1) ] |

953 | a= beta.pop() |

954 | beta.insert(Integer(0) ,a) |

955 | if d==r: |

956 | # This is the case of a 'simple' multiple point. |

957 | A= G.subs(c) / jacobian(H,X).subs(c).determinant().abs() |

958 | return A,Integer(1) ,A |

959 | elif r==Integer(1) : |

960 | # So 1 = r < d, and we have a smooth point. |

961 | return self._asymptotics_main_smooth(G,H[Integer(0) ],X,c,n,beta,N,numerical) |

962 | else: |

963 | # So 1 < r < d, and we have a non-smooth multiple point. |

964 | return self._asymptotics_main_multiple(G,H,X,c,n,beta,N,numerical) |

965 | #-------------------------------------------------------------------------------- |

966 | def _asymptotics_main_multiple(self,G,H,X,c,n,alpha,N,numerical): |

967 | r""" |

968 | This function is for internal use by _asymptotics_main(). |

969 | It calculates asymptotics in case they depend upon multiple points. |

970 | |

971 | Does not use `self`. |

972 | |

973 | INPUT: |

974 | |

975 | - ``G`` - A symbolic expression. |

976 | - ``H`` - A list of symbolic expressions. |

977 | - ``X`` - The list of variables occurring in `G` and `H`. |

978 | Call its length `d`. |

979 | - ``c`` - A dictionary with `X` as keys and numbers as values. |

980 | - ``n`` - The variable of the asymptotic expansion. |

981 | - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic |

982 | variables if `c` is a smooth point. |

983 | - ``N`` - A positive integer. |

984 | - ``numerical`` - Natural number. |

985 | If k=numerical > 0, then a numerical approximation of the coefficients |

986 | of `F_{n\alpha}` with k digits of precision will be returned. |

987 | Otherwise exact values will be returned. |

988 | |

989 | OUTPUT: |

990 | |

991 | The same as the function asymptotics(). |

992 | |

993 | NOTES: |

994 | |

995 | The formula used for computing the asymptotic expansion is |

996 | Theorem 3.4 of [RaWi2011]_. |

997 | |

998 | Currently this function cannot handle `c` with symbolic variable keys, |

999 | because _crit_cone_combo() crashes. |

1000 | |

1001 | REFERENCES: |

1002 | |

1003 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

1004 | coefficients of multivariate generating functions: improvements |

1005 | for smooth points", To appear. |

1006 | |

1007 | AUTHORS: |

1008 | |

1009 | - Alex Raichev (2008-10-01, 2010-09-28) |

1010 | """ |

1011 | I= sqrt(-Integer(1) ) |

1012 | d= len(X) # > r > 1 |

1013 | r= len(H) # > 1 |

1014 | C= copy(c) |

1015 | |

1016 | S= [var(self._new_var_name('s',X) + str(j)) for j in range(r-Integer(1) )] |

1017 | T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-Integer(1) )] |

1018 | Sstar= {} |

1019 | temp= self._crit_cone_combo(H,X,c,alpha) |

1020 | for j in range(r-Integer(1) ): |

1021 | Sstar[S[j]]= temp[j] |

1022 | thetastar= dict([(t,Integer(0) ) for t in T]) |

1023 | thetastar.update(Sstar) |

1024 | |

1025 | # Setup. |

1026 | print "Initializing auxiliary functions..." |

1027 | Hmul= mul(H) |

1028 | h= [function('h'+str(j),*tuple(X[:d-Integer(1) ])) for j in range(r)] # Implicit functions |

1029 | U = function('U',*tuple(X)) |

1030 | # All other functions are defined in terms of h, U, and explicit functions. |

1031 | Hcheck= mul([ X[d-Integer(1) ] -Integer(1) /h[j] for j in range(r)] ) |

1032 | Gcheck= -G/U *mul( [-h[j]/X[d-Integer(1) ] for j in range(r)] ) |

1033 | A= [(-Integer(1) )**(r-Integer(1) ) *X[d-Integer(1) ]**(-r+j)*diff(Gcheck.subs({X[d-Integer(1) ]:Integer(1) /X[d-Integer(1) ]}),X[d-Integer(1) ],j) for j in range(r)] |

1034 | e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-Integer(1) )]) |

1035 | ht= [hh.subs(e) for hh in h] |

1036 | Ht= [H[j].subs(e).subs({X[d-Integer(1) ]:Integer(1) /ht[j]}) for j in range(r)] |

1037 | hsumt= add([S[j]*ht[j] for j in range(r-Integer(1) )]) +(Integer(1) -add(S))*ht[r-Integer(1) ] |

1038 | At= [AA.subs(e).subs({X[d-Integer(1) ]:hsumt}) for AA in A] |

1039 | Phit = -log(C[X[d-Integer(1) ]] *hsumt)+ I*add([alpha[i]/alpha[d-Integer(1) ] *T[i] for i in range(d-Integer(1) )]) |

1040 | # atC Stores h and U and all their derivatives evaluated at C. |

1041 | atC = copy(C) |

1042 | atC.update(dict( [(hh.subs(C),Integer(1) /C[X[d-Integer(1) ]]) for hh in h ])) |

1043 | |

1044 | # Compute the derivatives of h up to order 2*N and evaluate at C. |

1045 | hderivs1= {} # First derivatives of h. |

1046 | for (i,j) in mrange([d-Integer(1) ,r]): |

1047 | s= solve( diff(H[j].subs({X[d-Integer(1) ]:Integer(1) /h[j]}),X[i]), diff(h[j],X[i]) )[Integer(0) ].rhs()\ |

1048 | .simplify() |

1049 | hderivs1.update({diff(h[j],X[i]):s}) |

1050 | atC.update({diff(h[j],X[i]).subs(C):s.subs(C).subs(atC)}) |

1051 | hderivs = self._diff_all(h,X[Integer(0) :d-Integer(1) ],Integer(2) *N,sub=hderivs1,rekey=h) |

1052 | for k in hderivs.keys(): |

1053 | atC.update({k.subs(C):hderivs[k].subs(atC)}) |

1054 | |

1055 | # Compute the derivatives of U up to order 2*N-2+ min{r,N}-1 and evaluate at C. |

1056 | # To do this, differentiate H = U*Hcheck over and over, evaluate at C, |

1057 | # and solve for the derivatives of U at C. |

1058 | # Need the derivatives of H with short keys to pass on to diff_prod later. |

1059 | m= min(r,N) |

1060 | end= [X[d-Integer(1) ] for j in range(r)] |

1061 | Hmulderivs= self._diff_all(Hmul,X,Integer(2) *N-Integer(2) +r,ending=end,sub_final=C) |

1062 | print "Calculating derivatives of auxiliary functions..." |

1063 | atC.update({U.subs(C):diff(Hmul,X[d-Integer(1) ],r).subs(C)/factorial(r)}) |

1064 | Uderivs={} |

1065 | p= Hmul.polynomial(CC).degree()-r |

1066 | if p == Integer(0) : |

1067 | # Then, using a proposition at the end of [RaWi2011], we can |

1068 | # conclude that all higher derivatives of U are zero. |

1069 | for l in (ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )): |

1070 | for s in UnorderedTuples(X,l): |

1071 | Uderivs[diff(U,s).subs(C)] = Integer(0) |

1072 | elif p > Integer(0) and p < Integer(2) *N-Integer(2) +m-Integer(1) : |

1073 | all_zero= True |

1074 | Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,p)),end,Uderivs,atC) |

1075 | # Check for a nonzero U derivative. |

1076 | if Uderivs.values() != [Integer(0) for i in range(len(Uderivs))]: |

1077 | all_zero= False |

1078 | if all_zero: |

1079 | # Then, using a proposition at the end of [RaWi2011], we can |

1080 | # conclude that all higher derivatives of U are zero. |

1081 | for l in (ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )): |

1082 | for s in UnorderedTuples(X,l): |

1083 | Uderivs.update({diff(U,s).subs(C):Integer(0) }) |

1084 | else: |

1085 | # Have to compute the rest of the derivatives. |

1086 | Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )),end,Uderivs,atC) |

1087 | else: |

1088 | Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )),end,Uderivs,atC) |

1089 | atC.update(Uderivs) |

1090 | Phit1= jacobian(Phit,T+S).subs(hderivs1) |

1091 | a= jacobian(Phit1,T+S).subs(hderivs1).subs(thetastar).subs(atC) |

1092 | a_inv= a.inverse() |

1093 | Phitu= Phit -(Integer(1) /Integer(2) ) *matrix([T+S]) *a *transpose(matrix([T+S])) |

1094 | Phitu= Phitu[Integer(0) ][Integer(0) ] |

1095 | |

1096 | # Compute all partial derivatives of At and Phitu up to orders 2*N-2 |

1097 | # and 2*N, respectively. Take advantage of the fact that At and Phitu |

1098 | # are sufficiently differentiable functions so that mixed partials |

1099 | # are equal. Thus only need to compute representative partials. |

1100 | # Choose nondecreasing sequences as representative differentiation- |

1101 | # order sequences. |

1102 | # To speed up later computations, create symbolic functions AA and BB |

1103 | # to stand in for the expressions At and Phitu respectively. |

1104 | print "Calculating derivatives of more auxiliary functions..." |

1105 | AA= [function('A'+str(j),*tuple(T+S)) for j in range(r)] |

1106 | At_derivs= self._diff_all(At,T+S,Integer(2) *N-Integer(2) ,sub=hderivs1,sub_final=[thetastar,atC],rekey=AA) |

1107 | BB= function('BB',*tuple(T+S)) |

1108 | Phitu_derivs= self._diff_all(Phitu,T+S,Integer(2) *N,sub=hderivs1,sub_final=[thetastar,atC],rekey=BB,zero_order=Integer(3) ) |

1109 | AABB_derivs= At_derivs |

1110 | AABB_derivs.update(Phitu_derivs) |

1111 | for j in range(r): |

1112 | AABB_derivs[AA[j]] = At[j].subs(thetastar).subs(atC) |

1113 | AABB_derivs[BB] = Phitu.subs(thetastar).subs(atC) |

1114 | print "Calculating second-order differential operator actions..." |

1115 | DD= self._diff_op(AA,BB,AABB_derivs,T+S,a_inv,r,N) |

1116 | |

1117 | L={} |

1118 | for (j,k) in CartesianProduct((ellipsis_range(Integer(0) ,Ellipsis,min(r-Integer(1) ,N-Integer(1) ))), (ellipsis_range(max(Integer(0) ,N-Integer(1) -r),Ellipsis,N-Integer(1) ))): |

1119 | if j+k <= N-Integer(1) : |

1120 | L[(j,k)] = add([ \ |

1121 | DD[(j,k,l)] /( (-Integer(1) )**k *Integer(2) **(k+l) *factorial(l) *factorial(k+l) ) \ |

1122 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k))] ) |

1123 | # The next line's QQ coercion is a workaround for the Sage 4.6 bug reported |

1124 | # on http://trac.sagemath.org/sage_trac/ticket/8659. |

1125 | # Once the bug is fixed, the QQ can be removed. |

1126 | det= QQ(a.determinant())**(-Integer(1) /Integer(2) ) *(Integer(2) *pi)**((r-d)/Integer(2) ) |

1127 | chunk= det *add([ (alpha[d-Integer(1) ]*n)**((r-d)/Integer(2) -q) *add([ \ |

1128 | L[(j,k)] *binomial(r-Integer(1) ,j) *stirling_number1(r-j,r+k-q) *(-Integer(1) )**(q-j-k) \ |

1129 | for (j,k) in CartesianProduct((ellipsis_range(Integer(0) ,Ellipsis,min(r-Integer(1) ,q))), (ellipsis_range(max(Integer(0) ,q-r),Ellipsis,q))) if j+k <= q ]) \ |

1130 | for q in range(N)]) |

1131 | |

1132 | chunk= chunk.subs(C).simplify() |

1133 | coeffs= chunk.coefficients(n) |

1134 | coeffs.reverse() |

1135 | coeffs= coeffs[:N] |

1136 | if numerical: # If a numerical approximation is desired... |

1137 | subexp_part = add( [co[Integer(0) ].subs(c).n(digits=numerical)*n**co[Integer(1) ] \ |

1138 | for co in coeffs] ) |

1139 | exp_scale= (Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )) \ |

1140 | .n(digits=numerical) |

1141 | else: |

1142 | subexp_part = add( [co[Integer(0) ].subs(c)*n**co[Integer(1) ] for co in coeffs] ) |

1143 | exp_scale= Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] ) |

1144 | return exp_scale**n*subexp_part, exp_scale, subexp_part |

1145 | #-------------------------------------------------------------------------------- |

1146 | def _asymptotics_main_smooth(self,G,H,X,c,n,alpha,N,numerical): |

1147 | r""" |

1148 | This function is for internal use by _asymptotics_main(). |

1149 | It calculates asymptotics in case they depend upon smooth points. |

1150 | |

1151 | Does not use `self`. |

1152 | |

1153 | INPUT: |

1154 | |

1155 | - ``G`` - A symbolic expression. |

1156 | - ``H`` - A list of symbolic expressions. |

1157 | - ``X`` - The list of variables occurring in `G` and `H`. |

1158 | Call its length `d`. |

1159 | - ``c`` - A dictionary with `X` as keys and numbers as values. |

1160 | - ``n`` - The variable of the asymptotic expansion. |

1161 | - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic |

1162 | variables if `c` is a smooth point. |

1163 | - ``N`` - A positive integer. |

1164 | - ``numerical`` - Natural number. |

1165 | If k=numerical > 0, then a numerical approximation of the coefficients |

1166 | of `F_{n\alpha}` with k digits of precision will be returned. |

1167 | Otherwise exact values will be returned. |

1168 | |

1169 | OUTPUT: |

1170 | |

1171 | The same as the function asymptotics(). |

1172 | |

1173 | NOTES: |

1174 | |

1175 | The formulas used for computing the asymptotic expansions are |

1176 | Theorems 3.2 and 3.3 [RaWi2008a]_ with `p=1`. |

1177 | Theorem 3.2 is a specialization of Theorem 3.4 of [RaWi2011]_ |

1178 | with `r=1`. |

1179 | |

1180 | REFERENCES: |

1181 | |

1182 | .. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

1183 | coefficients of multivariate generating functions: improvements |

1184 | for smooth points", Electronic Journal of Combinatorics, Vol. 15 |

1185 | (2008), R89. |

1186 | |

1187 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

1188 | coefficients of multivariate generating functions: improvements |

1189 | for smooth points", To appear. |

1190 | |

1191 | AUTHORS: |

1192 | |

1193 | - Alex Raichev (2008-10-01, 2010-09-28) |

1194 | """ |

1195 | I= sqrt(-Integer(1) ) |

1196 | d= len(X) # > 1 |

1197 | |

1198 | # If c is a tuple of rationals, then compute with it directly. |

1199 | # Otherwise, compute symbolically and plug in c at the end. |

1200 | if vector(c.values()) in QQ**d: |

1201 | C= c |

1202 | else: |

1203 | Cs= [var('cs' +str(j)) for j in range(d)] |

1204 | C= dict( [(X[j],Cs[j]) for j in range(d)] ) |

1205 | c= dict( [(Cs[j],c[X[j]]) for j in range(d)] ) |

1206 | |

1207 | # Setup. |

1208 | print "Initializing auxiliary functions..." |

1209 | h= function('h',*tuple(X[:d-Integer(1) ])) # Implicit functions |

1210 | U = function('U',*tuple(X)) # |

1211 | # All other functions are defined in terms of h, U, and explicit functions. |

1212 | Gcheck = -G/U *(h/X[d-Integer(1) ]) |

1213 | A= Gcheck.subs({X[d-Integer(1) ]:Integer(1) /h})/h |

1214 | T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-Integer(1) )] |

1215 | e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-Integer(1) )]) |

1216 | ht= h.subs(e) |

1217 | Ht= H.subs(e).subs({X[d-Integer(1) ]:Integer(1) /ht}) |

1218 | At= A.subs(e) |

1219 | Phit = -log(C[X[d-Integer(1) ]] *ht)\ |

1220 | + I* add([alpha[i]/alpha[d-Integer(1) ] *T[i] for i in range(d-Integer(1) )]) |

1221 | Tstar= dict([(t,Integer(0) ) for t in T]) |

1222 | atC = copy(C) |

1223 | atC.update({h.subs(C):Integer(1) /C[X[d-Integer(1) ]]}) # Stores h and U and all their derivatives |

1224 | # evaluated at C. |

1225 | |

1226 | # Compute the derivatives of h up to order 2*N and evaluate at C and store |

1227 | # in atC. Keep a copy of unevaluated h derivatives for use in the case |

1228 | # d=2 and v > 2 below. |

1229 | hderivs1= {} # First derivatives of h. |

1230 | for i in range(d-Integer(1) ): |

1231 | s= solve( diff(H.subs({X[d-Integer(1) ]:Integer(1) /h}),X[i]), diff(h,X[i]) )[Integer(0) ].rhs()\ |

1232 | .simplify() |

1233 | hderivs1.update({diff(h,X[i]):s}) |

1234 | atC.update({diff(h,X[i]).subs(C):s.subs(C).subs(atC)}) |

1235 | hderivs = self._diff_all(h,X[Integer(0) :d-Integer(1) ],Integer(2) *N,sub=hderivs1,rekey=h) |

1236 | for k in hderivs.keys(): |

1237 | atC.update({k.subs(C):hderivs[k].subs(atC)}) |

1238 | |

1239 | # Compute the derivatives of U up to order 2*N and evaluate at C. |

1240 | # To do this, differentiate H = U*Hcheck over and over, evaluate at C, |

1241 | # and solve for the derivatives of U at C. |

1242 | # Need the derivatives of H with short keys to pass on to diff_prod later. |

1243 | Hderivs= self._diff_all(H,X,Integer(2) *N,ending=[X[d-Integer(1) ]],sub_final=C) |

1244 | print "Calculating derivatives of auxiallary functions..." |

1245 | # For convenience in checking if all the nontrivial derivatives of U at c |

1246 | # are zero a few line below, store the value of U(c) in atC instead of in |

1247 | # Uderivs. |

1248 | Uderivs={} |

1249 | atC.update({U.subs(C):diff(H,X[d-Integer(1) ]).subs(C)}) |

1250 | end= [X[d-Integer(1) ]] |

1251 | Hcheck= X[d-Integer(1) ] - Integer(1) /h |

1252 | p= H.polynomial(CC).degree()-Integer(1) |

1253 | if p == Integer(0) : |

1254 | # Then, using a proposition at the end of [RaWi2011], we can |

1255 | # conclude that all higher derivatives of U are zero. |

1256 | for l in (ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N)): |

1257 | for s in UnorderedTuples(X,l): |

1258 | Uderivs[diff(U,s).subs(C)] = Integer(0) |

1259 | elif p > Integer(0) and p < Integer(2) *N: |

1260 | all_zero= True |

1261 | Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,p)),end,Uderivs,atC) |

1262 | # Check for a nonzero U derivative. |

1263 | if Uderivs.values() != [Integer(0) for i in range(len(Uderivs))]: |

1264 | all_zero= False |

1265 | if all_zero: |

1266 | # Then, using a proposition at the end of [RaWi2011], we can |

1267 | # conclude that all higher derivatives of U are zero. |

1268 | for l in (ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N)): |

1269 | for s in UnorderedTuples(X,l): |

1270 | Uderivs.update({diff(U,s).subs(C):Integer(0) }) |

1271 | else: |

1272 | # Have to compute the rest of the derivatives. |

1273 | Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N)),end,Uderivs,atC) |

1274 | else: |

1275 | Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N)),end,Uderivs,atC) |

1276 | atC.update(Uderivs) |

1277 | |

1278 | # In general, this algorithm is not designed to handle the case of a |

1279 | # singular Phit''(Tstar). However, when d=2 the algorithm can cope. |

1280 | if d==Integer(2) : |

1281 | # Compute v, the order of vanishing at Tstar of Phit. It is at least 2. |

1282 | v=Integer(2) |

1283 | Phitderiv= diff(Phit,T[Integer(0) ],Integer(2) ) |

1284 | splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify() |

1285 | while splat==Integer(0) : |

1286 | v= v+Integer(1) |

1287 | if v > Integer(2) *N: # Then need to compute more derivatives of h for atC. |

1288 | hderivs.update({diff(h,X[Integer(0) ],v) \ |

1289 | :diff(hderivs[diff(h,X[Integer(0) ],v-Integer(1) )],X[Integer(0) ]).subs(hderivs1)}) |

1290 | atC.update({diff(h,X[Integer(0) ],v).subs(C) \ |

1291 | :hderivs[diff(h,X[Integer(0) ],v)].subs(atC)}) |

1292 | Phitderiv= diff(Phitderiv,T[Integer(0) ]) |

1293 | splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify() |

1294 | if d==Integer(2) and v>Integer(2) : |

1295 | t= T[Integer(0) ] # Simplify variable names. |

1296 | a= splat/factorial(v) |

1297 | Phitu= Phit -a*t**v |

1298 | |

1299 | # Compute all partial derivatives of At and Phitu up to orders 2*(N-1) |

1300 | # and 2*(N-1)+v, respectively, in case v is even. |

1301 | # Otherwise, compute up to orders N-1 and N-1+v, respectively. |

1302 | # To speed up later computations, create symbolic functions AA and BB |

1303 | # to stand in for the expressions At and Phitu, respectively. |

1304 | print "Calculating derivatives of more auxiliary functions..." |

1305 | AA= function('AA',t) |

1306 | BB= function('BB',t) |

1307 | if v.mod(Integer(2) )==Integer(0) : |

1308 | At_derivs= self._diff_all(At,T,Integer(2) *N-Integer(2) , \ |

1309 | sub=hderivs1,sub_final=[Tstar,atC],rekey=AA) |

1310 | Phitu_derivs= self._diff_all(Phitu,T,Integer(2) *N-Integer(2) +v, \ |

1311 | sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+Integer(1) ,rekey=BB) |

1312 | else: |

1313 | At_derivs= self._diff_all(At,T,N-Integer(1) ,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA) |

1314 | Phitu_derivs= self._diff_all(Phitu,T,N-Integer(1) +v,sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+Integer(1) ,rekey=BB) |

1315 | AABB_derivs= At_derivs |

1316 | AABB_derivs.update(Phitu_derivs) |

1317 | AABB_derivs[AA] = At.subs(Tstar).subs(atC) |

1318 | AABB_derivs[BB] = Phitu.subs(Tstar).subs(atC) |

1319 | print "Calculating actions of the second order differential operator..." |

1320 | DD= self._diff_op_simple(AA,BB,AABB_derivs,t,v,a,N) |

1321 | # Plug above into asymptotic formula. |

1322 | L = [] |

1323 | if v.mod(Integer(2) ) == Integer(0) : |

1324 | for k in range(N): |

1325 | L.append( add([ \ |

1326 | (-Integer(1) )**l *gamma((Integer(2) *k+v*l+Integer(1) )/v) \ |

1327 | / (factorial(l) *factorial(Integer(2) *k+v*l)) \ |

1328 | * DD[(k,l)] for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)) ]) ) |

1329 | chunk= a**(-Integer(1) /v) /(pi*v) *add([ alpha[d-Integer(1) ]**(-(Integer(2) *k+Integer(1) )/v) \ |

1330 | * L[k] *n**(-(Integer(2) *k+Integer(1) )/v) for k in range(N) ]) |

1331 | else: |

1332 | zeta= exp(I*pi/(Integer(2) *v)) |

1333 | for k in range(N): |

1334 | L.append( add([ \ |

1335 | (-Integer(1) )**l *gamma((k+v*l+Integer(1) )/v) \ |

1336 | / (factorial(l) *factorial(k+v*l)) \ |

1337 | * (zeta**(k+v*l+Integer(1) ) +(-Integer(1) )**(k+v*l)*zeta**(-(k+v*l+Integer(1) ))) \ |

1338 | * DD[(k,l)] for l in (ellipsis_range(Integer(0) ,Ellipsis,k)) ]) ) |

1339 | chunk= abs(a)**(-Integer(1) /v) /(Integer(2) *pi*v) *add([ alpha[d-Integer(1) ]**(-(k+Integer(1) )/v) \ |

1340 | * L[k] *n**(-(k+Integer(1) )/v) for k in range(N) ]) |

1341 | # Asymptotics for d>=2 case. A singular Phit''(Tstar) will cause a crash |

1342 | # in this case. |

1343 | else: |

1344 | Phit1= jacobian(Phit,T).subs(hderivs1) |

1345 | a= jacobian(Phit1,T).subs(hderivs1).subs(Tstar).subs(atC) |

1346 | a_inv= a.inverse() |

1347 | Phitu= Phit -(Integer(1) /Integer(2) ) *matrix([T]) *a *transpose(matrix([T])) |

1348 | Phitu= Phitu[Integer(0) ][Integer(0) ] |

1349 | # Compute all partial derivatives of At and Phitu up to orders 2*N-2 |

1350 | # and 2*N, respectively. Take advantage of the fact that At and Phitu |

1351 | # are sufficiently differentiable functions so that mixed partials |

1352 | # are equal. Thus only need to compute representative partials. |

1353 | # Choose nondecreasing sequences as representative differentiation- |

1354 | # order sequences. |

1355 | # To speed up later computations, create symbolic functions AA and BB |

1356 | # to stand in for the expressions At and Phitu respectively. |

1357 | print "Calculating derivatives of more auxiliary functions..." |

1358 | AA= function('AA',*tuple(T)) |

1359 | At_derivs= self._diff_all(At,T,Integer(2) *N-Integer(2) ,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA) |

1360 | BB= function('BB',*tuple(T)) |

1361 | Phitu_derivs= self._diff_all(Phitu,T,Integer(2) *N,sub=hderivs1,sub_final=[Tstar,atC],rekey=BB,zero_order=Integer(3) ) |

1362 | AABB_derivs= At_derivs |

1363 | AABB_derivs.update(Phitu_derivs) |

1364 | AABB_derivs[AA] = At.subs(Tstar).subs(atC) |

1365 | AABB_derivs[BB] = Phitu.subs(Tstar).subs(atC) |

1366 | print "Calculating actions of the second order differential operator..." |

1367 | DD= self._diff_op(AA,BB,AABB_derivs,T,a_inv,Integer(1) ,N) |

1368 | |

1369 | # Plug above into asymptotic formula. |

1370 | L=[] |

1371 | for k in range(N): |

1372 | L.append( add([ \ |

1373 | DD[(Integer(0) ,k,l)] / ( (-Integer(1) )**k *Integer(2) **(l+k) *factorial(l) *factorial(l+k) ) \ |

1374 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k))]) ) |

1375 | chunk= add([ (Integer(2) *pi)**((Integer(1) -d)/Integer(2) ) *a.determinant()**(-Integer(1) /Integer(2) ) \ |

1376 | *alpha[d-Integer(1) ]**((Integer(1) -d)/Integer(2) -k) *L[k] \ |

1377 | *n**((Integer(1) -d)/Integer(2) -k) for k in range(N) ]) |

1378 | |

1379 | chunk= chunk.subs(c).simplify() |

1380 | coeffs= chunk.coefficients(n) |

1381 | coeffs.reverse() |

1382 | coeffs= coeffs[:N] |

1383 | if numerical: # If a numerical approximation is desired... |

1384 | subexp_part = add( [co[Integer(0) ].subs(c).n(digits=numerical)*n**co[Integer(1) ] for co in coeffs] ) |

1385 | exp_scale=(Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )) \ |

1386 | .n(digits=numerical) |

1387 | else: |

1388 | subexp_part = add( [co[Integer(0) ].subs(c)*n**co[Integer(1) ] for co in coeffs] ) |

1389 | exp_scale= Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] ) |

1390 | return exp_scale**n*subexp_part, exp_scale, subexp_part |

1391 | #----------------------------------------------------------------------------- |

1392 | def _crit_cone_combo(self,fs,X,c,alpha): |

1393 | r""" |

1394 | This function returns an auxiliary point associated to the multiple |

1395 | point `c` of the factors `fs`. |

1396 | It is for internal use by _asymptotics_main_multiple(). |

1397 | |

1398 | INPUT: |

1399 | |

1400 | - ``fs`` - A list of expressions in the variables of `X`. |

1401 | - ``X`` - A list of variables. |

1402 | - ``c`` - A dictionary with keys `X` and values in some field. |

1403 | - ``alpha`` - A list of rationals. |

1404 | |

1405 | OUTPUT: |

1406 | |

1407 | A solution of the matrix equation `y Gamma = a` for `y`, |

1408 | where `Gamma` is the matrix whose `j`th row is |

1409 | _direction(_log_grad(fj,X,c)) where `fj` |

1410 | is the `j`th item in `fs` and where `a` is _direction(alpha). |

1411 | |

1412 | EXAMPLES:: |

1413 | |

1414 | sage: from sage.combinat.amgf import * |

1415 | sage: R.<x,y,z>= PolynomialRing(QQ) |

1416 | sage: G,H= 1,1 |

1417 | sage: F= QuasiRationalExpression(G,H) |

1418 | sage: fs= [x + 2*y + z - 4, 2*x + y + z - 4] |

1419 | sage: c= {x:1,y:1,z:1} |

1420 | sage: alpha= [2,1,1] |

1421 | sage: F._crit_cone_combo(fs,R.gens(),c,alpha) |

1422 | [0, 1] |

1423 | |

1424 | NOTES: |

1425 | |

1426 | Use this function only when `Gamma` is well-defined and |

1427 | there is a unique solution to the matrix equation `y Gamma = a`. |

1428 | Fails otherwise. |

1429 | |

1430 | AUTHORS: |

1431 | |

1432 | - Alex Raichev (2008-10-01, 2008-11-25, 2009-03-04, 2010-09-08, |

1433 | 2010-12-02) |

1434 | """ |

1435 | # Assuming here that each _log_grad(f) has nonzero final component. |

1436 | # Then 'direction' will not throw a division by zero error. |

1437 | d= len(X) |

1438 | r= len(fs) |

1439 | Gamma= matrix([self._direction(self._log_grad(f,X,c)) for f in fs]) |

1440 | # solve_left() fails when working in SR :-(. So use solve() instead. |

1441 | #s= Gamma.solve_left(vector(alpha)/alpha[d-1]) |

1442 | V= [var('sss'+str(i)) for i in range(r)] |

1443 | M= matrix(V)*Gamma |

1444 | eqns= [M[Integer(0) ][i]== alpha[i]/alpha[d-Integer(1) ] for i in range(d)] |

1445 | s= solve(eqns,V,solution_dict=True)[Integer(0) ] # Assuming a unique solution. |

1446 | return [s[v] for v in V] |

1447 | # B ========================================================================== |

1448 | # C ========================================================================== |

1449 | def cohomologous_integrand(self,alpha,asy_var=None): |

1450 | r""" |

1451 | This function takes a multivariate Cauchy type integral |

1452 | `\int F / x^{asy_var\alpha+1} dx`, where `F=self`, and breaks it up |

1453 | into a list of nicer Cauchy type integrals for the purposes of |

1454 | computing asymptotics of the original integral as `asy_var\to\infty`. |

1455 | The sum of the nicer integrals is de Rham cohomologous to the original |

1456 | integral. |

1457 | It assumes that algebraic varieties corresponding to the irreducible |

1458 | factors of `self._H` intersect transversely (see notes below). |

1459 | |

1460 | INPUT: |

1461 | |

1462 | - ``alpha`` - A list of positive integers or symbolic variables. |

1463 | - ``asy_var`` - A symbolic variable (default: None). |

1464 | Eventually set to `var('n')` if None is given. |

1465 | |

1466 | OUTPUT: |

1467 | |

1468 | A list of the form `[chunk_1,\ldots,chunk_r]`, where each |

1469 | `chunk_j` has the form `[P,[B_1,1],\ldots,[B_l,1]]`. |

1470 | Here `l \le d`, `P` is a symbolic expression in the indeterminates of |

1471 | `R` and `asy_var`, `\{B_1,\ldots,B_l\} \subseteq \{Q_1,\ldots,Q_m\}`, |

1472 | and `[P,[B_1,1],\ldots,[B_l,1]]` represents the integral |

1473 | `\int P/(B_1 \cdots B_l) dx`. |

1474 | |

1475 | EXAMPLES:: |

1476 | |

1477 | sage: from sage.combinat.amgf import * |

1478 | sage: R.<x,y>= PolynomialRing(QQ) |

1479 | sage: G= 9*exp(x+y) |

1480 | sage: H= (3-2*x-y)*(3-x-2*y) |

1481 | sage: F= QuasiRationalExpression(G,H) |

1482 | sage: alpha= [4,3] |

1483 | sage: F.cohomologous_integrand(alpha) |

1484 | [[9*e^(x + y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]] |

1485 | |

1486 | sage: from sage.combinat.amgf import * |

1487 | sage: R.<x,y>= PolynomialRing(QQ) |

1488 | sage: G= 9*exp(x+y) |

1489 | sage: H= (3-2*x-y)^2*(3-x-2*y) |

1490 | sage: F= QuasiRationalExpression(G,H) |

1491 | sage: alpha= [4,3] |

1492 | sage: F.cohomologous_integrand(alpha) |

1493 | [[-3*(3*x*e^x - 8*y*e^x)*n*e^y/(x*y) - 3*((x - 2)*y*e^x + x*e^x)*e^y/(x*y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]] |

1494 | |

1495 | sage: from sage.combinat.amgf import * |

1496 | sage: R.<x,y,z>= PolynomialRing(QQ) |

1497 | sage: G= 16 |

1498 | sage: H= (4-2*x-y-z)^2*(4-x-2*y-z) |

1499 | sage: F= QuasiRationalExpression(G,H) |

1500 | sage: alpha= [3,3,2] |

1501 | sage: F.cohomologous_integrand(alpha) |

1502 | [[16*(4*y - 3*z)*n/(y*z) + 16*(2*y - z)/(y*z), [x + 2*y + z - 4, 1], [2*x + y + z - 4, 1]]] |

1503 | |

1504 | NOTES: |

1505 | |

1506 | The varieties corresponding to `Q_1,\ldots,Q_m` |

1507 | __intersect transversely__ iff for each point `c` of their intersection |

1508 | and for all `k \le l`, the Jacobian matrix of any `k` polynomials from |

1509 | `\{Q_1,\ldots,Q_m\}` has rank equal to `\min\{k,d\}` when evaluated at |

1510 | `c`. |

1511 | |

1512 | ALGORITHM: |

1513 | |

1514 | Let `asy_var= n` and consider the integral around a polycirle centered |

1515 | at the origin of the `d`-variate differential form |

1516 | `\frac{G(x) dx_1 \wedge\cdots\wedge dx_d}{H(x) x^{n\alpha+1}}`, where |

1517 | `G=self._G` and `H=self._H`. |

1518 | |

1519 | (1) Decompose `G/H` into a sum of partial fractions `P_1 +\cdots+ P_r` |

1520 | so that each term of the sum has at most `d` irreducible factors of `H` |

1521 | in the denominator. |

1522 | |

1523 | (2) For each differential form `P_j dx_1 \wedge\cdots\wedge dx_d`, |

1524 | find an equivalent form `\omega_j` in de Rham cohomology with no |

1525 | repeated irreducible factors of `H` in its denominator. |

1526 | |

1527 | AUTHOR: |

1528 | |

1529 | - Alex Raichev (2010-09-22) |

1530 | """ |

1531 | if not asy_var: |

1532 | asy_var = var('n') |

1533 | |

1534 | # Create symbolic (non-ring) variables. |

1535 | G= self._G |

1536 | H= self._H |

1537 | R= self._R |

1538 | d= self._d |

1539 | X= self._variables |

1540 | # Prepare input for partial_fraction_decomposition() which only works |

1541 | # for functions in the field of fractions of R. |

1542 | if G in R: |

1543 | numer=Integer(1) |

1544 | F= G/H |

1545 | else: |

1546 | numer= G |

1547 | F= R(Integer(1) )/H |

1548 | nstuff= Integer(1) /mul([X[j]**(alpha[j]*asy_var+Integer(1) ) for j in range(d)]) |

1549 | |

1550 | # Do steps (2) and (1). |

1551 | integrands= [] |

1552 | whole,parts= partial_fraction_decomposition(F) |

1553 | for f in parts: |

1554 | a= format_quotient(f) |

1555 | integrands.append( [a[Integer(0) ]*numer*nstuff] + a[Integer(1) :] ) |

1556 | |

1557 | # Do step (3). |

1558 | ce= decompose_via_cohomology(integrands) |

1559 | ce_new= [] |

1560 | for a in ce: |

1561 | ce_new.append( [(a[Integer(0) ]/nstuff).simplify_full().collect(n)] + a[Integer(1) :] ) |

1562 | return ce_new |

1563 | #------------------------------------------------------------------------------- |

1564 | def critical_cone(self,c,coordinate=None): |

1565 | r""" |

1566 | Returns the critical cone of a convenient multiple point `c`. |

1567 | |

1568 | INPUT: |

1569 | |

1570 | - ``c`` - A dictionary with keys `self.variables()` and values |

1571 | in a field. |

1572 | - ``coordinate`` - (optional) A natural number. |

1573 | |

1574 | OUTPUT: |

1575 | |

1576 | A list of vectors that generate the critical cone of `c` and |

1577 | the cone itself if the values of `c` lie in QQ. |

1578 | |

1579 | EXAMPLES:: |

1580 | |

1581 | sage: from sage.combinat.amgf import * |

1582 | sage: R.<x,y,z>= PolynomialRing(QQ) |

1583 | sage: G= 1 |

1584 | sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) |

1585 | sage: F= QuasiRationalExpression(G,H) |

1586 | sage: c= {z: 4/3, y: 1, x: 1/2} |

1587 | sage: F.critical_cone(c) |

1588 | ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N) |

1589 | |

1590 | NOTES: |

1591 | |

1592 | The _critical cone_ of a convenient multiple point `c` with |

1593 | with `c_k \del_k H_j(c) \neq 0` for all `j=1,\ldots,r` is |

1594 | the conical hull of the vectors `\gamma_j(c) = |

1595 | \left(\frac{c_1 \del_1 H_j(c)}{c_k \del_k H_j(c)},\ldots, |

1596 | \frac{c_d \del_d H_j(c)}{c_k \del_k H_j(c)} \right)`. |

1597 | Here `H_1,\ldots,H_r` are the irreducible germs of `self._H` around `c`. |

1598 | For more details, see [RaWi2011]_. |

1599 | |

1600 | If this function's optional argument `coordinate` isn't given, then |

1601 | this function searches (from `d` down to 1) for the first index `k` |

1602 | such that for all `j=1,\ldots,r` we have `c_k \del_k H_j(c) \neq 0` |

1603 | and sets `coordinate = k`. |

1604 | Almost. |

1605 | Since this is Python, all the indices actually start at 0. |

1606 | |

1607 | REFERENCES: |

1608 | |

1609 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

1610 | coefficients of multivariate generating functions: improvements |

1611 | for smooth points", To appear. |

1612 | |

1613 | AUTHORS: |

1614 | |

1615 | - Alex Raichev (2010-08-25) |

1616 | """ |

1617 | Hs= [SR(h[Integer(0) ]) for h in self._Hfac] # irreducible factors of H |

1618 | X= self._variables |

1619 | d= self._d |

1620 | # Ensure the variables of `c` lie in SR |

1621 | cc= {} |

1622 | for x in c.keys(): |

1623 | cc[SR(x)] = c[x] |

1624 | lg= [self._log_grad(h,X,cc) for h in Hs] |

1625 | if not coordinate: |

1626 | # Search (from d-1 down to 0) for a coordinate k such that |

1627 | # for all h in Hs we have cc[k] * diff(h,X[k]) !=0. |

1628 | # One is guaranteed to exist in the case of a convenient multiple |

1629 | # point. |

1630 | for k in reversed(range(d)): |

1631 | if Integer(0) not in [v[k] for v in lg]: |

1632 | coordinate= k |

1633 | break |

1634 | gamma= [self._direction(v,coordinate) for v in lg] |

1635 | if [[gg in QQ for gg in g] for g in gamma] == \ |

1636 | [[True for gg in g] for g in gamma]: |

1637 | return gamma,Cone(gamma) # Cone() needs rational vectors |

1638 | else: |

1639 | return gamma |

1640 | # D ============================================================================ |

1641 | def _diff_all(self,f,V,n,ending=[],sub=None,sub_final=None,zero_order=Integer(0) ,rekey=None): |

1642 | r""" |

1643 | This function returns a dictionary of representative mixed partial |

1644 | derivatives of `f` from order 1 up to order `n` with respect to the |

1645 | variables in `V`. |

1646 | The default is to key the dictionary by all nondecreasing sequences |

1647 | in `V` of length 1 up to length `n`. |

1648 | For internal use. |

1649 | |

1650 | Does not use `self`. |

1651 | |

1652 | INPUT: |

1653 | |

1654 | - ``f`` - An individual or list of `\mathcal{C}^{n+1}` functions. |

1655 | - ``V`` - A list of variables occurring in `f`. |

1656 | - ``n`` - A natural number. |

1657 | - ``ending`` - A list of variables in `V`. |

1658 | - ``sub`` - An individual or list of dictionaries. |

1659 | - ``sub_final`` - An individual or list of dictionaries. |

1660 | - ``rekey`` - A callable symbolic function in `V` or list thereof. |

1661 | - ``zero_order`` - A natural number. |

1662 | |

1663 | OUTPUT: |

1664 | |

1665 | The dictionary `{s_1:deriv_1,...,s_r:deriv_r}`. |

1666 | Here `s_1,\ldots,s_r` is a listing of |

1667 | all nondecreasing sequences of length 1 up to length `n` over the |

1668 | alphabet `V`, where `w > v` in `X` iff `str(w) > str(v)`, and |

1669 | `deriv_j` is the derivative of `f` with respect to the derivative |

1670 | sequence `s_j` and simplified with respect to the substitutions in `sub` |

1671 | and evaluated at `sub_final`. |

1672 | Moreover, all derivatives with respect to sequences of length less than |

1673 | `zero_order` (derivatives of order less than `zero_order`) will be made |

1674 | zero. |

1675 | |

1676 | If `rekey` is nonempty, then `s_1,\ldots,s_r` will be replaced by the |

1677 | symbolic derivatives of the functions in `rekey`. |

1678 | |

1679 | If `ending` is nonempty, then every derivative sequence `s_j` will be |

1680 | suffixed by `ending`. |

1681 | |

1682 | EXAMPLES:: |

1683 | |

1684 | I'd like to print the entire dictionaries, but that doesn't yield |

1685 | consistent output order for doctesting. |

1686 | Order of keys changes. :: |

1687 | |

1688 | sage: from sage.combinat.amgf import * |

1689 | sage: R.<x> = PolynomialRing(QQ) |

1690 | sage: G,H = 1,1 |

1691 | sage: F= QuasiRationalExpression(G,H) |

1692 | sage: f= function('f',x) |

1693 | sage: dd= F._diff_all(f,[x],3) |

1694 | sage: dd[(x,x,x)] |

1695 | D[0, 0, 0](f)(x) |

1696 | |

1697 | :: |

1698 | |

1699 | sage: d1= {diff(f,x): 4*x^3} |

1700 | sage: dd= F._diff_all(f,[x],3,sub=d1) |

1701 | sage: dd[(x,x,x)] |

1702 | 24*x |

1703 | |

1704 | :: |

1705 | |

1706 | sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f) |

1707 | sage: dd[diff(f,x,3)] |

1708 | 24*x |

1709 | |

1710 | :: |

1711 | |

1712 | sage: a= {x:1} |

1713 | sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f,sub_final=a) |

1714 | sage: dd[diff(f,x,3)] |

1715 | 24 |

1716 | |

1717 | :: |

1718 | |

1719 | sage: X= var('x,y,z') |

1720 | sage: f= function('f',*X) |

1721 | sage: dd= F._diff_all(f,X,2,ending=[y,y,y]) |

1722 | sage: dd[(z,y,y,y)] |

1723 | D[1, 1, 1, 2](f)(x, y, z) |

1724 | |

1725 | :: |

1726 | |

1727 | sage: g= function('g',*X) |

1728 | sage: dd= F._diff_all([f,g],X,2) |

1729 | sage: dd[(0,y,z)] |

1730 | D[1, 2](f)(x, y, z) |

1731 | |

1732 | :: |

1733 | |

1734 | sage: dd[(1,z,z)] |

1735 | D[2, 2](g)(x, y, z) |

1736 | |

1737 | :: |

1738 | |

1739 | sage: f= exp(x*y*z) |

1740 | sage: ff= function('ff',*X) |

1741 | sage: dd= F._diff_all(f,X,2,rekey=ff) |

1742 | sage: dd[diff(ff,x,z)] |

1743 | x*y^2*z*e^(x*y*z) + y*e^(x*y*z) |

1744 | |

1745 | AUTHORS: |

1746 | |

1747 | - Alex Raichev (2008-10-01, 2009-04-01, 2010-02-01) |

1748 | """ |

1749 | singleton=False |

1750 | if not isinstance(f,list): |

1751 | f= [f] |

1752 | singleton=True |

1753 | |

1754 | # Build the dictionary of derivatives iteratively from a list of nondecreasing |

1755 | # derivative-order sequences. |

1756 | derivs= {} |

1757 | r= len(f) |

1758 | if ending: |

1759 | seeds = [ending] |

1760 | start = Integer(1) |

1761 | else: |

1762 | seeds = [[v] for v in V] |

1763 | start = Integer(2) |

1764 | if singleton: |

1765 | for s in seeds: |

1766 | derivs[tuple(s)] = self._subs_all(diff(f[Integer(0) ],s),sub) |

1767 | for l in (ellipsis_range(start,Ellipsis,n)): |

1768 | for t in UnorderedTuples(V,l): |

1769 | s= tuple(t + ending) |

1770 | derivs[s] = self._subs_all(diff(derivs[s[Integer(1) :]],s[Integer(0) ]),sub) |

1771 | else: |

1772 | # Make the dictionary keys of the form (j,sequence of variables), |

1773 | # where j in range(r). |

1774 | for s in seeds: |

1775 | value= self._subs_all([diff(f[j],s) for j in range(r)],sub) |

1776 | derivs.update(dict([(tuple([j]+s),value[j]) for j in range(r)])) |

1777 | for l in (ellipsis_range(start,Ellipsis,n)): |

1778 | for t in UnorderedTuples(V,l): |

1779 | s= tuple(t + ending) |

1780 | value= self._subs_all(\ |

1781 | [diff(derivs[(j,)+s[Integer(1) :]],s[Integer(0) ]) for j in range(r)],sub) |

1782 | derivs.update(dict([((j,)+s,value[j]) for j in range(r)])) |

1783 | if zero_order: |

1784 | # Zero out all the derivatives of order < zero_order |

1785 | if singleton: |

1786 | for k in derivs.keys(): |

1787 | if len(k) < zero_order: |

1788 | derivs[k]= Integer(0) |

1789 | else: |

1790 | # Ignore the first of element of k, which is an index. |

1791 | for k in derivs.keys(): |

1792 | if len(k)-Integer(1) < zero_order: |

1793 | derivs[k]= Integer(0) |

1794 | if sub_final: |

1795 | # Substitute sub_final into the values of derivs. |

1796 | for k in derivs.keys(): |

1797 | derivs[k] = self._subs_all(derivs[k],sub_final) |

1798 | if rekey: |

1799 | # Rekey the derivs dictionary by the value of rekey. |

1800 | F= rekey |

1801 | if singleton: |

1802 | # F must be a singleton. |

1803 | derivs= dict( [(diff(F,list(k)), derivs[k]) for k in derivs.keys()] ) |

1804 | else: |

1805 | # F must be a list. |

1806 | derivs= dict( [(diff(F[k[Integer(0) ]],list(k)[Integer(1) :]), derivs[k]) for k in derivs.keys()] ) |

1807 | return derivs |

1808 | #------------------------------------------------------------------------------- |

1809 | def _diff_op(self,A,B,AB_derivs,V,M,r,N): |

1810 | r""" |

1811 | This function computes the derivatives `DD^(l+k)(A[j] B^l)` evaluated at a |

1812 | point `p` for various natural numbers `j,k,l` which depend on `r` and `N`. |

1813 | Here `DD` is a specific second-order linear differential operator that depends |

1814 | on `M`, `A` is a list of symbolic functions, `B` is symbolic function, |

1815 | and `AB_derivs` contains all the derivatives of `A` and `B` evaluated at `p` |

1816 | that are necessary for the computation. |

1817 | For internal use by the functions _asymptotics_main_multiple() and |

1818 | _asymptotics_main_smooth(). |

1819 | |

1820 | Does not use `self`. |

1821 | |

1822 | INPUT: |

1823 | |

1824 | - ``A`` - A single or length `r` list of symbolic functions in the |

1825 | variables `V`. |

1826 | - ``B`` - A symbolic function in the variables `V`. |

1827 | - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of |

1828 | `A[0],\ldots,A[r-1]` up to order `2N-2` and |

1829 | the (symbolic) derivatives of `B` up to order `2N`. |

1830 | The values of the dictionary are complex numbers that are |

1831 | the keys evaluated at a common point `p`. |

1832 | - ``V`` - The variables of the `A[j]` and `B`. |

1833 | - ``M`` - A symmetric `l \times l` matrix, where `l` is the length of `V`. |

1834 | - ``r,N`` - Natural numbers. |

1835 | |

1836 | OUTPUT: |

1837 | |

1838 | A dictionary whose keys are natural number tuples of the form `(j,k,l)`, |

1839 | where `l \le 2k`, `j \le r-1`, and `j+k \le N-1`, and whose values are |

1840 | `DD^(l+k)(A[j] B^l)` evaluated at a point `p`, where `DD` is the linear |

1841 | second-order differential operator |

1842 | `-\sum_{i=0}^{l-1} \sum_{j=0}^{l-1} M[i][j] |

1843 | \partial^2 /(\partial V[j] \partial V[i])`. |

1844 | |

1845 | EXAMPLES:: |

1846 | |

1847 | sage: from sage.combinat.amgf import * |

1848 | sage: R.<x> = PolynomialRing(QQ) |

1849 | sage: G,H = 1,1 |

1850 | sage: F= QuasiRationalExpression(G,H) |

1851 | sage: T= var('x,y') |

1852 | sage: A= function('A',*tuple(T)) |

1853 | sage: B= function('B',*tuple(T)) |

1854 | sage: AB_derivs= {} |

1855 | sage: M= matrix([[1,2],[2,1]]) |

1856 | sage: DD= F._diff_op(A,B,AB_derivs,T,M,1,2) |

1857 | sage: DD.keys() |

1858 | [(0, 1, 2), (0, 1, 1), (0, 1, 0), (0, 0, 0)] |

1859 | sage: len(DD[(0,1,2)]) |

1860 | 246 |

1861 | |

1862 | AUTHORS: |

1863 | |

1864 | - Alex Raichev (2008-10-01, 2010-01-12) |

1865 | """ |

1866 | if not isinstance(A,list): |

1867 | A= [A] |

1868 | |

1869 | # First, compute the necessary product derivatives of A and B. |

1870 | product_derivs= {} |

1871 | for (j,k) in mrange([r,N]): |

1872 | if j+k <N: |

1873 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)): |

1874 | for s in UnorderedTuples(V,Integer(2) *(k+l)): |

1875 | product_derivs[tuple([j,k,l]+s)] = \ |

1876 | diff(A[j]*B**l,s).subs(AB_derivs) |

1877 | |

1878 | # Second, compute DD^(k+l)(A[j]*B^l)(p) and store values in dictionary. |

1879 | DD= {} |

1880 | rows= M.nrows() |

1881 | for (j,k) in mrange([r,N]): |

1882 | if j+k <N: |

1883 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)): |

1884 | # Take advantage of the symmetry of M by ignoring |

1885 | # the upper-diagonal entries of M and multiplying by |

1886 | # appropriate powers of 2. |

1887 | if k+l == Integer(0) : |

1888 | DD[(j,k,l)] = product_derivs[(j,k,l)] |

1889 | continue |

1890 | S= [(a,b) for (a,b) in mrange([rows,rows]) if b <= a] |

1891 | P= cartesian_product_iterator([S for i in range(k+l)]) |

1892 | diffo= Integer(0) |

1893 | for t in P: |

1894 | if product_derivs[(j,k,l)+self._diff_seq(V,t)] != Integer(0) : |

1895 | MM= Integer(1) |

1896 | for (a,b) in t: |

1897 | MM= MM * M[a][b] |

1898 | if a != b: |

1899 | MM= Integer(2) *MM |

1900 | diffo= diffo + MM * product_derivs[(j,k,l)+self._diff_seq(V,t)] |

1901 | DD[(j,k,l)] = (-Integer(1) )**(k+l)*diffo |

1902 | return DD |

1903 | #------------------------------------------------------------------------------- |

1904 | def _diff_op_simple(self,A,B,AB_derivs,x,v,a,N): |

1905 | r""" |

1906 | This function computes `DD^(e k + v l)(A B^l)` evaluated at a point `p` |

1907 | for various natural numbers `e,k,l` that depend on `v` and `N`. |

1908 | Here `DD` is a specific linear differential operator that depends |

1909 | on `a` and `v`, `A` and `B` are symbolic functions, and `AB_derivs` contains |

1910 | all the derivatives of `A` and `B` evaluated at `p` that are necessary for |

1911 | the computation. |

1912 | For internal use by the function _asymptotics_main_smooth(). |

1913 | |

1914 | Does not use `self`. |

1915 | |

1916 | INPUT: |

1917 | |

1918 | - ``A``,``B`` - Symbolic functions in the variable `x`. |

1919 | - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of |

1920 | `A` up to order `2N` if `v` is even or `N` if `v` is odd and |

1921 | the (symbolic) derivatives of `B` up to order `2N+v` if `v` is even |

1922 | or `N+v` if `v` is odd. |

1923 | The values of the dictionary are complex numbers that are |

1924 | the keys evaluated at a common point `p`. |

1925 | - ``x`` - Symbolic variable. |

1926 | - ``a`` - A complex number. |

1927 | - ``v``,``N`` - Natural numbers. |

1928 | |

1929 | OUTPUT: |

1930 | |

1931 | A dictionary whose keys are natural number pairs of the form `(k,l)`, |

1932 | where `k < N` and `l \le 2k` and whose values are |

1933 | `DD^(e k + v l)(A B^l)` evaluated at a point `p`. |

1934 | Here `e=2` if `v` is even, `e=1` if `v` is odd, and `DD` is a particular |

1935 | linear differential operator |

1936 | `(a^{-1/v} d/dt)' if `v` is even and `(|a|^{-1/v} i \sgn a d/dt)` |

1937 | if `v` is odd. |

1938 | |

1939 | EXAMPLES:: |

1940 | |

1941 | sage: from sage.combinat.amgf import * |

1942 | sage: R.<x> = PolynomialRing(QQ) |

1943 | sage: G,H = 1,1 |

1944 | sage: F= QuasiRationalExpression(G,H) |

1945 | sage: A= function('A',x) |

1946 | sage: B= function('B',x) |

1947 | sage: AB_derivs= {} |

1948 | sage: F._diff_op_simple(A,B,AB_derivs,x,3,2,2) |

1949 | {(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)} |

1950 | |

1951 | AUTHORS: |

1952 | |

1953 | - Alex Raichev (2010-01-15) |

1954 | """ |

1955 | I= sqrt(-Integer(1) ) |

1956 | DD= {} |

1957 | if v.mod(Integer(2) ) == Integer(0) : |

1958 | for k in range(N): |

1959 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)): |

1960 | DD[(k,l)] = (a**(-Integer(1) /v))**(Integer(2) *k+v*l) \ |

1961 | * diff(A*B**l,x,Integer(2) *k+v*l).subs(AB_derivs) |

1962 | else: |

1963 | for k in range(N): |

1964 | for l in (ellipsis_range(Integer(0) ,Ellipsis,k)): |

1965 | DD[(k,l)] = (abs(a)**(-Integer(1) /v)*I*a/abs(a))**(k+v*l) \ |

1966 | * diff(A*B**l,x,k+v*l).subs(AB_derivs) |

1967 | return DD |

1968 | #------------------------------------------------------------------------------- |

1969 | def _diff_prod(self,f_derivs,u,g,X,interval,end,uderivs,atc): |

1970 | r""" |

1971 | This function takes various derivatives of the equation `f=ug`, |

1972 | evaluates at a point `c`, and solves for the derivatives of `u`. |

1973 | For internal use by the function _asymptotics_main_multiple(). |

1974 | |

1975 | Does not use `self`. |

1976 | |

1977 | INPUT: |

1978 | |

1979 | - ``f_derivs`` - A dictionary whose keys are tuples \code{s + end} for all |

1980 | `X`-variable sequences `s` with length in `interval` and whose |

1981 | values are the derivatives of a function `f` evaluated at `c`. |

1982 | - ``u`` - A callable symbolic function. |

1983 | - ``g`` - An expression or callable symbolic function. |

1984 | - ``X`` - A list of symbolic variables. |

1985 | - ``interval`` - A list of positive integers. |

1986 | Call the first and last values `n` and `nn`, respectively. |

1987 | - ``end`` - A possibly empty list of `z`'s where `z` is the last element of |

1988 | `X`. |

1989 | - ``uderivs`` - A dictionary whose keys are the symbolic |

1990 | derivatives of order 0 to order `n-1` of `u` evaluated at `c` |

1991 | and whose values are the corresponding derivatives evaluated at `c`. |

1992 | - ``atc`` - A dictionary whose keys are the keys of `c` and all the symbolic |

1993 | derivatives of order 0 to order `nn` of `g` evaluated `c` and whose |

1994 | values are the corresponding derivatives evaluated at `c`. |

1995 | |

1996 | OUTPUT: |

1997 | |

1998 | A dictionary whose keys are the derivatives of `u` up to order |

1999 | `nn` and whose values are those derivatives evaluated at `c`. |

2000 | |

2001 | EXAMPLES:: |

2002 | |

2003 | I'd like to print out the entire dictionary, but that does not give |

2004 | consistent output for doctesting. |

2005 | Order of keys changes :: |

2006 | |

2007 | sage: from sage.combinat.amgf import * |

2008 | sage: R.<x> = PolynomialRing(QQ) |

2009 | sage: G,H = 1,1 |

2010 | sage: F= QuasiRationalExpression(G,H) |

2011 | sage: u= function('u',x) |

2012 | sage: g= function('g',x) |

2013 | sage: dd= F._diff_prod({(x,):1,(x,x):1},u,g,[x],[1,2],[],{u(x=2):1},{x:2,g(x=2):3,diff(g,x)(x=2):5, diff(g,x,x)(x=2):7}) |

2014 | sage: dd[diff(u,x,2)(x=2)] |

2015 | 22/9 |

2016 | |

2017 | NOTES: |

2018 | |

2019 | This function works by differentiating the equation `f=ug` with respect |

2020 | to the variable sequence \code{s+end}, for all tuples `s` of `X` of |

2021 | lengths in `interval`, evaluating at the point `c`, |

2022 | and solving for the remaining derivatives of `u`. |

2023 | This function assumes that `u` never appears in the differentiations of |

2024 | `f=ug` after evaluating at `c`. |

2025 | |

2026 | AUTHORS: |

2027 | |

2028 | - Alex Raichev (2009-05-14, 2010-01-21) |

2029 | """ |

2030 | for l in interval: |

2031 | D= {} |

2032 | rhs= [] |

2033 | lhs= [] |

2034 | for t in UnorderedTuples(X,l): |

2035 | s= t+end |

2036 | lhs.append(f_derivs[tuple(s)]) |

2037 | rhs.append(diff(u*g,s).subs(atc).subs(uderivs)) |

2038 | # Since Sage's solve command can't take derivatives as variable |

2039 | # names, i make new variables based on t to stand in for |

2040 | # diff(u,t) and store them in D. |

2041 | D[diff(u,t).subs(atc)] = self._make_var([var('zing')]+t) |

2042 | eqns=[ lhs[i] == rhs[i].subs(uderivs).subs(D) for i in range(len(lhs))] |

2043 | vars= D.values() |

2044 | sol= solve(eqns,*vars,solution_dict=True) |

2045 | uderivs.update(self._subs_all(D,sol[Integer(0) ])) |

2046 | return uderivs |

2047 | #------------------------------------------------------------------------------- |

2048 | def _diff_seq(self,V,s): |

2049 | r""" |

2050 | Given a list `s` of tuples of natural numbers, this function returns the |

2051 | list of elements of `V` with indices the elements of the elements of `s`. |

2052 | This function is for internal use by the function _diff_op(). |

2053 | |

2054 | Does not use `self`. |

2055 | |

2056 | INPUT: |

2057 | |

2058 | - ``V`` - A list. |

2059 | - ``s`` - A list of tuples of natural numbers in the interval |

2060 | \code{range(len(V))}. |

2061 | |

2062 | OUTPUT: |

2063 | |

2064 | The tuple \code{tuple([V[tt] for tt in sorted(t)])}, where `t` is the |

2065 | list of elements of the elements of `s`. |

2066 | |

2067 | EXAMPLES:: |

2068 | |

2069 | sage: from sage.combinat.amgf import * |

2070 | sage: R.<x> = PolynomialRing(QQ) |

2071 | sage: G,H = 1,1 |

2072 | sage: F= QuasiRationalExpression(G,H) |

2073 | sage: V= list(var('x,t,z')) |

2074 | sage: F._diff_seq(V,([0,1],[0,2,1],[0,0])) |

2075 | (x, x, x, x, t, t, z) |

2076 | |

2077 | AUTHORS: |

2078 | |

2079 | - Alex Raichev (2009.05.19) |

2080 | """ |

2081 | t= [] |

2082 | for ss in s: |

2083 | t.extend(ss) |

2084 | return tuple([V[tt] for tt in sorted(t)]) |

2085 | #------------------------------------------------------------------------------- |

2086 | def _direction(self,v,coordinate=None): |

2087 | r""" |

2088 | This function returns \code{[vv/v[coordinate] for vv in v]} where |

2089 | coordinate is the last index of v if not specified otherwise. |

2090 | |

2091 | Does not use `self`. |

2092 | |

2093 | INPUT: |

2094 | |

2095 | - ``v`` - A vector. |

2096 | - ``coordinate`` - An index for `v` (default: None). |

2097 | |

2098 | OUTPUT: |

2099 | |

2100 | This function returns \code{[vv/v[coordinate] for vv in v]} where |

2101 | coordinate is the last index of v if not specified otherwise. |

2102 | |

2103 | EXAMPLES:: |

2104 | |

2105 | sage: from sage.combinat.amgf import * |

2106 | sage: R.<x> = PolynomialRing(QQ) |

2107 | sage: G,H = 1,1 |

2108 | sage: F= QuasiRationalExpression(G,H) |

2109 | sage: F._direction([2,3,5]) |

2110 | (2/5, 3/5, 1) |

2111 | sage: F._direction([2,3,5],0) |

2112 | (1, 3/2, 5/2) |

2113 | |

2114 | AUTHORS: |

2115 | |

2116 | - Alex Raichev (2010-08-25) |

2117 | """ |

2118 | if coordinate == None: |

2119 | coordinate= len(v)-Integer(1) |

2120 | try: |

2121 | v[Integer(0) ].variables() |

2122 | return tuple([(vv/v[coordinate]).simplify_full() for vv in v]) |

2123 | except: |

2124 | return tuple([vv/v[coordinate] for vv in v]) |

2125 | # E ============================================================================ |

2126 | # F ============================================================================ |

2127 | # G ============================================================================ |

2128 | # H ============================================================================ |

2129 | # I ============================================================================ |

2130 | def is_cmp(self,points): |

2131 | r""" |

2132 | Checks if the points in the list `points` are convenient multiple |

2133 | points of `V= \{ x\in CC^d : H(x) = 0\}`, where `H=self._H`. |

2134 | |

2135 | INPUT: |

2136 | |

2137 | - ``points`` - An individual or list of dictionaries with keys |

2138 | `self._variables` and values in some superfield of |

2139 | `self._R.base_ring()`. |

2140 | |

2141 | OUTPUT: |

2142 | |

2143 | A list of tuples `(p,verdict,comment)`, one for each point |

2144 | `p` in `points`, where `verdict` is True if `p` is a convenient |

2145 | multiple point and False otherwise, and where `comment` is a string |

2146 | comment relating to `verdict`, such as 'not a transverse intersection'. |

2147 | |

2148 | EXAMPLES:: |

2149 | |

2150 | sage: from sage.combinat.amgf import * |

2151 | sage: R.<x,y,z>= PolynomialRing(QQ) |

2152 | sage: G= 16 |

2153 | sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) |

2154 | sage: F= QuasiRationalExpression(G,H) |

2155 | sage: points= [{x:1/2,y:1,z:4/3},{x:1/2,y:1,z:2}] |

2156 | sage: F.is_cmp(points) |

2157 | [({y: 1, z: 4/3, x: 1/2}, True, 'all good'), ({y: 1, z: 2, x: 1/2}, False, 'not a singular point')] |

2158 | |

2159 | NOTES: |

2160 | |

2161 | A point `c` of `V` is a __convenient multiple point__ if `V` is locally |

2162 | a union of complex manifolds that intersect transversely at `c`; |

2163 | see [RaWi2011]_. |

2164 | |

2165 | REFERENCES: |

2166 | |

2167 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

2168 | coefficients of multivariate generating functions: improvements |

2169 | for smooth points", To appear. |

2170 | |

2171 | AUTHORS: |

2172 | |

2173 | - Alex Raichev (2011-04-18) |

2174 | """ |

2175 | H= self._H |

2176 | d= self._d |

2177 | Hs= [SR(h[Integer(0) ]) for h in self._Hfac] # irreducible factors of H |

2178 | X= self._variables |

2179 | J= [tuple([diff(h,x) for x in X]) for h in Hs] |

2180 | verdicts= [] |

2181 | if not isinstance(points,list): |

2182 | points= [points] |

2183 | for p in points: |

2184 | # Ensure variables in points lie in SR. |

2185 | pp= {} |

2186 | for x in p.keys(): |

2187 | pp[SR(x)]= p[x] |

2188 | # Test 1: Is p a zero of all factors of H? |

2189 | if [h.subs(pp) for h in Hs] != [Integer(0) for h in Hs]: |

2190 | # Failed test 1. Move on to next point. |

2191 | verdicts.append((p,False,'not a singular point')) |

2192 | continue |

2193 | # Test 2: Are the factors of H smooth and |

2194 | # do theyintersect transversely at p? |

2195 | J= [tuple([f.subs(pp) for f in dh]) for dh in J] |

2196 | l= len(J) |

2197 | if Set(J).cardinality() < l: |

2198 | # Fail. Move on to next point. |

2199 | verdicts.append((p,False,'not a transverse intersection')) |

2200 | continue |

2201 | temp= True |

2202 | for S in list(Set(J).subsets())[Integer(1) :]: # Subsets of size >= 1 |

2203 | k= len(list(S)) |

2204 | M = matrix(list(S)) |

2205 | if rank(M) != min(k,d): |

2206 | # Fail. |

2207 | temp= False |

2208 | verdicts.append((p,False,'not a transvere intersection')) |

2209 | break |

2210 | if not temp: |

2211 | # Move on to next point |

2212 | continue |

2213 | # Test 3: Is p convenient? |

2214 | Jlog= matrix([self._log_grad(h,X,pp) for h in Hs]) |

2215 | if [Integer(0) in f for f in Jlog.columns()] ==\ |

2216 | [True for f in Jlog.columns()]: |

2217 | # Fail. Move on to next point. |

2218 | verdict[p] = (False,'multiple point but not convenient') |

2219 | continue |

2220 | verdicts.append((p,True,'all good')) |

2221 | return verdicts |

2222 | # J ============================================================================ |

2223 | # K ============================================================================ |

2224 | # L ============================================================================ |

2225 | def _log_grad(self,f,X,c): |

2226 | r""" |

2227 | This function returns the logarithmic gradient of `f` with respect to the |

2228 | variables of `X` evalutated at `c`. |

2229 | |

2230 | Does not use `self`. |

2231 | |

2232 | INPUT: |

2233 | |

2234 | - ``f`` - An expression in the variables of `X`. |

2235 | - ``X`` - A list of variables. |

2236 | - ``c`` - A dictionary with keys `X` and values in a field `K`. |

2237 | |

2238 | OUTPUT: |

2239 | |

2240 | \code{[c[x] * diff(f,x).subs(c) for x in X]}. |

2241 | |

2242 | EXAMPLES:: |

2243 | |

2244 | sage: from sage.combinat.amgf import * |

2245 | sage: R.<x> = PolynomialRing(QQ) |

2246 | sage: G,H = 1,1 |

2247 | sage: F= QuasiRationalExpression(G,H) |

2248 | sage: X= var('x,y,z') |

2249 | sage: f= x*y*z^2 |

2250 | sage: c= {x:1,y:2,z:3} |

2251 | sage: f.gradient() |

2252 | (y*z^2, x*z^2, 2*x*y*z) |

2253 | sage: F._log_grad(f,X,c) |

2254 | (18, 18, 36) |

2255 | |

2256 | :: |

2257 | |

2258 | sage: R.<x,y,z>= PolynomialRing(QQ) |

2259 | sage: f= x*y*z^2 |

2260 | sage: c= {x:1,y:2,z:3} |

2261 | sage: F._log_grad(f,R.gens(),c) |

2262 | (18, 18, 36) |

2263 | |

2264 | AUTHORS: |

2265 | |

2266 | - Alex Raichev (2009-03-06) |

2267 | """ |

2268 | return tuple([SR(c[x] * diff(f,x).subs(c)).simplify() for x in X]) |

2269 | # M ============================================================================ |

2270 | def _make_var(self,L): |

2271 | r""" |

2272 | This function converts the list `L` to a string (without commas) and returns |

2273 | the string as a variable. |

2274 | For internal use by the function _diff_op() |

2275 | |

2276 | Does not use `self`. |

2277 | |

2278 | INPUT: |

2279 | |

2280 | - ``L`` - A list. |

2281 | |

2282 | OUTPUT: |

2283 | |

2284 | A variable whose name is the concatenation of the variable names in `L`. |

2285 | |

2286 | EXAMPLES:: |

2287 | |

2288 | sage: from sage.combinat.amgf import * |

2289 | sage: R.<x> = PolynomialRing(QQ) |

2290 | sage: G,H = 1,1 |

2291 | sage: F= QuasiRationalExpression(G,H) |

2292 | sage: L= list(var('x,y,hello')) |

2293 | sage: v= F._make_var(L) |

2294 | sage: print v, type(v) |

2295 | xyhello <type 'sage.symbolic.expression.Expression'> |

2296 | |

2297 | AUTHORS: |

2298 | |

2299 | - Alex Raichev (2010-01-21) |

2300 | """ |

2301 | return var(''.join([str(v) for v in L])) |

2302 | # N ============================================================================ |

2303 | def _new_var_name(self,name,V): |

2304 | r""" |

2305 | This function returns the first string in the sequence `name`, `name+name`, |

2306 | `name+name+name`,... that does not appear in the list `V`. |

2307 | It is for internal use by the function _asymptotics_main_multiple(). |

2308 | |

2309 | Does not use `self`. |

2310 | |

2311 | INPUT: |

2312 | |

2313 | - ``name`` - A string. |

2314 | - ``V`` - A list of variables. |

2315 | |

2316 | OUTPUT: |

2317 | |

2318 | The first string in the sequence `name`, `name+name`, |

2319 | `name+name+name`,... that does not appear in the list \code{str(V)}. |

2320 | |

2321 | EXAMPLES:: |

2322 | |

2323 | sage: from sage.combinat.amgf import * |

2324 | sage: R.<x> = PolynomialRing(QQ) |

2325 | sage: G,H = 1,1 |

2326 | sage: F= QuasiRationalExpression(G,H) |

2327 | sage: X= var('x,xx,y,z') |

2328 | sage: F._new_var_name('x',X) |

2329 | 'xxx' |

2330 | |

2331 | AUTHORS: |

2332 | |

2333 | - Alex Raichev (2008-10-01) |

2334 | """ |

2335 | newname= name |

2336 | while newname in str(V): |

2337 | newname= newname +name |

2338 | return newname |

2339 | # O ============================================================================ |

2340 | # P ============================================================================ |

2341 | # Q ============================================================================ |

2342 | # R ============================================================================ |

2343 | def relative_error(self,approx,alpha,interval,exp_scale=Integer(1) ): |

2344 | r""" |

2345 | Returns the relative error between the values of the Maclaurin |

2346 | coefficients of `self` with multi-indices `m alpha` for `m` in |

2347 | `interval` and the values of the functions in `approx`. |

2348 | |

2349 | INPUT: |

2350 | |

2351 | - ``approx`` - An individual or list of symbolic expressions in |

2352 | one variable. |

2353 | - ``alpha`` - A list of positive integers. |

2354 | - ``interval`` - A list of positive integers. |

2355 | - ``exp_scale`` - (optional) A number. Default: 1. |

2356 | |

2357 | OUTPUT: |

2358 | |

2359 | A list whose entries are of the form |

2360 | `[m\alpha,a_m,b_m,err_m]` for `m \in interval`. |

2361 | Here `m\alpha` is a tuple; `a_m` is the `m alpha` (multi-index) |

2362 | coefficient of the Maclaurin series for `F` divided by `exp_scale^m`; |

2363 | `b_m` is a list of the values of the functions in `approx` evaluated at |

2364 | `m` and divided by `exp_scale^m`; `err_m` is the list of relative errors |

2365 | `(a_m-f)/a_m` for `f` in `b_m`. |

2366 | All outputs are decimal approximations. |

2367 | |

2368 | EXAMPLES:: |

2369 | |

2370 | sage: from sage.combinat.amgf import * |

2371 | sage: R.<x,y>= PolynomialRing(QQ) |

2372 | sage: G=1 |

2373 | sage: H= 1-x-y-x*y |

2374 | sage: F= QuasiRationalExpression(G,H) |

2375 | sage: alpha= [1,1] |

2376 | sage: var('n') |

2377 | n |

2378 | sage: f= (0.573/sqrt(n))*5.83^n |

2379 | sage: es= 5.83 |

2380 | sage: F.relative_error(f,alpha,[1,2,4,8],es) |

2381 | Calculating errors table in the form |

2382 | exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... |

2383 | [[(1, 1), 0.514579759862779, [0.573000000000000], [-0.113530000000000]], [(2, 2), 0.382477808931739, [0.405172185619892], [-0.0593351461396876]], [(4, 4), 0.277863059517142, [0.286500000000000], [-0.0310834426780842]], [(8, 8), 0.199108827584423, [0.202586092809946], [-0.0174641439443390]]] |

2384 | sage: g= (0.573/sqrt(n) - 0.0674/n^(3/2))*5.83^n |

2385 | sage: F.relative_error([f,g],alpha,[1,2,4,8],es) |

2386 | Calculating errors table in the form |

2387 | exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... |

2388 | [[(1, 1), 0.514579759862779, [0.573000000000000, 0.505600000000000], [-0.113530000000000, 0.0174506666666667]], [(2, 2), 0.382477808931739, [0.405172185619892, 0.381342687093905], [-0.0593351461396876, 0.00296781097184384]], [(4, 4), 0.277863059517142, [0.286500000000000, 0.278075000000000], [-0.0310834426780842, -0.000762751562681505]], [(8, 8), 0.199108827584423, [0.202586092809946, 0.199607405494198], [-0.0174641439443390, -0.00250404723800224]]] |

2389 | |

2390 | AUTHORS: |

2391 | |

2392 | - Alex Raichev (2009-05-18, 2011-04-18) |

2393 | """ |

2394 | |

2395 | if not isinstance(approx,list): |

2396 | approx= [approx] |

2397 | av= approx[Integer(0) ].variables()[Integer(0) ] |

2398 | |

2399 | print "Calculating errors table in the form" |

2400 | print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..." |

2401 | |

2402 | # Get Maclaurin coefficients of self and scale them. |

2403 | # Then compute errors. |

2404 | n= interval[-Integer(1) ] |

2405 | mc= self.maclaurin_coefficients(alpha,n) |

2406 | mca={} |

2407 | stats=[] |

2408 | for m in interval: |

2409 | beta= tuple([m*a for a in alpha]) |

2410 | mc[beta]= mc[beta]/exp_scale**m |

2411 | mca[beta]= [f.subs({av:m})/exp_scale**m for f in approx] |

2412 | stats_row= [beta, mc[beta].n(), [a.n() for a in mca[beta]]] |

2413 | if mc[beta]==Integer(0) : |

2414 | stats_row.extend([None for a in mca[beta]]) |

2415 | else: |

2416 | stats_row.append([((mc[beta]-a)/mc[beta]).n() for a in mca[beta]]) |

2417 | stats.append(stats_row) |

2418 | return stats |

2419 | # S ============================================================================ |

2420 | def singular_points(self): |

2421 | r""" |

2422 | This function returns a Groebner basis ideal whose variety is the |

2423 | set of singular points of the algebraic variety |

2424 | `V= \{x\in\CC^d : H(x)=0\}`, where `H=sef._H`. |

2425 | |

2426 | INPUT: |

2427 | |

2428 | OUTPUT: |

2429 | |

2430 | A Groebner basis ideal whose variety is the set of singular points of |

2431 | the algebraic variety `V= \{x\in\CC^d : H(x)=0\}`. |

2432 | |

2433 | EXAMPLES:: |

2434 | |

2435 | sage: from sage.combinat.amgf import * |

2436 | sage: R.<x,y,z>= PolynomialRing(QQ) |

2437 | sage: G= 1 |

2438 | sage: H= (4-2*x-y-z)*(4-x-2*y-z) |

2439 | sage: F= QuasiRationalExpression(G,H) |

2440 | sage: F.singular_points() |

2441 | Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field |

2442 | |

2443 | AUTHORS: |

2444 | |

2445 | - Alex Raichev (2008-10-01, 2008-11-20, 2010-12-03, 2011-04-18) |

2446 | """ |

2447 | H= self._H |

2448 | R= self._R |

2449 | f= R.ideal(H).radical().gens()[Integer(0) ] # Compute the reduction of H. |

2450 | J= R.ideal([f] + f.gradient()) |

2451 | return R.ideal(J.groebner_basis()) |

2452 | #------------------------------------------------------------------------------- |

2453 | def smooth_critical(self,alpha): |

2454 | r""" |

2455 | This function returns a Groebner basis ideal whose variety is the set |

2456 | of smooth critical points of the algebraic variety |

2457 | `V= \{x\in\CC^d : H(x)=0\} for the direction `\alpha` where `H=self._H`. |

2458 | |

2459 | INPUT: |

2460 | |

2461 | - ``alpha`` - A `d`-tuple of positive integers and/or symbolic entries. |

2462 | |

2463 | OUTPUT: |

2464 | |

2465 | A Groebner basis ideal of smooth critical points of `V` for `\alpha`. |

2466 | |

2467 | EXAMPLES:: |

2468 | |

2469 | sage: R.<x,y> = PolynomialRing(QQ) |

2470 | sage: G=1 |

2471 | sage: H = (1-x-y-x*y)^2 |

2472 | sage: F= QuasiRationalExpression(G,H) |

2473 | sage: var('a1,a2') |

2474 | (a1, a2) |

2475 | sage: F.smooth_critical([a1,a2]) |

2476 | 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 |

2477 | |

2478 | NOTES: |

2479 | |

2480 | A point `c` of `V` is a __smooth critical point for `alpha`__ |

2481 | if the gradient of `f` at `c` is not identically zero and `\alpha` is in |

2482 | the span of the logarithmic gradient vector |

2483 | `(c[0] \partial_1 f(c)),\ldots,c[d-1] \partial_d f(c))`; see [RaWi2008a]_. |

2484 | |

2485 | REFERENCES: |

2486 | |

2487 | .. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

2488 | coefficients of multivariate generating functions: improvements |

2489 | for smooth points", Electronic Journal of Combinatorics, Vol. 15 |

2490 | (2008), R89. |

2491 | |

2492 | AUTHORS: |

2493 | |

2494 | - Alex Raichev (2008-10-01, 2008-11-20, 2009-03-09, 2010-12-02, 2011-04-18) |

2495 | """ |

2496 | H= self._H |

2497 | R= H.parent() |

2498 | B= R.base_ring() |

2499 | d= R.ngens() |

2500 | vars= R.gens() |

2501 | f= R.ideal(H).radical().gens()[Integer(0) ] # Compute the reduction of H. |

2502 | |

2503 | # Expand B by the variables of alpha if there are any. |

2504 | indets= [] |

2505 | indets_ind= [] |

2506 | for a in alpha: |

2507 | if not ((a in ZZ) and (a>Integer(0) )): |

2508 | try: |

2509 | CC(a) |

2510 | except: |

2511 | indets.append(var(a)) |

2512 | indets_ind.append(alpha.index(a)) |

2513 | else: |

2514 | print "The components of", alpha, \ |

2515 | "must be positive integers or symbolic variables." |

2516 | return |

2517 | indets= list(Set(indets)) # Delete duplicates in indets. |

2518 | if indets != []: |

2519 | BB= FractionField(PolynomialRing(B,tuple(indets))) |

2520 | S= R.change_ring(BB) |

2521 | vars= S.gens() |

2522 | # Coerce alpha into BB. |

2523 | for i in range(len(alpha)): |

2524 | alpha[i] = BB(alpha[i]) |

2525 | else: |

2526 | S= R |

2527 | |

2528 | # Find smooth, critical points for alpha. |

2529 | f= S(f) |

2530 | J= S.ideal([f] +[ alpha[d-Integer(1) ]*vars[i]*diff(f,vars[i]) \ |

2531 | -alpha[i]*vars[d-Integer(1) ]*diff(f,vars[d-Integer(1) ]) for i in range(d-Integer(1) )]) |

2532 | return S.ideal(J.groebner_basis()) |

2533 | #------------------------------------------------------------------------------- |

2534 | def _subs_all(self,f,sub,simplify=False): |

2535 | r""" |

2536 | This function returns the items of `f` substituted by the dictionaries of |

2537 | `sub` in order of their appearance in `sub`. |

2538 | |

2539 | Does not use `self`. |

2540 | |

2541 | INPUT: |

2542 | |

2543 | - ``f`` - An individual or list of symbolic expressions or dictionaries |

2544 | - ``sub`` - An individual or list of dictionaries. |

2545 | - ``simplify`` - Boolean (default: False). |

2546 | |

2547 | OUTPUT: |

2548 | |

2549 | The items of `f` substituted by the dictionaries of `sub` in order of |

2550 | their appearance in `sub`. The subs() command is used. |

2551 | If simplify is True, then simplify() is used after substitution. |

2552 | |

2553 | EXAMPLES:: |

2554 | |

2555 | sage: from sage.combinat.amgf import * |

2556 | sage: R.<x> = PolynomialRing(QQ) |

2557 | sage: G,H = 1,1 |

2558 | sage: F= QuasiRationalExpression(G,H) |

2559 | sage: var('x,y,z') |

2560 | (x, y, z) |

2561 | sage: a= {x:1} |

2562 | sage: b= {y:2} |

2563 | sage: c= {z:3} |

2564 | sage: F._subs_all(x+y+z,a) |

2565 | y + z + 1 |

2566 | sage: F._subs_all(x+y+z,[c,a]) |

2567 | y + 4 |

2568 | sage: F._subs_all([x+y+z,y^2],b) |

2569 | [x + z + 2, 4] |

2570 | sage: F._subs_all([x+y+z,y^2],[b,c]) |

2571 | [x + 5, 4] |

2572 | |

2573 | :: |

2574 | |

2575 | sage: var('x,y') |

2576 | (x, y) |

2577 | sage: a= {'foo':x^2+y^2, 'bar':x-y} |

2578 | sage: b= {x:1,y:2} |

2579 | sage: F._subs_all(a,b) |

2580 | {'foo': 5, 'bar': -1} |

2581 | |

2582 | AUTHORS: |

2583 | |

2584 | - Alex Raichev (2009-05-05) |

2585 | """ |

2586 | singleton= False |

2587 | if not isinstance(f,list): |

2588 | f= [f] |

2589 | singleton= True |

2590 | if not isinstance(sub,list): |

2591 | sub= [sub] |

2592 | g= [] |

2593 | for ff in f: |

2594 | for D in sub: |

2595 | if isinstance(ff,dict): |

2596 | ff= dict( [(k,ff[k].subs(D)) for k in ff.keys()] ) |

2597 | else: |

2598 | ff= ff.subs(D) |

2599 | g.append(ff) |

2600 | if singleton and simplify: |

2601 | if isinstance(g[Integer(0) ],dict): |

2602 | return g[Integer(0) ] |

2603 | else: |

2604 | return g[Integer(0) ].simplify() |

2605 | elif singleton and not simplify: |

2606 | return g[Integer(0) ] |

2607 | elif not singleton and simplify: |

2608 | G= [] |

2609 | for gg in g: |

2610 | if isinstance(gg,dict): |

2611 | G.append(gg) |

2612 | else: |

2613 | G.append(gg.simplify()) |

2614 | return G |

2615 | else: |

2616 | return g |

2617 | # T ============================================================================ |

2618 | def maclaurin_coefficients(self,alpha,n): |

2619 | r""" |

2620 | Returns the Maclaurin coefficients of self that have multi-indices |

2621 | `alpha`, `2*alpha`,...,`n*alpha`. |

2622 | |

2623 | INPUT: |

2624 | |

2625 | - ``n`` - positive integer |

2626 | - ``alpha`` - tuple of positive integers representing a multi-index. |

2627 | |

2628 | OUTPUT: |

2629 | |

2630 | A dictionary of the form (beta, beta Maclaurin coefficient of self). |

2631 | |

2632 | AUTHORS: |

2633 | |

2634 | - Alex Raichev (2011-04-08) |

2635 | """ |

2636 | # Do all computations in the Symbolic Ring. |

2637 | F= SR(self._G/self._H) |

2638 | v= F.variables() |

2639 | p= {} |

2640 | for x in v: |

2641 | p[x]=Integer(0) |

2642 | d= len(v) |

2643 | coeffs={} |

2644 | # Initialize sequence of variables to differentiate with. |

2645 | s=[] |

2646 | for i in range(d): |

2647 | s.extend([v[i] for j in range(alpha[i])]) |

2648 | F_deriv= diff(F,s) |

2649 | coeffs[tuple(alpha)]= F_deriv.subs(p)/ mul([factorial(a) for a in alpha]) |

2650 | old_beta= alpha |

2651 | for k in (ellipsis_range(Integer(2) ,Ellipsis,n)): |

2652 | # Update variable sequence to differentiate with. |

2653 | beta= [k*a for a in alpha] |

2654 | delta= [beta[i]-old_beta[i] for i in range(d)] |

2655 | s= [] |

2656 | for i in range(d): |

2657 | s.extend([v[i] for j in range(delta[i])]) |

2658 | F_deriv= diff(F_deriv,s) |

2659 | coeffs[tuple(beta)]= F_deriv.subs(p)/ mul([factorial(b) for b in beta]) |

2660 | old_beta= copy(beta) |

2661 | return coeffs |

2662 | # U ============================================================================ |

2663 | # V ============================================================================ |

2664 | def variables(self): |

2665 | r""" |

2666 | Returns the tuple of variables of `self`. |

2667 | |

2668 | EXAMPLES:: |

2669 | |

2670 | sage: from sage.combinat.amgf import * |

2671 | sage: R.<x,y>= PolynomialRing(QQ) |

2672 | sage: G= exp(x) |

2673 | sage: H= 1-y |

2674 | sage: F= QuasiRationalExpression(G,H) |

2675 | sage: F.variables() |

2676 | (x, y) |

2677 | |

2678 | sage: R.<x,y>= PolynomialRing(QQ,order='invlex') |

2679 | sage: G= exp(x) |

2680 | sage: H= 1-y |

2681 | sage: F= QuasiRationalExpression(G,H) |

2682 | sage: F.variables() |

2683 | (y, x) |

2684 | |

2685 | AUTHORS: |

2686 | |

2687 | - Alex Raichev (2011-04-01) |

2688 | """ |

2689 | return self._variables |

2690 | # W ============================================================================ |

2691 | # X ============================================================================ |

2692 | # Y ============================================================================ |

2693 | # Z ============================================================================ |