# Ticket #10519: amgf.sage

File amgf.sage, 105.5 KB (added by , 6 years ago) |
---|

Line | |
---|---|

1 | r""" |

2 | This code relates to analytic combinatorics. |

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

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

5 | multivariate generating functions. |

6 | |

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

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

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

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

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

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

13 | smooth or multiple points. |

14 | |

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

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

17 | and [RaWi2011]_. |

18 | |

19 | REFERENCES: |

20 | |

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

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

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

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

25 | |

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

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

28 | |

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

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

31 | |

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

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

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

35 | |

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

37 | examples of asymptotics derived from multivariate generating |

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

39 | |

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

41 | coefficients of multivariate generating functions: improvements |

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

43 | (2008), R89. |

44 | |

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

46 | coefficients of multivariate generating functions: improvements |

47 | for smooth points", To appear. |

48 | |

49 | AUTHORS: |

50 | |

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

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

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

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

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

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

57 | fixed _crit_cone_combo() to work in SR |

58 | |

59 | EXAMPLES:: |

60 | |

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

62 | A smooth point example. :: |

63 | |

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

65 | sage: G= 1 |

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

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

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

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

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

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

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

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

74 | [{y: -1/2*sqrt(13) - 3/2, x: -1/3*sqrt(13) - 2/3}, {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3}] |

75 | sage: c= s[1] |

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

77 | Initializing auxiliary functions... |

78 | Calculating derivatives of auxiallary functions... |

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

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

81 | Initializing auxiliary functions... |

82 | Calculating derivatives of auxiallary functions... |

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

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

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

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

87 | Calculating errors table in the form |

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

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

90 | |

91 | Another smooth point example. |

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

93 | |

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

95 | sage: q= 1/2 |

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

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

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

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

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

101 | [4, 1] |

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

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

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

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

106 | sage: c= s[0] |

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

108 | Initializing auxiliary functions... |

109 | Calculating derivatives of auxiallary functions... |

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

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

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

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

114 | Calculating errors table in the form |

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

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

117 | |

118 | A multiple point example. :: |

119 | |

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

121 | sage: G= 1 |

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

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

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

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

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

127 | sage: F.is_cmp(c) |

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

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

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

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

132 | sage: alpha in cc[1] |

133 | True |

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

135 | Initializing auxiliary functions... |

136 | Calculating derivatives of auxiliary functions... |

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

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

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

140 | |

141 | Another multiple point example. :: |

142 | |

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

144 | sage: G= 16 |

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

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

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

148 | (a1, a2, a3) |

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

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

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

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

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

154 | sage: F.is_cmp(c) |

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

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

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

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

159 | sage: alpha in cc[1] |

160 | True |

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

162 | Initializing auxiliary functions... |

163 | Calculating derivatives of auxiliary functions... |

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

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

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

167 | |

168 | """ |

169 | #***************************************************************************** |

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

171 | # |

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

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

174 | #***************************************************************************** |

175 | # |

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

177 | def algebraic_dependence(fs): |

178 | r""" |

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

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

181 | Otherwise it returns 0. |

182 | |

183 | INPUT: |

184 | |

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

186 | ring. |

187 | |

188 | OUTPUT: |

189 | |

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

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

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

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

194 | indeterminates different from those of self. |

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

196 | is the zero polynomial. |

197 | |

198 | EXAMPLES:: |

199 | |

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

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

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

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

204 | sage: g(fs) |

205 | 0 |

206 | |

207 | :: |

208 | |

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

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

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

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

213 | sage: g(fs) |

214 | 0 |

215 | sage: g.parent() |

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

217 | |

218 | :: |

219 | |

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

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

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

223 | 0 |

224 | |

225 | AUTHORS: |

226 | |

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

228 | """ |

229 | r= len(fs) |

230 | R= fs[0].parent() |

231 | B= R.base_ring() |

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

233 | d= len(Xs) |

234 | |

235 | # Expand R by r new variables. |

236 | T= 'T' |

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

238 | T= T+'T' |

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

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

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

242 | Xs= Vs[0:d] |

243 | Ts= Vs[d:] |

244 | |

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

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

247 | JJ= J.elimination_ideal(Xs) |

248 | g= JJ.gens()[0] |

249 | |

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

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

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

253 | return RRR(g) |

254 | #------------------------------------------------------------------------------- |

255 | def partial_fraction_decomposition(f): |

256 | r""" |

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

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

259 | |

260 | INPUT: |

261 | |

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

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

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

265 | |

266 | OUTPUT: |

267 | |

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

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

270 | See the notes below for more details. |

271 | |

272 | EXAMPLES:: |

273 | |

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

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

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

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

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

279 | True |

280 | |

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

282 | |

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

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

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

286 | sage: whole, parts |

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

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

289 | True |

290 | |

291 | :: |

292 | |

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

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

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

296 | sage: whole, parts |

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

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

299 | True |

300 | |

301 | :: |

302 | |

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

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

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

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

307 | True |

308 | |

309 | NOTES: |

310 | |

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

312 | partial fraction decomposition function. |

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

314 | partial fraction decompositions. |

315 | |

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

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

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

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

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

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

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

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

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

325 | coefficient field of `R`). |

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

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

328 | |

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

330 | [Lein1978]_. |

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

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

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

334 | |

335 | REFERENCES: |

336 | |

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

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

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

340 | |

341 | AUTHORS: |

342 | |

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

344 | """ |

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

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

347 | if d==1: |

348 | return f.partial_fraction_decomposition() |

349 | Q= f.denominator() |

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

351 | parts= [format_quotient(1/Q)] |

352 | # Decompose via nullstellensatz trick |

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

354 | parts= decompose_via_nullstellensatz(parts) |

355 | # Decompose via algebraic dependence trick |

356 | parts= decompose_via_algebraic_dependence(parts) |

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

358 | new_parts=[] |

359 | for p in parts: |

360 | f= unformat_quotient(p) |

361 | if f.denominator() == 1: |

362 | whole= whole +f |

363 | else: |

364 | new_parts.append(f) |

365 | # Put P back in |

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

367 | return whole, new_parts |

368 | #------------------------------------------------------------------------------- |

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

370 | r""" |

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

372 | For use by partial_fraction_decomposition() above. |

373 | |

374 | EXAMPLES: |

375 | |

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

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

378 | sage: combine_like_terms(parts) |

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

380 | |

381 | :: |

382 | |

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

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

385 | sage: combine_like_terms(parts) |

386 | [[1]] |

387 | |

388 | AUTHORS: |

389 | |

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

391 | """ |

392 | |

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

394 | # Sort parts by denominators |

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

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

397 | newnew_parts=[] |

398 | left,right=0,1 |

399 | glom= new_parts[left][0] |

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

401 | if new_parts[left][1:] == new_parts[right][1:]: |

402 | glom= glom +new_parts[right][0] |

403 | else: |

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

405 | left= right |

406 | glom= new_parts[left][0] |

407 | right= right +1 |

408 | if glom != 0: |

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

410 | if rationomial: |

411 | # Reduce fractions in newnew_parts. |

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

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

414 | # calls factor(). |

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

416 | return newnew_parts |

417 | #------------------------------------------------------------------------------- |

418 | def decompose_via_algebraic_dependence(parts): |

419 | r""" |

420 | Returns a decomposition of parts. |

421 | Used by partial_fraction_decomposition() above. |

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

423 | Recursive. |

424 | |

425 | REFERENCES: |

426 | |

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

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

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

430 | |

431 | AUTHORS: |

432 | |

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

434 | """ |

435 | decomposing_done= True |

436 | new_parts= [] |

437 | for p in parts: |

438 | p_parts= [p] |

439 | P= p[0] |

440 | Qs= p[1:] |

441 | m= len(Qs) |

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

443 | if G: |

444 | # Then the denominator factors are algebraically dependent |

445 | # and so we can decompose p. |

446 | decomposing_done= False |

447 | P= p[0] |

448 | Qs= p[1:] |

449 | |

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

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

452 | F= algebraic_dependence([q^e for q,e in Qs]) |

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

454 | |

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

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

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

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

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

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

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

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

463 | |

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

465 | # reduce terms and combine like terms. |

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

467 | p_parts_temp= combine_like_terms(p_parts_temp) |

468 | |

469 | # Substitute Qs into new_p. |

470 | Qpowsub= dict([(new_vars[j],Qs[j][0]^Qs[j][1]) for j in range(m)]) |

471 | p_parts=[] |

472 | for x in p_parts_temp: |

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

474 | yy=[] |

475 | for xx in x[1:]: |

476 | if xx[0] in new_vars: |

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

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

479 | else: |

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

481 | yy.append(xx) |

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

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

484 | new_parts.extend(p_parts) |

485 | if decomposing_done: |

486 | return new_parts |

487 | else: |

488 | return decompose_via_algebraic_dependence(new_parts) |

489 | #------------------------------------------------------------------------------- |

490 | def decompose_via_cohomology(parts): |

491 | r""" |

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

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

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

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

496 | |

497 | INPUT: |

498 | |

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

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

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

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

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

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

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

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

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

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

509 | |

510 | OUTPUT: |

511 | |

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

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

514 | |

515 | EXAMPLES:: |

516 | |

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

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

519 | [] |

520 | |

521 | :: |

522 | |

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

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

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

526 | |

527 | NOTES: |

528 | |

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

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

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

532 | [AiYu1983]_. |

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

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

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

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

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

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

539 | `c`. |

540 | |

541 | REFERENCES: |

542 | |

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

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

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

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

547 | |

548 | AUTHORS: |

549 | |

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

551 | """ |

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

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

554 | decomposing_done= True |

555 | new_parts= [] |

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

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

558 | for p in parts: |

559 | p_parts= [p] |

560 | P= p[0] |

561 | Qs= p[1:] |

562 | m= len(Qs) |

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

564 | # Then we can decompose p |

565 | p_parts= [] |

566 | decomposing_done= False |

567 | dets= [] |

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

569 | for v in vars_subsets: |

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

571 | # declared is first in vars list. |

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

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

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

575 | # Get a Nullstellensatz certificate. |

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

577 | for j in range(m): |

578 | if L[j] != 0: |

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

580 | new_Qs = copy.deepcopy(Qs) |

581 | if new_Qs[j][1] > 1: |

582 | new_Qs[j][1]= new_Qs[j][1] -1 |

583 | else: |

584 | del new_Qs[j] |

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

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

587 | if L[m+k] != 0: |

588 | new_Qs = copy.deepcopy(Qs) |

589 | for j in range(m): |

590 | if new_Qs[j][1] > 1: |

591 | new_Qs[j][1]= new_Qs[j][1] -1 |

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

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

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

595 | det= jac.determinant() |

596 | if det != 0: |

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

598 | *(-1)^j/new_Qs[j][1] *det] +new_Qs) |

599 | break |

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

601 | new_parts.extend(p_parts) |

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

603 | if decomposing_done: |

604 | return new_parts |

605 | else: |

606 | return decompose_via_cohomology(new_parts) |

607 | #------------------------------------------------------------------------------- |

608 | def decompose_via_nullstellensatz(parts): |

609 | r""" |

610 | Returns a decomposition of parts. |

611 | Used by partial_fraction_decomposition() above. |

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

613 | Recursive. |

614 | |

615 | REFERENCES: |

616 | |

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

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

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

620 | |

621 | AUTHORS: |

622 | |

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

624 | """ |

625 | decomposing_done= True |

626 | new_parts= [] |

627 | R= parts[0][0].parent() |

628 | for p in parts: |

629 | p_parts= [p] |

630 | P= p[0] |

631 | Qs= p[1:] |

632 | m= len(Qs) |

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

634 | # Then we can decompose p. |

635 | decomposing_done= False |

636 | L= R(1).lift(R.ideal([q^e for q,e in Qs])) |

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

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

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

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

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

642 | new_parts.extend(p_parts) |

643 | if decomposing_done: |

644 | return new_parts |

645 | else: |

646 | return decompose_via_nullstellensatz(new_parts) |

647 | #------------------------------------------------------------------------------- |

648 | def format_quotient(f): |

649 | r""" |

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

651 | |

652 | AUTHORS: |

653 | |

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

655 | """ |

656 | P= f.numerator() |

657 | Q= f.denominator() |

658 | Qs= Q.factor() |

659 | if Qs.unit() != 1: |

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

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

662 | return [P]+Qs |

663 | #------------------------------------------------------------------------------- |

664 | def permutation_sign(v,vars): |

665 | r""" |

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

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

668 | For internal use by _cohom_equiv_main(). |

669 | |

670 | INPUT: |

671 | |

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

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

674 | |

675 | OUTPUT: |

676 | |

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

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

679 | elements of `v` removed. |

680 | |

681 | EXAMPLES:: |

682 | |

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

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

685 | sage: permutation_sign(v,vars) |

686 | -1 |

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

688 | sage: permutation_sign(v,vars) |

689 | 1 |

690 | |

691 | AUTHORS: |

692 | |

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

694 | """ |

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

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

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

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

699 | C.sort() |

700 | P= Permutation(B+C) |

701 | return P.signature() |

702 | #------------------------------------------------------------------------------- |

703 | def unformat_quotient(part): |

704 | r""" |

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

706 | Inverse of format_quotient() above. |

707 | |

708 | AUTHORS: |

709 | |

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

711 | """ |

712 | P= part[0] |

713 | Qs= part[1:] |

714 | Q= prod([q^e for q,e in Qs]) |

715 | return P/Q |

716 | #=============================================================================== |

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

718 | class QuasiRationalExpression(object): |

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

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

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

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

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

724 | self._G = G |

725 | self._H = H |

726 | R= H.parent() |

727 | self._R = R |

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

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

730 | |

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

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

733 | try: |

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

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

736 | except: |

737 | Gv= Set([]) |

738 | try: |

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

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

741 | except: |

742 | Hv= Set([]) |

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

744 | # depends upon. |

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

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

747 | #------------------------------------------------------------------------------- |

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

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

750 | r""" |

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

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

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

754 | where `F = self`. |

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

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

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

758 | |

759 | INPUT: |

760 | |

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

762 | point, possibly of symbolic variables. |

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

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

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

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

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

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

769 | Otherwise exact values will be returned. |

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

771 | The variable of the asymptotic expansion. |

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

773 | |

774 | OUTPUT: |

775 | |

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

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

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

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

780 | `asy`. |

781 | |

782 | EXAMPLES:: |

783 | |

784 | A smooth point example :: |

785 | |

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

787 | sage: G= 1 |

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

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

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

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

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

793 | Initializing auxiliary functions... |

794 | Calculating derivatives of auxiallary functions... |

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

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

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

798 | |

799 | A multiple point example :: |

800 | |

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

802 | sage: G= 1 |

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

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

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

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

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

808 | Initializing auxiliary functions... |

809 | Calculating derivatives of auxiliary functions... |

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

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

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

813 | |

814 | NOTES: |

815 | |

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

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

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

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

820 | see the documentation for _asymptotics_main_smooth() and |

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

822 | work. |

823 | |

824 | ALGORITHM: |

825 | |

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

827 | |

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

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

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

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

832 | |

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

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

835 | in the denominator. |

836 | |

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

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

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

840 | |

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

842 | |

843 | (5) Add the expansions. |

844 | |

845 | REFERENCES: |

846 | |

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

848 | coefficients of multivariate generating functions: improvements |

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

850 | (2008), R89. |

851 | |

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

853 | coefficients of multivariate generating functions: improvements |

854 | for smooth points", submitted. |

855 | |

856 | AUTHORS: |

857 | |

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

859 | """ |

860 | # The variable for asymptotic expansions. |

861 | if not asy_var: |

862 | asy_var= var('n') |

863 | |

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

865 | R= self._R |

866 | d= self._d |

867 | X= list(self._variables) |

868 | |

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

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

871 | |

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

873 | # calculations are no longer needed. |

874 | # Calculate asymptotics. |

875 | cc={} |

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

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

878 | c= cc |

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

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

881 | subexp_parts= [] |

882 | for chunk in new_integrands: |

883 | # Convert chunk into Symbolic Ring |

884 | GG= SR(chunk[0]) |

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

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

887 | subexp_parts.append(asy[2]) |

888 | exp_scale= asy[1] # same for all chunk in new_integrands |

889 | |

890 | # Do step (5). |

891 | subexp_part= add(subexp_parts) |

892 | return exp_scale^asy_var *subexp_part, exp_scale, subexp_part |

893 | #-------------------------------------------------------------------------------- |

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

895 | r""" |

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

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

898 | whether to call _asymptotics_main_smooth() or |

899 | _asymptotics_main_multiple(). |

900 | |

901 | Does not use `self`. |

902 | |

903 | INPUT: |

904 | |

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

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

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

908 | Call its length `d`. |

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

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

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

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

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

914 | - ``numerical`` - Natural number. |

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

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

917 | Otherwise exact values will be returned. |

918 | |

919 | OUTPUT: |

920 | |

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

922 | |

923 | AUTHORS: |

924 | |

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

926 | """ |

927 | d= len(X) |

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

929 | |

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

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

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

933 | # convenient smooth or multiple points. |

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

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

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

937 | # in the same way. |

938 | x= X[d-1] |

939 | beta= copy(alpha) |

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

941 | X.pop() |

942 | X.insert(0,x) |

943 | x= X[d-1] |

944 | a= beta.pop() |

945 | beta.insert(0,a) |

946 | if d==r: |

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

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

949 | return A,1,A |

950 | elif r==1: |

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

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

953 | else: |

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

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

956 | #-------------------------------------------------------------------------------- |

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

958 | r""" |

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

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

961 | |

962 | Does not use `self`. |

963 | |

964 | INPUT: |

965 | |

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

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

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

969 | Call its length `d`. |

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

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

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

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

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

975 | - ``numerical`` - Natural number. |

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

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

978 | Otherwise exact values will be returned. |

979 | |

980 | OUTPUT: |

981 | |

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

983 | |

984 | NOTES: |

985 | |

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

987 | Theorem 3.4 of [RaWi2011]_. |

988 | |

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

990 | because _crit_cone_combo() crashes. |

991 | |

992 | REFERENCES: |

993 | |

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

995 | coefficients of multivariate generating functions: improvements |

996 | for smooth points", submitted. |

997 | |

998 | AUTHORS: |

999 | |

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

1001 | """ |

1002 | I= sqrt(-1) |

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

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

1005 | C= copy(c) |

1006 | |

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

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

1009 | Sstar= {} |

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

1011 | for j in range(r-1): |

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

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

1014 | thetastar.update(Sstar) |

1015 | |

1016 | # Setup. |

1017 | print "Initializing auxiliary functions..." |

1018 | Hmul= mul(H) |

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

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

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

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

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

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

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

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

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

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

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

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

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

1032 | atC = copy(C) |

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

1034 | |

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

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

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

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

1039 | .simplify() |

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

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

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

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

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

1045 | |

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

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

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

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

1050 | m= min(r,N) |

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

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

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

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

1055 | Uderivs={} |

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

1057 | if p == 0: |

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

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

1060 | for l in [1..2*N-2+m-1]: |

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

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

1063 | elif p > 0 and p < 2*N-2+m-1: |

1064 | all_zero= True |

1065 | Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,[1..p],end,Uderivs,atC) |

1066 | # Check for a nonzero U derivative. |

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

1068 | all_zero= False |

1069 | if all_zero: |

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

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

1072 | for l in [p+1..2*N-2+m-1]: |

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

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

1075 | else: |

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

1077 | Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,[p+1..2*N-2+m-1],end,Uderivs,atC) |

1078 | else: |

1079 | Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,[1..2*N-2+m-1],end,Uderivs,atC) |

1080 | atC.update(Uderivs) |

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

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

1083 | a_inv= a.inverse() |

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

1085 | Phitu= Phitu[0][0] |

1086 | |

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

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

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

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

1091 | # Choose nondecreasing sequences as representative differentiation- |

1092 | # order sequences. |

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

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

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

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

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

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

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

1100 | AABB_derivs= At_derivs |

1101 | AABB_derivs.update(Phitu_derivs) |

1102 | for j in range(r): |

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

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

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

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

1107 | |

1108 | L={} |

1109 | for (j,k) in CartesianProduct([0..min(r-1,N-1)], [max(0,N-1-r)..N-1]): |

1110 | if j+k <= N-1: |

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

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

1113 | for l in [0..2*k]] ) |

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

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

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

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

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

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

1120 | for (j,k) in CartesianProduct([0..min(r-1,q)], [max(0,q-r)..q]) if j+k <= q ]) \ |

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

1122 | |

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

1124 | coeffs= chunk.coefficients(n) |

1125 | coeffs.reverse() |

1126 | coeffs= coeffs[:N] |

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

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

1129 | for co in coeffs] ) |

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

1131 | .n(digits=numerical) |

1132 | else: |

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

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

1135 | return exp_scale^n*subexp_part, exp_scale, subexp_part |

1136 | #-------------------------------------------------------------------------------- |

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

1138 | r""" |

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

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

1141 | |

1142 | Does not use `self`. |

1143 | |

1144 | INPUT: |

1145 | |

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

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

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

1149 | Call its length `d`. |

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

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

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

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

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

1155 | - ``numerical`` - Natural number. |

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

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

1158 | Otherwise exact values will be returned. |

1159 | |

1160 | OUTPUT: |

1161 | |

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

1163 | |

1164 | NOTES: |

1165 | |

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

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

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

1169 | with `r=1`. |

1170 | |

1171 | REFERENCES: |

1172 | |

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

1174 | coefficients of multivariate generating functions: improvements |

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

1176 | (2008), R89. |

1177 | |

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

1179 | coefficients of multivariate generating functions: improvements |

1180 | for smooth points", submitted. |

1181 | |

1182 | AUTHORS: |

1183 | |

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

1185 | """ |

1186 | I= sqrt(-1) |

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

1188 | |

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

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

1191 | if vector(c.values()) in QQ^d: |

1192 | C= c |

1193 | else: |

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

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

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

1197 | |

1198 | # Setup. |

1199 | print "Initializing auxiliary functions..." |

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

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

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

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

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

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

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

1207 | ht= h.subs(e) |

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

1209 | At= A.subs(e) |

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

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

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

1213 | atC = copy(C) |

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

1215 | # evaluated at C. |

1216 | |

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

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

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

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

1221 | for i in range(d-1): |

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

1223 | .simplify() |

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

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

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

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

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

1229 | |

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

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

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

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

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

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

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

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

1238 | # Uderivs. |

1239 | Uderivs={} |

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

1241 | end= [X[d-1]] |

1242 | Hcheck= X[d-1] - 1/h |

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

1244 | if p == 0: |

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

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

1247 | for l in [1..2*N]: |

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

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

1250 | elif p > 0 and p < 2*N: |

1251 | all_zero= True |

1252 | Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,[1..p],end,Uderivs,atC) |

1253 | # Check for a nonzero U derivative. |

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

1255 | all_zero= False |

1256 | if all_zero: |

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

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

1259 | for l in [p+1..2*N]: |

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

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

1262 | else: |

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

1264 | Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,[p+1..2*N],end,Uderivs,atC) |

1265 | else: |

1266 | Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,[1..2*N],end,Uderivs,atC) |

1267 | atC.update(Uderivs) |

1268 | |

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

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

1271 | if d==2: |

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

1273 | v=2 |

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

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

1276 | while splat==0: |

1277 | v= v+1 |

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

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

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

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

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

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

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

1285 | if d==2 and v>2: |

1286 | t= T[0] # Simplify variable names. |

1287 | a= splat/factorial(v) |

1288 | Phitu= Phit -a*t^v |

1289 | |

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

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

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

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

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

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

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

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

1298 | if v.mod(2)==0: |

1299 | At_derivs= self._diff_all(At,T,2*N-2, \ |

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

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

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

1303 | else: |

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

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

1306 | AABB_derivs= At_derivs |

1307 | AABB_derivs.update(Phitu_derivs) |

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

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

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

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

1312 | # Plug above into asymptotic formula. |

1313 | L = [] |

1314 | if v.mod(2) == 0: |

1315 | for k in range(N): |

1316 | L.append( add([ \ |

1317 | (-1)^l *gamma((2*k+v*l+1)/v) \ |

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

1319 | * DD[(k,l)] for l in [0..2*k] ]) ) |

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

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

1322 | else: |

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

1324 | for k in range(N): |

1325 | L.append( add([ \ |

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

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

1328 | * (zeta^(k+v*l+1) +(-1)^(k+v*l)*zeta^(-(k+v*l+1))) \ |

1329 | * DD[(k,l)] for l in [0..k] ]) ) |

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

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

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

1333 | # in this case. |

1334 | else: |

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

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

1337 | a_inv= a.inverse() |

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

1339 | Phitu= Phitu[0][0] |

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

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

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

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

1344 | # Choose nondecreasing sequences as representative differentiation- |

1345 | # order sequences. |

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

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

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

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

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

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

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

1353 | AABB_derivs= At_derivs |

1354 | AABB_derivs.update(Phitu_derivs) |

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

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

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

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

1359 | |

1360 | # Plug above into asymptotic formula. |

1361 | L=[] |

1362 | for k in range(N): |

1363 | L.append( add([ \ |

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

1365 | for l in [0..2*k]]) ) |

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

1367 | *alpha[d-1]^((1-d)/2 -k) *L[k] \ |

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

1369 | |

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

1371 | coeffs= chunk.coefficients(n) |

1372 | coeffs.reverse() |

1373 | coeffs= coeffs[:N] |

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

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

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

1377 | .n(digits=numerical) |

1378 | else: |

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

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

1381 | return exp_scale^n*subexp_part, exp_scale, subexp_part |

1382 | #----------------------------------------------------------------------------- |

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

1384 | r""" |

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

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

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

1388 | |

1389 | INPUT: |

1390 | |

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

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

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

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

1395 | |

1396 | OUTPUT: |

1397 | |

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

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

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

1401 | is the `j`th item in `fs` and where `a` is _direction(alpha). |

1402 | |

1403 | EXAMPLES:: |

1404 | |

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

1406 | sage: G,H= 1,1 |

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

1408 | sage: fs= [x + 2*y + z - 4, 2*x + y + z - 4] |

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

1410 | sage: alpha= [2,1,1] |

1411 | sage: F._crit_cone_combo(fs,R.gens(),c,alpha) |

1412 | [0, 1] |

1413 | |

1414 | NOTES: |

1415 | |

1416 | Use this function only when `Gamma` is well-defined and |

1417 | there is a unique solution to the matrix equation `y Gamma = a`. |

1418 | Fails otherwise. |

1419 | |

1420 | AUTHORS: |

1421 | |

1422 | - Alex Raichev (2008-10-01, 2008-11-25, 2009-03-04, 2010-09-08, |

1423 | 2010-12-02) |

1424 | """ |

1425 | # Assuming here that each _log_grad(f) has nonzero final component. |

1426 | # Then 'direction' will not throw a division by zero error. |

1427 | d= len(X) |

1428 | r= len(fs) |

1429 | Gamma= matrix([self._direction(self._log_grad(f,X,c)) for f in fs]) |

1430 | # solve_left() fails when working in SR :-(. So use solve() instead. |

1431 | #s= Gamma.solve_left(vector(alpha)/alpha[d-1]) |

1432 | V= [var('sss'+str(i)) for i in range(r)] |

1433 | M= matrix(V)*Gamma |

1434 | eqns= [M[0][i]== alpha[i]/alpha[d-1] for i in range(d)] |

1435 | s= solve(eqns,V,solution_dict=True)[0] # Assuming a unique solution. |

1436 | return [s[v] for v in V] |

1437 | # B ========================================================================== |

1438 | # C ========================================================================== |

1439 | def cohomologous_integrand(self,alpha,asy_var=None): |

1440 | r""" |

1441 | This function takes a multivariate Cauchy type integral |

1442 | `\int F / x^{asy_var\alpha+1} dx`, where `F=self`, and breaks it up |

1443 | into a list of nicer Cauchy type integrals for the purposes of |

1444 | computing asymptotics of the original integral as `asy_var\to\infty`. |

1445 | The sum of the nicer integrals is de Rham cohomologous to the original |

1446 | integral. |

1447 | It assumes that algebraic varieties corresponding to the irreducible |

1448 | factors of `self._H` intersect transversely (see notes below). |

1449 | |

1450 | INPUT: |

1451 | |

1452 | - ``alpha`` - A list of positive integers or symbolic variables. |

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

1454 | Eventually set to `var('n')` if None is given. |

1455 | |

1456 | OUTPUT: |

1457 | |

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

1459 | `chunk_j` has the form `[P,[B_1,1],\ldots,[B_l,1]]`. |

1460 | Here `l \le d`, `P` is a symbolic expression in the indeterminates of |

1461 | `R` and `asy_var`, `\{B_1,\ldots,B_l\} \subseteq \{Q_1,\ldots,Q_m\}`, |

1462 | and `[P,[B_1,1],\ldots,[B_l,1]]` represents the integral |

1463 | `\int P/(B_1 \cdots B_l) dx`. |

1464 | |

1465 | EXAMPLES:: |

1466 | |

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

1468 | sage: G= 9*exp(x+y) |

1469 | sage: H= (3-2*x-y)*(3-x-2*y) |

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

1471 | sage: alpha= [4,3] |

1472 | sage: F.cohomologous_integrand(alpha) |

1473 | [[9*e^(x + y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]] |

1474 | |

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

1476 | sage: G= 9*exp(x+y) |

1477 | sage: H= (3-2*x-y)^2*(3-x-2*y) |

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

1479 | sage: alpha= [4,3] |

1480 | sage: F.cohomologous_integrand(alpha) |

1481 | [[-3*(3*x*e^x - 8*y*e^x)*n*e^y/(x*y) - 3*((x - 2)*y*e^x + x*e^x)*e^y/(x*y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]] |

1482 | |

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

1484 | sage: G= 16 |

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

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

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

1488 | sage: F.cohomologous_integrand(alpha) |

1489 | [[16*(4*y - 3*z)*n/(y*z) + 16*(2*y - z)/(y*z), [x + 2*y + z - 4, 1], [2*x + y + z - 4, 1]]] |

1490 | |

1491 | NOTES: |

1492 | |

1493 | The varieties corresponding to `Q_1,\ldots,Q_m` |

1494 | __intersect transversely__ iff for each point `c` of their intersection |

1495 | and for all `k \le l`, the Jacobian matrix of any `k` polynomials from |

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

1497 | `c`. |

1498 | |

1499 | ALGORITHM: |

1500 | |

1501 | Let `asy_var= n` and consider the integral around a polycirle centered |

1502 | at the origin of the `d`-variate differential form |

1503 | `\frac{G(x) dx_1 \wedge\cdots\wedge dx_d}{H(x) x^{n\alpha+1}}`, where |

1504 | `G=self._G` and `H=self._H`. |

1505 | |

1506 | (1) Decompose `G/H` into a sum of partial fractions `P_1 +\cdots+ P_r` |

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

1508 | in the denominator. |

1509 | |

1510 | (2) For each differential form `P_j dx_1 \wedge\cdots\wedge dx_d`, |

1511 | find an equivalent form `\omega_j` in de Rham cohomology with no |

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

1513 | |

1514 | AUTHOR: |

1515 | |

1516 | - Alex Raichev (2010-09-22) |

1517 | """ |

1518 | if not asy_var: |

1519 | asy_var = var('n') |

1520 | |

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

1522 | G= self._G |

1523 | H= self._H |

1524 | R= self._R |

1525 | d= self._d |

1526 | X= self._variables |

1527 | # Prepare input for partial_fraction_decomposition() which only works |

1528 | # for functions in the field of fractions of R. |

1529 | if G in R: |

1530 | numer=1 |

1531 | F= G/H |

1532 | else: |

1533 | numer= G |

1534 | F= R(1)/H |

1535 | nstuff= 1/mul([X[j]^(alpha[j]*asy_var+1) for j in range(d)]) |

1536 | |

1537 | # Do steps (2) and (1). |

1538 | integrands= [] |

1539 | whole,parts= partial_fraction_decomposition(F) |

1540 | for f in parts: |

1541 | a= format_quotient(f) |

1542 | integrands.append( [a[0]*numer*nstuff] + a[1:] ) |

1543 | |

1544 | # Do step (3). |

1545 | ce= decompose_via_cohomology(integrands) |

1546 | ce_new= [] |

1547 | for a in ce: |

1548 | ce_new.append( [(a[0]/nstuff).simplify_full().collect(n)] + a[1:] ) |

1549 | return ce_new |

1550 | #------------------------------------------------------------------------------- |

1551 | def critical_cone(self,c,coordinate=None): |

1552 | r""" |

1553 | Returns the critical cone of a convenient multiple point `c`. |

1554 | |

1555 | INPUT: |

1556 | |

1557 | - ``c`` - A dictionary with keys `self.variables()` and values |

1558 | in a field. |

1559 | - ``coordinate`` - (optional) A natural number. |

1560 | |

1561 | OUTPUT: |

1562 | |

1563 | A list of vectors that generate the critical cone of `c` and |

1564 | the cone itself if the values of `c` lie in QQ. |

1565 | |

1566 | EXAMPLES:: |

1567 | |

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

1569 | sage: G= 1 |

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

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

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

1573 | sage: F.critical_cone(c) |

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

1575 | |

1576 | NOTES: |

1577 | |

1578 | The _critical cone_ of a convenient multiple point `c` with |

1579 | with `c_k \del_k H_j(c) \neq 0` for all `j=1,\ldots,r` is |

1580 | the conical hull of the vectors `\gamma_j(c) = |

1581 | \left(\frac{c_1 \del_1 H_j(c)}{c_k \del_k H_j(c)},\ldots, |

1582 | \frac{c_d \del_d H_j(c)}{c_k \del_k H_j(c)} \right)`. |

1583 | Here `H_1,\ldots,H_r` are the irreducible germs of `self._H` around `c`. |

1584 | For more details, see [RaWi2011]_. |

1585 | |

1586 | If this function's optional argument `coordinate` isn't given, then |

1587 | this function searches (from `d` down to 1) for the first index `k` |

1588 | such that for all `j=1,\ldots,r` we have `c_k \del_k H_j(c) \neq 0` |

1589 | and sets `coordinate = k`. |

1590 | Almost. |

1591 | Since this is Python, all the indices actually start at 0. |

1592 | |

1593 | REFERENCES: |

1594 | |

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

1596 | coefficients of multivariate generating functions: improvements |

1597 | for smooth points", submitted. |

1598 | |

1599 | AUTHORS: |

1600 | |

1601 | - Alex Raichev (2010-08-25) |

1602 | """ |

1603 | Hs= [SR(h[0]) for h in self._Hfac] # irreducible factors of H |

1604 | X= self._variables |

1605 | d= self._d |

1606 | # Ensure the variables of `c` lie in SR |

1607 | cc= {} |

1608 | for x in c.keys(): |

1609 | cc[SR(x)] = c[x] |

1610 | lg= [self._log_grad(h,X,cc) for h in Hs] |

1611 | if not coordinate: |

1612 | # Search (from d-1 down to 0) for a coordinate k such that |

1613 | # for all h in Hs we have cc[k] * diff(h,X[k]) !=0. |

1614 | # One is guaranteed to exist in the case of a convenient multiple |

1615 | # point. |

1616 | for k in reversed(range(d)): |

1617 | if 0 not in [v[k] for v in lg]: |

1618 | coordinate= k |

1619 | break |

1620 | gamma= [self._direction(v,coordinate) for v in lg] |

1621 | if [[gg in QQ for gg in g] for g in gamma] == \ |

1622 | [[True for gg in g] for g in gamma]: |

1623 | return gamma,Cone(gamma) # Cone() needs rational vectors |

1624 | else: |

1625 | return gamma |

1626 | # D ============================================================================ |

1627 | def _diff_all(self,f,V,n,ending=[],sub=None,sub_final=None,zero_order=0,rekey=None): |

1628 | r""" |

1629 | This function returns a dictionary of representative mixed partial |

1630 | derivatives of `f` from order 1 up to order `n` with respect to the |

1631 | variables in `V`. |

1632 | The default is to key the dictionary by all nondecreasing sequences |

1633 | in `V` of length 1 up to length `n`. |

1634 | For internal use. |

1635 | |

1636 | Does not use `self`. |

1637 | |

1638 | INPUT: |

1639 | |

1640 | - ``f`` - An individual or list of `\mathcal{C}^{n+1}` functions. |

1641 | - ``V`` - A list of variables occurring in `f`. |

1642 | - ``n`` - A natural number. |

1643 | - ``ending`` - A list of variables in `V`. |

1644 | - ``sub`` - An individual or list of dictionaries. |

1645 | - ``sub_final`` - An individual or list of dictionaries. |

1646 | - ``rekey`` - A callable symbolic function in `V` or list thereof. |

1647 | - ``zero_order`` - A natural number. |

1648 | |

1649 | OUTPUT: |

1650 | |

1651 | The dictionary `{s_1:deriv_1,...,s_r:deriv_r}`. |

1652 | Here `s_1,\ldots,s_r` is a listing of |

1653 | all nondecreasing sequences of length 1 up to length `n` over the |

1654 | alphabet `V`, where `w > v` in `X` iff `str(w) > str(v)`, and |

1655 | `deriv_j` is the derivative of `f` with respect to the derivative |

1656 | sequence `s_j` and simplified with respect to the substitutions in `sub` |

1657 | and evaluated at `sub_final`. |

1658 | Moreover, all derivatives with respect to sequences of length less than |

1659 | `zero_order` (derivatives of order less than `zero_order`) will be made |

1660 | zero. |

1661 | |

1662 | If `rekey` is nonempty, then `s_1,\ldots,s_r` will be replaced by the |

1663 | symbolic derivatives of the functions in `rekey`. |

1664 | |

1665 | If `ending` is nonempty, then every derivative sequence `s_j` will be |

1666 | suffixed by `ending`. |

1667 | |

1668 | EXAMPLES:: |

1669 | |

1670 | I'd like to print the entire dictionaries, but that doesn't yield |

1671 | consistent output order for doctesting. |

1672 | Order of keys changes. :: |

1673 | |

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

1675 | sage: G,H = 1,1 |

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

1677 | sage: f= function('f',x) |

1678 | sage: dd= F._diff_all(f,[x],3) |

1679 | sage: dd[(x,x,x)] |

1680 | D[0, 0, 0](f)(x) |

1681 | |

1682 | :: |

1683 | |

1684 | sage: d1= {diff(f,x): 4*x^3} |

1685 | sage: dd= F._diff_all(f,[x],3,sub=d1) |

1686 | sage: dd[(x,x,x)] |

1687 | 24*x |

1688 | |

1689 | :: |

1690 | |

1691 | sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f) |

1692 | sage: dd[diff(f,x,3)] |

1693 | 24*x |

1694 | |

1695 | :: |

1696 | |

1697 | sage: a= {x:1} |

1698 | sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f,sub_final=a) |

1699 | sage: dd[diff(f,x,3)] |

1700 | 24 |

1701 | |

1702 | :: |

1703 | |

1704 | sage: X= var('x,y,z') |

1705 | sage: f= function('f',*X) |

1706 | sage: dd= F._diff_all(f,X,2,ending=[y,y,y]) |

1707 | sage: dd[(z,y,y,y)] |

1708 | D[1, 1, 1, 2](f)(x, y, z) |

1709 | |

1710 | :: |

1711 | |

1712 | sage: g= function('g',*X) |

1713 | sage: dd= F._diff_all([f,g],X,2) |

1714 | sage: dd[(0,y,z)] |

1715 | D[1, 2](f)(x, y, z) |

1716 | |

1717 | :: |

1718 | |

1719 | sage: dd[(1,z,z)] |

1720 | D[2, 2](g)(x, y, z) |

1721 | |

1722 | :: |

1723 | |

1724 | sage: f= exp(x*y*z) |

1725 | sage: ff= function('ff',*X) |

1726 | sage: dd= F._diff_all(f,X,2,rekey=ff) |

1727 | sage: dd[diff(ff,x,z)] |

1728 | x*y^2*z*e^(x*y*z) + y*e^(x*y*z) |

1729 | |

1730 | AUTHORS: |

1731 | |

1732 | - Alex Raichev (2008-10-01, 2009-04-01, 2010-02-01) |

1733 | """ |

1734 | singleton=False |

1735 | if not isinstance(f,list): |

1736 | f= [f] |

1737 | singleton=True |

1738 | |

1739 | # Build the dictionary of derivatives iteratively from a list of nondecreasing |

1740 | # derivative-order sequences. |

1741 | derivs= {} |

1742 | r= len(f) |

1743 | if ending: |

1744 | seeds = [ending] |

1745 | start = 1 |

1746 | else: |

1747 | seeds = [[v] for v in V] |

1748 | start = 2 |

1749 | if singleton: |

1750 | for s in seeds: |

1751 | derivs[tuple(s)] = self._subs_all(diff(f[0],s),sub) |

1752 | for l in [start..n]: |

1753 | for t in UnorderedTuples(V,l): |

1754 | s= tuple(t + ending) |

1755 | derivs[s] = self._subs_all(diff(derivs[s[1:]],s[0]),sub) |

1756 | else: |

1757 | # Make the dictionary keys of the form (j,sequence of variables), |

1758 | # where j in range(r). |

1759 | for s in seeds: |

1760 | value= self._subs_all([diff(f[j],s) for j in range(r)],sub) |

1761 | derivs.update(dict([(tuple([j]+s),value[j]) for j in range(r)])) |

1762 | for l in [start..n]: |

1763 | for t in UnorderedTuples(V,l): |

1764 | s= tuple(t + ending) |

1765 | value= self._subs_all(\ |

1766 | [diff(derivs[(j,)+s[1:]],s[0]) for j in range(r)],sub) |

1767 | derivs.update(dict([((j,)+s,value[j]) for j in range(r)])) |

1768 | if zero_order: |

1769 | # Zero out all the derivatives of order < zero_order |

1770 | if singleton: |

1771 | for k in derivs.keys(): |

1772 | if len(k) < zero_order: |

1773 | derivs[k]= 0 |

1774 | else: |

1775 | # Ignore the first of element of k, which is an index. |

1776 | for k in derivs.keys(): |

1777 | if len(k)-1 < zero_order: |

1778 | derivs[k]= 0 |

1779 | if sub_final: |

1780 | # Substitute sub_final into the values of derivs. |

1781 | for k in derivs.keys(): |

1782 | derivs[k] = self._subs_all(derivs[k],sub_final) |

1783 | if rekey: |

1784 | # Rekey the derivs dictionary by the value of rekey. |

1785 | F= rekey |

1786 | if singleton: |

1787 | # F must be a singleton. |

1788 | derivs= dict( [(diff(F,list(k)), derivs[k]) for k in derivs.keys()] ) |

1789 | else: |

1790 | # F must be a list. |

1791 | derivs= dict( [(diff(F[k[0]],list(k)[1:]), derivs[k]) for k in derivs.keys()] ) |

1792 | return derivs |

1793 | #------------------------------------------------------------------------------- |

1794 | def _diff_op(self,A,B,AB_derivs,V,M,r,N): |

1795 | r""" |

1796 | This function computes the derivatives `DD^(l+k)(A[j] B^l)` evaluated at a |

1797 | point `p` for various natural numbers `j,k,l` which depend on `r` and `N`. |

1798 | Here `DD` is a specific second-order linear differential operator that depends |

1799 | on `M`, `A` is a list of symbolic functions, `B` is symbolic function, |

1800 | and `AB_derivs` contains all the derivatives of `A` and `B` evaluated at `p` |

1801 | that are necessary for the computation. |

1802 | For internal use by the functions _asymptotics_main_multiple() and |

1803 | _asymptotics_main_smooth(). |

1804 | |

1805 | Does not use `self`. |

1806 | |

1807 | INPUT: |

1808 | |

1809 | - ``A`` - A single or length `r` list of symbolic functions in the |

1810 | variables `V`. |

1811 | - ``B`` - A symbolic function in the variables `V`. |

1812 | - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of |

1813 | `A[0],\ldots,A[r-1]` up to order `2N-2` and |

1814 | the (symbolic) derivatives of `B` up to order `2N`. |

1815 | The values of the dictionary are complex numbers that are |

1816 | the keys evaluated at a common point `p`. |

1817 | - ``V`` - The variables of the `A[j]` and `B`. |

1818 | - ``M`` - A symmetric `l \times l` matrix, where `l` is the length of `V`. |

1819 | - ``r,N`` - Natural numbers. |

1820 | |

1821 | OUTPUT: |

1822 | |

1823 | A dictionary whose keys are natural number tuples of the form `(j,k,l)`, |

1824 | where `l \le 2k`, `j \le r-1`, and `j+k \le N-1`, and whose values are |

1825 | `DD^(l+k)(A[j] B^l)` evaluated at a point `p`, where `DD` is the linear |

1826 | second-order differential operator |

1827 | `-\sum_{i=0}^{l-1} \sum_{j=0}^{l-1} M[i][j] |

1828 | \partial^2 /(\partial V[j] \partial V[i])`. |

1829 | |

1830 | EXAMPLES:: |

1831 | |

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

1833 | sage: G,H = 1,1 |

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

1835 | sage: T= var('x,y') |

1836 | sage: A= function('A',*tuple(T)) |

1837 | sage: B= function('B',*tuple(T)) |

1838 | sage: AB_derivs= {} |

1839 | sage: M= matrix([[1,2],[2,1]]) |

1840 | sage: DD= F._diff_op(A,B,AB_derivs,T,M,1,2) |

1841 | sage: DD.keys() |

1842 | [(0, 1, 2), (0, 1, 1), (0, 1, 0), (0, 0, 0)] |

1843 | sage: len(DD[(0,1,2)]) |

1844 | 246 |

1845 | |

1846 | AUTHORS: |

1847 | |

1848 | - Alex Raichev (2008-10-01, 2010-01-12) |

1849 | """ |

1850 | if not isinstance(A,list): |

1851 | A= [A] |

1852 | |

1853 | # First, compute the necessary product derivatives of A and B. |

1854 | product_derivs= {} |

1855 | for (j,k) in mrange([r,N]): |

1856 | if j+k <N: |

1857 | for l in [0..2*k]: |

1858 | for s in UnorderedTuples(V,2*(k+l)): |

1859 | product_derivs[tuple([j,k,l]+s)] = \ |

1860 | diff(A[j]*B^l,s).subs(AB_derivs) |

1861 | |

1862 | # Second, compute DD^(k+l)(A[j]*B^l)(p) and store values in dictionary. |

1863 | DD= {} |

1864 | rows= M.nrows() |

1865 | for (j,k) in mrange([r,N]): |

1866 | if j+k <N: |

1867 | for l in [0..2*k]: |

1868 | # Take advantage of the symmetry of M by ignoring |

1869 | # the upper-diagonal entries of M and multiplying by |

1870 | # appropriate powers of 2. |

1871 | if k+l == 0: |

1872 | DD[(j,k,l)] = product_derivs[(j,k,l)] |

1873 | continue |

1874 | S= [(a,b) for (a,b) in mrange([rows,rows]) if b <= a] |

1875 | P= cartesian_product_iterator([S for i in range(k+l)]) |

1876 | diffo= 0 |

1877 | for t in P: |

1878 | if product_derivs[(j,k,l)+self._diff_seq(V,t)] != 0: |

1879 | MM= 1 |

1880 | for (a,b) in t: |

1881 | MM= MM * M[a][b] |

1882 | if a != b: |

1883 | MM= 2*MM |

1884 | diffo= diffo + MM * product_derivs[(j,k,l)+self._diff_seq(V,t)] |

1885 | DD[(j,k,l)] = (-1)^(k+l)*diffo |

1886 | return DD |

1887 | #------------------------------------------------------------------------------- |

1888 | def _diff_op_simple(self,A,B,AB_derivs,x,v,a,N): |

1889 | r""" |

1890 | This function computes `DD^(e k + v l)(A B^l)` evaluated at a point `p` |

1891 | for various natural numbers `e,k,l` that depend on `v` and `N`. |

1892 | Here `DD` is a specific linear differential operator that depends |

1893 | on `a` and `v`, `A` and `B` are symbolic functions, and `AB_derivs` contains |

1894 | all the derivatives of `A` and `B` evaluated at `p` that are necessary for |

1895 | the computation. |

1896 | For internal use by the function _asymptotics_main_smooth(). |

1897 | |

1898 | Does not use `self`. |

1899 | |

1900 | INPUT: |

1901 | |

1902 | - ``A``,``B`` - Symbolic functions in the variable `x`. |

1903 | - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of |

1904 | `A` up to order `2N` if `v` is even or `N` if `v` is odd and |

1905 | the (symbolic) derivatives of `B` up to order `2N+v` if `v` is even |

1906 | or `N+v` if `v` is odd. |

1907 | The values of the dictionary are complex numbers that are |

1908 | the keys evaluated at a common point `p`. |

1909 | - ``x`` - Symbolic variable. |

1910 | - ``a`` - A complex number. |

1911 | - ``v``,``N`` - Natural numbers. |

1912 | |

1913 | OUTPUT: |

1914 | |

1915 | A dictionary whose keys are natural number pairs of the form `(k,l)`, |

1916 | where `k < N` and `l \le 2k` and whose values are |

1917 | `DD^(e k + v l)(A B^l)` evaluated at a point `p`. |

1918 | Here `e=2` if `v` is even, `e=1` if `v` is odd, and `DD` is a particular |

1919 | linear differential operator |

1920 | `(a^{-1/v} d/dt)' if `v` is even and `(|a|^{-1/v} i \sgn a d/dt)` |

1921 | if `v` is odd. |

1922 | |

1923 | EXAMPLES:: |

1924 | |

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

1926 | sage: G,H = 1,1 |

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

1928 | sage: A= function('A',x) |

1929 | sage: B= function('B',x) |

1930 | sage: AB_derivs= {} |

1931 | sage: F._diff_op_simple(A,B,AB_derivs,x,3,2,2) |

1932 | {(1, 0): 1/2*I*2^(2/3)*D[0](A)(x), (0, 0): A(x), (1, 1): 1/4*(A(x)*D[0, 0, 0, 0](B)(x) + B(x)*D[0, 0, 0, 0](A)(x) + 4*D[0](A)(x)*D[0, 0, 0](B)(x) + 4*D[0](B)(x)*D[0, 0, 0](A)(x) + 6*D[0, 0](A)(x)*D[0, 0](B)(x))*2^(2/3)} |

1933 | |

1934 | AUTHORS: |

1935 | |

1936 | - Alex Raichev (2010-01-15) |

1937 | """ |

1938 | I= sqrt(-1) |

1939 | DD= {} |

1940 | if v.mod(2) == 0: |

1941 | for k in range(N): |

1942 | for l in [0..2*k]: |

1943 | DD[(k,l)] = (a^(-1/v))^(2*k+v*l) \ |

1944 | * diff(A*B^l,x,2*k+v*l).subs(AB_derivs) |

1945 | else: |

1946 | for k in range(N): |

1947 | for l in [0..k]: |

1948 | DD[(k,l)] = (abs(a)^(-1/v)*I*a/abs(a))^(k+v*l) \ |

1949 | * diff(A*B^l,x,k+v*l).subs(AB_derivs) |

1950 | return DD |

1951 | #------------------------------------------------------------------------------- |

1952 | def _diff_prod(self,f_derivs,u,g,X,interval,end,uderivs,atc): |

1953 | r""" |

1954 | This function takes various derivatives of the equation `f=ug`, |

1955 | evaluates at a point `c`, and solves for the derivatives of `u`. |

1956 | For internal use by the function _asymptotics_main_multiple(). |

1957 | |

1958 | Does not use `self`. |

1959 | |

1960 | INPUT: |

1961 | |

1962 | - ``f_derivs`` - A dictionary whose keys are tuples \code{s + end} for all |

1963 | `X`-variable sequences `s` with length in `interval` and whose |

1964 | values are the derivatives of a function `f` evaluated at `c`. |

1965 | - ``u`` - A callable symbolic function. |

1966 | - ``g`` - An expression or callable symbolic function. |

1967 | - ``X`` - A list of symbolic variables. |

1968 | - ``interval`` - A list of positive integers. |

1969 | Call the first and last values `n` and `nn`, respectively. |

1970 | - ``end`` - A possibly empty list of `z`'s where `z` is the last element of |

1971 | `X`. |

1972 | - ``uderivs`` - A dictionary whose keys are the symbolic |

1973 | derivatives of order 0 to order `n-1` of `u` evaluated at `c` |

1974 | and whose values are the corresponding derivatives evaluated at `c`. |

1975 | - ``atc`` - A dictionary whose keys are the keys of `c` and all the symbolic |

1976 | derivatives of order 0 to order `nn` of `g` evaluated `c` and whose |

1977 | values are the corresponding derivatives evaluated at `c`. |

1978 | |

1979 | OUTPUT: |

1980 | |

1981 | A dictionary whose keys are the derivatives of `u` up to order |

1982 | `nn` and whose values are those derivatives evaluated at `c`. |

1983 | |

1984 | EXAMPLES:: |

1985 | |

1986 | I'd like to print out the entire dictionary, but that does not give |

1987 | consistent output for doctesting. |

1988 | Order of keys changes :: |

1989 | |

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

1991 | sage: G,H = 1,1 |

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

1993 | sage: u= function('u',x) |

1994 | sage: g= function('g',x) |

1995 | sage: dd= F._diff_prod({(x,):1,(x,x):1},u,g,[x],[1,2],[],{u(x=2):1},{x:2,g(x=2):3,diff(g,x)(x=2):5, diff(g,x,x)(x=2):7}) |

1996 | sage: dd[diff(u,x,2)(x=2)] |

1997 | 22/9 |

1998 | |

1999 | NOTES: |

2000 | |

2001 | This function works by differentiating the equation `f=ug` with respect |

2002 | to the variable sequence \code{s+end}, for all tuples `s` of `X` of |

2003 | lengths in `interval`, evaluating at the point `c`, |

2004 | and solving for the remaining derivatives of `u`. |

2005 | This function assumes that `u` never appears in the differentiations of |

2006 | `f=ug` after evaluating at `c`. |

2007 | |

2008 | AUTHORS: |

2009 | |

2010 | - Alex Raichev (2009-05-14, 2010-01-21) |

2011 | """ |

2012 | for l in interval: |

2013 | D= {} |

2014 | rhs= [] |

2015 | lhs= [] |

2016 | for t in UnorderedTuples(X,l): |

2017 | s= t+end |

2018 | lhs.append(f_derivs[tuple(s)]) |

2019 | rhs.append(diff(u*g,s).subs(atc).subs(uderivs)) |

2020 | # Since Sage's solve command can't take derivatives as variable |

2021 | # names, i make new variables based on t to stand in for |

2022 | # diff(u,t) and store them in D. |

2023 | D[diff(u,t).subs(atc)] = self._make_var([var('zing')]+t) |

2024 | eqns=[ lhs[i] == rhs[i].subs(uderivs).subs(D) for i in range(len(lhs))] |

2025 | vars= D.values() |

2026 | sol= solve(eqns,*vars,solution_dict=True) |

2027 | uderivs.update(self._subs_all(D,sol[0])) |

2028 | return uderivs |

2029 | #------------------------------------------------------------------------------- |

2030 | def _diff_seq(self,V,s): |

2031 | r""" |

2032 | Given a list `s` of tuples of natural numbers, this function returns the |

2033 | list of elements of `V` with indices the elements of the elements of `s`. |

2034 | This function is for internal use by the function _diff_op(). |

2035 | |

2036 | Does not use `self`. |

2037 | |

2038 | INPUT: |

2039 | |

2040 | - ``V`` - A list. |

2041 | - ``s`` - A list of tuples of natural numbers in the interval |

2042 | \code{range(len(V))}. |

2043 | |

2044 | OUTPUT: |

2045 | |

2046 | The tuple \code{tuple([V[tt] for tt in sorted(t)])}, where `t` is the |

2047 | list of elements of the elements of `s`. |

2048 | |

2049 | EXAMPLES:: |

2050 | |

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

2052 | sage: G,H = 1,1 |

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

2054 | sage: V= list(var('x,t,z')) |

2055 | sage: F._diff_seq(V,([0,1],[0,2,1],[0,0])) |

2056 | (x, x, x, x, t, t, z) |

2057 | |

2058 | AUTHORS: |

2059 | |

2060 | - Alex Raichev (2009.05.19) |

2061 | """ |

2062 | t= [] |

2063 | for ss in s: |

2064 | t.extend(ss) |

2065 | return tuple([V[tt] for tt in sorted(t)]) |

2066 | #------------------------------------------------------------------------------- |

2067 | def _direction(self,v,coordinate=None): |

2068 | r""" |

2069 | This function returns \code{[vv/v[coordinate] for vv in v]} where |

2070 | coordinate is the last index of v if not specified otherwise. |

2071 | |

2072 | Does not use `self`. |

2073 | |

2074 | INPUT: |

2075 | |

2076 | - ``v`` - A vector. |

2077 | - ``coordinate`` - An index for `v` (default: None). |

2078 | |

2079 | OUTPUT: |

2080 | |

2081 | This function returns \code{[vv/v[coordinate] for vv in v]} where |

2082 | coordinate is the last index of v if not specified otherwise. |

2083 | |

2084 | EXAMPLES:: |

2085 | |

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

2087 | sage: G,H = 1,1 |

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

2089 | sage: F._direction([2,3,5]) |

2090 | (2/5, 3/5, 1) |

2091 | sage: F._direction([2,3,5],0) |

2092 | (1, 3/2, 5/2) |

2093 | |

2094 | AUTHORS: |

2095 | |

2096 | - Alex Raichev (2010-08-25) |

2097 | """ |

2098 | if coordinate == None: |

2099 | coordinate= len(v)-1 |

2100 | try: |

2101 | v[0].variables() |

2102 | return tuple([(vv/v[coordinate]).simplify_full() for vv in v]) |

2103 | except: |

2104 | return tuple([vv/v[coordinate] for vv in v]) |

2105 | # E ============================================================================ |

2106 | # F ============================================================================ |

2107 | # G ============================================================================ |

2108 | # H ============================================================================ |

2109 | # I ============================================================================ |

2110 | def is_cmp(self,points): |

2111 | r""" |

2112 | Checks if the points in the list `points` are convenient multiple |

2113 | points of `V= \{ x\in CC^d : H(x) = 0\}`, where `H=self._H`. |

2114 | |

2115 | INPUT: |

2116 | |

2117 | - ``points`` - An individual or list of dictionaries with keys |

2118 | `self._variables` and values in some superfield of |

2119 | `self._R.base_ring()`. |

2120 | |

2121 | OUTPUT: |

2122 | |

2123 | A list of tuples `(p,verdict,comment)`, one for each point |

2124 | `p` in `points`, where `verdict` is True if `p` is a convenient |

2125 | multiple point and False otherwise, and where `comment` is a string |

2126 | comment relating to `verdict`, such as 'not a transverse intersection'. |

2127 | |

2128 | EXAMPLES:: |

2129 | |

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

2131 | sage: G= 16 |

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

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

2134 | sage: points= [{x:1/2,y:1,z:4/3},{x:1/2,y:1,z:2}] |

2135 | sage: F.is_cmp(points) |

2136 | [({y: 1, z: 4/3, x: 1/2}, True, 'all good'), ({y: 1, z: 2, x: 1/2}, False, 'not a singular point')] |

2137 | |

2138 | NOTES: |

2139 | |

2140 | A point `c` of `V` is a __convenient multiple point__ if `V` is locally |

2141 | a union of complex manifolds that intersect transversely at `c`; |

2142 | see [RaWi2011]_. |

2143 | |

2144 | REFERENCES: |

2145 | |

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

2147 | coefficients of multivariate generating functions: improvements |

2148 | for smooth points", submitted. |

2149 | |

2150 | AUTHORS: |

2151 | |

2152 | - Alex Raichev (2011-04-18) |

2153 | """ |

2154 | H= self._H |

2155 | d= self._d |

2156 | Hs= [SR(h[0]) for h in self._Hfac] # irreducible factors of H |

2157 | X= self._variables |

2158 | J= [tuple([diff(h,x) for x in X]) for h in Hs] |

2159 | verdicts= [] |

2160 | if not isinstance(points,list): |

2161 | points= [points] |

2162 | for p in points: |

2163 | # Ensure variables in points lie in SR. |

2164 | pp= {} |

2165 | for x in p.keys(): |

2166 | pp[SR(x)]= p[x] |

2167 | # Test 1: Is p a zero of all factors of H? |

2168 | if [h.subs(pp) for h in Hs] != [0 for h in Hs]: |

2169 | # Failed test 1. Move on to next point. |

2170 | verdicts.append((p,False,'not a singular point')) |

2171 | continue |

2172 | # Test 2: Are the factors of H smooth and |

2173 | # do theyintersect transversely at p? |

2174 | J= [tuple([f.subs(pp) for f in dh]) for dh in J] |

2175 | l= len(J) |

2176 | if Set(J).cardinality() < l: |

2177 | # Fail. Move on to next point. |

2178 | verdicts.append((p,False,'not a transverse intersection')) |

2179 | continue |

2180 | temp= True |

2181 | for S in list(Set(J).subsets())[1:]: # Subsets of size >= 1 |

2182 | k= len(list(S)) |

2183 | M = matrix(list(S)) |

2184 | if rank(M) != min(k,d): |

2185 | # Fail. |

2186 | temp= False |

2187 | verdicts.append((p,False,'not a transvere intersection')) |

2188 | break |

2189 | if not temp: |

2190 | # Move on to next point |

2191 | continue |

2192 | # Test 3: Is p convenient? |

2193 | Jlog= matrix([self._log_grad(h,X,pp) for h in Hs]) |

2194 | if [0 in f for f in Jlog.columns()] ==\ |

2195 | [True for f in Jlog.columns()]: |

2196 | # Fail. Move on to next point. |

2197 | verdict[p] = (False,'multiple point but not convenient') |

2198 | continue |

2199 | verdicts.append((p,True,'all good')) |

2200 | return verdicts |

2201 | # J ============================================================================ |

2202 | # K ============================================================================ |

2203 | # L ============================================================================ |

2204 | def _log_grad(self,f,X,c): |

2205 | r""" |

2206 | This function returns the logarithmic gradient of `f` with respect to the |

2207 | variables of `X` evalutated at `c`. |

2208 | |

2209 | Does not use `self`. |

2210 | |

2211 | INPUT: |

2212 | |

2213 | - ``f`` - An expression in the variables of `X`. |

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

2215 | - ``c`` - A dictionary with keys `X` and values in a field `K`. |

2216 | |

2217 | OUTPUT: |

2218 | |

2219 | \code{[c[x] * diff(f,x).subs(c) for x in X]}. |

2220 | |

2221 | EXAMPLES:: |

2222 | |

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

2224 | sage: G,H = 1,1 |

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

2226 | sage: X= var('x,y,z') |

2227 | sage: f= x*y*z^2 |

2228 | sage: c= {x:1,y:2,z:3} |

2229 | sage: f.gradient() |

2230 | (y*z^2, x*z^2, 2*x*y*z) |

2231 | sage: F._log_grad(f,X,c) |

2232 | (18, 18, 36) |

2233 | |

2234 | :: |

2235 | |

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

2237 | sage: f= x*y*z^2 |

2238 | sage: c= {x:1,y:2,z:3} |

2239 | sage: F._log_grad(f,R.gens(),c) |

2240 | (18, 18, 36) |

2241 | |

2242 | AUTHORS: |

2243 | |

2244 | - Alex Raichev (2009-03-06) |

2245 | """ |

2246 | return tuple([SR(c[x] * diff(f,x).subs(c)).simplify() for x in X]) |

2247 | # M ============================================================================ |

2248 | def _make_var(self,L): |

2249 | r""" |

2250 | This function converts the list `L` to a string (without commas) and returns |

2251 | the string as a variable. |

2252 | For internal use by the function _diff_op() |

2253 | |

2254 | Does not use `self`. |

2255 | |

2256 | INPUT: |

2257 | |

2258 | - ``L`` - A list. |

2259 | |

2260 | OUTPUT: |

2261 | |

2262 | A variable whose name is the concatenation of the variable names in `L`. |

2263 | |

2264 | EXAMPLES:: |

2265 | |

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

2267 | sage: G,H = 1,1 |

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

2269 | sage: L= list(var('x,y,hello')) |

2270 | sage: v= F._make_var(L) |

2271 | sage: print v, type(v) |

2272 | xyhello <type 'sage.symbolic.expression.Expression'> |

2273 | |

2274 | AUTHORS: |

2275 | |

2276 | - Alex Raichev (2010-01-21) |

2277 | """ |

2278 | return var(''.join([str(v) for v in L])) |

2279 | # N ============================================================================ |

2280 | def _new_var_name(self,name,V): |

2281 | r""" |

2282 | This function returns the first string in the sequence `name`, `name+name`, |

2283 | `name+name+name`,... that does not appear in the list `V`. |

2284 | It is for internal use by the function _asymptotics_main_multiple(). |

2285 | |

2286 | Does not use `self`. |

2287 | |

2288 | INPUT: |

2289 | |

2290 | - ``name`` - A string. |

2291 | - ``V`` - A list of variables. |

2292 | |

2293 | OUTPUT: |

2294 | |

2295 | The first string in the sequence `name`, `name+name`, |

2296 | `name+name+name`,... that does not appear in the list \code{str(V)}. |

2297 | |

2298 | EXAMPLES:: |

2299 | |

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

2301 | sage: G,H = 1,1 |

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

2303 | sage: X= var('x,xx,y,z') |

2304 | sage: F._new_var_name('x',X) |

2305 | 'xxx' |

2306 | |

2307 | AUTHORS: |

2308 | |

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

2310 | """ |

2311 | newname= name |

2312 | while newname in str(V): |

2313 | newname= newname +name |

2314 | return newname |

2315 | # O ============================================================================ |

2316 | # P ============================================================================ |

2317 | # Q ============================================================================ |

2318 | # R ============================================================================ |

2319 | def relative_error(self,approx,alpha,interval,exp_scale=1): |

2320 | r""" |

2321 | Returns the relative error between the values of the Maclaurin |

2322 | coefficients of `self` with multi-indices `m alpha` for `m` in |

2323 | `interval` and the values of the functions in `approx`. |

2324 | |

2325 | INPUT: |

2326 | |

2327 | - ``approx`` - An individual or list of symbolic expressions in |

2328 | one variable. |

2329 | - ``alpha`` - A list of positive integers. |

2330 | - ``interval`` - A list of positive integers. |

2331 | - ``exp_scale`` - (optional) A number. Default: 1. |

2332 | |

2333 | OUTPUT: |

2334 | |

2335 | A list whose entries are of the form |

2336 | `[m\alpha,a_m,b_m,err_m]` for `m \in interval`. |

2337 | Here `m\alpha` is a tuple; `a_m` is the `m alpha` (multi-index) |

2338 | coefficient of the Maclaurin series for `F` divided by `exp_scale^m`; |

2339 | `b_m` is a list of the values of the functions in `approx` evaluated at |

2340 | `m` and divided by `exp_scale^m`; `err_m` is the list of relative errors |

2341 | `(a_m-f)/a_m` for `f` in `b_m`. |

2342 | All outputs are decimal approximations. |

2343 | |

2344 | EXAMPLES:: |

2345 | |

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

2347 | sage: G=1 |

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

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

2350 | sage: alpha= [1,1] |

2351 | sage: var('n') |

2352 | n |

2353 | sage: f= (0.573/sqrt(n))*5.83^n |

2354 | sage: es= 5.83 |

2355 | sage: F.relative_error(f,alpha,[1,2,4,8],es) |

2356 | Calculating errors table in the form |

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

2358 | [[(1, 1), 0.514579759862779, [0.573000000000000], [-0.113530000000000]], [(2, 2), 0.382477808931739, [0.405172185619892], [-0.0593351461396876]], [(4, 4), 0.277863059517142, [0.286500000000000], [-0.0310834426780842]], [(8, 8), 0.199108827584423, [0.202586092809946], [-0.0174641439443390]]] |

2359 | sage: g= (0.573/sqrt(n) - 0.0674/n^(3/2))*5.83^n |

2360 | sage: F.relative_error([f,g],alpha,[1,2,4,8],es) |

2361 | Calculating errors table in the form |

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

2363 | [[(1, 1), 0.514579759862779, [0.573000000000000, 0.505600000000000], [-0.113530000000000, 0.0174506666666667]], [(2, 2), 0.382477808931739, [0.405172185619892, 0.381342687093905], [-0.0593351461396876, 0.00296781097184384]], [(4, 4), 0.277863059517142, [0.286500000000000, 0.278075000000000], [-0.0310834426780842, -0.000762751562681505]], [(8, 8), 0.199108827584423, [0.202586092809946, 0.199607405494198], [-0.0174641439443390, -0.00250404723800224]]] |

2364 | |

2365 | AUTHORS: |

2366 | |

2367 | - Alex Raichev (2009-05-18, 2011-04-18) |

2368 | """ |

2369 | |

2370 | if not isinstance(approx,list): |

2371 | approx= [approx] |

2372 | av= approx[0].variables()[0] |

2373 | |

2374 | print "Calculating errors table in the form" |

2375 | print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..." |

2376 | |

2377 | # Get Maclaurin coefficients of self and scale them. |

2378 | # Then compute errors. |

2379 | n= interval[-1] |

2380 | mc= self.maclaurin_coefficients(alpha,n) |

2381 | mca={} |

2382 | stats=[] |

2383 | for m in interval: |

2384 | beta= tuple([m*a for a in alpha]) |

2385 | mc[beta]= mc[beta]/exp_scale^m |

2386 | mca[beta]= [f.subs({av:m})/exp_scale^m for f in approx] |

2387 | stats_row= [beta, mc[beta].n(), [a.n() for a in mca[beta]]] |

2388 | if mc[beta]==0: |

2389 | stats_row.extend([None for a in mca[beta]]) |

2390 | else: |

2391 | stats_row.append([((mc[beta]-a)/mc[beta]).n() for a in mca[beta]]) |

2392 | stats.append(stats_row) |

2393 | return stats |

2394 | # S ============================================================================ |

2395 | def singular_points(self): |

2396 | r""" |

2397 | This function returns a Groebner basis ideal whose variety is the |

2398 | set of singular points of the algebraic variety |

2399 | `V= \{x\in\CC^d : H(x)=0\}`, where `H=sef._H`. |

2400 | |

2401 | INPUT: |

2402 | |

2403 | OUTPUT: |

2404 | |

2405 | A Groebner basis ideal whose variety is the set of singular points of |

2406 | the algebraic variety `V= \{x\in\CC^d : H(x)=0\}`. |

2407 | |

2408 | EXAMPLES:: |

2409 | |

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

2411 | sage: G= 1 |

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

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

2414 | sage: F.singular_points() |

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

2416 | |

2417 | AUTHORS: |

2418 | |

2419 | - Alex Raichev (2008-10-01, 2008-11-20, 2010-12-03, 2011-04-18) |

2420 | """ |

2421 | H= self._H |

2422 | R= self._R |

2423 | f= R.ideal(H).radical().gens()[0] # Compute the reduction of H. |

2424 | J= R.ideal([f] + f.gradient()) |

2425 | return R.ideal(J.groebner_basis()) |

2426 | #------------------------------------------------------------------------------- |

2427 | def smooth_critical(self,alpha): |

2428 | r""" |

2429 | This function returns a Groebner basis ideal whose variety is the set |

2430 | of smooth critical points of the algebraic variety |

2431 | `V= \{x\in\CC^d : H(x)=0\} for the direction `\alpha` where `H=self._H`. |

2432 | |

2433 | INPUT: |

2434 | |

2435 | - ``alpha`` - A `d`-tuple of positive integers and/or symbolic entries. |

2436 | |

2437 | OUTPUT: |

2438 | |

2439 | A Groebner basis ideal of smooth critical points of `V` for `\alpha`. |

2440 | |

2441 | EXAMPLES:: |

2442 | |

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

2444 | sage: G=1 |

2445 | sage: H = (1-x-y-x*y)^2 |

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

2447 | sage: var('a1,a2') |

2448 | (a1, a2) |

2449 | sage: F.smooth_critical([a1,a2]) |

2450 | Ideal (y^2 + 2*a1/a2*y - 1, x + (a2/(-a1))*y + (-a2 + a1)/(-a1)) of Multivariate Polynomial Ring in x, y over Fraction Field of Multivariate Polynomial Ring in a2, a1 over Rational Field |

2451 | |

2452 | NOTES: |

2453 | |

2454 | A point `c` of `V` is a __smooth critical point for `alpha`__ |

2455 | if the gradient of `f` at `c` is not identically zero and `\alpha` is in |

2456 | the span of the logarithmic gradient vector |

2457 | `(c[0] \partial_1 f(c)),\ldots,c[d-1] \partial_d f(c))`; see [RaWi2008a]_. |

2458 | |

2459 | REFERENCES: |

2460 | |

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

2462 | coefficients of multivariate generating functions: improvements |

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

2464 | (2008), R89. |

2465 | |

2466 | AUTHORS: |

2467 | |

2468 | - Alex Raichev (2008-10-01, 2008-11-20, 2009-03-09, 2010-12-02, 2011-04-18) |

2469 | """ |

2470 | H= self._H |

2471 | R= H.parent() |

2472 | B= R.base_ring() |

2473 | d= R.ngens() |

2474 | vars= R.gens() |

2475 | f= R.ideal(H).radical().gens()[0] # Compute the reduction of H. |

2476 | |

2477 | # Expand B by the variables of alpha if there are any. |

2478 | indets= [] |

2479 | indets_ind= [] |

2480 | for a in alpha: |

2481 | if not ((a in ZZ) and (a>0)): |

2482 | try: |

2483 | CC(a) |

2484 | except: |

2485 | indets.append(var(a)) |

2486 | indets_ind.append(alpha.index(a)) |

2487 | else: |

2488 | print "The components of", alpha, \ |

2489 | "must be positive integers or symbolic variables." |

2490 | return |

2491 | indets= list(Set(indets)) # Delete duplicates in indets. |

2492 | if indets != []: |

2493 | BB= FractionField(PolynomialRing(B,tuple(indets))) |

2494 | S= R.change_ring(BB) |

2495 | vars= S.gens() |

2496 | # Coerce alpha into BB. |

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

2498 | alpha[i] = BB(alpha[i]) |

2499 | else: |

2500 | S= R |

2501 | |

2502 | # Find smooth, critical points for alpha. |

2503 | f= S(f) |

2504 | J= S.ideal([f] +[ alpha[d-1]*vars[i]*diff(f,vars[i]) \ |

2505 | -alpha[i]*vars[d-1]*diff(f,vars[d-1]) for i in range(d-1)]) |

2506 | return S.ideal(J.groebner_basis()) |

2507 | #------------------------------------------------------------------------------- |

2508 | def _subs_all(self,f,sub,simplify=False): |

2509 | r""" |

2510 | This function returns the items of `f` substituted by the dictionaries of |

2511 | `sub` in order of their appearance in `sub`. |

2512 | |

2513 | Does not use `self`. |

2514 | |

2515 | INPUT: |

2516 | |

2517 | - ``f`` - An individual or list of symbolic expressions or dictionaries |

2518 | - ``sub`` - An individual or list of dictionaries. |

2519 | - ``simplify`` - Boolean (default: False). |

2520 | |

2521 | OUTPUT: |

2522 | |

2523 | The items of `f` substituted by the dictionaries of `sub` in order of |

2524 | their appearance in `sub`. The subs() command is used. |

2525 | If simplify is True, then simplify() is used after substitution. |

2526 | |

2527 | EXAMPLES:: |

2528 | |

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

2530 | sage: G,H = 1,1 |

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

2532 | sage: var('x,y,z') |

2533 | (x, y, z) |

2534 | sage: a= {x:1} |

2535 | sage: b= {y:2} |

2536 | sage: c= {z:3} |

2537 | sage: F._subs_all(x+y+z,a) |

2538 | y + z + 1 |

2539 | sage: F._subs_all(x+y+z,[c,a]) |

2540 | y + 4 |

2541 | sage: F._subs_all([x+y+z,y^2],b) |

2542 | [x + z + 2, 4] |

2543 | sage: F._subs_all([x+y+z,y^2],[b,c]) |

2544 | [x + 5, 4] |

2545 | |

2546 | :: |

2547 | |

2548 | sage: var('x,y') |

2549 | (x, y) |

2550 | sage: a= {'foo':x^2+y^2, 'bar':x-y} |

2551 | sage: b= {x:1,y:2} |

2552 | sage: F._subs_all(a,b) |

2553 | {'foo': 5, 'bar': -1} |

2554 | |

2555 | AUTHORS: |

2556 | |

2557 | - Alex Raichev (2009-05-05) |

2558 | """ |

2559 | singleton= False |

2560 | if not isinstance(f,list): |

2561 | f= [f] |

2562 | singleton= True |

2563 | if not isinstance(sub,list): |

2564 | sub= [sub] |

2565 | g= [] |

2566 | for ff in f: |

2567 | for D in sub: |

2568 | if isinstance(ff,dict): |

2569 | ff= dict( [(k,ff[k].subs(D)) for k in ff.keys()] ) |

2570 | else: |

2571 | ff= ff.subs(D) |

2572 | g.append(ff) |

2573 | if singleton and simplify: |

2574 | if isinstance(g[0],dict): |

2575 | return g[0] |

2576 | else: |

2577 | return g[0].simplify() |

2578 | elif singleton and not simplify: |

2579 | return g[0] |

2580 | elif not singleton and simplify: |

2581 | G= [] |

2582 | for gg in g: |

2583 | if isinstance(gg,dict): |

2584 | G.append(gg) |

2585 | else: |

2586 | G.append(gg.simplify()) |

2587 | return G |

2588 | else: |

2589 | return g |

2590 | # T ============================================================================ |

2591 | def maclaurin_coefficients(self,alpha,n): |

2592 | r""" |

2593 | Returns the Maclaurin coefficients of self that have multi-indices |

2594 | `alpha`, `2*alpha`,...,`n*alpha`. |

2595 | |

2596 | INPUT: |

2597 | |

2598 | - ``n`` - positive integer |

2599 | - ``alpha`` - tuple of positive integers representing a multi-index. |

2600 | |

2601 | OUTPUT: |

2602 | |

2603 | A dictionary of the form (beta, beta Maclaurin coefficient of self). |

2604 | |

2605 | AUTHORS: |

2606 | |

2607 | - Alex Raichev (2011-04-08) |

2608 | """ |

2609 | # Do all computations in the Symbolic Ring. |

2610 | F= SR(self._G/self._H) |

2611 | v= F.variables() |

2612 | p= {} |

2613 | for x in v: |

2614 | p[x]=0 |

2615 | d= len(v) |

2616 | coeffs={} |

2617 | # Initialize sequence of variables to differentiate with. |

2618 | s=[] |

2619 | for i in range(d): |

2620 | s.extend([v[i] for j in range(alpha[i])]) |

2621 | F_deriv= diff(F,s) |

2622 | coeffs[tuple(alpha)]= F_deriv.subs(p)/ mul([factorial(a) for a in alpha]) |

2623 | old_beta= alpha |

2624 | for k in [2..n]: |

2625 | # Update variable sequence to differentiate with. |

2626 | beta= [k*a for a in alpha] |

2627 | delta= [beta[i]-old_beta[i] for i in range(d)] |

2628 | s= [] |

2629 | for i in range(d): |

2630 | s.extend([v[i] for j in range(delta[i])]) |

2631 | F_deriv= diff(F_deriv,s) |

2632 | coeffs[tuple(beta)]= F_deriv.subs(p)/ mul([factorial(b) for b in beta]) |

2633 | old_beta= copy(beta) |

2634 | return coeffs |

2635 | # U ============================================================================ |

2636 | # V ============================================================================ |

2637 | def variables(self): |

2638 | r""" |

2639 | Returns the tuple of variables of `self`. |

2640 | |

2641 | EXAMPLES:: |

2642 | |

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

2644 | sage: G= exp(x) |

2645 | sage: H= 1-y |

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

2647 | sage: F.variables() |

2648 | (x, y) |

2649 | |

2650 | sage: R.<x,y>= PolynomialRing(QQ,order='invlex') |

2651 | sage: G= exp(x) |

2652 | sage: H= 1-y |

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

2654 | sage: F.variables() |

2655 | (y, x) |

2656 | |

2657 | AUTHORS: |

2658 | |

2659 | - Alex Raichev (2011-04-01) |

2660 | """ |

2661 | return self._variables |

2662 | # W ============================================================================ |

2663 | # X ============================================================================ |

2664 | # Y ============================================================================ |

2665 | # Z ============================================================================ |