# Ticket #10519: trac_10519.patch

File trac_10519.patch, 111.5 KB (added by , 8 years ago) |
---|

Line | |
---|---|

1 | Detected SAGE64 flag |

2 | Building Sage on OS X in 64-bit mode |

3 | # HG changeset patch |

4 | # User Alex Raichev <tortoise.said@gmail.com> |

5 | # Date 1308279234 -43200 |

6 | # Node ID 14c95b266433e0c817c09748d5077a8f459de251 |

7 | # Parent ce324e28c3334398d3552640e2cb1520d22465a3 |

8 | # 10519: initial version |

9 | |

10 | r""" |

11 | This code relates to analytic combinatorics. |

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

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

14 | multivariate generating functions. |

15 | |

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

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

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

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

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

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

22 | smooth or multiple points. |

23 | |

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

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

26 | and [RaWi2011]_. |

27 | |

28 | REFERENCES: |

29 | |

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

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

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

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

34 | |

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

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

37 | |

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

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

40 | |

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

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

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

44 | |

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

46 | examples of asymptotics derived from multivariate generating |

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

48 | |

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

50 | coefficients of multivariate generating functions: improvements |

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

52 | (2008), R89. |

53 | |

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

55 | coefficients of multivariate generating functions: improvements |

56 | for smooth points", To appear. |

57 | |

58 | AUTHORS: |

59 | |

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

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

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

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

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

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

66 | fixed _crit_cone_combo() to work in SR |

67 | |

68 | EXAMPLES:: |

69 | |

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

71 | A smooth point example. :: |

72 | |

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

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

75 | sage: G= 1 |

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

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

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

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

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

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

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

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

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

85 | sage: c= s[1] |

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

87 | Initializing auxiliary functions... |

88 | Calculating derivatives of auxiallary functions... |

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

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

91 | Initializing auxiliary functions... |

92 | Calculating derivatives of auxiallary functions... |

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

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

95 | [(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))] |

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

97 | Calculating errors table in the form |

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

99 | [[(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]]] |

100 | |

101 | Another smooth point example. |

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

103 | |

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

105 | sage: q= 1/2 |

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

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

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

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

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

111 | [4, 1] |

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

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

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

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

116 | sage: c= s[0] |

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

118 | Initializing auxiliary functions... |

119 | Calculating derivatives of auxiallary functions... |

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

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

122 | (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))) |

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

124 | Calculating errors table in the form |

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

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

127 | |

128 | A multiple point example. :: |

129 | |

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

131 | sage: G= 1 |

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

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

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

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

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

137 | sage: F.is_cmp(c) |

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

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

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

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

142 | sage: alpha in cc[1] |

143 | True |

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

145 | Initializing auxiliary functions... |

146 | Calculating derivatives of auxiliary functions... |

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

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

149 | (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))) |

150 | |

151 | Another multiple point example. :: |

152 | |

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

154 | sage: G= 16 |

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

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

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

158 | (a1, a2, a3) |

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

160 | [[-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]]] |

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

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

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

164 | sage: F.is_cmp(c) |

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

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

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

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

169 | sage: alpha in cc[1] |

170 | True |

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

172 | Initializing auxiliary functions... |

173 | Calculating derivatives of auxiliary functions... |

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

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

176 | (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))) |

177 | |

178 | """ |

179 | #***************************************************************************** |

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

181 | # |

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

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

184 | #***************************************************************************** |

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

186 | |

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

188 | def algebraic_dependence(fs): |

189 | r""" |

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

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

192 | Otherwise it returns 0. |

193 | |

194 | INPUT: |

195 | |

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

197 | ring. |

198 | |

199 | OUTPUT: |

200 | |

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

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

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

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

205 | indeterminates different from those of self. |

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

207 | is the zero polynomial. |

208 | |

209 | EXAMPLES:: |

210 | |

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

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

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

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

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

216 | sage: g(fs) |

217 | 0 |

218 | |

219 | :: |

220 | |

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

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

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

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

225 | sage: g(fs) |

226 | 0 |

227 | sage: g.parent() |

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

229 | |

230 | :: |

231 | |

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

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

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

235 | 0 |

236 | |

237 | AUTHORS: |

238 | |

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

240 | """ |

241 | r= len(fs) |

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

243 | B= R.base_ring() |

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

245 | d= len(Xs) |

246 | |

247 | # Expand R by r new variables. |

248 | T= 'T' |

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

250 | T= T+'T' |

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

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

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

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

255 | Ts= Vs[d:] |

256 | |

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

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

259 | JJ= J.elimination_ideal(Xs) |

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

261 | |

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

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

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

265 | return RRR(g) |

266 | #------------------------------------------------------------------------------- |

267 | def partial_fraction_decomposition(f): |

268 | r""" |

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

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

271 | |

272 | INPUT: |

273 | |

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

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

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

277 | |

278 | OUTPUT: |

279 | |

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

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

282 | See the notes below for more details. |

283 | |

284 | EXAMPLES:: |

285 | |

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

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

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

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

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

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

292 | True |

293 | |

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

295 | |

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

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

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

299 | sage: whole, parts |

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

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

302 | True |

303 | |

304 | :: |

305 | |

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

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

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

309 | sage: whole, parts |

310 | (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)]) |

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

312 | True |

313 | |

314 | :: |

315 | |

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

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

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

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

320 | True |

321 | |

322 | NOTES: |

323 | |

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

325 | partial fraction decomposition function. |

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

327 | partial fraction decompositions. |

328 | |

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

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

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

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

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

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

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

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

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

338 | coefficient field of `R`). |

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

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

341 | |

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

343 | [Lein1978]_. |

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

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

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

347 | |

348 | REFERENCES: |

349 | |

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

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

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

353 | |

354 | AUTHORS: |

355 | |

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

357 | """ |

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

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

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

361 | return f.partial_fraction_decomposition() |

362 | Q= f.denominator() |

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

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

365 | # Decompose via nullstellensatz trick |

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

367 | parts= decompose_via_nullstellensatz(parts) |

368 | # Decompose via algebraic dependence trick |

369 | parts= decompose_via_algebraic_dependence(parts) |

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

371 | new_parts=[] |

372 | for p in parts: |

373 | f= unformat_quotient(p) |

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

375 | whole= whole +f |

376 | else: |

377 | new_parts.append(f) |

378 | # Put P back in |

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

380 | return whole, new_parts |

381 | #------------------------------------------------------------------------------- |

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

383 | r""" |

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

385 | For use by partial_fraction_decomposition() above. |

386 | |

387 | EXAMPLES: |

388 | |

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

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

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

392 | sage: combine_like_terms(parts) |

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

394 | |

395 | :: |

396 | |

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

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

399 | sage: combine_like_terms(parts) |

400 | [[1]] |

401 | |

402 | AUTHORS: |

403 | |

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

405 | """ |

406 | |

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

408 | # Sort parts by denominators |

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

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

411 | newnew_parts=[] |

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

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

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

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

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

417 | else: |

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

419 | left= right |

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

421 | right= right +Integer(1) |

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

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

424 | if rationomial: |

425 | # Reduce fractions in newnew_parts. |

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

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

428 | # calls factor(). |

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

430 | return newnew_parts |

431 | #------------------------------------------------------------------------------- |

432 | def decompose_via_algebraic_dependence(parts): |

433 | r""" |

434 | Returns a decomposition of parts. |

435 | Used by partial_fraction_decomposition() above. |

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

437 | Recursive. |

438 | |

439 | REFERENCES: |

440 | |

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

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

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

444 | |

445 | AUTHORS: |

446 | |

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

448 | """ |

449 | decomposing_done= True |

450 | new_parts= [] |

451 | for p in parts: |

452 | p_parts= [p] |

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

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

455 | m= len(Qs) |

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

457 | if G: |

458 | # Then the denominator factors are algebraically dependent |

459 | # and so we can decompose p. |

460 | decomposing_done= False |

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

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

463 | |

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

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

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

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

468 | |

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

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

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

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

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

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

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

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

477 | |

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

479 | # reduce terms and combine like terms. |

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

481 | p_parts_temp= combine_like_terms(p_parts_temp) |

482 | |

483 | # Substitute Qs into new_p. |

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

485 | p_parts=[] |

486 | for x in p_parts_temp: |

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

488 | yy=[] |

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

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

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

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

493 | else: |

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

495 | yy.append(xx) |

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

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

498 | new_parts.extend(p_parts) |

499 | if decomposing_done: |

500 | return new_parts |

501 | else: |

502 | return decompose_via_algebraic_dependence(new_parts) |

503 | #------------------------------------------------------------------------------- |

504 | def decompose_via_cohomology(parts): |

505 | r""" |

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

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

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

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

510 | |

511 | INPUT: |

512 | |

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

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

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

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

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

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

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

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

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

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

523 | |

524 | OUTPUT: |

525 | |

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

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

528 | |

529 | EXAMPLES:: |

530 | |

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

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

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

534 | [] |

535 | |

536 | :: |

537 | |

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

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

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

541 | |

542 | NOTES: |

543 | |

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

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

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

547 | [AiYu1983]_. |

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

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

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

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

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

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

554 | `c`. |

555 | |

556 | REFERENCES: |

557 | |

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

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

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

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

562 | |

563 | AUTHORS: |

564 | |

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

566 | """ |

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

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

569 | decomposing_done= True |

570 | new_parts= [] |

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

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

573 | for p in parts: |

574 | p_parts= [p] |

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

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

577 | m= len(Qs) |

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

579 | # Then we can decompose p |

580 | p_parts= [] |

581 | decomposing_done= False |

582 | dets= [] |

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

584 | for v in vars_subsets: |

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

586 | # declared is first in vars list. |

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

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

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

590 | # Get a Nullstellensatz certificate. |

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

592 | for j in range(m): |

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

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

595 | new_Qs = copy.deepcopy(Qs) |

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

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

598 | else: |

599 | del new_Qs[j] |

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

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

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

603 | new_Qs = copy.deepcopy(Qs) |

604 | for j in range(m): |

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

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

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

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

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

610 | det= jac.determinant() |

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

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

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

614 | break |

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

616 | new_parts.extend(p_parts) |

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

618 | if decomposing_done: |

619 | return new_parts |

620 | else: |

621 | return decompose_via_cohomology(new_parts) |

622 | #------------------------------------------------------------------------------- |

623 | def decompose_via_nullstellensatz(parts): |

624 | r""" |

625 | Returns a decomposition of parts. |

626 | Used by partial_fraction_decomposition() above. |

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

628 | Recursive. |

629 | |

630 | REFERENCES: |

631 | |

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

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

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

635 | |

636 | AUTHORS: |

637 | |

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

639 | """ |

640 | decomposing_done= True |

641 | new_parts= [] |

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

643 | for p in parts: |

644 | p_parts= [p] |

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

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

647 | m= len(Qs) |

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

649 | # Then we can decompose p. |

650 | decomposing_done= False |

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

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

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

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

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

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

657 | new_parts.extend(p_parts) |

658 | if decomposing_done: |

659 | return new_parts |

660 | else: |

661 | return decompose_via_nullstellensatz(new_parts) |

662 | #------------------------------------------------------------------------------- |

663 | def format_quotient(f): |

664 | r""" |

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

666 | |

667 | AUTHORS: |

668 | |

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

670 | """ |

671 | P= f.numerator() |

672 | Q= f.denominator() |

673 | Qs= Q.factor() |

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

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

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

677 | return [P]+Qs |

678 | #------------------------------------------------------------------------------- |

679 | def permutation_sign(v,vars): |

680 | r""" |

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

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

683 | For internal use by _cohom_equiv_main(). |

684 | |

685 | INPUT: |

686 | |

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

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

689 | |

690 | OUTPUT: |

691 | |

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

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

694 | elements of `v` removed. |

695 | |

696 | EXAMPLES:: |

697 | |

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

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

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

701 | sage: permutation_sign(v,vars) |

702 | -1 |

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

704 | sage: permutation_sign(v,vars) |

705 | 1 |

706 | |

707 | AUTHORS: |

708 | |

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

710 | """ |

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

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

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

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

715 | C.sort() |

716 | P= Permutation(B+C) |

717 | return P.signature() |

718 | #------------------------------------------------------------------------------- |

719 | def unformat_quotient(part): |

720 | r""" |

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

722 | Inverse of format_quotient() above. |

723 | |

724 | AUTHORS: |

725 | |

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

727 | """ |

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

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

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

731 | return P/Q |

732 | #=============================================================================== |

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

734 | class QuasiRationalExpression(object): |

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

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

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

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

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

740 | self._G = G |

741 | self._H = H |

742 | R= H.parent() |

743 | self._R = R |

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

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

746 | |

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

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

749 | try: |

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

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

752 | except: |

753 | Gv= Set([]) |

754 | try: |

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

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

757 | except: |

758 | Hv= Set([]) |

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

760 | # depends upon. |

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

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

763 | #------------------------------------------------------------------------------- |

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

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

766 | r""" |

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

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

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

770 | where `F = self`. |

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

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

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

774 | |

775 | INPUT: |

776 | |

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

778 | point, possibly of symbolic variables. |

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

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

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

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

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

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

785 | Otherwise exact values will be returned. |

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

787 | The variable of the asymptotic expansion. |

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

789 | |

790 | OUTPUT: |

791 | |

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

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

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

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

796 | `asy`. |

797 | |

798 | EXAMPLES:: |

799 | |

800 | A smooth point example :: |

801 | |

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

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

804 | sage: G= 1 |

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

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

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

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

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

810 | Initializing auxiliary functions... |

811 | Calculating derivatives of auxiallary functions... |

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

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

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

815 | |

816 | A multiple point example :: |

817 | |

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

819 | sage: G= 1 |

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

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

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

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

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

825 | Initializing auxiliary functions... |

826 | Calculating derivatives of auxiliary functions... |

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

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

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

830 | |

831 | NOTES: |

832 | |

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

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

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

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

837 | see the documentation for _asymptotics_main_smooth() and |

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

839 | work. |

840 | |

841 | ALGORITHM: |

842 | |

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

844 | |

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

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

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

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

849 | |

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

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

852 | in the denominator. |

853 | |

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

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

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

857 | |

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

859 | |

860 | (5) Add the expansions. |

861 | |

862 | REFERENCES: |

863 | |

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

865 | coefficients of multivariate generating functions: improvements |

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

867 | (2008), R89. |

868 | |

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

870 | coefficients of multivariate generating functions: improvements |

871 | for smooth points", To appear. |

872 | |

873 | AUTHORS: |

874 | |

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

876 | """ |

877 | # The variable for asymptotic expansions. |

878 | if not asy_var: |

879 | asy_var= var('n') |

880 | |

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

882 | R= self._R |

883 | d= self._d |

884 | X= list(self._variables) |

885 | |

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

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

888 | |

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

890 | # calculations are no longer needed. |

891 | # Calculate asymptotics. |

892 | cc={} |

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

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

895 | c= cc |

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

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

898 | subexp_parts= [] |

899 | for chunk in new_integrands: |

900 | # Convert chunk into Symbolic Ring |

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

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

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

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

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

906 | |

907 | # Do step (5). |

908 | subexp_part= add(subexp_parts) |

909 | return exp_scale**asy_var *subexp_part, exp_scale, subexp_part |

910 | #-------------------------------------------------------------------------------- |

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

912 | r""" |

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

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

915 | whether to call _asymptotics_main_smooth() or |

916 | _asymptotics_main_multiple(). |

917 | |

918 | Does not use `self`. |

919 | |

920 | INPUT: |

921 | |

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

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

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

925 | Call its length `d`. |

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

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

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

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

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

931 | - ``numerical`` - Natural number. |

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

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

934 | Otherwise exact values will be returned. |

935 | |

936 | OUTPUT: |

937 | |

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

939 | |

940 | AUTHORS: |

941 | |

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

943 | """ |

944 | d= len(X) |

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

946 | |

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

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

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

950 | # convenient smooth or multiple points. |

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

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

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

954 | # in the same way. |

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

956 | beta= copy(alpha) |

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

958 | X.pop() |

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

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

961 | a= beta.pop() |

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

963 | if d==r: |

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

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

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

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

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

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

970 | else: |

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

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

973 | #-------------------------------------------------------------------------------- |

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

975 | r""" |

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

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

978 | |

979 | Does not use `self`. |

980 | |

981 | INPUT: |

982 | |

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

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

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

986 | Call its length `d`. |

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

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

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

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

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

992 | - ``numerical`` - Natural number. |

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

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

995 | Otherwise exact values will be returned. |

996 | |

997 | OUTPUT: |

998 | |

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

1000 | |

1001 | NOTES: |

1002 | |

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

1004 | Theorem 3.4 of [RaWi2011]_. |

1005 | |

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

1007 | because _crit_cone_combo() crashes. |

1008 | |

1009 | REFERENCES: |

1010 | |

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

1012 | coefficients of multivariate generating functions: improvements |

1013 | for smooth points", To appear. |

1014 | |

1015 | AUTHORS: |

1016 | |

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

1018 | """ |

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

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

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

1022 | C= copy(c) |

1023 | |

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

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

1026 | Sstar= {} |

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

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

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

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

1031 | thetastar.update(Sstar) |

1032 | |

1033 | # Setup. |

1034 | print "Initializing auxiliary functions..." |

1035 | Hmul= mul(H) |

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

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

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

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

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

1041 | 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)] |

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

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

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

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

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

1047 | 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) )]) |

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

1049 | atC = copy(C) |

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

1051 | |

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

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

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

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

1056 | .simplify() |

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

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

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

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

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

1062 | |

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

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

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

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

1067 | m= min(r,N) |

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

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

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

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

1072 | Uderivs={} |

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

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

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

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

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

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

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

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

1081 | all_zero= True |

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

1083 | # Check for a nonzero U derivative. |

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

1085 | all_zero= False |

1086 | if all_zero: |

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

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

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

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

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

1092 | else: |

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

1094 | 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) |

1095 | else: |

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

1097 | atC.update(Uderivs) |

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

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

1100 | a_inv= a.inverse() |

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

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

1103 | |

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

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

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

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

1108 | # Choose nondecreasing sequences as representative differentiation- |

1109 | # order sequences. |

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

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

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

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

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

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

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

1117 | AABB_derivs= At_derivs |

1118 | AABB_derivs.update(Phitu_derivs) |

1119 | for j in range(r): |

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

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

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

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

1124 | |

1125 | L={} |

1126 | 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) ))): |

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

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

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

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

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

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

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

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

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

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

1137 | 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 ]) \ |

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

1139 | |

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

1141 | coeffs= chunk.coefficients(n) |

1142 | coeffs.reverse() |

1143 | coeffs= coeffs[:N] |

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

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

1146 | for co in coeffs] ) |

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

1148 | .n(digits=numerical) |

1149 | else: |

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

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

1152 | return exp_scale**n*subexp_part, exp_scale, subexp_part |

1153 | #-------------------------------------------------------------------------------- |

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

1155 | r""" |

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

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

1158 | |

1159 | Does not use `self`. |

1160 | |

1161 | INPUT: |

1162 | |

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

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

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

1166 | Call its length `d`. |

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

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

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

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

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

1172 | - ``numerical`` - Natural number. |

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

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

1175 | Otherwise exact values will be returned. |

1176 | |

1177 | OUTPUT: |

1178 | |

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

1180 | |

1181 | NOTES: |

1182 | |

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

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

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

1186 | with `r=1`. |

1187 | |

1188 | REFERENCES: |

1189 | |

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

1191 | coefficients of multivariate generating functions: improvements |

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

1193 | (2008), R89. |

1194 | |

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

1196 | coefficients of multivariate generating functions: improvements |

1197 | for smooth points", To appear. |

1198 | |

1199 | AUTHORS: |

1200 | |

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

1202 | """ |

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

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

1205 | |

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

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

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

1209 | C= c |

1210 | else: |

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

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

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

1214 | |

1215 | # Setup. |

1216 | print "Initializing auxiliary functions..." |

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

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

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

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

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

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

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

1224 | ht= h.subs(e) |

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

1226 | At= A.subs(e) |

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

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

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

1230 | atC = copy(C) |

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

1232 | # evaluated at C. |

1233 | |

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

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

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

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

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

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

1240 | .simplify() |

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

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

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

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

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

1246 | |

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

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

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

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

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

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

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

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

1255 | # Uderivs. |

1256 | Uderivs={} |

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

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

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

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

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

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

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

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

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

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

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

1268 | all_zero= True |

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

1270 | # Check for a nonzero U derivative. |

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

1272 | all_zero= False |

1273 | if all_zero: |

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

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

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

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

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

1279 | else: |

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

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

1282 | else: |

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

1284 | atC.update(Uderivs) |

1285 | |

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

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

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

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

1290 | v=Integer(2) |

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

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

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

1294 | v= v+Integer(1) |

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

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

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

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

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

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

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

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

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

1304 | a= splat/factorial(v) |

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

1306 | |

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

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

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

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

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

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

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

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

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

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

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

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

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

1320 | else: |

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

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

1323 | AABB_derivs= At_derivs |

1324 | AABB_derivs.update(Phitu_derivs) |

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

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

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

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

1329 | # Plug above into asymptotic formula. |

1330 | L = [] |

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

1332 | for k in range(N): |

1333 | L.append( add([ \ |

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

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

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

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

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

1339 | else: |

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

1341 | for k in range(N): |

1342 | L.append( add([ \ |

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

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

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

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

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

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

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

1350 | # in this case. |

1351 | else: |

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

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

1354 | a_inv= a.inverse() |

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

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

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

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

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

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

1361 | # Choose nondecreasing sequences as representative differentiation- |

1362 | # order sequences. |

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

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

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

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

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

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

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

1370 | AABB_derivs= At_derivs |

1371 | AABB_derivs.update(Phitu_derivs) |

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

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

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

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

1376 | |

1377 | # Plug above into asymptotic formula. |

1378 | L=[] |

1379 | for k in range(N): |

1380 | L.append( add([ \ |

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

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

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

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

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

1386 | |

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

1388 | coeffs= chunk.coefficients(n) |

1389 | coeffs.reverse() |

1390 | coeffs= coeffs[:N] |

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

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

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

1394 | .n(digits=numerical) |

1395 | else: |

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

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

1398 | return exp_scale**n*subexp_part, exp_scale, subexp_part |

1399 | #----------------------------------------------------------------------------- |

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

1401 | r""" |

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

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

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

1405 | |

1406 | INPUT: |

1407 | |

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

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

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

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

1412 | |

1413 | OUTPUT: |

1414 | |

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

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

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

1418 | is the `j`th item in `fs` and where `a` is _direction(alpha). |

1419 | |

1420 | EXAMPLES:: |

1421 | |

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

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

1424 | sage: G,H= 1,1 |

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

1426 | sage: fs= [x + 2*y + z - 4, 2*x + y + z - 4] |

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

1428 | sage: alpha= [2,1,1] |

1429 | sage: F._crit_cone_combo(fs,R.gens(),c,alpha) |

1430 | [0, 1] |

1431 | |

1432 | NOTES: |

1433 | |

1434 | Use this function only when `Gamma` is well-defined and |

1435 | there is a unique solution to the matrix equation `y Gamma = a`. |

1436 | Fails otherwise. |

1437 | |

1438 | AUTHORS: |

1439 | |

1440 | - Alex Raichev (2008-10-01, 2008-11-25, 2009-03-04, 2010-09-08, |

1441 | 2010-12-02) |

1442 | """ |

1443 | # Assuming here that each _log_grad(f) has nonzero final component. |

1444 | # Then 'direction' will not throw a division by zero error. |

1445 | d= len(X) |

1446 | r= len(fs) |

1447 | Gamma= matrix([self._direction(self._log_grad(f,X,c)) for f in fs]) |

1448 | # solve_left() fails when working in SR :-(. So use solve() instead. |

1449 | #s= Gamma.solve_left(vector(alpha)/alpha[d-1]) |

1450 | V= [var('sss'+str(i)) for i in range(r)] |

1451 | M= matrix(V)*Gamma |

1452 | eqns= [M[Integer(0) ][i]== alpha[i]/alpha[d-Integer(1) ] for i in range(d)] |

1453 | s= solve(eqns,V,solution_dict=True)[Integer(0) ] # Assuming a unique solution. |

1454 | return [s[v] for v in V] |

1455 | # B ========================================================================== |

1456 | # C ========================================================================== |

1457 | def cohomologous_integrand(self,alpha,asy_var=None): |

1458 | r""" |

1459 | This function takes a multivariate Cauchy type integral |

1460 | `\int F / x^{asy_var\alpha+1} dx`, where `F=self`, and breaks it up |

1461 | into a list of nicer Cauchy type integrals for the purposes of |

1462 | computing asymptotics of the original integral as `asy_var\to\infty`. |

1463 | The sum of the nicer integrals is de Rham cohomologous to the original |

1464 | integral. |

1465 | It assumes that algebraic varieties corresponding to the irreducible |

1466 | factors of `self._H` intersect transversely (see notes below). |

1467 | |

1468 | INPUT: |

1469 | |

1470 | - ``alpha`` - A list of positive integers or symbolic variables. |

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

1472 | Eventually set to `var('n')` if None is given. |

1473 | |

1474 | OUTPUT: |

1475 | |

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

1477 | `chunk_j` has the form `[P,[B_1,1],\ldots,[B_l,1]]`. |

1478 | Here `l \le d`, `P` is a symbolic expression in the indeterminates of |

1479 | `R` and `asy_var`, `\{B_1,\ldots,B_l\} \subseteq \{Q_1,\ldots,Q_m\}`, |

1480 | and `[P,[B_1,1],\ldots,[B_l,1]]` represents the integral |

1481 | `\int P/(B_1 \cdots B_l) dx`. |

1482 | |

1483 | EXAMPLES:: |

1484 | |

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

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

1487 | sage: G= 9*exp(x+y) |

1488 | sage: H= (3-2*x-y)*(3-x-2*y) |

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

1490 | sage: alpha= [4,3] |

1491 | sage: F.cohomologous_integrand(alpha) |

1492 | [[9*e^(x + y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]] |

1493 | |

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

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

1496 | sage: G= 9*exp(x+y) |

1497 | sage: H= (3-2*x-y)^2*(3-x-2*y) |

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

1499 | sage: alpha= [4,3] |

1500 | sage: F.cohomologous_integrand(alpha) |

1501 | [[-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]]] |

1502 | |

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

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

1505 | sage: G= 16 |

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

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

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

1509 | sage: F.cohomologous_integrand(alpha) |

1510 | [[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]]] |

1511 | |

1512 | NOTES: |

1513 | |

1514 | The varieties corresponding to `Q_1,\ldots,Q_m` |

1515 | __intersect transversely__ iff for each point `c` of their intersection |

1516 | and for all `k \le l`, the Jacobian matrix of any `k` polynomials from |

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

1518 | `c`. |

1519 | |

1520 | ALGORITHM: |

1521 | |

1522 | Let `asy_var= n` and consider the integral around a polycirle centered |

1523 | at the origin of the `d`-variate differential form |

1524 | `\frac{G(x) dx_1 \wedge\cdots\wedge dx_d}{H(x) x^{n\alpha+1}}`, where |

1525 | `G=self._G` and `H=self._H`. |

1526 | |

1527 | (1) Decompose `G/H` into a sum of partial fractions `P_1 +\cdots+ P_r` |

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

1529 | in the denominator. |

1530 | |

1531 | (2) For each differential form `P_j dx_1 \wedge\cdots\wedge dx_d`, |

1532 | find an equivalent form `\omega_j` in de Rham cohomology with no |

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

1534 | |

1535 | AUTHOR: |

1536 | |

1537 | - Alex Raichev (2010-09-22) |

1538 | """ |

1539 | if not asy_var: |

1540 | asy_var = var('n') |

1541 | |

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

1543 | G= self._G |

1544 | H= self._H |

1545 | R= self._R |

1546 | d= self._d |

1547 | X= self._variables |

1548 | # Prepare input for partial_fraction_decomposition() which only works |

1549 | # for functions in the field of fractions of R. |

1550 | if G in R: |

1551 | numer=Integer(1) |

1552 | F= G/H |

1553 | else: |

1554 | numer= G |

1555 | F= R(Integer(1) )/H |

1556 | nstuff= Integer(1) /mul([X[j]**(alpha[j]*asy_var+Integer(1) ) for j in range(d)]) |

1557 | |

1558 | # Do steps (2) and (1). |

1559 | integrands= [] |

1560 | whole,parts= partial_fraction_decomposition(F) |

1561 | for f in parts: |

1562 | a= format_quotient(f) |

1563 | integrands.append( [a[Integer(0) ]*numer*nstuff] + a[Integer(1) :] ) |

1564 | |

1565 | # Do step (3). |

1566 | ce= decompose_via_cohomology(integrands) |

1567 | ce_new= [] |

1568 | for a in ce: |

1569 | ce_new.append( [(a[Integer(0) ]/nstuff).simplify_full().collect(n)] + a[Integer(1) :] ) |

1570 | return ce_new |

1571 | #------------------------------------------------------------------------------- |

1572 | def critical_cone(self,c,coordinate=None): |

1573 | r""" |

1574 | Returns the critical cone of a convenient multiple point `c`. |

1575 | |

1576 | INPUT: |

1577 | |

1578 | - ``c`` - A dictionary with keys `self.variables()` and values |

1579 | in a field. |

1580 | - ``coordinate`` - (optional) A natural number. |

1581 | |

1582 | OUTPUT: |

1583 | |

1584 | A list of vectors that generate the critical cone of `c` and |

1585 | the cone itself if the values of `c` lie in QQ. |

1586 | |

1587 | EXAMPLES:: |

1588 | |

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

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

1591 | sage: G= 1 |

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

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

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

1595 | sage: F.critical_cone(c) |

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

1597 | |

1598 | NOTES: |

1599 | |

1600 | The _critical cone_ of a convenient multiple point `c` with |

1601 | with `c_k \del_k H_j(c) \neq 0` for all `j=1,\ldots,r` is |

1602 | the conical hull of the vectors `\gamma_j(c) = |

1603 | \left(\frac{c_1 \del_1 H_j(c)}{c_k \del_k H_j(c)},\ldots, |

1604 | \frac{c_d \del_d H_j(c)}{c_k \del_k H_j(c)} \right)`. |

1605 | Here `H_1,\ldots,H_r` are the irreducible germs of `self._H` around `c`. |

1606 | For more details, see [RaWi2011]_. |

1607 | |

1608 | If this function's optional argument `coordinate` isn't given, then |

1609 | this function searches (from `d` down to 1) for the first index `k` |

1610 | such that for all `j=1,\ldots,r` we have `c_k \del_k H_j(c) \neq 0` |

1611 | and sets `coordinate = k`. |

1612 | Almost. |

1613 | Since this is Python, all the indices actually start at 0. |

1614 | |

1615 | REFERENCES: |

1616 | |

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

1618 | coefficients of multivariate generating functions: improvements |

1619 | for smooth points", To appear. |

1620 | |

1621 | AUTHORS: |

1622 | |

1623 | - Alex Raichev (2010-08-25) |

1624 | """ |

1625 | Hs= [SR(h[Integer(0) ]) for h in self._Hfac] # irreducible factors of H |

1626 | X= self._variables |

1627 | d= self._d |

1628 | # Ensure the variables of `c` lie in SR |

1629 | cc= {} |

1630 | for x in c.keys(): |

1631 | cc[SR(x)] = c[x] |

1632 | lg= [self._log_grad(h,X,cc) for h in Hs] |

1633 | if not coordinate: |

1634 | # Search (from d-1 down to 0) for a coordinate k such that |

1635 | # for all h in Hs we have cc[k] * diff(h,X[k]) !=0. |

1636 | # One is guaranteed to exist in the case of a convenient multiple |

1637 | # point. |

1638 | for k in reversed(range(d)): |

1639 | if Integer(0) not in [v[k] for v in lg]: |

1640 | coordinate= k |

1641 | break |

1642 | gamma= [self._direction(v,coordinate) for v in lg] |

1643 | if [[gg in QQ for gg in g] for g in gamma] == \ |

1644 | [[True for gg in g] for g in gamma]: |

1645 | return gamma,Cone(gamma) # Cone() needs rational vectors |

1646 | else: |

1647 | return gamma |

1648 | # D ============================================================================ |

1649 | def _diff_all(self,f,V,n,ending=[],sub=None,sub_final=None,zero_order=Integer(0) ,rekey=None): |

1650 | r""" |

1651 | This function returns a dictionary of representative mixed partial |

1652 | derivatives of `f` from order 1 up to order `n` with respect to the |

1653 | variables in `V`. |

1654 | The default is to key the dictionary by all nondecreasing sequences |

1655 | in `V` of length 1 up to length `n`. |

1656 | For internal use. |

1657 | |

1658 | Does not use `self`. |

1659 | |

1660 | INPUT: |

1661 | |

1662 | - ``f`` - An individual or list of `\mathcal{C}^{n+1}` functions. |

1663 | - ``V`` - A list of variables occurring in `f`. |

1664 | - ``n`` - A natural number. |

1665 | - ``ending`` - A list of variables in `V`. |

1666 | - ``sub`` - An individual or list of dictionaries. |

1667 | - ``sub_final`` - An individual or list of dictionaries. |

1668 | - ``rekey`` - A callable symbolic function in `V` or list thereof. |

1669 | - ``zero_order`` - A natural number. |

1670 | |

1671 | OUTPUT: |

1672 | |

1673 | The dictionary `{s_1:deriv_1,...,s_r:deriv_r}`. |

1674 | Here `s_1,\ldots,s_r` is a listing of |

1675 | all nondecreasing sequences of length 1 up to length `n` over the |

1676 | alphabet `V`, where `w > v` in `X` iff `str(w) > str(v)`, and |

1677 | `deriv_j` is the derivative of `f` with respect to the derivative |

1678 | sequence `s_j` and simplified with respect to the substitutions in `sub` |

1679 | and evaluated at `sub_final`. |

1680 | Moreover, all derivatives with respect to sequences of length less than |

1681 | `zero_order` (derivatives of order less than `zero_order`) will be made |

1682 | zero. |

1683 | |

1684 | If `rekey` is nonempty, then `s_1,\ldots,s_r` will be replaced by the |

1685 | symbolic derivatives of the functions in `rekey`. |

1686 | |

1687 | If `ending` is nonempty, then every derivative sequence `s_j` will be |

1688 | suffixed by `ending`. |

1689 | |

1690 | EXAMPLES:: |

1691 | |

1692 | I'd like to print the entire dictionaries, but that doesn't yield |

1693 | consistent output order for doctesting. |

1694 | Order of keys changes. :: |

1695 | |

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

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

1698 | sage: G,H = 1,1 |

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

1700 | sage: f= function('f',x) |

1701 | sage: dd= F._diff_all(f,[x],3) |

1702 | sage: dd[(x,x,x)] |

1703 | D[0, 0, 0](f)(x) |

1704 | |

1705 | :: |

1706 | |

1707 | sage: d1= {diff(f,x): 4*x^3} |

1708 | sage: dd= F._diff_all(f,[x],3,sub=d1) |

1709 | sage: dd[(x,x,x)] |

1710 | 24*x |

1711 | |

1712 | :: |

1713 | |

1714 | sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f) |

1715 | sage: dd[diff(f,x,3)] |

1716 | 24*x |

1717 | |

1718 | :: |

1719 | |

1720 | sage: a= {x:1} |

1721 | sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f,sub_final=a) |

1722 | sage: dd[diff(f,x,3)] |

1723 | 24 |

1724 | |

1725 | :: |

1726 | |

1727 | sage: X= var('x,y,z') |

1728 | sage: f= function('f',*X) |

1729 | sage: dd= F._diff_all(f,X,2,ending=[y,y,y]) |

1730 | sage: dd[(z,y,y,y)] |

1731 | D[1, 1, 1, 2](f)(x, y, z) |

1732 | |

1733 | :: |

1734 | |

1735 | sage: g= function('g',*X) |

1736 | sage: dd= F._diff_all([f,g],X,2) |

1737 | sage: dd[(0,y,z)] |

1738 | D[1, 2](f)(x, y, z) |

1739 | |

1740 | :: |

1741 | |

1742 | sage: dd[(1,z,z)] |

1743 | D[2, 2](g)(x, y, z) |

1744 | |

1745 | :: |

1746 | |

1747 | sage: f= exp(x*y*z) |

1748 | sage: ff= function('ff',*X) |

1749 | sage: dd= F._diff_all(f,X,2,rekey=ff) |

1750 | sage: dd[diff(ff,x,z)] |

1751 | x*y^2*z*e^(x*y*z) + y*e^(x*y*z) |

1752 | |

1753 | AUTHORS: |

1754 | |

1755 | - Alex Raichev (2008-10-01, 2009-04-01, 2010-02-01) |

1756 | """ |

1757 | singleton=False |

1758 | if not isinstance(f,list): |

1759 | f= [f] |

1760 | singleton=True |

1761 | |

1762 | # Build the dictionary of derivatives iteratively from a list of nondecreasing |

1763 | # derivative-order sequences. |

1764 | derivs= {} |

1765 | r= len(f) |

1766 | if ending: |

1767 | seeds = [ending] |

1768 | start = Integer(1) |

1769 | else: |

1770 | seeds = [[v] for v in V] |

1771 | start = Integer(2) |

1772 | if singleton: |

1773 | for s in seeds: |

1774 | derivs[tuple(s)] = self._subs_all(diff(f[Integer(0) ],s),sub) |

1775 | for l in (ellipsis_range(start,Ellipsis,n)): |

1776 | for t in UnorderedTuples(V,l): |

1777 | s= tuple(t + ending) |

1778 | derivs[s] = self._subs_all(diff(derivs[s[Integer(1) :]],s[Integer(0) ]),sub) |

1779 | else: |

1780 | # Make the dictionary keys of the form (j,sequence of variables), |

1781 | # where j in range(r). |

1782 | for s in seeds: |

1783 | value= self._subs_all([diff(f[j],s) for j in range(r)],sub) |

1784 | derivs.update(dict([(tuple([j]+s),value[j]) for j in range(r)])) |

1785 | for l in (ellipsis_range(start,Ellipsis,n)): |

1786 | for t in UnorderedTuples(V,l): |

1787 | s= tuple(t + ending) |

1788 | value= self._subs_all(\ |

1789 | [diff(derivs[(j,)+s[Integer(1) :]],s[Integer(0) ]) for j in range(r)],sub) |

1790 | derivs.update(dict([((j,)+s,value[j]) for j in range(r)])) |

1791 | if zero_order: |

1792 | # Zero out all the derivatives of order < zero_order |

1793 | if singleton: |

1794 | for k in derivs.keys(): |

1795 | if len(k) < zero_order: |

1796 | derivs[k]= Integer(0) |

1797 | else: |

1798 | # Ignore the first of element of k, which is an index. |

1799 | for k in derivs.keys(): |

1800 | if len(k)-Integer(1) < zero_order: |

1801 | derivs[k]= Integer(0) |

1802 | if sub_final: |

1803 | # Substitute sub_final into the values of derivs. |

1804 | for k in derivs.keys(): |

1805 | derivs[k] = self._subs_all(derivs[k],sub_final) |

1806 | if rekey: |

1807 | # Rekey the derivs dictionary by the value of rekey. |

1808 | F= rekey |

1809 | if singleton: |

1810 | # F must be a singleton. |

1811 | derivs= dict( [(diff(F,list(k)), derivs[k]) for k in derivs.keys()] ) |

1812 | else: |

1813 | # F must be a list. |

1814 | derivs= dict( [(diff(F[k[Integer(0) ]],list(k)[Integer(1) :]), derivs[k]) for k in derivs.keys()] ) |

1815 | return derivs |

1816 | #------------------------------------------------------------------------------- |

1817 | def _diff_op(self,A,B,AB_derivs,V,M,r,N): |

1818 | r""" |

1819 | This function computes the derivatives `DD^(l+k)(A[j] B^l)` evaluated at a |

1820 | point `p` for various natural numbers `j,k,l` which depend on `r` and `N`. |

1821 | Here `DD` is a specific second-order linear differential operator that depends |

1822 | on `M`, `A` is a list of symbolic functions, `B` is symbolic function, |

1823 | and `AB_derivs` contains all the derivatives of `A` and `B` evaluated at `p` |

1824 | that are necessary for the computation. |

1825 | For internal use by the functions _asymptotics_main_multiple() and |

1826 | _asymptotics_main_smooth(). |

1827 | |

1828 | Does not use `self`. |

1829 | |

1830 | INPUT: |

1831 | |

1832 | - ``A`` - A single or length `r` list of symbolic functions in the |

1833 | variables `V`. |

1834 | - ``B`` - A symbolic function in the variables `V`. |

1835 | - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of |

1836 | `A[0],\ldots,A[r-1]` up to order `2N-2` and |

1837 | the (symbolic) derivatives of `B` up to order `2N`. |

1838 | The values of the dictionary are complex numbers that are |

1839 | the keys evaluated at a common point `p`. |

1840 | - ``V`` - The variables of the `A[j]` and `B`. |

1841 | - ``M`` - A symmetric `l \times l` matrix, where `l` is the length of `V`. |

1842 | - ``r,N`` - Natural numbers. |

1843 | |

1844 | OUTPUT: |

1845 | |

1846 | A dictionary whose keys are natural number tuples of the form `(j,k,l)`, |

1847 | where `l \le 2k`, `j \le r-1`, and `j+k \le N-1`, and whose values are |

1848 | `DD^(l+k)(A[j] B^l)` evaluated at a point `p`, where `DD` is the linear |

1849 | second-order differential operator |

1850 | `-\sum_{i=0}^{l-1} \sum_{j=0}^{l-1} M[i][j] |

1851 | \partial^2 /(\partial V[j] \partial V[i])`. |

1852 | |

1853 | EXAMPLES:: |

1854 | |

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

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

1857 | sage: G,H = 1,1 |

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

1859 | sage: T= var('x,y') |

1860 | sage: A= function('A',*tuple(T)) |

1861 | sage: B= function('B',*tuple(T)) |

1862 | sage: AB_derivs= {} |

1863 | sage: M= matrix([[1,2],[2,1]]) |

1864 | sage: DD= F._diff_op(A,B,AB_derivs,T,M,1,2) |

1865 | sage: DD.keys() |

1866 | [(0, 1, 2), (0, 1, 1), (0, 1, 0), (0, 0, 0)] |

1867 | sage: len(DD[(0,1,2)]) |

1868 | 246 |

1869 | |

1870 | AUTHORS: |

1871 | |

1872 | - Alex Raichev (2008-10-01, 2010-01-12) |

1873 | """ |

1874 | if not isinstance(A,list): |

1875 | A= [A] |

1876 | |

1877 | # First, compute the necessary product derivatives of A and B. |

1878 | product_derivs= {} |

1879 | for (j,k) in mrange([r,N]): |

1880 | if j+k <N: |

1881 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)): |

1882 | for s in UnorderedTuples(V,Integer(2) *(k+l)): |

1883 | product_derivs[tuple([j,k,l]+s)] = \ |

1884 | diff(A[j]*B**l,s).subs(AB_derivs) |

1885 | |

1886 | # Second, compute DD^(k+l)(A[j]*B^l)(p) and store values in dictionary. |

1887 | DD= {} |

1888 | rows= M.nrows() |

1889 | for (j,k) in mrange([r,N]): |

1890 | if j+k <N: |

1891 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)): |

1892 | # Take advantage of the symmetry of M by ignoring |

1893 | # the upper-diagonal entries of M and multiplying by |

1894 | # appropriate powers of 2. |

1895 | if k+l == Integer(0) : |

1896 | DD[(j,k,l)] = product_derivs[(j,k,l)] |

1897 | continue |

1898 | S= [(a,b) for (a,b) in mrange([rows,rows]) if b <= a] |

1899 | P= cartesian_product_iterator([S for i in range(k+l)]) |

1900 | diffo= Integer(0) |

1901 | for t in P: |

1902 | if product_derivs[(j,k,l)+self._diff_seq(V,t)] != Integer(0) : |

1903 | MM= Integer(1) |

1904 | for (a,b) in t: |

1905 | MM= MM * M[a][b] |

1906 | if a != b: |

1907 | MM= Integer(2) *MM |

1908 | diffo= diffo + MM * product_derivs[(j,k,l)+self._diff_seq(V,t)] |

1909 | DD[(j,k,l)] = (-Integer(1) )**(k+l)*diffo |

1910 | return DD |

1911 | #------------------------------------------------------------------------------- |

1912 | def _diff_op_simple(self,A,B,AB_derivs,x,v,a,N): |

1913 | r""" |

1914 | This function computes `DD^(e k + v l)(A B^l)` evaluated at a point `p` |

1915 | for various natural numbers `e,k,l` that depend on `v` and `N`. |

1916 | Here `DD` is a specific linear differential operator that depends |

1917 | on `a` and `v`, `A` and `B` are symbolic functions, and `AB_derivs` contains |

1918 | all the derivatives of `A` and `B` evaluated at `p` that are necessary for |

1919 | the computation. |

1920 | For internal use by the function _asymptotics_main_smooth(). |

1921 | |

1922 | Does not use `self`. |

1923 | |

1924 | INPUT: |

1925 | |

1926 | - ``A``,``B`` - Symbolic functions in the variable `x`. |

1927 | - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of |

1928 | `A` up to order `2N` if `v` is even or `N` if `v` is odd and |

1929 | the (symbolic) derivatives of `B` up to order `2N+v` if `v` is even |

1930 | or `N+v` if `v` is odd. |

1931 | The values of the dictionary are complex numbers that are |

1932 | the keys evaluated at a common point `p`. |

1933 | - ``x`` - Symbolic variable. |

1934 | - ``a`` - A complex number. |

1935 | - ``v``,``N`` - Natural numbers. |

1936 | |

1937 | OUTPUT: |

1938 | |

1939 | A dictionary whose keys are natural number pairs of the form `(k,l)`, |

1940 | where `k < N` and `l \le 2k` and whose values are |

1941 | `DD^(e k + v l)(A B^l)` evaluated at a point `p`. |

1942 | Here `e=2` if `v` is even, `e=1` if `v` is odd, and `DD` is a particular |

1943 | linear differential operator |

1944 | `(a^{-1/v} d/dt)' if `v` is even and `(|a|^{-1/v} i \sgn a d/dt)` |

1945 | if `v` is odd. |

1946 | |

1947 | EXAMPLES:: |

1948 | |

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

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

1951 | sage: G,H = 1,1 |

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

1953 | sage: A= function('A',x) |

1954 | sage: B= function('B',x) |

1955 | sage: AB_derivs= {} |

1956 | sage: F._diff_op_simple(A,B,AB_derivs,x,3,2,2) |

1957 | {(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)} |

1958 | |

1959 | AUTHORS: |

1960 | |

1961 | - Alex Raichev (2010-01-15) |

1962 | """ |

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

1964 | DD= {} |

1965 | if v.mod(Integer(2) ) == Integer(0) : |

1966 | for k in range(N): |

1967 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)): |

1968 | DD[(k,l)] = (a**(-Integer(1) /v))**(Integer(2) *k+v*l) \ |

1969 | * diff(A*B**l,x,Integer(2) *k+v*l).subs(AB_derivs) |

1970 | else: |

1971 | for k in range(N): |

1972 | for l in (ellipsis_range(Integer(0) ,Ellipsis,k)): |

1973 | DD[(k,l)] = (abs(a)**(-Integer(1) /v)*I*a/abs(a))**(k+v*l) \ |

1974 | * diff(A*B**l,x,k+v*l).subs(AB_derivs) |

1975 | return DD |

1976 | #------------------------------------------------------------------------------- |

1977 | def _diff_prod(self,f_derivs,u,g,X,interval,end,uderivs,atc): |

1978 | r""" |

1979 | This function takes various derivatives of the equation `f=ug`, |

1980 | evaluates at a point `c`, and solves for the derivatives of `u`. |

1981 | For internal use by the function _asymptotics_main_multiple(). |

1982 | |

1983 | Does not use `self`. |

1984 | |

1985 | INPUT: |

1986 | |

1987 | - ``f_derivs`` - A dictionary whose keys are tuples \code{s + end} for all |

1988 | `X`-variable sequences `s` with length in `interval` and whose |

1989 | values are the derivatives of a function `f` evaluated at `c`. |

1990 | - ``u`` - A callable symbolic function. |

1991 | - ``g`` - An expression or callable symbolic function. |

1992 | - ``X`` - A list of symbolic variables. |

1993 | - ``interval`` - A list of positive integers. |

1994 | Call the first and last values `n` and `nn`, respectively. |

1995 | - ``end`` - A possibly empty list of `z`'s where `z` is the last element of |

1996 | `X`. |

1997 | - ``uderivs`` - A dictionary whose keys are the symbolic |

1998 | derivatives of order 0 to order `n-1` of `u` evaluated at `c` |

1999 | and whose values are the corresponding derivatives evaluated at `c`. |

2000 | - ``atc`` - A dictionary whose keys are the keys of `c` and all the symbolic |

2001 | derivatives of order 0 to order `nn` of `g` evaluated `c` and whose |

2002 | values are the corresponding derivatives evaluated at `c`. |

2003 | |

2004 | OUTPUT: |

2005 | |

2006 | A dictionary whose keys are the derivatives of `u` up to order |

2007 | `nn` and whose values are those derivatives evaluated at `c`. |

2008 | |

2009 | EXAMPLES:: |

2010 | |

2011 | I'd like to print out the entire dictionary, but that does not give |

2012 | consistent output for doctesting. |

2013 | Order of keys changes :: |

2014 | |

2015 | sage: from sage.combinat.amgf import * |

2016 | sage: R.<x> = PolynomialRing(QQ) |

2017 | sage: G,H = 1,1 |

2018 | sage: F= QuasiRationalExpression(G,H) |

2019 | sage: u= function('u',x) |

2020 | sage: g= function('g',x) |

2021 | 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}) |

2022 | sage: dd[diff(u,x,2)(x=2)] |

2023 | 22/9 |

2024 | |

2025 | NOTES: |

2026 | |

2027 | This function works by differentiating the equation `f=ug` with respect |

2028 | to the variable sequence \code{s+end}, for all tuples `s` of `X` of |

2029 | lengths in `interval`, evaluating at the point `c`, |

2030 | and solving for the remaining derivatives of `u`. |

2031 | This function assumes that `u` never appears in the differentiations of |

2032 | `f=ug` after evaluating at `c`. |

2033 | |

2034 | AUTHORS: |

2035 | |

2036 | - Alex Raichev (2009-05-14, 2010-01-21) |

2037 | """ |

2038 | for l in interval: |

2039 | D= {} |

2040 | rhs= [] |

2041 | lhs= [] |

2042 | for t in UnorderedTuples(X,l): |

2043 | s= t+end |

2044 | lhs.append(f_derivs[tuple(s)]) |

2045 | rhs.append(diff(u*g,s).subs(atc).subs(uderivs)) |

2046 | # Since Sage's solve command can't take derivatives as variable |

2047 | # names, i make new variables based on t to stand in for |

2048 | # diff(u,t) and store them in D. |

2049 | D[diff(u,t).subs(atc)] = self._make_var([var('zing')]+t) |

2050 | eqns=[ lhs[i] == rhs[i].subs(uderivs).subs(D) for i in range(len(lhs))] |

2051 | vars= D.values() |

2052 | sol= solve(eqns,*vars,solution_dict=True) |

2053 | uderivs.update(self._subs_all(D,sol[Integer(0) ])) |

2054 | return uderivs |

2055 | #------------------------------------------------------------------------------- |

2056 | def _diff_seq(self,V,s): |

2057 | r""" |

2058 | Given a list `s` of tuples of natural numbers, this function returns the |

2059 | list of elements of `V` with indices the elements of the elements of `s`. |

2060 | This function is for internal use by the function _diff_op(). |

2061 | |

2062 | Does not use `self`. |

2063 | |

2064 | INPUT: |

2065 | |

2066 | - ``V`` - A list. |

2067 | - ``s`` - A list of tuples of natural numbers in the interval |

2068 | \code{range(len(V))}. |

2069 | |

2070 | OUTPUT: |

2071 | |

2072 | The tuple \code{tuple([V[tt] for tt in sorted(t)])}, where `t` is the |

2073 | list of elements of the elements of `s`. |

2074 | |

2075 | EXAMPLES:: |

2076 | |

2077 | sage: from sage.combinat.amgf import * |

2078 | sage: R.<x> = PolynomialRing(QQ) |

2079 | sage: G,H = 1,1 |

2080 | sage: F= QuasiRationalExpression(G,H) |

2081 | sage: V= list(var('x,t,z')) |

2082 | sage: F._diff_seq(V,([0,1],[0,2,1],[0,0])) |

2083 | (x, x, x, x, t, t, z) |

2084 | |

2085 | AUTHORS: |

2086 | |

2087 | - Alex Raichev (2009.05.19) |

2088 | """ |

2089 | t= [] |

2090 | for ss in s: |

2091 | t.extend(ss) |

2092 | return tuple([V[tt] for tt in sorted(t)]) |

2093 | #------------------------------------------------------------------------------- |

2094 | def _direction(self,v,coordinate=None): |

2095 | r""" |

2096 | This function returns \code{[vv/v[coordinate] for vv in v]} where |

2097 | coordinate is the last index of v if not specified otherwise. |

2098 | |

2099 | Does not use `self`. |

2100 | |

2101 | INPUT: |

2102 | |

2103 | - ``v`` - A vector. |

2104 | - ``coordinate`` - An index for `v` (default: None). |

2105 | |

2106 | OUTPUT: |

2107 | |

2108 | This function returns \code{[vv/v[coordinate] for vv in v]} where |

2109 | coordinate is the last index of v if not specified otherwise. |

2110 | |

2111 | EXAMPLES:: |

2112 | |

2113 | sage: from sage.combinat.amgf import * |

2114 | sage: R.<x> = PolynomialRing(QQ) |

2115 | sage: G,H = 1,1 |

2116 | sage: F= QuasiRationalExpression(G,H) |

2117 | sage: F._direction([2,3,5]) |

2118 | (2/5, 3/5, 1) |

2119 | sage: F._direction([2,3,5],0) |

2120 | (1, 3/2, 5/2) |

2121 | |

2122 | AUTHORS: |

2123 | |

2124 | - Alex Raichev (2010-08-25) |

2125 | """ |

2126 | if coordinate == None: |

2127 | coordinate= len(v)-Integer(1) |

2128 | try: |

2129 | v[Integer(0) ].variables() |

2130 | return tuple([(vv/v[coordinate]).simplify_full() for vv in v]) |

2131 | except: |

2132 | return tuple([vv/v[coordinate] for vv in v]) |

2133 | # E ============================================================================ |

2134 | # F ============================================================================ |

2135 | # G ============================================================================ |

2136 | # H ============================================================================ |

2137 | # I ============================================================================ |

2138 | def is_cmp(self,points): |

2139 | r""" |

2140 | Checks if the points in the list `points` are convenient multiple |

2141 | points of `V= \{ x\in CC^d : H(x) = 0\}`, where `H=self._H`. |

2142 | |

2143 | INPUT: |

2144 | |

2145 | - ``points`` - An individual or list of dictionaries with keys |

2146 | `self._variables` and values in some superfield of |

2147 | `self._R.base_ring()`. |

2148 | |

2149 | OUTPUT: |

2150 | |

2151 | A list of tuples `(p,verdict,comment)`, one for each point |

2152 | `p` in `points`, where `verdict` is True if `p` is a convenient |

2153 | multiple point and False otherwise, and where `comment` is a string |

2154 | comment relating to `verdict`, such as 'not a transverse intersection'. |

2155 | |

2156 | EXAMPLES:: |

2157 | |

2158 | sage: from sage.combinat.amgf import * |

2159 | sage: R.<x,y,z>= PolynomialRing(QQ) |

2160 | sage: G= 16 |

2161 | sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) |

2162 | sage: F= QuasiRationalExpression(G,H) |

2163 | sage: points= [{x:1/2,y:1,z:4/3},{x:1/2,y:1,z:2}] |

2164 | sage: F.is_cmp(points) |

2165 | [({y: 1, z: 4/3, x: 1/2}, True, 'all good'), ({y: 1, z: 2, x: 1/2}, False, 'not a singular point')] |

2166 | |

2167 | NOTES: |

2168 | |

2169 | A point `c` of `V` is a __convenient multiple point__ if `V` is locally |

2170 | a union of complex manifolds that intersect transversely at `c`; |

2171 | see [RaWi2011]_. |

2172 | |

2173 | REFERENCES: |

2174 | |

2175 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

2176 | coefficients of multivariate generating functions: improvements |

2177 | for smooth points", To appear. |

2178 | |

2179 | AUTHORS: |

2180 | |

2181 | - Alex Raichev (2011-04-18) |

2182 | """ |

2183 | H= self._H |

2184 | d= self._d |

2185 | Hs= [SR(h[Integer(0) ]) for h in self._Hfac] # irreducible factors of H |

2186 | X= self._variables |

2187 | J= [tuple([diff(h,x) for x in X]) for h in Hs] |

2188 | verdicts= [] |

2189 | if not isinstance(points,list): |

2190 | points= [points] |

2191 | for p in points: |

2192 | # Ensure variables in points lie in SR. |

2193 | pp= {} |

2194 | for x in p.keys(): |

2195 | pp[SR(x)]= p[x] |

2196 | # Test 1: Is p a zero of all factors of H? |

2197 | if [h.subs(pp) for h in Hs] != [Integer(0) for h in Hs]: |

2198 | # Failed test 1. Move on to next point. |

2199 | verdicts.append((p,False,'not a singular point')) |

2200 | continue |

2201 | # Test 2: Are the factors of H smooth and |

2202 | # do theyintersect transversely at p? |

2203 | J= [tuple([f.subs(pp) for f in dh]) for dh in J] |

2204 | l= len(J) |

2205 | if Set(J).cardinality() < l: |

2206 | # Fail. Move on to next point. |

2207 | verdicts.append((p,False,'not a transverse intersection')) |

2208 | continue |

2209 | temp= True |

2210 | for S in list(Set(J).subsets())[Integer(1) :]: # Subsets of size >= 1 |

2211 | k= len(list(S)) |

2212 | M = matrix(list(S)) |

2213 | if rank(M) != min(k,d): |

2214 | # Fail. |

2215 | temp= False |

2216 | verdicts.append((p,False,'not a transvere intersection')) |

2217 | break |

2218 | if not temp: |

2219 | # Move on to next point |

2220 | continue |

2221 | # Test 3: Is p convenient? |

2222 | Jlog= matrix([self._log_grad(h,X,pp) for h in Hs]) |

2223 | if [Integer(0) in f for f in Jlog.columns()] ==\ |

2224 | [True for f in Jlog.columns()]: |

2225 | # Fail. Move on to next point. |

2226 | verdict[p] = (False,'multiple point but not convenient') |

2227 | continue |

2228 | verdicts.append((p,True,'all good')) |

2229 | return verdicts |

2230 | # J ============================================================================ |

2231 | # K ============================================================================ |

2232 | # L ============================================================================ |

2233 | def _log_grad(self,f,X,c): |

2234 | r""" |

2235 | This function returns the logarithmic gradient of `f` with respect to the |

2236 | variables of `X` evalutated at `c`. |

2237 | |

2238 | Does not use `self`. |

2239 | |

2240 | INPUT: |

2241 | |

2242 | - ``f`` - An expression in the variables of `X`. |

2243 | - ``X`` - A list of variables. |

2244 | - ``c`` - A dictionary with keys `X` and values in a field `K`. |

2245 | |

2246 | OUTPUT: |

2247 | |

2248 | \code{[c[x] * diff(f,x).subs(c) for x in X]}. |

2249 | |

2250 | EXAMPLES:: |

2251 | |

2252 | sage: from sage.combinat.amgf import * |

2253 | sage: R.<x> = PolynomialRing(QQ) |

2254 | sage: G,H = 1,1 |

2255 | sage: F= QuasiRationalExpression(G,H) |

2256 | sage: X= var('x,y,z') |

2257 | sage: f= x*y*z^2 |

2258 | sage: c= {x:1,y:2,z:3} |

2259 | sage: f.gradient() |

2260 | (y*z^2, x*z^2, 2*x*y*z) |

2261 | sage: F._log_grad(f,X,c) |

2262 | (18, 18, 36) |

2263 | |

2264 | :: |

2265 | |

2266 | sage: R.<x,y,z>= PolynomialRing(QQ) |

2267 | sage: f= x*y*z^2 |

2268 | sage: c= {x:1,y:2,z:3} |

2269 | sage: F._log_grad(f,R.gens(),c) |

2270 | (18, 18, 36) |

2271 | |

2272 | AUTHORS: |

2273 | |

2274 | - Alex Raichev (2009-03-06) |

2275 | """ |

2276 | return tuple([SR(c[x] * diff(f,x).subs(c)).simplify() for x in X]) |

2277 | # M ============================================================================ |

2278 | def _make_var(self,L): |

2279 | r""" |

2280 | This function converts the list `L` to a string (without commas) and returns |

2281 | the string as a variable. |

2282 | For internal use by the function _diff_op() |

2283 | |

2284 | Does not use `self`. |

2285 | |

2286 | INPUT: |

2287 | |

2288 | - ``L`` - A list. |

2289 | |

2290 | OUTPUT: |

2291 | |

2292 | A variable whose name is the concatenation of the variable names in `L`. |

2293 | |

2294 | EXAMPLES:: |

2295 | |

2296 | sage: from sage.combinat.amgf import * |

2297 | sage: R.<x> = PolynomialRing(QQ) |

2298 | sage: G,H = 1,1 |

2299 | sage: F= QuasiRationalExpression(G,H) |

2300 | sage: L= list(var('x,y,hello')) |

2301 | sage: v= F._make_var(L) |

2302 | sage: print v, type(v) |

2303 | xyhello <type 'sage.symbolic.expression.Expression'> |

2304 | |

2305 | AUTHORS: |

2306 | |

2307 | - Alex Raichev (2010-01-21) |

2308 | """ |

2309 | return var(''.join([str(v) for v in L])) |

2310 | # N ============================================================================ |

2311 | def _new_var_name(self,name,V): |

2312 | r""" |

2313 | This function returns the first string in the sequence `name`, `name+name`, |

2314 | `name+name+name`,... that does not appear in the list `V`. |

2315 | It is for internal use by the function _asymptotics_main_multiple(). |

2316 | |

2317 | Does not use `self`. |

2318 | |

2319 | INPUT: |

2320 | |

2321 | - ``name`` - A string. |

2322 | - ``V`` - A list of variables. |

2323 | |

2324 | OUTPUT: |

2325 | |

2326 | The first string in the sequence `name`, `name+name`, |

2327 | `name+name+name`,... that does not appear in the list \code{str(V)}. |

2328 | |

2329 | EXAMPLES:: |

2330 | |

2331 | sage: from sage.combinat.amgf import * |

2332 | sage: R.<x> = PolynomialRing(QQ) |

2333 | sage: G,H = 1,1 |

2334 | sage: F= QuasiRationalExpression(G,H) |

2335 | sage: X= var('x,xx,y,z') |

2336 | sage: F._new_var_name('x',X) |

2337 | 'xxx' |

2338 | |

2339 | AUTHORS: |

2340 | |

2341 | - Alex Raichev (2008-10-01) |

2342 | """ |

2343 | newname= name |

2344 | while newname in str(V): |

2345 | newname= newname +name |

2346 | return newname |

2347 | # O ============================================================================ |

2348 | # P ============================================================================ |

2349 | # Q ============================================================================ |

2350 | # R ============================================================================ |

2351 | def relative_error(self,approx,alpha,interval,exp_scale=Integer(1) ): |

2352 | r""" |

2353 | Returns the relative error between the values of the Maclaurin |

2354 | coefficients of `self` with multi-indices `m alpha` for `m` in |

2355 | `interval` and the values of the functions in `approx`. |

2356 | |

2357 | INPUT: |

2358 | |

2359 | - ``approx`` - An individual or list of symbolic expressions in |

2360 | one variable. |

2361 | - ``alpha`` - A list of positive integers. |

2362 | - ``interval`` - A list of positive integers. |

2363 | - ``exp_scale`` - (optional) A number. Default: 1. |

2364 | |

2365 | OUTPUT: |

2366 | |

2367 | A list whose entries are of the form |

2368 | `[m\alpha,a_m,b_m,err_m]` for `m \in interval`. |

2369 | Here `m\alpha` is a tuple; `a_m` is the `m alpha` (multi-index) |

2370 | coefficient of the Maclaurin series for `F` divided by `exp_scale^m`; |

2371 | `b_m` is a list of the values of the functions in `approx` evaluated at |

2372 | `m` and divided by `exp_scale^m`; `err_m` is the list of relative errors |

2373 | `(a_m-f)/a_m` for `f` in `b_m`. |

2374 | All outputs are decimal approximations. |

2375 | |

2376 | EXAMPLES:: |

2377 | |

2378 | sage: from sage.combinat.amgf import * |

2379 | sage: R.<x,y>= PolynomialRing(QQ) |

2380 | sage: G=1 |

2381 | sage: H= 1-x-y-x*y |

2382 | sage: F= QuasiRationalExpression(G,H) |

2383 | sage: alpha= [1,1] |

2384 | sage: var('n') |

2385 | n |

2386 | sage: f= (0.573/sqrt(n))*5.83^n |

2387 | sage: es= 5.83 |

2388 | sage: F.relative_error(f,alpha,[1,2,4,8],es) |

2389 | Calculating errors table in the form |

2390 | exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... |

2391 | [[(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]]] |

2392 | sage: g= (0.573/sqrt(n) - 0.0674/n^(3/2))*5.83^n |

2393 | sage: F.relative_error([f,g],alpha,[1,2,4,8],es) |

2394 | Calculating errors table in the form |

2395 | exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... |

2396 | [[(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]]] |

2397 | |

2398 | AUTHORS: |

2399 | |

2400 | - Alex Raichev (2009-05-18, 2011-04-18) |

2401 | """ |

2402 | |

2403 | if not isinstance(approx,list): |

2404 | approx= [approx] |

2405 | av= approx[Integer(0) ].variables()[Integer(0) ] |

2406 | |

2407 | print "Calculating errors table in the form" |

2408 | print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..." |

2409 | |

2410 | # Get Maclaurin coefficients of self and scale them. |

2411 | # Then compute errors. |

2412 | n= interval[-Integer(1) ] |

2413 | mc= self.maclaurin_coefficients(alpha,n) |

2414 | mca={} |

2415 | stats=[] |

2416 | for m in interval: |

2417 | beta= tuple([m*a for a in alpha]) |

2418 | mc[beta]= mc[beta]/exp_scale**m |

2419 | mca[beta]= [f.subs({av:m})/exp_scale**m for f in approx] |

2420 | stats_row= [beta, mc[beta].n(), [a.n() for a in mca[beta]]] |

2421 | if mc[beta]==Integer(0) : |

2422 | stats_row.extend([None for a in mca[beta]]) |

2423 | else: |

2424 | stats_row.append([((mc[beta]-a)/mc[beta]).n() for a in mca[beta]]) |

2425 | stats.append(stats_row) |

2426 | return stats |

2427 | # S ============================================================================ |

2428 | def singular_points(self): |

2429 | r""" |

2430 | This function returns a Groebner basis ideal whose variety is the |

2431 | set of singular points of the algebraic variety |

2432 | `V= \{x\in\CC^d : H(x)=0\}`, where `H=sef._H`. |

2433 | |

2434 | INPUT: |

2435 | |

2436 | OUTPUT: |

2437 | |

2438 | A Groebner basis ideal whose variety is the set of singular points of |

2439 | the algebraic variety `V= \{x\in\CC^d : H(x)=0\}`. |

2440 | |

2441 | EXAMPLES:: |

2442 | |

2443 | sage: from sage.combinat.amgf import * |

2444 | sage: R.<x,y,z>= PolynomialRing(QQ) |

2445 | sage: G= 1 |

2446 | sage: H= (4-2*x-y-z)*(4-x-2*y-z) |

2447 | sage: F= QuasiRationalExpression(G,H) |

2448 | sage: F.singular_points() |

2449 | Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field |

2450 | |

2451 | AUTHORS: |

2452 | |

2453 | - Alex Raichev (2008-10-01, 2008-11-20, 2010-12-03, 2011-04-18) |

2454 | """ |

2455 | H= self._H |

2456 | R= self._R |

2457 | f= R.ideal(H).radical().gens()[Integer(0) ] # Compute the reduction of H. |

2458 | J= R.ideal([f] + f.gradient()) |

2459 | return R.ideal(J.groebner_basis()) |

2460 | #------------------------------------------------------------------------------- |

2461 | def smooth_critical(self,alpha): |

2462 | r""" |

2463 | This function returns a Groebner basis ideal whose variety is the set |

2464 | of smooth critical points of the algebraic variety |

2465 | `V= \{x\in\CC^d : H(x)=0\} for the direction `\alpha` where `H=self._H`. |

2466 | |

2467 | INPUT: |

2468 | |

2469 | - ``alpha`` - A `d`-tuple of positive integers and/or symbolic entries. |

2470 | |

2471 | OUTPUT: |

2472 | |

2473 | A Groebner basis ideal of smooth critical points of `V` for `\alpha`. |

2474 | |

2475 | EXAMPLES:: |

2476 | |

2477 | sage: from sage.combinat.amgf import * |

2478 | sage: R.<x,y> = PolynomialRing(QQ) |

2479 | sage: G=1 |

2480 | sage: H = (1-x-y-x*y)^2 |

2481 | sage: F= QuasiRationalExpression(G,H) |

2482 | sage: var('a1,a2') |

2483 | (a1, a2) |

2484 | sage: F.smooth_critical([a1,a2]) |

2485 | 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 |

2486 | |

2487 | NOTES: |

2488 | |

2489 | A point `c` of `V` is a __smooth critical point for `alpha`__ |

2490 | if the gradient of `f` at `c` is not identically zero and `\alpha` is in |

2491 | the span of the logarithmic gradient vector |

2492 | `(c[0] \partial_1 f(c)),\ldots,c[d-1] \partial_d f(c))`; see [RaWi2008a]_. |

2493 | |

2494 | REFERENCES: |

2495 | |

2496 | .. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of |

2497 | coefficients of multivariate generating functions: improvements |

2498 | for smooth points", Electronic Journal of Combinatorics, Vol. 15 |

2499 | (2008), R89. |

2500 | |

2501 | AUTHORS: |

2502 | |

2503 | - Alex Raichev (2008-10-01, 2008-11-20, 2009-03-09, 2010-12-02, 2011-04-18) |

2504 | """ |

2505 | H= self._H |

2506 | R= H.parent() |

2507 | B= R.base_ring() |

2508 | d= R.ngens() |

2509 | vars= R.gens() |

2510 | f= R.ideal(H).radical().gens()[Integer(0) ] # Compute the reduction of H. |

2511 | |

2512 | # Expand B by the variables of alpha if there are any. |

2513 | indets= [] |

2514 | indets_ind= [] |

2515 | for a in alpha: |

2516 | if not ((a in ZZ) and (a>Integer(0) )): |

2517 | try: |

2518 | CC(a) |

2519 | except: |

2520 | indets.append(var(a)) |

2521 | indets_ind.append(alpha.index(a)) |

2522 | else: |

2523 | print "The components of", alpha, \ |

2524 | "must be positive integers or symbolic variables." |

2525 | return |

2526 | indets= list(Set(indets)) # Delete duplicates in indets. |

2527 | if indets != []: |

2528 | BB= FractionField(PolynomialRing(B,tuple(indets))) |

2529 | S= R.change_ring(BB) |

2530 | vars= S.gens() |

2531 | # Coerce alpha into BB. |

2532 | for i in range(len(alpha)): |

2533 | alpha[i] = BB(alpha[i]) |

2534 | else: |

2535 | S= R |

2536 | |

2537 | # Find smooth, critical points for alpha. |

2538 | f= S(f) |

2539 | J= S.ideal([f] +[ alpha[d-Integer(1) ]*vars[i]*diff(f,vars[i]) \ |

2540 | -alpha[i]*vars[d-Integer(1) ]*diff(f,vars[d-Integer(1) ]) for i in range(d-Integer(1) )]) |

2541 | return S.ideal(J.groebner_basis()) |

2542 | #------------------------------------------------------------------------------- |

2543 | def _subs_all(self,f,sub,simplify=False): |

2544 | r""" |

2545 | This function returns the items of `f` substituted by the dictionaries of |

2546 | `sub` in order of their appearance in `sub`. |

2547 | |

2548 | Does not use `self`. |

2549 | |

2550 | INPUT: |

2551 | |

2552 | - ``f`` - An individual or list of symbolic expressions or dictionaries |

2553 | - ``sub`` - An individual or list of dictionaries. |

2554 | - ``simplify`` - Boolean (default: False). |

2555 | |

2556 | OUTPUT: |

2557 | |

2558 | The items of `f` substituted by the dictionaries of `sub` in order of |

2559 | their appearance in `sub`. The subs() command is used. |

2560 | If simplify is True, then simplify() is used after substitution. |

2561 | |

2562 | EXAMPLES:: |

2563 | |

2564 | sage: from sage.combinat.amgf import * |

2565 | sage: R.<x> = PolynomialRing(QQ) |

2566 | sage: G,H = 1,1 |

2567 | sage: F= QuasiRationalExpression(G,H) |

2568 | sage: var('x,y,z') |

2569 | (x, y, z) |

2570 | sage: a= {x:1} |

2571 | sage: b= {y:2} |

2572 | sage: c= {z:3} |

2573 | sage: F._subs_all(x+y+z,a) |

2574 | y + z + 1 |

2575 | sage: F._subs_all(x+y+z,[c,a]) |

2576 | y + 4 |

2577 | sage: F._subs_all([x+y+z,y^2],b) |

2578 | [x + z + 2, 4] |

2579 | sage: F._subs_all([x+y+z,y^2],[b,c]) |

2580 | [x + 5, 4] |

2581 | |

2582 | :: |

2583 | |

2584 | sage: var('x,y') |

2585 | (x, y) |

2586 | sage: a= {'foo':x^2+y^2, 'bar':x-y} |

2587 | sage: b= {x:1,y:2} |

2588 | sage: F._subs_all(a,b) |

2589 | {'foo': 5, 'bar': -1} |

2590 | |

2591 | AUTHORS: |

2592 | |

2593 | - Alex Raichev (2009-05-05) |

2594 | """ |

2595 | singleton= False |

2596 | if not isinstance(f,list): |

2597 | f= [f] |

2598 | singleton= True |

2599 | if not isinstance(sub,list): |

2600 | sub= [sub] |

2601 | g= [] |

2602 | for ff in f: |

2603 | for D in sub: |

2604 | if isinstance(ff,dict): |

2605 | ff= dict( [(k,ff[k].subs(D)) for k in ff.keys()] ) |

2606 | else: |

2607 | ff= ff.subs(D) |

2608 | g.append(ff) |

2609 | if singleton and simplify: |

2610 | if isinstance(g[Integer(0) ],dict): |

2611 | return g[Integer(0) ] |

2612 | else: |

2613 | return g[Integer(0) ].simplify() |

2614 | elif singleton and not simplify: |

2615 | return g[Integer(0) ] |

2616 | elif not singleton and simplify: |

2617 | G= [] |

2618 | for gg in g: |

2619 | if isinstance(gg,dict): |

2620 | G.append(gg) |

2621 | else: |

2622 | G.append(gg.simplify()) |

2623 | return G |

2624 | else: |

2625 | return g |

2626 | # T ============================================================================ |

2627 | def maclaurin_coefficients(self,alpha,n): |

2628 | r""" |

2629 | Returns the Maclaurin coefficients of self that have multi-indices |

2630 | `alpha`, `2*alpha`,...,`n*alpha`. |

2631 | |

2632 | INPUT: |

2633 | |

2634 | - ``n`` - positive integer |

2635 | - ``alpha`` - tuple of positive integers representing a multi-index. |

2636 | |

2637 | OUTPUT: |

2638 | |

2639 | A dictionary of the form (beta, beta Maclaurin coefficient of self). |

2640 | |

2641 | AUTHORS: |

2642 | |

2643 | - Alex Raichev (2011-04-08) |

2644 | """ |

2645 | # Do all computations in the Symbolic Ring. |

2646 | F= SR(self._G/self._H) |

2647 | v= F.variables() |

2648 | p= {} |

2649 | for x in v: |

2650 | p[x]=Integer(0) |

2651 | d= len(v) |

2652 | coeffs={} |

2653 | # Initialize sequence of variables to differentiate with. |

2654 | s=[] |

2655 | for i in range(d): |

2656 | s.extend([v[i] for j in range(alpha[i])]) |

2657 | F_deriv= diff(F,s) |

2658 | coeffs[tuple(alpha)]= F_deriv.subs(p)/ mul([factorial(a) for a in alpha]) |

2659 | old_beta= alpha |

2660 | for k in (ellipsis_range(Integer(2) ,Ellipsis,n)): |

2661 | # Update variable sequence to differentiate with. |

2662 | beta= [k*a for a in alpha] |

2663 | delta= [beta[i]-old_beta[i] for i in range(d)] |

2664 | s= [] |

2665 | for i in range(d): |

2666 | s.extend([v[i] for j in range(delta[i])]) |

2667 | F_deriv= diff(F_deriv,s) |

2668 | coeffs[tuple(beta)]= F_deriv.subs(p)/ mul([factorial(b) for b in beta]) |

2669 | old_beta= copy(beta) |

2670 | return coeffs |

2671 | # U ============================================================================ |

2672 | # V ============================================================================ |

2673 | def variables(self): |

2674 | r""" |

2675 | Returns the tuple of variables of `self`. |

2676 | |

2677 | EXAMPLES:: |

2678 | |

2679 | sage: from sage.combinat.amgf import * |

2680 | sage: R.<x,y>= PolynomialRing(QQ) |

2681 | sage: G= exp(x) |

2682 | sage: H= 1-y |

2683 | sage: F= QuasiRationalExpression(G,H) |

2684 | sage: F.variables() |

2685 | (x, y) |

2686 | |

2687 | sage: R.<x,y>= PolynomialRing(QQ,order='invlex') |

2688 | sage: G= exp(x) |

2689 | sage: H= 1-y |

2690 | sage: F= QuasiRationalExpression(G,H) |

2691 | sage: F.variables() |

2692 | (y, x) |

2693 | |

2694 | AUTHORS: |

2695 | |

2696 | - Alex Raichev (2011-04-01) |

2697 | """ |

2698 | return self._variables |

2699 | # W ============================================================================ |

2700 | # X ============================================================================ |

2701 | # Y ============================================================================ |

2702 | # Z ============================================================================ |