| 1 | r""" |
| 2 | Asymptotics of Maclaurin coefficients of generating functions |
| 3 | |
| 4 | This code relates to analytic combinatorics. More specifically, it is a |
| 5 | collection of functions designed to compute asymptotics of Maclaurin |
| 6 | coefficients of certain classes of multivariate generating functions. |
| 7 | |
| 8 | The main function asymptotics() returns the first `N` terms of |
| 9 | the asymptotic expansion of the Maclaurin coefficients `F_{n\alpha}` |
| 10 | of the multivariate meromorphic function `F=G/H` as `n\to\infty`. |
| 11 | It assumes that `F` is holomorphic in a neighborhood of the origin, |
| 12 | that `H` is a polynomial, and that asymptotics in the direction of |
| 13 | `\alpha` (a tuple of positive integers) are controlled by convenient |
| 14 | smooth or multiple points. |
| 15 | |
| 16 | For an introduction to this subject, see [PeWi2008]_. |
| 17 | The main algorithms and formulas implemented here come from [RaWi2008a]_ |
| 18 | and [RaWi2011]_. |
| 19 | |
| 20 | REFERENCES: |
| 21 | |
| 22 | .. [AiYu1983] I. A. A\u\izenberg and A. P. Yuzhakov, "Integral |
| 23 | representations and residues in multidimensional complex analysis", |
| 24 | Translations of Mathematical Monographs, 58. American Mathematical |
| 25 | Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X. |
| 26 | |
| 27 | .. [DeLo2006] Wolfram Decker and Christoph Lossen, "Computing in |
| 28 | algebraic geometry", Chapter 7.1, Springer-Verlag, 2006. |
| 29 | |
| 30 | .. [DiEm2005] Alicia Dickenstein and Ioannis Z. Emiris (editors), |
| 31 | "Solving polynomial equations", Chapter 9.0, Springer-Verlag, 2005. |
| 32 | |
| 33 | .. [Lein1978] E. K. Leinartas, "On expansion of rational functions of |
| 34 | several variables into partial fractions", Soviet Math. (Iz. VUZ) |
| 35 | 22 (1978), no. 10, 35--38. |
| 36 | |
| 37 | .. [PeWi2008] Robin Pemantle and Mark C. Wilson, "Twenty combinatorial |
| 38 | examples of asymptotics derived from multivariate generating |
| 39 | functions", SIAM Rev. (2008) vol. 50 (2) pp. 199-272. |
| 40 | |
| 41 | .. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of |
| 42 | coefficients of multivariate generating functions: improvements |
| 43 | for smooth points", Electronic Journal of Combinatorics, Vol. 15 |
| 44 | (2008), R89. |
| 45 | |
| 46 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |
| 47 | coefficients of multivariate generating functions: improvements |
| 48 | for smooth points", To appear. |
| 49 | |
| 50 | AUTHORS: |
| 51 | |
| 52 | - Alex Raichev (2008-10-01) : Initial version |
| 53 | - Alex Raichev (2010-09-28) : Corrected many functions |
| 54 | - Alex Raichev (2010-12-15) : Updated documentation |
| 55 | - Alex Raichev (2011-03-09) : Fixed a division by zero bug in relative_error() |
| 56 | - Alex Raichev (2011-04-26) : Rewrote in object-oreinted style |
| 57 | - Alex Raichev (2011-05-06) : Fixed bug in cohomologous_integrand() and |
| 58 | fixed _crit_cone_combo() to work in SR |
| 59 | |
| 60 | EXAMPLES: |
| 61 | |
| 62 | These come from [RaWi2008a]_ and [RaWi2011]_. |
| 63 | A smooth point example:: |
| 64 | |
| 65 | sage: R.<x,y>= PolynomialRing(QQ) |
| 66 | sage: G= 1 |
| 67 | sage: H= 1-x-y-x*y |
| 68 | sage: F= QuasiRationalExpression(G,H) |
| 69 | sage: alpha= [3,2] |
| 70 | sage: F.maclaurin_coefficients(alpha,8) |
| 71 | {(21, 14): 1279919346549, (6, 4): 1289, (3, 2): 25, (18, 12): 19403906105, (24, 16): 85275509086721, (15, 10): 298199265, (9, 6): 75517, (12, 8): 4673345} |
| 72 | sage: I= F.smooth_critical(alpha); I |
| 73 | Ideal (y^2 + 3*y - 1, x - 2/3*y - 1/3) of Multivariate Polynomial Ring in x, y over Rational Field |
| 74 | sage: s= solve(I.gens(),F.variables(),solution_dict=true); s |
| 75 | [{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}] |
| 76 | sage: c= s[1] |
| 77 | sage: asys= [F.asymptotics(c,alpha,n,numerical=3) for n in [1..2]]; asys |
| 78 | Initializing auxiliary functions... |
| 79 | Calculating derivatives of auxiallary functions... |
| 80 | Calculating derivatives of more auxiliary functions... |
| 81 | Calculating actions of the second order differential operator... |
| 82 | Initializing auxiliary functions... |
| 83 | Calculating derivatives of auxiallary functions... |
| 84 | Calculating derivatives of more auxiliary functions... |
| 85 | Calculating actions of the second order differential operator... |
| 86 | [(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))] |
| 87 | sage: F.relative_error([f[0] for f in asys],alpha,[1,2,4,8,16],asys[0][1]) |
| 88 | Calculating errors table in the form |
| 89 | exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... |
| 90 | [[(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]]] |
| 91 | |
| 92 | Another smooth point example. Here it turns out that the terms 2--4 of the |
| 93 | asymptotic expansion are 0. |
| 94 | |
| 95 | :: |
| 96 | |
| 97 | sage: R.<x,y> = PolynomialRing(QQ) |
| 98 | sage: q= 1/2 |
| 99 | sage: qq= q.denominator() |
| 100 | sage: G= 1-q*x |
| 101 | sage: H= 1-q*x +q*x*y -x^2*y |
| 102 | sage: F= QuasiRationalExpression(G,H) |
| 103 | sage: alpha= list(qq*vector([2,1-q])); alpha |
| 104 | [4, 1] |
| 105 | sage: I= F.smooth_critical(alpha); I |
| 106 | Ideal (y^2 - 2*y + 1, x + 1/4*y - 5/4) of Multivariate Polynomial Ring in x, y over Rational Field |
| 107 | sage: s= solve(I.gens(),F.variables(),solution_dict=true); s |
| 108 | [{y: 1, x: 1}] |
| 109 | sage: c= s[0] |
| 110 | sage: asy= F.asymptotics(c,alpha,5); asy |
| 111 | Initializing auxiliary functions... |
| 112 | Calculating derivatives of auxiallary functions... |
| 113 | Calculating derivatives of more auxiliary functions... |
| 114 | Calculating actions of the second order differential operator... |
| 115 | (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))) |
| 116 | sage: F.relative_error(asy[0],alpha,[1,2],asy[1]) |
| 117 | Calculating errors table in the form |
| 118 | exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... |
| 119 | [[(4, 1), 0.187500000000000, [0.185581424449526], [0.0102324029358597]], [(8, 2), 0.152343750000000, [0.151986595969759], [0.00234439568568301]]] |
| 120 | |
| 121 | A multiple point example. |
| 122 | |
| 123 | :: |
| 124 | |
| 125 | sage: R.<x,y,z>= PolynomialRing(QQ) |
| 126 | sage: G= 1 |
| 127 | sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) |
| 128 | sage: F= QuasiRationalExpression(G,H) |
| 129 | sage: I= F.singular_points(); I |
| 130 | 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 |
| 131 | sage: c= {x: 1/2, y: 1, z: 4/3} |
| 132 | sage: F.is_cmp(c) |
| 133 | [({y: 1, z: 4/3, x: 1/2}, True, 'all good')] |
| 134 | sage: cc= F.critical_cone(c); cc |
| 135 | ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N) |
| 136 | sage: alpha= [8,3,3] |
| 137 | sage: alpha in cc[1] |
| 138 | True |
| 139 | sage: asy= F.asymptotics(c,alpha,2); asy |
| 140 | Initializing auxiliary functions... |
| 141 | Calculating derivatives of auxiliary functions... |
| 142 | Calculating derivatives of more auxiliary functions... |
| 143 | Calculating second-order differential operator actions... |
| 144 | (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))) |
| 145 | |
| 146 | Another multiple point example. |
| 147 | |
| 148 | :: |
| 149 | |
| 150 | sage: R.<x,y,z>= PolynomialRing(QQ) |
| 151 | sage: G= 16 |
| 152 | sage: H= (4-2*x-y-z)^2*(4-x-2*y-z) |
| 153 | sage: F= QuasiRationalExpression(G,H) |
| 154 | sage: var('a1,a2,a3') |
| 155 | (a1, a2, a3) |
| 156 | sage: F.cohomologous_integrand([a1,a2,a3]) |
| 157 | [[-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]]] |
| 158 | sage: I= F.singular_points(); I |
| 159 | Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field |
| 160 | sage: c= {x:1, y:1, z:1} |
| 161 | sage: F.is_cmp(c) |
| 162 | [({y: 1, z: 1, x: 1}, True, 'all good')] |
| 163 | sage: alpha= [3,3,2] |
| 164 | sage: cc= F.critical_cone(c); cc |
| 165 | ([(1, 2, 1), (2, 1, 1)], 2-d cone in 3-d lattice N) |
| 166 | sage: alpha in cc[1] |
| 167 | True |
| 168 | sage: asy= F.asymptotics(c,alpha,2); asy |
| 169 | Initializing auxiliary functions... |
| 170 | Calculating derivatives of auxiliary functions... |
| 171 | Calculating derivatives of more auxiliary functions... |
| 172 | Calculating second-order differential operator actions... |
| 173 | (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))) |
| 174 | |
| 175 | """ |
| 176 | #***************************************************************************** |
| 177 | # Copyright (C) 2008 Alexander Raichev <tortoise.said@gmail.com> |
| 178 | # |
| 179 | # Distributed under the terms of the GNU General Public License (GPL) |
| 180 | # http://www.gnu.org/licenses/ |
| 181 | #***************************************************************************** |
| 182 | |
| 183 | from copy import copy, deepcopy |
| 184 | |
| 185 | # Sage libraries |
| 186 | from sage.calculus.functional import diff |
| 187 | from sage.calculus.functions import jacobian |
| 188 | from sage.calculus.var import function, var |
| 189 | from sage.combinat.cartesian_product import CartesianProduct |
| 190 | from sage.combinat.combinat import stirling_number1 |
| 191 | from sage.combinat.permutation import Permutation |
| 192 | from sage.combinat.tuple import UnorderedTuples |
| 193 | from sage.functions.log import exp, log |
| 194 | from sage.functions.other import factorial, gamma, sqrt |
| 195 | from sage.geometry.cone import Cone |
| 196 | from sage.matrix.constructor import matrix |
| 197 | from sage.misc.misc import add, ellipsis_range, mul |
| 198 | from sage.misc.misc_c import prod |
| 199 | from sage.misc.mrange import cartesian_product_iterator, mrange |
| 200 | from sage.modules.free_module_element import vector |
| 201 | from sage.rings.arith import binomial |
| 202 | from sage.rings.all import CC |
| 203 | from sage.rings.fraction_field import FractionField |
| 204 | from sage.rings.integer import Integer |
| 205 | from sage.rings.integer_ring import ZZ |
| 206 | from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing |
| 207 | from sage.rings.rational_field import QQ |
| 208 | from sage.sets.set import Set |
| 209 | from sage.structure.sage_object import SageObject |
| 210 | from sage.symbolic.constants import pi |
| 211 | from sage.symbolic.relation import solve |
| 212 | from sage.symbolic.ring import SR |
| 213 | |
| 214 | |
| 215 | # Functions to incorporate later into existing Sage classes. |
| 216 | def algebraic_dependence(fs): |
| 217 | r""" |
| 218 | This function returns an irreducible annihilating polynomial for the |
| 219 | polynomials in `fs`, if those polynomials are algebraically dependent. |
| 220 | Otherwise it returns 0. |
| 221 | |
| 222 | INPUT: |
| 223 | |
| 224 | - ``fs`` - A list of polynomials `f_1,\ldots,f_r' from a common polynomial |
| 225 | ring. |
| 226 | |
| 227 | OUTPUT: |
| 228 | |
| 229 | If the polynomials in `fs` are algebraically dependent, then the |
| 230 | output is an irreducible polynomial `g` in `K[T_1,\ldots,T_r]` such that |
| 231 | `g(f_1,\ldots,f_r) = 0`. |
| 232 | Here `K` is the coefficient ring of self and `T_1,\ldots,T_r` are |
| 233 | indeterminates different from those of self. |
| 234 | If the polynomials in `fs` are algebraically independent, then the output |
| 235 | is the zero polynomial. |
| 236 | |
| 237 | EXAMPLES:: |
| 238 | |
| 239 | sage: from sage.combinat.amgf import * |
| 240 | sage: R.<x>= PolynomialRing(QQ) |
| 241 | sage: fs= [x-3, x^2+1] |
| 242 | sage: g= algebraic_dependence(fs); g |
| 243 | 10 + 6*T0 - T1 + T0^2 |
| 244 | sage: g(fs) |
| 245 | 0 |
| 246 | |
| 247 | :: |
| 248 | |
| 249 | sage: R.<w,z>= PolynomialRing(QQ) |
| 250 | sage: fs= [w, (w^2 + z^2 - 1)^2, w*z - 2] |
| 251 | sage: g= algebraic_dependence(fs); g |
| 252 | 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 |
| 253 | sage: g(fs) |
| 254 | 0 |
| 255 | sage: g.parent() |
| 256 | Multivariate Polynomial Ring in T0, T1, T2 over Rational Field |
| 257 | |
| 258 | :: |
| 259 | |
| 260 | sage: R.<x,y,z>= PolynomialRing(GF(7)) |
| 261 | sage: fs= [x*y,y*z] |
| 262 | sage: g= algebraic_dependence(fs); g |
| 263 | 0 |
| 264 | |
| 265 | AUTHORS: |
| 266 | |
| 267 | - Alex Raichev (2011-01-06) |
| 268 | """ |
| 269 | r= len(fs) |
| 270 | R= fs[Integer(0) ].parent() |
| 271 | B= R.base_ring() |
| 272 | Xs= list(R.gens()) |
| 273 | d= len(Xs) |
| 274 | |
| 275 | # Expand R by r new variables. |
| 276 | T= 'T' |
| 277 | while T in [str(x) for x in Xs]: |
| 278 | T= T+'T' |
| 279 | Ts= [T + str(j) for j in range(r)] |
| 280 | RR= PolynomialRing(B,d+r,tuple(Xs+Ts)) |
| 281 | Vs= list(RR.gens()) |
| 282 | Xs= Vs[Integer(0) :d] |
| 283 | Ts= Vs[d:] |
| 284 | |
| 285 | # Find an irreducible annihilating polynomial g for the fs if there is one. |
| 286 | J= RR.ideal([ Ts[j] -RR(fs[j]) for j in range(r)]) |
| 287 | JJ= J.elimination_ideal(Xs) |
| 288 | g= JJ.gens()[Integer(0) ] |
| 289 | |
| 290 | # Shrink the ambient ring of g to include only Ts. |
| 291 | # I chose the negdeglex order simply because i find it useful in my work. |
| 292 | RRR= PolynomialRing(B,r,tuple(Ts),order='negdeglex') |
| 293 | return RRR(g) |
| 294 | #------------------------------------------------------------------------------- |
| 295 | def partial_fraction_decomposition(f): |
| 296 | r""" |
| 297 | Return a partial fraction decomposition of the rational expression `f`. |
| 298 | Works for univariate and mulitivariate `f`. |
| 299 | |
| 300 | INPUT: |
| 301 | |
| 302 | - ``f`` - An element of the field of fractions `F` of a polynomial ring |
| 303 | `R` whose coefficients ring is a field. |
| 304 | In the univariate case, the coefficient ring doesn't have to be a field. |
| 305 | |
| 306 | OUTPUT: |
| 307 | |
| 308 | A tuple `whole,parts` where `whole \in R` and `parts` is the list of |
| 309 | terms (in `F`) of a partial fraction decomposition of `f - whole`. |
| 310 | See the notes below for more details. |
| 311 | |
| 312 | EXAMPLES:: |
| 313 | |
| 314 | sage: from sage.combinat.amgf import * |
| 315 | sage: S.<t> = QQ[] |
| 316 | sage: q = 1/(t+1) + 2/(t+2) + 3/(t-3) |
| 317 | sage: whole, parts = partial_fraction_decomposition(q); parts |
| 318 | [3/(t - 3), 1/(t + 1), 2/(t + 2)] |
| 319 | sage: whole +sum(parts) == q |
| 320 | True |
| 321 | |
| 322 | Notice that there is a whole part below despite the appearance of q :: |
| 323 | |
| 324 | sage: R.<x,y>= PolynomialRing(QQ) |
| 325 | sage: q= 1/( x *y *(x*y-1) ) |
| 326 | sage: whole,parts= partial_fraction_decomposition(q) |
| 327 | sage: whole, parts |
| 328 | (-1, [x*y/(x*y - 1), (-1)/(x*y)]) |
| 329 | sage: q == whole +sum(parts) |
| 330 | True |
| 331 | |
| 332 | :: |
| 333 | |
| 334 | sage: R.<x,y>= PolynomialRing(QQ) |
| 335 | sage: q= x +1/x +1/(x*y-2) +1/(x^2+y^2-1) |
| 336 | sage: whole,parts= partial_fraction_decomposition(q) |
| 337 | sage: whole, parts |
| 338 | (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)]) |
| 339 | sage: whole +sum(parts)==q |
| 340 | True |
| 341 | |
| 342 | :: |
| 343 | |
| 344 | sage: R.<x,y,z>= PolynomialRing(QQ) |
| 345 | sage: q= 1/x +1/(x*y-z-2)^2 +1/(x^2+y^2 +z^2-1)^3 +1/(x*y-3) |
| 346 | sage: whole,parts= partial_fraction_decomposition(q) |
| 347 | sage: whole +sum(parts)==q |
| 348 | True |
| 349 | |
| 350 | NOTES: |
| 351 | |
| 352 | In the case of univariate `f` this function calls the old univariate |
| 353 | partial fraction decomposition function. |
| 354 | In the multivariate case, it uses a different notion of and algorithm for |
| 355 | partial fraction decompositions. |
| 356 | |
| 357 | Let `f= P/Q` where `P,Q \in R`, let `Q_1^{e_1} \cdots Q_m^{e_m}` be the |
| 358 | unique factorization of `Q` in `R` into irreducible factors, and let `d` be |
| 359 | the number of indeterminates of `R`. |
| 360 | Then `f` can be written as a sum `\sum P_A/\prod_{j\in A} Q_j^{b_j}`, |
| 361 | where the `b_j \le e_j` are positive integers, the `P_A` are in `R`, and |
| 362 | the sum is taken over all size `\le m` subsets `A \subseteq \{ 1, \ldots, d \}` |
| 363 | such that `S:= \{ Q_j : j\in A \}` is algebraically independent and the |
| 364 | ideal generated by `S` is not all of `R` (that is, by the Nullstellensatz, |
| 365 | the polynomials of `S` have no common roots in the algebraic closure of the |
| 366 | coefficient field of `R`). |
| 367 | Following [Lein1978]_, i call any such decomposition of `f` a |
| 368 | `\emph{partial fraction decomposition}`. |
| 369 | |
| 370 | The algorithm used below comes from Theorem 1, Lemma 2, and Lemma 3 of |
| 371 | [Lein1978]_. |
| 372 | By the way, that paper has a typo in equation (c) on the |
| 373 | second page and equation (b) on the third page. |
| 374 | The right sides of (c) and (b) should be multiplied by `P`. |
| 375 | |
| 376 | REFERENCES: |
| 377 | |
| 378 | .. [Lein1978] E. K. Leinartas, "On expansion of rational functions of |
| 379 | several variables into partial fractions", Soviet Math. (Iz. VUZ) |
| 380 | 22 (1978), no. 10, 35--38. |
| 381 | |
| 382 | AUTHORS: |
| 383 | |
| 384 | - Alex Raichev (2011-01-10) |
| 385 | """ |
| 386 | R= f.denominator().parent() |
| 387 | d= len(R.gens()) |
| 388 | if d==Integer(1) : |
| 389 | return f.partial_fraction_decomposition() |
| 390 | Q= f.denominator() |
| 391 | whole,P= f.numerator().quo_rem(Q) |
| 392 | parts= [format_quotient(Integer(1) /Q)] |
| 393 | # Decompose via nullstellensatz trick |
| 394 | # (which is faster than the algebraic dependence trick) |
| 395 | parts= decompose_via_nullstellensatz(parts) |
| 396 | # Decompose via algebraic dependence trick |
| 397 | parts= decompose_via_algebraic_dependence(parts) |
| 398 | # Rewrite parts back in terms of rational expressions |
| 399 | new_parts=[] |
| 400 | for p in parts: |
| 401 | f= unformat_quotient(p) |
| 402 | if f.denominator() == Integer(1) : |
| 403 | whole= whole +f |
| 404 | else: |
| 405 | new_parts.append(f) |
| 406 | # Put P back in |
| 407 | new_parts= [P*f for f in new_parts] |
| 408 | return whole, new_parts |
| 409 | #------------------------------------------------------------------------------- |
| 410 | def combine_like_terms(parts,rationomial=True): |
| 411 | r""" |
| 412 | Combines like terms of the fractions represented by parts. |
| 413 | For use by partial_fraction_decomposition() above. |
| 414 | |
| 415 | EXAMPLES: |
| 416 | |
| 417 | sage: from sage.combinat.amgf import * |
| 418 | sage: R.<x,y>= PolynomialRing(QQ) |
| 419 | sage: parts =[[1, [x*y, 1]], [x, [x*y, 1]]] |
| 420 | sage: combine_like_terms(parts) |
| 421 | [[x + 1, [y, 1], [x, 1]]] |
| 422 | |
| 423 | :: |
| 424 | |
| 425 | sage: R.<x>= PolynomialRing(QQ) |
| 426 | sage: parts =[[1, [x, 1]], [x-1, [x, 1]]] |
| 427 | sage: combine_like_terms(parts) |
| 428 | [[1]] |
| 429 | |
| 430 | AUTHORS: |
| 431 | |
| 432 | - Alex Raichev (2011-01-10) |
| 433 | """ |
| 434 | |
| 435 | if parts == []: return parts |
| 436 | # Sort parts by denominators |
| 437 | new_parts= sorted(parts, key= lambda p:p[Integer(1) :]) |
| 438 | # Combine items of parts with same denominators. |
| 439 | newnew_parts=[] |
| 440 | left,right=Integer(0) ,Integer(1) |
| 441 | glom= new_parts[left][Integer(0) ] |
| 442 | while right <= len(new_parts)-Integer(1) : |
| 443 | if new_parts[left][Integer(1) :] == new_parts[right][Integer(1) :]: |
| 444 | glom= glom +new_parts[right][Integer(0) ] |
| 445 | else: |
| 446 | newnew_parts.append([glom]+new_parts[left][Integer(1) :]) |
| 447 | left= right |
| 448 | glom= new_parts[left][Integer(0) ] |
| 449 | right= right +Integer(1) |
| 450 | if glom != Integer(0) : |
| 451 | newnew_parts.append([glom]+new_parts[left][Integer(1) :]) |
| 452 | if rationomial: |
| 453 | # Reduce fractions in newnew_parts. |
| 454 | # Todo: speed up below by working directly with newnew_parts and |
| 455 | # thereby make fewer calls to format_quotient() which in turn |
| 456 | # calls factor(). |
| 457 | newnew_parts= [format_quotient(unformat_quotient(part)) for part in newnew_parts] |
| 458 | return newnew_parts |
| 459 | #------------------------------------------------------------------------------- |
| 460 | def decompose_via_algebraic_dependence(parts): |
| 461 | r""" |
| 462 | Returns a decomposition of parts. |
| 463 | Used by partial_fraction_decomposition() above. |
| 464 | Implements Lemma 2 of [Lein1978]_. |
| 465 | Recursive. |
| 466 | |
| 467 | REFERENCES: |
| 468 | |
| 469 | .. [Lein1978] E. K. Leinartas, "On expansion of rational functions of |
| 470 | several variables into partial fractions", Soviet Math. (Iz. VUZ) |
| 471 | 22 (1978), no. 10, 35--38. |
| 472 | |
| 473 | AUTHORS: |
| 474 | |
| 475 | - Alex Raichev (2011-01-10) |
| 476 | """ |
| 477 | decomposing_done= True |
| 478 | new_parts= [] |
| 479 | for p in parts: |
| 480 | p_parts= [p] |
| 481 | P= p[Integer(0) ] |
| 482 | Qs= p[Integer(1) :] |
| 483 | m= len(Qs) |
| 484 | G= algebraic_dependence([q for q,e in Qs]) |
| 485 | if G: |
| 486 | # Then the denominator factors are algebraically dependent |
| 487 | # and so we can decompose p. |
| 488 | decomposing_done= False |
| 489 | P= p[Integer(0) ] |
| 490 | Qs= p[Integer(1) :] |
| 491 | |
| 492 | # Todo: speed up step below by using |
| 493 | # G to calculate F. See [Lein1978]_ Lemma 1. |
| 494 | F= algebraic_dependence([q**e for q,e in Qs]) |
| 495 | new_vars= F.parent().gens() |
| 496 | |
| 497 | # Note that new_vars[j] corresponds to Qs[j] so that |
| 498 | # F([q^e for q,e in Qs]) = 0. |
| 499 | # Assuming below that F.parent() has negdeglex term order |
| 500 | # so that F.lt() is indeed the monomial we want. |
| 501 | FF= (F.lt() -F)/(F.lc()) |
| 502 | numers= map(mul,zip(FF.coefficients(),FF.monomials())) |
| 503 | e= list(F.lt().exponents())[Integer(0) :m] |
| 504 | denom= [[new_vars[j], e[Integer(0) ][j]+Integer(1) ] for j in range(m)] |
| 505 | |
| 506 | # Before making things messy by substituting in Qs, |
| 507 | # reduce terms and combine like terms. |
| 508 | p_parts_temp= [format_quotient(unformat_quotient([a]+denom)) for a in numers] |
| 509 | p_parts_temp= combine_like_terms(p_parts_temp) |
| 510 | |
| 511 | # Substitute Qs into new_p. |
| 512 | Qpowsub= dict([(new_vars[j],Qs[j][Integer(0) ]**Qs[j][Integer(1) ]) for j in range(m)]) |
| 513 | p_parts=[] |
| 514 | for x in p_parts_temp: |
| 515 | y= P*F.parent()(x[Integer(0) ]).subs(Qpowsub) |
| 516 | yy=[] |
| 517 | for xx in x[Integer(1) :]: |
| 518 | if xx[Integer(0) ] in new_vars: |
| 519 | j= new_vars.index(xx[Integer(0) ]) |
| 520 | yy.append([Qs[j][Integer(0) ],Qs[j][Integer(1) ]*xx[Integer(1) ]]) |
| 521 | else: |
| 522 | # Occasionally xx[0] is an integer. |
| 523 | yy.append(xx) |
| 524 | p_parts.append([y]+yy) |
| 525 | # Done one step of decomposing p. Add it to new_parts. |
| 526 | new_parts.extend(p_parts) |
| 527 | if decomposing_done: |
| 528 | return new_parts |
| 529 | else: |
| 530 | return decompose_via_algebraic_dependence(new_parts) |
| 531 | #------------------------------------------------------------------------------- |
| 532 | def decompose_via_cohomology(parts): |
| 533 | r""" |
| 534 | Given each nice (described below) differential form |
| 535 | `(P/Q) dx_1 \wedge\cdots\wedge dx_d` in `parts`, |
| 536 | this function returns a differential form equivalent in De Rham |
| 537 | cohomology that has no repeated factors in the denominator. |
| 538 | |
| 539 | INPUT: |
| 540 | |
| 541 | - ``parts`` - A list of the form `[chunk_1,\ldots,chunk_r]`, where each |
| 542 | `chunk_j` has the form `[P,[Q_1,e_1],\ldots,[Q_m,e_m]]`, |
| 543 | `Q_1,\ldots,Q_m` are irreducible elements of a common polynomial |
| 544 | ring `R` such that their corresponding algebraic varieties |
| 545 | `\{x\in F^d : B_j(x)=0\}` intersect transversely (where `F` is the |
| 546 | algebraic closure of the field of coefficients of `R`), |
| 547 | `e_1,\ldots,e_m` are positive integers, `m \le d`, and |
| 548 | `P` is a symbolic expression in some of the indeterminates of `R`. |
| 549 | Here `[P,[Q_1,e_1],\ldots,[Q_m,e_m]]` represents the fraction |
| 550 | `P/(Q_1^e_1 \cdots Q_m^e_m)`. |
| 551 | |
| 552 | OUTPUT: |
| 553 | |
| 554 | A list of the form `[chunky_1,\ldots,chunky_s]`, where each |
| 555 | `chunky_j` has the form `[P,[Q_1,1],\ldots,[Q_m,1]]`. |
| 556 | |
| 557 | EXAMPLES:: |
| 558 | |
| 559 | sage: from sage.combinat.amgf import * |
| 560 | sage: R.<x,y>= PolynomialRing(QQ) |
| 561 | sage: decompose_via_cohomology([[ 1, [x,3] ]]) |
| 562 | [] |
| 563 | |
| 564 | :: |
| 565 | |
| 566 | sage: R.<x,y>= PolynomialRing(QQ) |
| 567 | sage: decompose_via_cohomology([[ 1, [x,3], [x*y-1,2] ]]) |
| 568 | [[-3*y^2, [x, 1], [x*y - 1, 1]], [-5*y^3, [x*y - 1, 1]]] |
| 569 | |
| 570 | NOTES: |
| 571 | |
| 572 | This is a recursive function thats stops calling itself when all the |
| 573 | `e_j` equal 1 or `parts == []`. |
| 574 | The algorithm used here is that of Theorem 17.4 of |
| 575 | [AiYu1983]_. |
| 576 | The algebraic varieties `\{x\in F^d : Q_j(x)=0\}` |
| 577 | (where `F` is the algebraic closure of the field of coefficients of `R`) |
| 578 | corresponding to the `Q_j` __intersect transversely__ iff for each |
| 579 | point `c` of their intersection and for all `k \le m`, |
| 580 | the Jacobian matrix of any `k` polynomials from |
| 581 | `\{Q_1,\ldots,Q_m\}` has rank equal to `\min\{k,d\}` when evaluated at |
| 582 | `c`. |
| 583 | |
| 584 | REFERENCES: |
| 585 | |
| 586 | .. [AiYu1983] I. A. A\u\izenberg and A. P. Yuzhakov, "Integral |
| 587 | representations and residues in multidimensional complex analysis", |
| 588 | Translations of Mathematical Monographs, 58. American Mathematical |
| 589 | Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X. |
| 590 | |
| 591 | AUTHORS: |
| 592 | |
| 593 | - Alex Raichev (2008-10-01, 2011-01-15) |
| 594 | """ |
| 595 | if parts == []: return parts |
| 596 | decomposing_done= True |
| 597 | new_parts= [] |
| 598 | R= parts[Integer(0) ][Integer(1) ][Integer(0) ].parent() |
| 599 | V= list(R.gens()) |
| 600 | for p in parts: |
| 601 | p_parts= [p] |
| 602 | P= p[Integer(0) ] |
| 603 | Qs= p[Integer(1) :] |
| 604 | m= len(Qs) |
| 605 | if sum([e for q,e in Qs]) > m: |
| 606 | # Then we can decompose p |
| 607 | p_parts= [] |
| 608 | decomposing_done=False |
| 609 | dets= [] |
| 610 | vars_subsets= Set(V).subsets(m) |
| 611 | for v in vars_subsets: |
| 612 | # Sort variables so that first polynomial ring indeterminate |
| 613 | # declared is first in vars list. |
| 614 | v= sorted(v,reverse=True) |
| 615 | jac= jacobian([q for q,e in Qs],v) |
| 616 | dets.append(jac.determinant()) |
| 617 | # Get a Nullstellensatz certificate. |
| 618 | L= R(Integer(1) ).lift(R.ideal([q for q,e in Qs] +dets)) |
| 619 | for j in range(m): |
| 620 | if L[j] != Integer(0) : |
| 621 | # Make a copy of (and not a reference to) the nested list Qs. |
| 622 | new_Qs = deepcopy(Qs) |
| 623 | if new_Qs[j][Integer(1) ] > Integer(1) : |
| 624 | new_Qs[j][Integer(1) ]= new_Qs[j][Integer(1) ] -Integer(1) |
| 625 | else: |
| 626 | del new_Qs[j] |
| 627 | p_parts.append( [P*L[j]] +new_Qs ) |
| 628 | for k in range(vars_subsets.cardinality()): |
| 629 | if L[m+k] != Integer(0) : |
| 630 | new_Qs = deepcopy(Qs) |
| 631 | for j in range(m): |
| 632 | if new_Qs[j][Integer(1) ] > Integer(1) : |
| 633 | new_Qs[j][Integer(1) ]= new_Qs[j][Integer(1) ] -Integer(1) |
| 634 | v=sorted(vars_subsets[k],reverse=True) |
| 635 | jac= jacobian([SR(P*L[m+k])] +[ SR(Qs[jj][Integer(0) ]) for \ |
| 636 | jj in range(m) if jj !=j], [SR(vv) for vv in v]) |
| 637 | det= jac.determinant() |
| 638 | if det != Integer(0) : |
| 639 | p_parts.append([permutation_sign(v,V) \ |
| 640 | *(-Integer(1) )**j/new_Qs[j][Integer(1) ] *det] +new_Qs) |
| 641 | break |
| 642 | # Done one step of decomposing p. Add it to new_parts. |
| 643 | new_parts.extend(p_parts) |
| 644 | new_parts= combine_like_terms(new_parts,rationomial=False) |
| 645 | if decomposing_done: |
| 646 | return new_parts |
| 647 | else: |
| 648 | return decompose_via_cohomology(new_parts) |
| 649 | #------------------------------------------------------------------------------- |
| 650 | def decompose_via_nullstellensatz(parts): |
| 651 | r""" |
| 652 | Returns a decomposition of parts. |
| 653 | Used by partial_fraction_decomposition() above. |
| 654 | Implements Lemma 3 of [Lein1978]_. |
| 655 | Recursive. |
| 656 | |
| 657 | REFERENCES: |
| 658 | |
| 659 | .. [Lein1978] E. K. Leinartas, "On expansion of rational functions of |
| 660 | several variables into partial fractions", Soviet Math. (Iz. VUZ) |
| 661 | 22 (1978), no. 10, 35--38. |
| 662 | |
| 663 | AUTHORS: |
| 664 | |
| 665 | - Alex Raichev (2011-01-10) |
| 666 | """ |
| 667 | decomposing_done= True |
| 668 | new_parts= [] |
| 669 | R= parts[Integer(0) ][Integer(0) ].parent() |
| 670 | for p in parts: |
| 671 | p_parts= [p] |
| 672 | P= p[Integer(0) ] |
| 673 | Qs= p[Integer(1) :] |
| 674 | m= len(Qs) |
| 675 | if R(Integer(1) ) in R.ideal([q for q,e in Qs]): |
| 676 | # Then we can decompose p. |
| 677 | decomposing_done= False |
| 678 | L= R(Integer(1) ).lift(R.ideal([q**e for q,e in Qs])) |
| 679 | p_parts= [ [P*L[i]] + \ |
| 680 | [[Qs[j][Integer(0) ],Qs[j][Integer(1) ]] for j in range(m) if j != i] \ |
| 681 | for i in range(m) if L[i]!=Integer(0) ] |
| 682 | # The procedure above yields no like terms to combine. |
| 683 | # Done one step of decomposing p. Add it to new_parts. |
| 684 | new_parts.extend(p_parts) |
| 685 | if decomposing_done: |
| 686 | return new_parts |
| 687 | else: |
| 688 | return decompose_via_nullstellensatz(new_parts) |
| 689 | #------------------------------------------------------------------------------- |
| 690 | def format_quotient(f): |
| 691 | r""" |
| 692 | Formats `f` for use by partial_fraction_decomposition() above. |
| 693 | |
| 694 | AUTHORS: |
| 695 | |
| 696 | - Alex Raichev (2011-01-10) |
| 697 | """ |
| 698 | P= f.numerator() |
| 699 | Q= f.denominator() |
| 700 | Qs= Q.factor() |
| 701 | if Qs.unit() != Integer(1) : |
| 702 | P= P/Qs.unit() |
| 703 | Qs= sorted([[q,e] for q,e in Qs]) # sorting for future bookkeeping |
| 704 | return [P]+Qs |
| 705 | #------------------------------------------------------------------------------- |
| 706 | def permutation_sign(v,vars): |
| 707 | r""" |
| 708 | This function returns the sign of the permutation on `1,\ldots,len(vars)` |
| 709 | that is induced by the sublist `v` of `vars`. |
| 710 | For internal use by _cohom_equiv_main(). |
| 711 | |
| 712 | INPUT: |
| 713 | |
| 714 | - ``v`` - A sublist of `vars`. |
| 715 | - ``vars`` - A list of predefined variables or numbers. |
| 716 | |
| 717 | OUTPUT: |
| 718 | |
| 719 | The sign of the permutation obtained by taking indices (and adding 1) |
| 720 | within `vars` of the list `v,w`, where `w` is the list `vars` with the |
| 721 | elements of `v` removed. |
| 722 | |
| 723 | EXAMPLES:: |
| 724 | |
| 725 | sage: from sage.combinat.amgf import * |
| 726 | sage: vars= ['a','b','c','d','e'] |
| 727 | sage: v= ['b','d'] |
| 728 | sage: permutation_sign(v,vars) |
| 729 | -1 |
| 730 | sage: v= ['d','b'] |
| 731 | sage: permutation_sign(v,vars) |
| 732 | 1 |
| 733 | |
| 734 | AUTHORS: |
| 735 | |
| 736 | - Alex Raichev (2008-10-01) |
| 737 | """ |
| 738 | # Convert variable lists to lists of numbers in {1,...,len(vars)} |
| 739 | A= [x+Integer(1) for x in range(len(vars))] |
| 740 | B= [vars.index(x)+Integer(1) for x in v] |
| 741 | C= list(Set(A).difference(Set(B))) |
| 742 | C.sort() |
| 743 | P= Permutation(B+C) |
| 744 | return P.signature() |
| 745 | #------------------------------------------------------------------------------- |
| 746 | def unformat_quotient(part): |
| 747 | r""" |
| 748 | Unformats `f` for use by partial_fraction_decomposition() above. |
| 749 | Inverse of format_quotient() above. |
| 750 | |
| 751 | AUTHORS: |
| 752 | |
| 753 | - Alex Raichev (2011-01-10) |
| 754 | """ |
| 755 | P= part[Integer(0) ] |
| 756 | Qs= part[Integer(1) :] |
| 757 | Q= prod([q**e for q,e in Qs]) |
| 758 | return P/Q |
| 759 | #=============================================================================== |
| 760 | # Class for calculation of asymptotics of multivariate generating functions. |
| 761 | class QuasiRationalExpression(SageObject): |
| 762 | "Store an expression G/H, where H comes from a polynomial ring R and \ |
| 763 | G comes from R or the Symbolic Ring." |
| 764 | def __init__(self, G, H): |
| 765 | # Store important information about object as attributes of self. |
| 766 | # G, H, H's ring, ring dimension, H's factorization. |
| 767 | self._G = G |
| 768 | self._H = H |
| 769 | R= H.parent() |
| 770 | self._R = R |
| 771 | self._d = len(R.gens()) |
| 772 | self._Hfac= list(H.factor()) |
| 773 | |
| 774 | # Variables of self as elements of the SR. |
| 775 | # Remember that G might be in SR and not in R. |
| 776 | try: |
| 777 | # This fails if G is a constant, for example. |
| 778 | Gv= Set([R(x) for x in G.variables()]) |
| 779 | except: |
| 780 | Gv= Set([]) |
| 781 | try: |
| 782 | # This fails if H is a constant, for example. |
| 783 | Hv= Set(H.variables()) |
| 784 | except: |
| 785 | Hv= Set([]) |
| 786 | # Preserve the ring ordering of the variables which some methods below |
| 787 | # depends upon. |
| 788 | V= sorted(list(Gv.union(Hv)),reverse=True) |
| 789 | self._variables = tuple([SR(x) for x in V]) |
| 790 | #------------------------------------------------------------------------------- |
| 791 | # Keeping methods in alphabetical order (ignoring initial single underscores) |
| 792 | def asymptotics(self,c,alpha,N,numerical=Integer(0) ,asy_var=None): |
| 793 | r""" |
| 794 | This function returns the first `N` terms of the asymptotic expansion |
| 795 | of the Maclaurin coefficients `F_{n\alpha}` of the |
| 796 | multivariate meromorphic function `F=G/H` as `n\to\infty`, |
| 797 | where `F = self`. |
| 798 | It assumes that `F` is holomorphic in a neighborhood of the origin, |
| 799 | that `H` is a polynomial, and that `c` is a convenient strictly minimal |
| 800 | smooth or multiple of `F` that is critical for `\alpha`. |
| 801 | |
| 802 | INPUT: |
| 803 | |
| 804 | - ``alpha`` - A `d`-tuple of positive integers or, if `c` is a smooth |
| 805 | point, possibly of symbolic variables. |
| 806 | - ``c`` - A dictionary with keys `self._variables` and values from a |
| 807 | superfield of the field of `self._R.base_ring()`. |
| 808 | - ``N`` - A positive integer. |
| 809 | - ``numerical`` - A natural number (default: 0). |
| 810 | If k=numerical > 0, then a numerical approximation of the coefficients |
| 811 | of `F_{n\alpha}` with k digits of precision will be returned. |
| 812 | Otherwise exact values will be returned. |
| 813 | - ``asy_var`` - A symbolic variable (default: None). |
| 814 | The variable of the asymptotic expansion. |
| 815 | If none is given, `var('n')` will be assigned. |
| 816 | |
| 817 | OUTPUT: |
| 818 | |
| 819 | The tuple `(asy,exp_part,subexp_part)`, where `asy` is first `N` terms |
| 820 | of the asymptotic expansion of the Maclaurin coefficients `F_{n\alpha}` |
| 821 | of the function `F=self` as `n\to\infty`, `exp_part^n` is the exponential |
| 822 | factor of `asy`, and `subexp_part` is the subexponential factor of |
| 823 | `asy`. |
| 824 | |
| 825 | EXAMPLES:: |
| 826 | |
| 827 | A smooth point example :: |
| 828 | |
| 829 | sage: R.<x,y>= PolynomialRing(QQ) |
| 830 | sage: G= 1 |
| 831 | sage: H= 1-x-y-x*y |
| 832 | sage: F= QuasiRationalExpression(G,H) |
| 833 | sage: alpha= [3,2] |
| 834 | sage: c= {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3} |
| 835 | sage: F.asymptotics(c,alpha,2,numerical=3) |
| 836 | Initializing auxiliary functions... |
| 837 | Calculating derivatives of auxiallary functions... |
| 838 | Calculating derivatives of more auxiliary functions... |
| 839 | Calculating actions of the second order differential operator... |
| 840 | ((0.369/sqrt(n) - 0.0186/n^(3/2))*71.2^n, 71.2, 0.369/sqrt(n) - 0.0186/n^(3/2)) |
| 841 | |
| 842 | A multiple point example :: |
| 843 | |
| 844 | sage: R.<x,y,z>= PolynomialRing(QQ) |
| 845 | sage: G= 1 |
| 846 | sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) |
| 847 | sage: F= QuasiRationalExpression(G,H) |
| 848 | sage: c= {z: 4/3, y: 1, x: 1/2} |
| 849 | sage: alpha= [8,3,3] |
| 850 | sage: F.asymptotics(c,alpha,1) |
| 851 | Initializing auxiliary functions... |
| 852 | Calculating derivatives of auxiliary functions... |
| 853 | Calculating derivatives of more auxiliary functions... |
| 854 | Calculating second-order differential operator actions... |
| 855 | (1/7*sqrt(3)*sqrt(7)*108^n/(sqrt(pi)*sqrt(n)), 108, 1/7*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(n))) |
| 856 | |
| 857 | NOTES: |
| 858 | |
| 859 | A zero `c` of `H` is __strictly minimal__ if there is no zero `x` of `H` |
| 860 | such that `x_j < c_j` for all `0 \le j < d`. |
| 861 | For definitions of the terms "smooth critical point for `\alpha`" and |
| 862 | "multiple critical point for `\alpha`", |
| 863 | see the documentation for _asymptotics_main_smooth() and |
| 864 | _asymptotics_main_multiple(), which are the functions that do most of the |
| 865 | work. |
| 866 | |
| 867 | ALGORITHM: |
| 868 | |
| 869 | The algorithm used here comes from [RaWi2011]_. |
| 870 | |
| 871 | (1) Use Cauchy's multivariate integral formula to write `F_{n\alpha}` as |
| 872 | an integral around a polycirle centered at the origin of the |
| 873 | differential form `\frac{G(x) dx[0] \wedge\cdots\wedge |
| 874 | dx[d-1]}{H(x)x^\alpha}`. |
| 875 | |
| 876 | (2) Decompose `G/H` into a sum of partial fractions `P[0] +\cdots+ P[r]` |
| 877 | so that each term of the sum has at most `d` irreducible factors of `H` |
| 878 | in the denominator. |
| 879 | |
| 880 | (3) For each differential form `P[j] dx[0] \wedge\cdots\wedge dx[d-1]`, |
| 881 | find an equivalent form `\omega[j]` in de Rham cohomology with no |
| 882 | repeated irreducible factors of `H` in its denominator. |
| 883 | |
| 884 | (4) Compute an asymptotic expansion for each integral `\omega[j]`. |
| 885 | |
| 886 | (5) Add the expansions. |
| 887 | |
| 888 | REFERENCES: |
| 889 | |
| 890 | .. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of |
| 891 | coefficients of multivariate generating functions: improvements |
| 892 | for smooth points", Electronic Journal of Combinatorics, Vol. 15 |
| 893 | (2008), R89. |
| 894 | |
| 895 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |
| 896 | coefficients of multivariate generating functions: improvements |
| 897 | for smooth points", To appear. |
| 898 | |
| 899 | AUTHORS: |
| 900 | |
| 901 | - Alex Raichev (2008-10-01, 2010-09-28, 2011-04-27) |
| 902 | """ |
| 903 | # The variable for asymptotic expansions. |
| 904 | if not asy_var: |
| 905 | asy_var= var('n') |
| 906 | |
| 907 | # Create symbolic (non-ring) variables. |
| 908 | R= self._R |
| 909 | d= self._d |
| 910 | X= list(self._variables) |
| 911 | |
| 912 | # Do steps (1)--(3) |
| 913 | new_integrands= self.cohomologous_integrand(alpha,asy_var) |
| 914 | |
| 915 | # Coerce everything into the Symbolic Ring, as polynomial ring |
| 916 | # calculations are no longer needed. |
| 917 | # Calculate asymptotics. |
| 918 | cc={} |
| 919 | for k in c.keys(): |
| 920 | cc[SR(k)] = SR(c[k]) |
| 921 | c= cc |
| 922 | for i in range(len(alpha)): |
| 923 | alpha[i] = SR(alpha[i]) |
| 924 | subexp_parts= [] |
| 925 | for chunk in new_integrands: |
| 926 | # Convert chunk into Symbolic Ring |
| 927 | GG= SR(chunk[Integer(0) ]) |
| 928 | HH= [SR(f) for (f,e) in chunk[Integer(1) :]] |
| 929 | asy= self._asymptotics_main(GG,HH,X,c,asy_var,alpha,N,numerical) |
| 930 | subexp_parts.append(asy[Integer(2) ]) |
| 931 | exp_scale= asy[Integer(1) ] # same for all chunk in new_integrands |
| 932 | |
| 933 | # Do step (5). |
| 934 | subexp_part= add(subexp_parts) |
| 935 | return exp_scale**asy_var *subexp_part, exp_scale, subexp_part |
| 936 | #-------------------------------------------------------------------------------- |
| 937 | def _asymptotics_main(self,G,H,X,c,n,alpha,N,numerical): |
| 938 | r""" |
| 939 | This function is for internal use by asymptotics(). |
| 940 | It finds a variable in `X` to use to calculate asymptotics and decides |
| 941 | whether to call _asymptotics_main_smooth() or |
| 942 | _asymptotics_main_multiple(). |
| 943 | |
| 944 | Does not use `self`. |
| 945 | |
| 946 | INPUT: |
| 947 | |
| 948 | - ``G`` - A symbolic expression. |
| 949 | - ``H`` - A list of symbolic expressions. |
| 950 | - ``X`` - The list of variables occurring in `G` and `H`. |
| 951 | Call its length `d`. |
| 952 | - ``c`` - A dictionary with `X` as keys and numbers as values. |
| 953 | - ``n`` - The variable of the asymptotic expansion. |
| 954 | - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic |
| 955 | variables if `c` is a smooth point. |
| 956 | - ``N`` - A positive integer. |
| 957 | - ``numerical`` - Natural number. |
| 958 | If k=numerical > 0, then a numerical approximation of the coefficients |
| 959 | of `F_{n\alpha}` with k digits of precision will be returned. |
| 960 | Otherwise exact values will be returned. |
| 961 | |
| 962 | OUTPUT: |
| 963 | |
| 964 | The same as the function asymptotics(). |
| 965 | |
| 966 | AUTHORS: |
| 967 | |
| 968 | - Alex Raichev (2008-10-01, 2010-09-28) |
| 969 | """ |
| 970 | d= len(X) |
| 971 | r= len(H) # We know 1 <= r <= d. |
| 972 | |
| 973 | # Find a good variable x in X to do asymptotics calculations with, that is, |
| 974 | # a variable x in X such that for all h in H, diff(h,x).subs(c) != 0. |
| 975 | # A good variable is guaranteed to exist since we are dealing with |
| 976 | # convenient smooth or multiple points. |
| 977 | # Search for good x in X from back to front (to be consistent with |
| 978 | # [RaWi2008a]_ [RaWi2011]_ which uses X[d-1] as a good variable). |
| 979 | # Put each not good x found at the beginning of X and reshuffle alpha |
| 980 | # in the same way. |
| 981 | x= X[d-Integer(1) ] |
| 982 | beta= copy(alpha) |
| 983 | while Integer(0) in [diff(h,x).subs(c) for h in H]: |
| 984 | X.pop() |
| 985 | X.insert(Integer(0) ,x) |
| 986 | x= X[d-Integer(1) ] |
| 987 | a= beta.pop() |
| 988 | beta.insert(Integer(0) ,a) |
| 989 | if d==r: |
| 990 | # This is the case of a 'simple' multiple point. |
| 991 | A= G.subs(c) / jacobian(H,X).subs(c).determinant().abs() |
| 992 | return A,Integer(1) ,A |
| 993 | elif r==Integer(1) : |
| 994 | # So 1 = r < d, and we have a smooth point. |
| 995 | return self._asymptotics_main_smooth(G,H[Integer(0) ],X,c,n,beta,N,numerical) |
| 996 | else: |
| 997 | # So 1 < r < d, and we have a non-smooth multiple point. |
| 998 | return self._asymptotics_main_multiple(G,H,X,c,n,beta,N,numerical) |
| 999 | #-------------------------------------------------------------------------------- |
| 1000 | def _asymptotics_main_multiple(self,G,H,X,c,n,alpha,N,numerical): |
| 1001 | r""" |
| 1002 | This function is for internal use by _asymptotics_main(). |
| 1003 | It calculates asymptotics in case they depend upon multiple points. |
| 1004 | |
| 1005 | Does not use `self`. |
| 1006 | |
| 1007 | INPUT: |
| 1008 | |
| 1009 | - ``G`` - A symbolic expression. |
| 1010 | - ``H`` - A list of symbolic expressions. |
| 1011 | - ``X`` - The list of variables occurring in `G` and `H`. |
| 1012 | Call its length `d`. |
| 1013 | - ``c`` - A dictionary with `X` as keys and numbers as values. |
| 1014 | - ``n`` - The variable of the asymptotic expansion. |
| 1015 | - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic |
| 1016 | variables if `c` is a smooth point. |
| 1017 | - ``N`` - A positive integer. |
| 1018 | - ``numerical`` - Natural number. |
| 1019 | If k=numerical > 0, then a numerical approximation of the coefficients |
| 1020 | of `F_{n\alpha}` with k digits of precision will be returned. |
| 1021 | Otherwise exact values will be returned. |
| 1022 | |
| 1023 | OUTPUT: |
| 1024 | |
| 1025 | The same as the function asymptotics(). |
| 1026 | |
| 1027 | NOTES: |
| 1028 | |
| 1029 | The formula used for computing the asymptotic expansion is |
| 1030 | Theorem 3.4 of [RaWi2011]_. |
| 1031 | |
| 1032 | Currently this function cannot handle `c` with symbolic variable keys, |
| 1033 | because _crit_cone_combo() crashes. |
| 1034 | |
| 1035 | REFERENCES: |
| 1036 | |
| 1037 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |
| 1038 | coefficients of multivariate generating functions: improvements |
| 1039 | for smooth points", To appear. |
| 1040 | |
| 1041 | AUTHORS: |
| 1042 | |
| 1043 | - Alex Raichev (2008-10-01, 2010-09-28) |
| 1044 | """ |
| 1045 | I= sqrt(-Integer(1) ) |
| 1046 | d= len(X) # > r > 1 |
| 1047 | r= len(H) # > 1 |
| 1048 | C= copy(c) |
| 1049 | |
| 1050 | S= [var(self._new_var_name('s',X) + str(j)) for j in range(r-Integer(1) )] |
| 1051 | T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-Integer(1) )] |
| 1052 | Sstar= {} |
| 1053 | temp= self._crit_cone_combo(H,X,c,alpha) |
| 1054 | for j in range(r-Integer(1) ): |
| 1055 | Sstar[S[j]]= temp[j] |
| 1056 | thetastar= dict([(t,Integer(0) ) for t in T]) |
| 1057 | thetastar.update(Sstar) |
| 1058 | |
| 1059 | # Setup. |
| 1060 | print "Initializing auxiliary functions..." |
| 1061 | Hmul= mul(H) |
| 1062 | h= [function('h'+str(j),*tuple(X[:d-Integer(1) ])) for j in range(r)] # Implicit functions |
| 1063 | U = function('U',*tuple(X)) |
| 1064 | # All other functions are defined in terms of h, U, and explicit functions. |
| 1065 | Hcheck= mul([ X[d-Integer(1) ] -Integer(1) /h[j] for j in range(r)] ) |
| 1066 | Gcheck= -G/U *mul( [-h[j]/X[d-Integer(1) ] for j in range(r)] ) |
| 1067 | A= [(-Integer(1) )**(r-Integer(1) ) *X[d-Integer(1) ]**(-r+j)*diff(Gcheck.subs({X[d-Integer(1) ]:Integer(1) /X[d-Integer(1) ]}),X[d-Integer(1) ],j) for j in range(r)] |
| 1068 | e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-Integer(1) )]) |
| 1069 | ht= [hh.subs(e) for hh in h] |
| 1070 | Ht= [H[j].subs(e).subs({X[d-Integer(1) ]:Integer(1) /ht[j]}) for j in range(r)] |
| 1071 | hsumt= add([S[j]*ht[j] for j in range(r-Integer(1) )]) +(Integer(1) -add(S))*ht[r-Integer(1) ] |
| 1072 | At= [AA.subs(e).subs({X[d-Integer(1) ]:hsumt}) for AA in A] |
| 1073 | Phit = -log(C[X[d-Integer(1) ]] *hsumt)+ I*add([alpha[i]/alpha[d-Integer(1) ] *T[i] for i in range(d-Integer(1) )]) |
| 1074 | # atC Stores h and U and all their derivatives evaluated at C. |
| 1075 | atC = copy(C) |
| 1076 | atC.update(dict( [(hh.subs(C),Integer(1) /C[X[d-Integer(1) ]]) for hh in h ])) |
| 1077 | |
| 1078 | # Compute the derivatives of h up to order 2*N and evaluate at C. |
| 1079 | hderivs1= {} # First derivatives of h. |
| 1080 | for (i,j) in mrange([d-Integer(1) ,r]): |
| 1081 | s= solve( diff(H[j].subs({X[d-Integer(1) ]:Integer(1) /h[j]}),X[i]), diff(h[j],X[i]) )[Integer(0) ].rhs()\ |
| 1082 | .simplify() |
| 1083 | hderivs1.update({diff(h[j],X[i]):s}) |
| 1084 | atC.update({diff(h[j],X[i]).subs(C):s.subs(C).subs(atC)}) |
| 1085 | hderivs = self._diff_all(h,X[Integer(0) :d-Integer(1) ],Integer(2) *N,sub=hderivs1,rekey=h) |
| 1086 | for k in hderivs.keys(): |
| 1087 | atC.update({k.subs(C):hderivs[k].subs(atC)}) |
| 1088 | |
| 1089 | # Compute the derivatives of U up to order 2*N-2+ min{r,N}-1 and evaluate at C. |
| 1090 | # To do this, differentiate H = U*Hcheck over and over, evaluate at C, |
| 1091 | # and solve for the derivatives of U at C. |
| 1092 | # Need the derivatives of H with short keys to pass on to diff_prod later. |
| 1093 | m= min(r,N) |
| 1094 | end= [X[d-Integer(1) ] for j in range(r)] |
| 1095 | Hmulderivs= self._diff_all(Hmul,X,Integer(2) *N-Integer(2) +r,ending=end,sub_final=C) |
| 1096 | print "Calculating derivatives of auxiliary functions..." |
| 1097 | atC.update({U.subs(C):diff(Hmul,X[d-Integer(1) ],r).subs(C)/factorial(r)}) |
| 1098 | Uderivs={} |
| 1099 | p= Hmul.polynomial(CC).degree()-r |
| 1100 | if p == Integer(0) : |
| 1101 | # Then, using a proposition at the end of [RaWi2011], we can |
| 1102 | # conclude that all higher derivatives of U are zero. |
| 1103 | for l in (ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )): |
| 1104 | for s in UnorderedTuples(X,l): |
| 1105 | Uderivs[diff(U,s).subs(C)] = Integer(0) |
| 1106 | elif p > Integer(0) and p < Integer(2) *N-Integer(2) +m-Integer(1) : |
| 1107 | all_zero= True |
| 1108 | Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,p)),end,Uderivs,atC) |
| 1109 | # Check for a nonzero U derivative. |
| 1110 | if Uderivs.values() != [Integer(0) for i in range(len(Uderivs))]: |
| 1111 | all_zero= False |
| 1112 | if all_zero: |
| 1113 | # Then, using a proposition at the end of [RaWi2011], we can |
| 1114 | # conclude that all higher derivatives of U are zero. |
| 1115 | for l in (ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )): |
| 1116 | for s in UnorderedTuples(X,l): |
| 1117 | Uderivs.update({diff(U,s).subs(C):Integer(0) }) |
| 1118 | else: |
| 1119 | # Have to compute the rest of the derivatives. |
| 1120 | Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )),end,Uderivs,atC) |
| 1121 | else: |
| 1122 | Uderivs= self._diff_prod(Hmulderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N-Integer(2) +m-Integer(1) )),end,Uderivs,atC) |
| 1123 | atC.update(Uderivs) |
| 1124 | Phit1= jacobian(Phit,T+S).subs(hderivs1) |
| 1125 | a= jacobian(Phit1,T+S).subs(hderivs1).subs(thetastar).subs(atC) |
| 1126 | a_inv= a.inverse() |
| 1127 | Phitu= Phit -(Integer(1) /Integer(2) ) *matrix([T+S]) *a *matrix([T+S]).transpose() |
| 1128 | Phitu= Phitu[Integer(0) ][Integer(0) ] |
| 1129 | |
| 1130 | # Compute all partial derivatives of At and Phitu up to orders 2*N-2 |
| 1131 | # and 2*N, respectively. Take advantage of the fact that At and Phitu |
| 1132 | # are sufficiently differentiable functions so that mixed partials |
| 1133 | # are equal. Thus only need to compute representative partials. |
| 1134 | # Choose nondecreasing sequences as representative differentiation- |
| 1135 | # order sequences. |
| 1136 | # To speed up later computations, create symbolic functions AA and BB |
| 1137 | # to stand in for the expressions At and Phitu respectively. |
| 1138 | print "Calculating derivatives of more auxiliary functions..." |
| 1139 | AA= [function('A'+str(j),*tuple(T+S)) for j in range(r)] |
| 1140 | At_derivs= self._diff_all(At,T+S,Integer(2) *N-Integer(2) ,sub=hderivs1,sub_final=[thetastar,atC],rekey=AA) |
| 1141 | BB= function('BB',*tuple(T+S)) |
| 1142 | Phitu_derivs= self._diff_all(Phitu,T+S,Integer(2) *N,sub=hderivs1,sub_final=[thetastar,atC],rekey=BB,zero_order=Integer(3) ) |
| 1143 | AABB_derivs= At_derivs |
| 1144 | AABB_derivs.update(Phitu_derivs) |
| 1145 | for j in range(r): |
| 1146 | AABB_derivs[AA[j]] = At[j].subs(thetastar).subs(atC) |
| 1147 | AABB_derivs[BB] = Phitu.subs(thetastar).subs(atC) |
| 1148 | print "Calculating second-order differential operator actions..." |
| 1149 | DD= self._diff_op(AA,BB,AABB_derivs,T+S,a_inv,r,N) |
| 1150 | |
| 1151 | L={} |
| 1152 | for (j,k) in CartesianProduct((ellipsis_range(Integer(0) ,Ellipsis,min(r-Integer(1) ,N-Integer(1) ))), (ellipsis_range(max(Integer(0) ,N-Integer(1) -r),Ellipsis,N-Integer(1) ))): |
| 1153 | if j+k <= N-Integer(1) : |
| 1154 | L[(j,k)] = add([ \ |
| 1155 | DD[(j,k,l)] /( (-Integer(1) )**k *Integer(2) **(k+l) *factorial(l) *factorial(k+l) ) \ |
| 1156 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k))] ) |
| 1157 | # The next line's QQ coercion is a workaround for the Sage 4.6 bug reported |
| 1158 | # on http://trac.sagemath.org/sage_trac/ticket/8659. |
| 1159 | # Once the bug is fixed, the QQ can be removed. |
| 1160 | det= QQ(a.determinant())**(-Integer(1) /Integer(2) ) *(Integer(2) *pi)**((r-d)/Integer(2) ) |
| 1161 | chunk= det *add([ (alpha[d-Integer(1) ]*n)**((r-d)/Integer(2) -q) *add([ \ |
| 1162 | L[(j,k)] *binomial(r-Integer(1) ,j) *stirling_number1(r-j,r+k-q) *(-Integer(1) )**(q-j-k) \ |
| 1163 | for (j,k) in CartesianProduct((ellipsis_range(Integer(0) ,Ellipsis,min(r-Integer(1) ,q))), (ellipsis_range(max(Integer(0) ,q-r),Ellipsis,q))) if j+k <= q ]) \ |
| 1164 | for q in range(N)]) |
| 1165 | |
| 1166 | chunk= chunk.subs(C).simplify() |
| 1167 | coeffs= chunk.coefficients(n) |
| 1168 | coeffs.reverse() |
| 1169 | coeffs= coeffs[:N] |
| 1170 | if numerical: # If a numerical approximation is desired... |
| 1171 | subexp_part = add( [co[Integer(0) ].subs(c).n(digits=numerical)*n**co[Integer(1) ] \ |
| 1172 | for co in coeffs] ) |
| 1173 | exp_scale= (Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )) \ |
| 1174 | .n(digits=numerical) |
| 1175 | else: |
| 1176 | subexp_part = add( [co[Integer(0) ].subs(c)*n**co[Integer(1) ] for co in coeffs] ) |
| 1177 | exp_scale= Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] ) |
| 1178 | return exp_scale**n*subexp_part, exp_scale, subexp_part |
| 1179 | #-------------------------------------------------------------------------------- |
| 1180 | def _asymptotics_main_smooth(self,G,H,X,c,n,alpha,N,numerical): |
| 1181 | r""" |
| 1182 | This function is for internal use by _asymptotics_main(). |
| 1183 | It calculates asymptotics in case they depend upon smooth points. |
| 1184 | |
| 1185 | Does not use `self`. |
| 1186 | |
| 1187 | INPUT: |
| 1188 | |
| 1189 | - ``G`` - A symbolic expression. |
| 1190 | - ``H`` - A list of symbolic expressions. |
| 1191 | - ``X`` - The list of variables occurring in `G` and `H`. |
| 1192 | Call its length `d`. |
| 1193 | - ``c`` - A dictionary with `X` as keys and numbers as values. |
| 1194 | - ``n`` - The variable of the asymptotic expansion. |
| 1195 | - ``alpha`` - A `d`-tuple of positive natural numbers or possibly of symbolic |
| 1196 | variables if `c` is a smooth point. |
| 1197 | - ``N`` - A positive integer. |
| 1198 | - ``numerical`` - Natural number. |
| 1199 | If k=numerical > 0, then a numerical approximation of the coefficients |
| 1200 | of `F_{n\alpha}` with k digits of precision will be returned. |
| 1201 | Otherwise exact values will be returned. |
| 1202 | |
| 1203 | OUTPUT: |
| 1204 | |
| 1205 | The same as the function asymptotics(). |
| 1206 | |
| 1207 | NOTES: |
| 1208 | |
| 1209 | The formulas used for computing the asymptotic expansions are |
| 1210 | Theorems 3.2 and 3.3 [RaWi2008a]_ with `p=1`. |
| 1211 | Theorem 3.2 is a specialization of Theorem 3.4 of [RaWi2011]_ |
| 1212 | with `r=1`. |
| 1213 | |
| 1214 | REFERENCES: |
| 1215 | |
| 1216 | .. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of |
| 1217 | coefficients of multivariate generating functions: improvements |
| 1218 | for smooth points", Electronic Journal of Combinatorics, Vol. 15 |
| 1219 | (2008), R89. |
| 1220 | |
| 1221 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |
| 1222 | coefficients of multivariate generating functions: improvements |
| 1223 | for smooth points", To appear. |
| 1224 | |
| 1225 | AUTHORS: |
| 1226 | |
| 1227 | - Alex Raichev (2008-10-01, 2010-09-28) |
| 1228 | """ |
| 1229 | I= sqrt(-Integer(1) ) |
| 1230 | d= len(X) # > 1 |
| 1231 | |
| 1232 | # If c is a tuple of rationals, then compute with it directly. |
| 1233 | # Otherwise, compute symbolically and plug in c at the end. |
| 1234 | if vector(c.values()) in QQ**d: |
| 1235 | C= c |
| 1236 | else: |
| 1237 | Cs= [var('cs' +str(j)) for j in range(d)] |
| 1238 | C= dict( [(X[j],Cs[j]) for j in range(d)] ) |
| 1239 | c= dict( [(Cs[j],c[X[j]]) for j in range(d)] ) |
| 1240 | |
| 1241 | # Setup. |
| 1242 | print "Initializing auxiliary functions..." |
| 1243 | h= function('h',*tuple(X[:d-Integer(1) ])) # Implicit functions |
| 1244 | U = function('U',*tuple(X)) # |
| 1245 | # All other functions are defined in terms of h, U, and explicit functions. |
| 1246 | Gcheck = -G/U *(h/X[d-Integer(1) ]) |
| 1247 | A= Gcheck.subs({X[d-Integer(1) ]:Integer(1) /h})/h |
| 1248 | T= [var(self._new_var_name('t',X) + str(i)) for i in range(d-Integer(1) )] |
| 1249 | e= dict([(X[i],C[X[i]]*exp(I*T[i])) for i in range(d-Integer(1) )]) |
| 1250 | ht= h.subs(e) |
| 1251 | Ht= H.subs(e).subs({X[d-Integer(1) ]:Integer(1) /ht}) |
| 1252 | At= A.subs(e) |
| 1253 | Phit = -log(C[X[d-Integer(1) ]] *ht)\ |
| 1254 | + I* add([alpha[i]/alpha[d-Integer(1) ] *T[i] for i in range(d-Integer(1) )]) |
| 1255 | Tstar= dict([(t,Integer(0) ) for t in T]) |
| 1256 | atC = copy(C) |
| 1257 | atC.update({h.subs(C):Integer(1) /C[X[d-Integer(1) ]]}) # Stores h and U and all their derivatives |
| 1258 | # evaluated at C. |
| 1259 | |
| 1260 | # Compute the derivatives of h up to order 2*N and evaluate at C and store |
| 1261 | # in atC. Keep a copy of unevaluated h derivatives for use in the case |
| 1262 | # d=2 and v > 2 below. |
| 1263 | hderivs1= {} # First derivatives of h. |
| 1264 | for i in range(d-Integer(1) ): |
| 1265 | s= solve( diff(H.subs({X[d-Integer(1) ]:Integer(1) /h}),X[i]), diff(h,X[i]) )[Integer(0) ].rhs()\ |
| 1266 | .simplify() |
| 1267 | hderivs1.update({diff(h,X[i]):s}) |
| 1268 | atC.update({diff(h,X[i]).subs(C):s.subs(C).subs(atC)}) |
| 1269 | hderivs = self._diff_all(h,X[Integer(0) :d-Integer(1) ],Integer(2) *N,sub=hderivs1,rekey=h) |
| 1270 | for k in hderivs.keys(): |
| 1271 | atC.update({k.subs(C):hderivs[k].subs(atC)}) |
| 1272 | |
| 1273 | # Compute the derivatives of U up to order 2*N and evaluate at C. |
| 1274 | # To do this, differentiate H = U*Hcheck over and over, evaluate at C, |
| 1275 | # and solve for the derivatives of U at C. |
| 1276 | # Need the derivatives of H with short keys to pass on to diff_prod later. |
| 1277 | Hderivs= self._diff_all(H,X,Integer(2) *N,ending=[X[d-Integer(1) ]],sub_final=C) |
| 1278 | print "Calculating derivatives of auxiallary functions..." |
| 1279 | # For convenience in checking if all the nontrivial derivatives of U at c |
| 1280 | # are zero a few line below, store the value of U(c) in atC instead of in |
| 1281 | # Uderivs. |
| 1282 | Uderivs={} |
| 1283 | atC.update({U.subs(C):diff(H,X[d-Integer(1) ]).subs(C)}) |
| 1284 | end= [X[d-Integer(1) ]] |
| 1285 | Hcheck= X[d-Integer(1) ] - Integer(1) /h |
| 1286 | p= H.polynomial(CC).degree()-Integer(1) |
| 1287 | if p == Integer(0) : |
| 1288 | # Then, using a proposition at the end of [RaWi2011], we can |
| 1289 | # conclude that all higher derivatives of U are zero. |
| 1290 | for l in (ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N)): |
| 1291 | for s in UnorderedTuples(X,l): |
| 1292 | Uderivs[diff(U,s).subs(C)] = Integer(0) |
| 1293 | elif p > Integer(0) and p < Integer(2) *N: |
| 1294 | all_zero= True |
| 1295 | Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,p)),end,Uderivs,atC) |
| 1296 | # Check for a nonzero U derivative. |
| 1297 | if Uderivs.values() != [Integer(0) for i in range(len(Uderivs))]: |
| 1298 | all_zero= False |
| 1299 | if all_zero: |
| 1300 | # Then, using a proposition at the end of [RaWi2011], we can |
| 1301 | # conclude that all higher derivatives of U are zero. |
| 1302 | for l in (ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N)): |
| 1303 | for s in UnorderedTuples(X,l): |
| 1304 | Uderivs.update({diff(U,s).subs(C):Integer(0) }) |
| 1305 | else: |
| 1306 | # Have to compute the rest of the derivatives. |
| 1307 | Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(p+Integer(1) ,Ellipsis,Integer(2) *N)),end,Uderivs,atC) |
| 1308 | else: |
| 1309 | Uderivs= self._diff_prod(Hderivs,U,Hcheck,X,(ellipsis_range(Integer(1) ,Ellipsis,Integer(2) *N)),end,Uderivs,atC) |
| 1310 | atC.update(Uderivs) |
| 1311 | |
| 1312 | # In general, this algorithm is not designed to handle the case of a |
| 1313 | # singular Phit''(Tstar). However, when d=2 the algorithm can cope. |
| 1314 | if d==Integer(2) : |
| 1315 | # Compute v, the order of vanishing at Tstar of Phit. It is at least 2. |
| 1316 | v=Integer(2) |
| 1317 | Phitderiv= diff(Phit,T[Integer(0) ],Integer(2) ) |
| 1318 | splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify() |
| 1319 | while splat==Integer(0) : |
| 1320 | v= v+Integer(1) |
| 1321 | if v > Integer(2) *N: # Then need to compute more derivatives of h for atC. |
| 1322 | hderivs.update({diff(h,X[Integer(0) ],v) \ |
| 1323 | :diff(hderivs[diff(h,X[Integer(0) ],v-Integer(1) )],X[Integer(0) ]).subs(hderivs1)}) |
| 1324 | atC.update({diff(h,X[Integer(0) ],v).subs(C) \ |
| 1325 | :hderivs[diff(h,X[Integer(0) ],v)].subs(atC)}) |
| 1326 | Phitderiv= diff(Phitderiv,T[Integer(0) ]) |
| 1327 | splat= Phitderiv.subs(Tstar).subs(atC).subs(c).simplify() |
| 1328 | if d==Integer(2) and v>Integer(2) : |
| 1329 | t= T[Integer(0) ] # Simplify variable names. |
| 1330 | a= splat/factorial(v) |
| 1331 | Phitu= Phit -a*t**v |
| 1332 | |
| 1333 | # Compute all partial derivatives of At and Phitu up to orders 2*(N-1) |
| 1334 | # and 2*(N-1)+v, respectively, in case v is even. |
| 1335 | # Otherwise, compute up to orders N-1 and N-1+v, respectively. |
| 1336 | # To speed up later computations, create symbolic functions AA and BB |
| 1337 | # to stand in for the expressions At and Phitu, respectively. |
| 1338 | print "Calculating derivatives of more auxiliary functions..." |
| 1339 | AA= function('AA',t) |
| 1340 | BB= function('BB',t) |
| 1341 | if v.mod(Integer(2) )==Integer(0) : |
| 1342 | At_derivs= self._diff_all(At,T,Integer(2) *N-Integer(2) , \ |
| 1343 | sub=hderivs1,sub_final=[Tstar,atC],rekey=AA) |
| 1344 | Phitu_derivs= self._diff_all(Phitu,T,Integer(2) *N-Integer(2) +v, \ |
| 1345 | sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+Integer(1) ,rekey=BB) |
| 1346 | else: |
| 1347 | At_derivs= self._diff_all(At,T,N-Integer(1) ,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA) |
| 1348 | Phitu_derivs= self._diff_all(Phitu,T,N-Integer(1) +v,sub=hderivs1,sub_final=[Tstar,atC],zero_order=v+Integer(1) ,rekey=BB) |
| 1349 | AABB_derivs= At_derivs |
| 1350 | AABB_derivs.update(Phitu_derivs) |
| 1351 | AABB_derivs[AA] = At.subs(Tstar).subs(atC) |
| 1352 | AABB_derivs[BB] = Phitu.subs(Tstar).subs(atC) |
| 1353 | print "Calculating actions of the second order differential operator..." |
| 1354 | DD= self._diff_op_simple(AA,BB,AABB_derivs,t,v,a,N) |
| 1355 | # Plug above into asymptotic formula. |
| 1356 | L = [] |
| 1357 | if v.mod(Integer(2) ) == Integer(0) : |
| 1358 | for k in range(N): |
| 1359 | L.append( add([ \ |
| 1360 | (-Integer(1) )**l *gamma((Integer(2) *k+v*l+Integer(1) )/v) \ |
| 1361 | / (factorial(l) *factorial(Integer(2) *k+v*l)) \ |
| 1362 | * DD[(k,l)] for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)) ]) ) |
| 1363 | chunk= a**(-Integer(1) /v) /(pi*v) *add([ alpha[d-Integer(1) ]**(-(Integer(2) *k+Integer(1) )/v) \ |
| 1364 | * L[k] *n**(-(Integer(2) *k+Integer(1) )/v) for k in range(N) ]) |
| 1365 | else: |
| 1366 | zeta= exp(I*pi/(Integer(2) *v)) |
| 1367 | for k in range(N): |
| 1368 | L.append( add([ \ |
| 1369 | (-Integer(1) )**l *gamma((k+v*l+Integer(1) )/v) \ |
| 1370 | / (factorial(l) *factorial(k+v*l)) \ |
| 1371 | * (zeta**(k+v*l+Integer(1) ) +(-Integer(1) )**(k+v*l)*zeta**(-(k+v*l+Integer(1) ))) \ |
| 1372 | * DD[(k,l)] for l in (ellipsis_range(Integer(0) ,Ellipsis,k)) ]) ) |
| 1373 | chunk= abs(a)**(-Integer(1) /v) /(Integer(2) *pi*v) *add([ alpha[d-Integer(1) ]**(-(k+Integer(1) )/v) \ |
| 1374 | * L[k] *n**(-(k+Integer(1) )/v) for k in range(N) ]) |
| 1375 | # Asymptotics for d>=2 case. A singular Phit''(Tstar) will cause a crash |
| 1376 | # in this case. |
| 1377 | else: |
| 1378 | Phit1= jacobian(Phit,T).subs(hderivs1) |
| 1379 | a= jacobian(Phit1,T).subs(hderivs1).subs(Tstar).subs(atC) |
| 1380 | a_inv= a.inverse() |
| 1381 | Phitu= Phit -(Integer(1) /Integer(2) ) *matrix([T]) *a *matrix([T]).transpose() |
| 1382 | Phitu= Phitu[Integer(0) ][Integer(0) ] |
| 1383 | # Compute all partial derivatives of At and Phitu up to orders 2*N-2 |
| 1384 | # and 2*N, respectively. Take advantage of the fact that At and Phitu |
| 1385 | # are sufficiently differentiable functions so that mixed partials |
| 1386 | # are equal. Thus only need to compute representative partials. |
| 1387 | # Choose nondecreasing sequences as representative differentiation- |
| 1388 | # order sequences. |
| 1389 | # To speed up later computations, create symbolic functions AA and BB |
| 1390 | # to stand in for the expressions At and Phitu respectively. |
| 1391 | print "Calculating derivatives of more auxiliary functions..." |
| 1392 | AA= function('AA',*tuple(T)) |
| 1393 | At_derivs= self._diff_all(At,T,Integer(2) *N-Integer(2) ,sub=hderivs1,sub_final=[Tstar,atC],rekey=AA) |
| 1394 | BB= function('BB',*tuple(T)) |
| 1395 | Phitu_derivs= self._diff_all(Phitu,T,Integer(2) *N,sub=hderivs1,sub_final=[Tstar,atC],rekey=BB,zero_order=Integer(3) ) |
| 1396 | AABB_derivs= At_derivs |
| 1397 | AABB_derivs.update(Phitu_derivs) |
| 1398 | AABB_derivs[AA] = At.subs(Tstar).subs(atC) |
| 1399 | AABB_derivs[BB] = Phitu.subs(Tstar).subs(atC) |
| 1400 | print "Calculating actions of the second order differential operator..." |
| 1401 | DD= self._diff_op(AA,BB,AABB_derivs,T,a_inv,Integer(1) ,N) |
| 1402 | |
| 1403 | # Plug above into asymptotic formula. |
| 1404 | L=[] |
| 1405 | for k in range(N): |
| 1406 | L.append( add([ \ |
| 1407 | DD[(Integer(0) ,k,l)] / ( (-Integer(1) )**k *Integer(2) **(l+k) *factorial(l) *factorial(l+k) ) \ |
| 1408 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k))]) ) |
| 1409 | chunk= add([ (Integer(2) *pi)**((Integer(1) -d)/Integer(2) ) *a.determinant()**(-Integer(1) /Integer(2) ) \ |
| 1410 | *alpha[d-Integer(1) ]**((Integer(1) -d)/Integer(2) -k) *L[k] \ |
| 1411 | *n**((Integer(1) -d)/Integer(2) -k) for k in range(N) ]) |
| 1412 | |
| 1413 | chunk= chunk.subs(c).simplify() |
| 1414 | coeffs= chunk.coefficients(n) |
| 1415 | coeffs.reverse() |
| 1416 | coeffs= coeffs[:N] |
| 1417 | if numerical: # If a numerical approximation is desired... |
| 1418 | subexp_part = add( [co[Integer(0) ].subs(c).n(digits=numerical)*n**co[Integer(1) ] for co in coeffs] ) |
| 1419 | exp_scale=(Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] )) \ |
| 1420 | .n(digits=numerical) |
| 1421 | else: |
| 1422 | subexp_part = add( [co[Integer(0) ].subs(c)*n**co[Integer(1) ] for co in coeffs] ) |
| 1423 | exp_scale= Integer(1) /mul( [(C[X[i]]**alpha[i]).subs(c) for i in range(d)] ) |
| 1424 | return exp_scale**n*subexp_part, exp_scale, subexp_part |
| 1425 | #----------------------------------------------------------------------------- |
| 1426 | def _crit_cone_combo(self,fs,X,c,alpha): |
| 1427 | r""" |
| 1428 | This function returns an auxiliary point associated to the multiple |
| 1429 | point `c` of the factors `fs`. |
| 1430 | It is for internal use by _asymptotics_main_multiple(). |
| 1431 | |
| 1432 | INPUT: |
| 1433 | |
| 1434 | - ``fs`` - A list of expressions in the variables of `X`. |
| 1435 | - ``X`` - A list of variables. |
| 1436 | - ``c`` - A dictionary with keys `X` and values in some field. |
| 1437 | - ``alpha`` - A list of rationals. |
| 1438 | |
| 1439 | OUTPUT: |
| 1440 | |
| 1441 | A solution of the matrix equation `y Gamma = a` for `y`, |
| 1442 | where `Gamma` is the matrix whose `j`th row is |
| 1443 | _direction(_log_grad(fj,X,c)) where `fj` |
| 1444 | is the `j`th item in `fs` and where `a` is _direction(alpha). |
| 1445 | |
| 1446 | EXAMPLES:: |
| 1447 | |
| 1448 | sage: R.<x,y,z>= PolynomialRing(QQ) |
| 1449 | sage: G,H= 1,1 |
| 1450 | sage: F= QuasiRationalExpression(G,H) |
| 1451 | sage: fs= [x + 2*y + z - 4, 2*x + y + z - 4] |
| 1452 | sage: c= {x:1,y:1,z:1} |
| 1453 | sage: alpha= [2,1,1] |
| 1454 | sage: F._crit_cone_combo(fs,R.gens(),c,alpha) |
| 1455 | [0, 1] |
| 1456 | |
| 1457 | NOTES: |
| 1458 | |
| 1459 | Use this function only when `Gamma` is well-defined and |
| 1460 | there is a unique solution to the matrix equation `y Gamma = a`. |
| 1461 | Fails otherwise. |
| 1462 | |
| 1463 | AUTHORS: |
| 1464 | |
| 1465 | - Alex Raichev (2008-10-01, 2008-11-25, 2009-03-04, 2010-09-08, |
| 1466 | 2010-12-02) |
| 1467 | """ |
| 1468 | # Assuming here that each _log_grad(f) has nonzero final component. |
| 1469 | # Then 'direction' will not throw a division by zero error. |
| 1470 | d= len(X) |
| 1471 | r= len(fs) |
| 1472 | Gamma= matrix([self._direction(self._log_grad(f,X,c)) for f in fs]) |
| 1473 | # solve_left() fails when working in SR :-(. So use solve() instead. |
| 1474 | #s= Gamma.solve_left(vector(alpha)/alpha[d-1]) |
| 1475 | V= [var('sss'+str(i)) for i in range(r)] |
| 1476 | M= matrix(V)*Gamma |
| 1477 | eqns= [M[Integer(0) ][i]== alpha[i]/alpha[d-Integer(1) ] for i in range(d)] |
| 1478 | s= solve(eqns,V,solution_dict=True)[Integer(0) ] # Assuming a unique solution. |
| 1479 | return [s[v] for v in V] |
| 1480 | # B ========================================================================== |
| 1481 | # C ========================================================================== |
| 1482 | def cohomologous_integrand(self,alpha,asy_var=None): |
| 1483 | r""" |
| 1484 | This function takes a multivariate Cauchy type integral |
| 1485 | `\int F / x^{asy_var\alpha+1} dx`, where `F=self`, and breaks it up |
| 1486 | into a list of nicer Cauchy type integrals for the purposes of |
| 1487 | computing asymptotics of the original integral as `asy_var\to\infty`. |
| 1488 | The sum of the nicer integrals is de Rham cohomologous to the original |
| 1489 | integral. |
| 1490 | It assumes that algebraic varieties corresponding to the irreducible |
| 1491 | factors of `self._H` intersect transversely (see notes below). |
| 1492 | |
| 1493 | INPUT: |
| 1494 | |
| 1495 | - ``alpha`` - A list of positive integers or symbolic variables. |
| 1496 | - ``asy_var`` - A symbolic variable (default: None). |
| 1497 | Eventually set to `var('n')` if None is given. |
| 1498 | |
| 1499 | OUTPUT: |
| 1500 | |
| 1501 | A list of the form `[chunk_1,\ldots,chunk_r]`, where each |
| 1502 | `chunk_j` has the form `[P,[B_1,1],\ldots,[B_l,1]]`. |
| 1503 | Here `l \le d`, `P` is a symbolic expression in the indeterminates of |
| 1504 | `R` and `asy_var`, `\{B_1,\ldots,B_l\} \subseteq \{Q_1,\ldots,Q_m\}`, |
| 1505 | and `[P,[B_1,1],\ldots,[B_l,1]]` represents the integral |
| 1506 | `\int P/(B_1 \cdots B_l) dx`. |
| 1507 | |
| 1508 | EXAMPLES:: |
| 1509 | |
| 1510 | sage: R.<x,y>= PolynomialRing(QQ) |
| 1511 | sage: G= 9*exp(x+y) |
| 1512 | sage: H= (3-2*x-y)*(3-x-2*y) |
| 1513 | sage: F= QuasiRationalExpression(G,H) |
| 1514 | sage: alpha= [4,3] |
| 1515 | sage: F.cohomologous_integrand(alpha) |
| 1516 | [[9*e^(x + y), [x + 2*y - 3, 1], [2*x + y - 3, 1]]] |
| 1517 | |
| 1518 | sage: R.<x,y>= PolynomialRing(QQ) |
| 1519 | sage: G= 9*exp(x+y) |
| 1520 | sage: H= (3-2*x-y)^2*(3-x-2*y) |
| 1521 | sage: F= QuasiRationalExpression(G,H) |
| 1522 | sage: alpha= [4,3] |
| 1523 | sage: F.cohomologous_integrand(alpha) |
| 1524 | [[-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]]] |
| 1525 | |
| 1526 | sage: R.<x,y,z>= PolynomialRing(QQ) |
| 1527 | sage: G= 16 |
| 1528 | sage: H= (4-2*x-y-z)^2*(4-x-2*y-z) |
| 1529 | sage: F= QuasiRationalExpression(G,H) |
| 1530 | sage: alpha= [3,3,2] |
| 1531 | sage: F.cohomologous_integrand(alpha) |
| 1532 | [[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]]] |
| 1533 | |
| 1534 | NOTES: |
| 1535 | |
| 1536 | The varieties corresponding to `Q_1,\ldots,Q_m` |
| 1537 | __intersect transversely__ iff for each point `c` of their intersection |
| 1538 | and for all `k \le l`, the Jacobian matrix of any `k` polynomials from |
| 1539 | `\{Q_1,\ldots,Q_m\}` has rank equal to `\min\{k,d\}` when evaluated at |
| 1540 | `c`. |
| 1541 | |
| 1542 | ALGORITHM: |
| 1543 | |
| 1544 | Let `asy_var= n` and consider the integral around a polycirle centered |
| 1545 | at the origin of the `d`-variate differential form |
| 1546 | `\frac{G(x) dx_1 \wedge\cdots\wedge dx_d}{H(x) x^{n\alpha+1}}`, where |
| 1547 | `G=self._G` and `H=self._H`. |
| 1548 | |
| 1549 | (1) Decompose `G/H` into a sum of partial fractions `P_1 +\cdots+ P_r` |
| 1550 | so that each term of the sum has at most `d` irreducible factors of `H` |
| 1551 | in the denominator. |
| 1552 | |
| 1553 | (2) For each differential form `P_j dx_1 \wedge\cdots\wedge dx_d`, |
| 1554 | find an equivalent form `\omega_j` in de Rham cohomology with no |
| 1555 | repeated irreducible factors of `H` in its denominator. |
| 1556 | |
| 1557 | AUTHOR: |
| 1558 | |
| 1559 | - Alex Raichev (2010-09-22) |
| 1560 | """ |
| 1561 | if not asy_var: |
| 1562 | asy_var = var('n') |
| 1563 | |
| 1564 | # Create symbolic (non-ring) variables. |
| 1565 | G= self._G |
| 1566 | H= self._H |
| 1567 | R= self._R |
| 1568 | d= self._d |
| 1569 | X= self._variables |
| 1570 | # Prepare input for partial_fraction_decomposition() which only works |
| 1571 | # for functions in the field of fractions of R. |
| 1572 | if G in R: |
| 1573 | numer=Integer(1) |
| 1574 | F= G/H |
| 1575 | else: |
| 1576 | numer= G |
| 1577 | F= R(Integer(1) )/H |
| 1578 | nstuff= Integer(1) /mul([X[j]**(alpha[j]*asy_var+Integer(1) ) for j in range(d)]) |
| 1579 | |
| 1580 | # Do steps (2) and (1). |
| 1581 | integrands= [] |
| 1582 | whole,parts= partial_fraction_decomposition(F) |
| 1583 | for f in parts: |
| 1584 | a= format_quotient(f) |
| 1585 | integrands.append( [a[Integer(0) ]*numer*nstuff] + a[Integer(1) :] ) |
| 1586 | |
| 1587 | # Do step (3). |
| 1588 | ce= decompose_via_cohomology(integrands) |
| 1589 | ce_new= [] |
| 1590 | for a in ce: |
| 1591 | ce_new.append( [(a[Integer(0) ]/nstuff).simplify_full().collect(n)] + a[Integer(1) :] ) |
| 1592 | return ce_new |
| 1593 | #------------------------------------------------------------------------------- |
| 1594 | def critical_cone(self,c,coordinate=None): |
| 1595 | r""" |
| 1596 | Returns the critical cone of a convenient multiple point `c`. |
| 1597 | |
| 1598 | INPUT: |
| 1599 | |
| 1600 | - ``c`` - A dictionary with keys `self.variables()` and values |
| 1601 | in a field. |
| 1602 | - ``coordinate`` - (optional) A natural number. |
| 1603 | |
| 1604 | OUTPUT: |
| 1605 | |
| 1606 | A list of vectors that generate the critical cone of `c` and |
| 1607 | the cone itself if the values of `c` lie in QQ. |
| 1608 | |
| 1609 | EXAMPLES:: |
| 1610 | |
| 1611 | sage: R.<x,y,z>= PolynomialRing(QQ) |
| 1612 | sage: G= 1 |
| 1613 | sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) |
| 1614 | sage: F= QuasiRationalExpression(G,H) |
| 1615 | sage: c= {z: 4/3, y: 1, x: 1/2} |
| 1616 | sage: F.critical_cone(c) |
| 1617 | ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N) |
| 1618 | |
| 1619 | NOTES: |
| 1620 | |
| 1621 | The _critical cone_ of a convenient multiple point `c` with |
| 1622 | with `c_k \del_k H_j(c) \neq 0` for all `j=1,\ldots,r` is |
| 1623 | the conical hull of the vectors `\gamma_j(c) = |
| 1624 | \left(\frac{c_1 \del_1 H_j(c)}{c_k \del_k H_j(c)},\ldots, |
| 1625 | \frac{c_d \del_d H_j(c)}{c_k \del_k H_j(c)} \right)`. |
| 1626 | Here `H_1,\ldots,H_r` are the irreducible germs of `self._H` around `c`. |
| 1627 | For more details, see [RaWi2011]_. |
| 1628 | |
| 1629 | If this function's optional argument `coordinate` isn't given, then |
| 1630 | this function searches (from `d` down to 1) for the first index `k` |
| 1631 | such that for all `j=1,\ldots,r` we have `c_k \del_k H_j(c) \neq 0` |
| 1632 | and sets `coordinate = k`. |
| 1633 | Almost. |
| 1634 | Since this is Python, all the indices actually start at 0. |
| 1635 | |
| 1636 | REFERENCES: |
| 1637 | |
| 1638 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |
| 1639 | coefficients of multivariate generating functions: improvements |
| 1640 | for smooth points", To appear. |
| 1641 | |
| 1642 | AUTHORS: |
| 1643 | |
| 1644 | - Alex Raichev (2010-08-25) |
| 1645 | """ |
| 1646 | Hs= [SR(h[Integer(0) ]) for h in self._Hfac] # irreducible factors of H |
| 1647 | X= self._variables |
| 1648 | d= self._d |
| 1649 | # Ensure the variables of `c` lie in SR |
| 1650 | cc= {} |
| 1651 | for x in c.keys(): |
| 1652 | cc[SR(x)] = c[x] |
| 1653 | lg= [self._log_grad(h,X,cc) for h in Hs] |
| 1654 | if not coordinate: |
| 1655 | # Search (from d-1 down to 0) for a coordinate k such that |
| 1656 | # for all h in Hs we have cc[k] * diff(h,X[k]) !=0. |
| 1657 | # One is guaranteed to exist in the case of a convenient multiple |
| 1658 | # point. |
| 1659 | for k in reversed(range(d)): |
| 1660 | if Integer(0) not in [v[k] for v in lg]: |
| 1661 | coordinate= k |
| 1662 | break |
| 1663 | gamma= [self._direction(v,coordinate) for v in lg] |
| 1664 | if [[gg in QQ for gg in g] for g in gamma] == \ |
| 1665 | [[True for gg in g] for g in gamma]: |
| 1666 | return gamma,Cone(gamma) # Cone() needs rational vectors |
| 1667 | else: |
| 1668 | return gamma |
| 1669 | # D ============================================================================ |
| 1670 | def _diff_all(self,f,V,n,ending=[],sub=None,sub_final=None,zero_order=Integer(0) ,rekey=None): |
| 1671 | r""" |
| 1672 | This function returns a dictionary of representative mixed partial |
| 1673 | derivatives of `f` from order 1 up to order `n` with respect to the |
| 1674 | variables in `V`. |
| 1675 | The default is to key the dictionary by all nondecreasing sequences |
| 1676 | in `V` of length 1 up to length `n`. |
| 1677 | For internal use. |
| 1678 | |
| 1679 | Does not use `self`. |
| 1680 | |
| 1681 | INPUT: |
| 1682 | |
| 1683 | - ``f`` - An individual or list of `\mathcal{C}^{n+1}` functions. |
| 1684 | - ``V`` - A list of variables occurring in `f`. |
| 1685 | - ``n`` - A natural number. |
| 1686 | - ``ending`` - A list of variables in `V`. |
| 1687 | - ``sub`` - An individual or list of dictionaries. |
| 1688 | - ``sub_final`` - An individual or list of dictionaries. |
| 1689 | - ``rekey`` - A callable symbolic function in `V` or list thereof. |
| 1690 | - ``zero_order`` - A natural number. |
| 1691 | |
| 1692 | OUTPUT: |
| 1693 | |
| 1694 | The dictionary `{s_1:deriv_1,...,s_r:deriv_r}`. |
| 1695 | Here `s_1,\ldots,s_r` is a listing of |
| 1696 | all nondecreasing sequences of length 1 up to length `n` over the |
| 1697 | alphabet `V`, where `w > v` in `X` iff `str(w) > str(v)`, and |
| 1698 | `deriv_j` is the derivative of `f` with respect to the derivative |
| 1699 | sequence `s_j` and simplified with respect to the substitutions in `sub` |
| 1700 | and evaluated at `sub_final`. |
| 1701 | Moreover, all derivatives with respect to sequences of length less than |
| 1702 | `zero_order` (derivatives of order less than `zero_order`) will be made |
| 1703 | zero. |
| 1704 | |
| 1705 | If `rekey` is nonempty, then `s_1,\ldots,s_r` will be replaced by the |
| 1706 | symbolic derivatives of the functions in `rekey`. |
| 1707 | |
| 1708 | If `ending` is nonempty, then every derivative sequence `s_j` will be |
| 1709 | suffixed by `ending`. |
| 1710 | |
| 1711 | EXAMPLES: |
| 1712 | |
| 1713 | I'd like to print the entire dictionaries, but that doesn't yield |
| 1714 | consistent output order for doctesting. |
| 1715 | Order of keys changes. |
| 1716 | |
| 1717 | :: |
| 1718 | |
| 1719 | sage: R.<x> = PolynomialRing(QQ) |
| 1720 | sage: G,H = 1,1 |
| 1721 | sage: F= QuasiRationalExpression(G,H) |
| 1722 | sage: f= function('f',x) |
| 1723 | sage: dd= F._diff_all(f,[x],3) |
| 1724 | sage: dd[(x,x,x)] |
| 1725 | D[0, 0, 0](f)(x) |
| 1726 | |
| 1727 | :: |
| 1728 | |
| 1729 | sage: d1= {diff(f,x): 4*x^3} |
| 1730 | sage: dd= F._diff_all(f,[x],3,sub=d1) |
| 1731 | sage: dd[(x,x,x)] |
| 1732 | 24*x |
| 1733 | |
| 1734 | :: |
| 1735 | |
| 1736 | sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f) |
| 1737 | sage: dd[diff(f,x,3)] |
| 1738 | 24*x |
| 1739 | |
| 1740 | :: |
| 1741 | |
| 1742 | sage: a= {x:1} |
| 1743 | sage: dd= F._diff_all(f,[x],3,sub=d1,rekey=f,sub_final=a) |
| 1744 | sage: dd[diff(f,x,3)] |
| 1745 | 24 |
| 1746 | |
| 1747 | :: |
| 1748 | |
| 1749 | sage: X= var('x,y,z') |
| 1750 | sage: f= function('f',*X) |
| 1751 | sage: dd= F._diff_all(f,X,2,ending=[y,y,y]) |
| 1752 | sage: dd[(z,y,y,y)] |
| 1753 | D[1, 1, 1, 2](f)(x, y, z) |
| 1754 | |
| 1755 | :: |
| 1756 | |
| 1757 | sage: g= function('g',*X) |
| 1758 | sage: dd= F._diff_all([f,g],X,2) |
| 1759 | sage: dd[(0,y,z)] |
| 1760 | D[1, 2](f)(x, y, z) |
| 1761 | |
| 1762 | :: |
| 1763 | |
| 1764 | sage: dd[(1,z,z)] |
| 1765 | D[2, 2](g)(x, y, z) |
| 1766 | |
| 1767 | :: |
| 1768 | |
| 1769 | sage: f= exp(x*y*z) |
| 1770 | sage: ff= function('ff',*X) |
| 1771 | sage: dd= F._diff_all(f,X,2,rekey=ff) |
| 1772 | sage: dd[diff(ff,x,z)] |
| 1773 | x*y^2*z*e^(x*y*z) + y*e^(x*y*z) |
| 1774 | |
| 1775 | AUTHORS: |
| 1776 | |
| 1777 | - Alex Raichev (2008-10-01, 2009-04-01, 2010-02-01) |
| 1778 | """ |
| 1779 | singleton=False |
| 1780 | if not isinstance(f,list): |
| 1781 | f= [f] |
| 1782 | singleton=True |
| 1783 | |
| 1784 | # Build the dictionary of derivatives iteratively from a list of nondecreasing |
| 1785 | # derivative-order sequences. |
| 1786 | derivs= {} |
| 1787 | r= len(f) |
| 1788 | if ending: |
| 1789 | seeds = [ending] |
| 1790 | start = Integer(1) |
| 1791 | else: |
| 1792 | seeds = [[v] for v in V] |
| 1793 | start = Integer(2) |
| 1794 | if singleton: |
| 1795 | for s in seeds: |
| 1796 | derivs[tuple(s)] = self._subs_all(diff(f[Integer(0) ],s),sub) |
| 1797 | for l in (ellipsis_range(start,Ellipsis,n)): |
| 1798 | for t in UnorderedTuples(V,l): |
| 1799 | s= tuple(t + ending) |
| 1800 | derivs[s] = self._subs_all(diff(derivs[s[Integer(1) :]],s[Integer(0) ]),sub) |
| 1801 | else: |
| 1802 | # Make the dictionary keys of the form (j,sequence of variables), |
| 1803 | # where j in range(r). |
| 1804 | for s in seeds: |
| 1805 | value= self._subs_all([diff(f[j],s) for j in range(r)],sub) |
| 1806 | derivs.update(dict([(tuple([j]+s),value[j]) for j in range(r)])) |
| 1807 | for l in (ellipsis_range(start,Ellipsis,n)): |
| 1808 | for t in UnorderedTuples(V,l): |
| 1809 | s= tuple(t + ending) |
| 1810 | value= self._subs_all(\ |
| 1811 | [diff(derivs[(j,)+s[Integer(1) :]],s[Integer(0) ]) for j in range(r)],sub) |
| 1812 | derivs.update(dict([((j,)+s,value[j]) for j in range(r)])) |
| 1813 | if zero_order: |
| 1814 | # Zero out all the derivatives of order < zero_order |
| 1815 | if singleton: |
| 1816 | for k in derivs.keys(): |
| 1817 | if len(k) < zero_order: |
| 1818 | derivs[k]= Integer(0) |
| 1819 | else: |
| 1820 | # Ignore the first of element of k, which is an index. |
| 1821 | for k in derivs.keys(): |
| 1822 | if len(k)-Integer(1) < zero_order: |
| 1823 | derivs[k]= Integer(0) |
| 1824 | if sub_final: |
| 1825 | # Substitute sub_final into the values of derivs. |
| 1826 | for k in derivs.keys(): |
| 1827 | derivs[k] = self._subs_all(derivs[k],sub_final) |
| 1828 | if rekey: |
| 1829 | # Rekey the derivs dictionary by the value of rekey. |
| 1830 | F= rekey |
| 1831 | if singleton: |
| 1832 | # F must be a singleton. |
| 1833 | derivs= dict( [(diff(F,list(k)), derivs[k]) for k in derivs.keys()] ) |
| 1834 | else: |
| 1835 | # F must be a list. |
| 1836 | derivs= dict( [(diff(F[k[Integer(0) ]],list(k)[Integer(1) :]), derivs[k]) for k in derivs.keys()] ) |
| 1837 | return derivs |
| 1838 | #------------------------------------------------------------------------------- |
| 1839 | def _diff_op(self,A,B,AB_derivs,V,M,r,N): |
| 1840 | r""" |
| 1841 | This function computes the derivatives `DD^(l+k)(A[j] B^l)` evaluated at a |
| 1842 | point `p` for various natural numbers `j,k,l` which depend on `r` and `N`. |
| 1843 | Here `DD` is a specific second-order linear differential operator that depends |
| 1844 | on `M`, `A` is a list of symbolic functions, `B` is symbolic function, |
| 1845 | and `AB_derivs` contains all the derivatives of `A` and `B` evaluated at `p` |
| 1846 | that are necessary for the computation. |
| 1847 | For internal use by the functions _asymptotics_main_multiple() and |
| 1848 | _asymptotics_main_smooth(). |
| 1849 | |
| 1850 | Does not use `self`. |
| 1851 | |
| 1852 | INPUT: |
| 1853 | |
| 1854 | - ``A`` - A single or length `r` list of symbolic functions in the |
| 1855 | variables `V`. |
| 1856 | - ``B`` - A symbolic function in the variables `V`. |
| 1857 | - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of |
| 1858 | `A[0],\ldots,A[r-1]` up to order `2N-2` and |
| 1859 | the (symbolic) derivatives of `B` up to order `2N`. |
| 1860 | The values of the dictionary are complex numbers that are |
| 1861 | the keys evaluated at a common point `p`. |
| 1862 | - ``V`` - The variables of the `A[j]` and `B`. |
| 1863 | - ``M`` - A symmetric `l \times l` matrix, where `l` is the length of `V`. |
| 1864 | - ``r,N`` - Natural numbers. |
| 1865 | |
| 1866 | OUTPUT: |
| 1867 | |
| 1868 | A dictionary whose keys are natural number tuples of the form `(j,k,l)`, |
| 1869 | where `l \le 2k`, `j \le r-1`, and `j+k \le N-1`, and whose values are |
| 1870 | `DD^(l+k)(A[j] B^l)` evaluated at a point `p`, where `DD` is the linear |
| 1871 | second-order differential operator |
| 1872 | `-\sum_{i=0}^{l-1} \sum_{j=0}^{l-1} M[i][j] |
| 1873 | \partial^2 /(\partial V[j] \partial V[i])`. |
| 1874 | |
| 1875 | EXAMPLES:: |
| 1876 | |
| 1877 | sage: R.<x> = PolynomialRing(QQ) |
| 1878 | sage: G,H = 1,1 |
| 1879 | sage: F= QuasiRationalExpression(G,H) |
| 1880 | sage: T= var('x,y') |
| 1881 | sage: A= function('A',*tuple(T)) |
| 1882 | sage: B= function('B',*tuple(T)) |
| 1883 | sage: AB_derivs= {} |
| 1884 | sage: M= matrix([[1,2],[2,1]]) |
| 1885 | sage: DD= F._diff_op(A,B,AB_derivs,T,M,1,2) |
| 1886 | sage: DD.keys() |
| 1887 | [(0, 1, 2), (0, 1, 1), (0, 1, 0), (0, 0, 0)] |
| 1888 | sage: len(DD[(0,1,2)]) |
| 1889 | 246 |
| 1890 | |
| 1891 | AUTHORS: |
| 1892 | |
| 1893 | - Alex Raichev (2008-10-01, 2010-01-12) |
| 1894 | """ |
| 1895 | if not isinstance(A,list): |
| 1896 | A= [A] |
| 1897 | |
| 1898 | # First, compute the necessary product derivatives of A and B. |
| 1899 | product_derivs= {} |
| 1900 | for (j,k) in mrange([r,N]): |
| 1901 | if j+k <N: |
| 1902 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)): |
| 1903 | for s in UnorderedTuples(V,Integer(2) *(k+l)): |
| 1904 | product_derivs[tuple([j,k,l]+s)] = \ |
| 1905 | diff(A[j]*B**l,s).subs(AB_derivs) |
| 1906 | |
| 1907 | # Second, compute DD^(k+l)(A[j]*B^l)(p) and store values in dictionary. |
| 1908 | DD= {} |
| 1909 | rows= M.nrows() |
| 1910 | for (j,k) in mrange([r,N]): |
| 1911 | if j+k <N: |
| 1912 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)): |
| 1913 | # Take advantage of the symmetry of M by ignoring |
| 1914 | # the upper-diagonal entries of M and multiplying by |
| 1915 | # appropriate powers of 2. |
| 1916 | if k+l == Integer(0) : |
| 1917 | DD[(j,k,l)] = product_derivs[(j,k,l)] |
| 1918 | continue |
| 1919 | S= [(a,b) for (a,b) in mrange([rows,rows]) if b <= a] |
| 1920 | P= cartesian_product_iterator([S for i in range(k+l)]) |
| 1921 | diffo= Integer(0) |
| 1922 | for t in P: |
| 1923 | if product_derivs[(j,k,l)+self._diff_seq(V,t)] != Integer(0) : |
| 1924 | MM= Integer(1) |
| 1925 | for (a,b) in t: |
| 1926 | MM= MM * M[a][b] |
| 1927 | if a != b: |
| 1928 | MM= Integer(2) *MM |
| 1929 | diffo= diffo + MM * product_derivs[(j,k,l)+self._diff_seq(V,t)] |
| 1930 | DD[(j,k,l)] = (-Integer(1) )**(k+l)*diffo |
| 1931 | return DD |
| 1932 | #------------------------------------------------------------------------------- |
| 1933 | def _diff_op_simple(self,A,B,AB_derivs,x,v,a,N): |
| 1934 | r""" |
| 1935 | This function computes `DD^(e k + v l)(A B^l)` evaluated at a point `p` |
| 1936 | for various natural numbers `e,k,l` that depend on `v` and `N`. |
| 1937 | Here `DD` is a specific linear differential operator that depends |
| 1938 | on `a` and `v`, `A` and `B` are symbolic functions, and `AB_derivs` contains |
| 1939 | all the derivatives of `A` and `B` evaluated at `p` that are necessary for |
| 1940 | the computation. |
| 1941 | For internal use by the function _asymptotics_main_smooth(). |
| 1942 | |
| 1943 | Does not use `self`. |
| 1944 | |
| 1945 | INPUT: |
| 1946 | |
| 1947 | - ``A``,``B`` - Symbolic functions in the variable `x`. |
| 1948 | - ``AB_derivs`` - A dictionary whose keys are the (symbolic) derivatives of |
| 1949 | `A` up to order `2N` if `v` is even or `N` if `v` is odd and |
| 1950 | the (symbolic) derivatives of `B` up to order `2N+v` if `v` is even |
| 1951 | or `N+v` if `v` is odd. |
| 1952 | The values of the dictionary are complex numbers that are |
| 1953 | the keys evaluated at a common point `p`. |
| 1954 | - ``x`` - Symbolic variable. |
| 1955 | - ``a`` - A complex number. |
| 1956 | - ``v``,``N`` - Natural numbers. |
| 1957 | |
| 1958 | OUTPUT: |
| 1959 | |
| 1960 | A dictionary whose keys are natural number pairs of the form `(k,l)`, |
| 1961 | where `k < N` and `l \le 2k` and whose values are |
| 1962 | `DD^(e k + v l)(A B^l)` evaluated at a point `p`. |
| 1963 | Here `e=2` if `v` is even, `e=1` if `v` is odd, and `DD` is a particular |
| 1964 | linear differential operator |
| 1965 | `(a^{-1/v} d/dt)' if `v` is even and `(|a|^{-1/v} i \sgn a d/dt)` |
| 1966 | if `v` is odd. |
| 1967 | |
| 1968 | EXAMPLES:: |
| 1969 | |
| 1970 | sage: R.<x> = PolynomialRing(QQ) |
| 1971 | sage: G,H = 1,1 |
| 1972 | sage: F= QuasiRationalExpression(G,H) |
| 1973 | sage: A= function('A',x) |
| 1974 | sage: B= function('B',x) |
| 1975 | sage: AB_derivs= {} |
| 1976 | sage: F._diff_op_simple(A,B,AB_derivs,x,3,2,2) |
| 1977 | {(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)} |
| 1978 | |
| 1979 | AUTHORS: |
| 1980 | |
| 1981 | - Alex Raichev (2010-01-15) |
| 1982 | """ |
| 1983 | I= sqrt(-Integer(1) ) |
| 1984 | DD= {} |
| 1985 | if v.mod(Integer(2) ) == Integer(0) : |
| 1986 | for k in range(N): |
| 1987 | for l in (ellipsis_range(Integer(0) ,Ellipsis,Integer(2) *k)): |
| 1988 | DD[(k,l)] = (a**(-Integer(1) /v))**(Integer(2) *k+v*l) \ |
| 1989 | * diff(A*B**l,x,Integer(2) *k+v*l).subs(AB_derivs) |
| 1990 | else: |
| 1991 | for k in range(N): |
| 1992 | for l in (ellipsis_range(Integer(0) ,Ellipsis,k)): |
| 1993 | DD[(k,l)] = (abs(a)**(-Integer(1) /v)*I*a/abs(a))**(k+v*l) \ |
| 1994 | * diff(A*B**l,x,k+v*l).subs(AB_derivs) |
| 1995 | return DD |
| 1996 | #------------------------------------------------------------------------------- |
| 1997 | def _diff_prod(self,f_derivs,u,g,X,interval,end,uderivs,atc): |
| 1998 | r""" |
| 1999 | This function takes various derivatives of the equation `f=ug`, |
| 2000 | evaluates at a point `c`, and solves for the derivatives of `u`. |
| 2001 | For internal use by the function _asymptotics_main_multiple(). |
| 2002 | |
| 2003 | Does not use `self`. |
| 2004 | |
| 2005 | INPUT: |
| 2006 | |
| 2007 | - ``f_derivs`` - A dictionary whose keys are tuples \code{s + end} for all |
| 2008 | `X`-variable sequences `s` with length in `interval` and whose |
| 2009 | values are the derivatives of a function `f` evaluated at `c`. |
| 2010 | - ``u`` - A callable symbolic function. |
| 2011 | - ``g`` - An expression or callable symbolic function. |
| 2012 | - ``X`` - A list of symbolic variables. |
| 2013 | - ``interval`` - A list of positive integers. |
| 2014 | Call the first and last values `n` and `nn`, respectively. |
| 2015 | - ``end`` - A possibly empty list of `z`'s where `z` is the last element of |
| 2016 | `X`. |
| 2017 | - ``uderivs`` - A dictionary whose keys are the symbolic |
| 2018 | derivatives of order 0 to order `n-1` of `u` evaluated at `c` |
| 2019 | and whose values are the corresponding derivatives evaluated at `c`. |
| 2020 | - ``atc`` - A dictionary whose keys are the keys of `c` and all the symbolic |
| 2021 | derivatives of order 0 to order `nn` of `g` evaluated `c` and whose |
| 2022 | values are the corresponding derivatives evaluated at `c`. |
| 2023 | |
| 2024 | OUTPUT: |
| 2025 | |
| 2026 | A dictionary whose keys are the derivatives of `u` up to order |
| 2027 | `nn` and whose values are those derivatives evaluated at `c`. |
| 2028 | |
| 2029 | EXAMPLES:: |
| 2030 | |
| 2031 | I'd like to print out the entire dictionary, but that does not give |
| 2032 | consistent output for doctesting. |
| 2033 | Order of keys changes :: |
| 2034 | |
| 2035 | sage: R.<x> = PolynomialRing(QQ) |
| 2036 | sage: G,H = 1,1 |
| 2037 | sage: F= QuasiRationalExpression(G,H) |
| 2038 | sage: u= function('u',x) |
| 2039 | sage: g= function('g',x) |
| 2040 | 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}) |
| 2041 | sage: dd[diff(u,x,2)(x=2)] |
| 2042 | 22/9 |
| 2043 | |
| 2044 | NOTES: |
| 2045 | |
| 2046 | This function works by differentiating the equation `f=ug` with respect |
| 2047 | to the variable sequence \code{s+end}, for all tuples `s` of `X` of |
| 2048 | lengths in `interval`, evaluating at the point `c`, |
| 2049 | and solving for the remaining derivatives of `u`. |
| 2050 | This function assumes that `u` never appears in the differentiations of |
| 2051 | `f=ug` after evaluating at `c`. |
| 2052 | |
| 2053 | AUTHORS: |
| 2054 | |
| 2055 | - Alex Raichev (2009-05-14, 2010-01-21) |
| 2056 | """ |
| 2057 | for l in interval: |
| 2058 | D= {} |
| 2059 | rhs= [] |
| 2060 | lhs= [] |
| 2061 | for t in UnorderedTuples(X,l): |
| 2062 | s= t+end |
| 2063 | lhs.append(f_derivs[tuple(s)]) |
| 2064 | rhs.append(diff(u*g,s).subs(atc).subs(uderivs)) |
| 2065 | # Since Sage's solve command can't take derivatives as variable |
| 2066 | # names, i make new variables based on t to stand in for |
| 2067 | # diff(u,t) and store them in D. |
| 2068 | D[diff(u,t).subs(atc)] = self._make_var([var('zing')]+t) |
| 2069 | eqns=[ lhs[i] == rhs[i].subs(uderivs).subs(D) for i in range(len(lhs))] |
| 2070 | vars= D.values() |
| 2071 | sol= solve(eqns,*vars,solution_dict=True) |
| 2072 | uderivs.update(self._subs_all(D,sol[Integer(0) ])) |
| 2073 | return uderivs |
| 2074 | #------------------------------------------------------------------------------- |
| 2075 | def _diff_seq(self,V,s): |
| 2076 | r""" |
| 2077 | Given a list `s` of tuples of natural numbers, this function returns the |
| 2078 | list of elements of `V` with indices the elements of the elements of `s`. |
| 2079 | This function is for internal use by the function _diff_op(). |
| 2080 | |
| 2081 | Does not use `self`. |
| 2082 | |
| 2083 | INPUT: |
| 2084 | |
| 2085 | - ``V`` - A list. |
| 2086 | - ``s`` - A list of tuples of natural numbers in the interval |
| 2087 | \code{range(len(V))}. |
| 2088 | |
| 2089 | OUTPUT: |
| 2090 | |
| 2091 | The tuple \code{tuple([V[tt] for tt in sorted(t)])}, where `t` is the |
| 2092 | list of elements of the elements of `s`. |
| 2093 | |
| 2094 | EXAMPLES:: |
| 2095 | |
| 2096 | sage: R.<x> = PolynomialRing(QQ) |
| 2097 | sage: G,H = 1,1 |
| 2098 | sage: F= QuasiRationalExpression(G,H) |
| 2099 | sage: V= list(var('x,t,z')) |
| 2100 | sage: F._diff_seq(V,([0,1],[0,2,1],[0,0])) |
| 2101 | (x, x, x, x, t, t, z) |
| 2102 | |
| 2103 | AUTHORS: |
| 2104 | |
| 2105 | - Alex Raichev (2009.05.19) |
| 2106 | """ |
| 2107 | t= [] |
| 2108 | for ss in s: |
| 2109 | t.extend(ss) |
| 2110 | return tuple([V[tt] for tt in sorted(t)]) |
| 2111 | #------------------------------------------------------------------------------- |
| 2112 | def _direction(self,v,coordinate=None): |
| 2113 | r""" |
| 2114 | This function returns \code{[vv/v[coordinate] for vv in v]} where |
| 2115 | coordinate is the last index of v if not specified otherwise. |
| 2116 | |
| 2117 | Does not use `self`. |
| 2118 | |
| 2119 | INPUT: |
| 2120 | |
| 2121 | - ``v`` - A vector. |
| 2122 | - ``coordinate`` - An index for `v` (default: None). |
| 2123 | |
| 2124 | OUTPUT: |
| 2125 | |
| 2126 | This function returns \code{[vv/v[coordinate] for vv in v]} where |
| 2127 | coordinate is the last index of v if not specified otherwise. |
| 2128 | |
| 2129 | EXAMPLES:: |
| 2130 | |
| 2131 | sage: R.<x> = PolynomialRing(QQ) |
| 2132 | sage: G,H = 1,1 |
| 2133 | sage: F= QuasiRationalExpression(G,H) |
| 2134 | sage: F._direction([2,3,5]) |
| 2135 | (2/5, 3/5, 1) |
| 2136 | sage: F._direction([2,3,5],0) |
| 2137 | (1, 3/2, 5/2) |
| 2138 | |
| 2139 | AUTHORS: |
| 2140 | |
| 2141 | - Alex Raichev (2010-08-25) |
| 2142 | """ |
| 2143 | if coordinate == None: |
| 2144 | coordinate= len(v)-Integer(1) |
| 2145 | try: |
| 2146 | v[Integer(0) ].variables() |
| 2147 | return tuple([(vv/v[coordinate]).simplify_full() for vv in v]) |
| 2148 | except: |
| 2149 | return tuple([vv/v[coordinate] for vv in v]) |
| 2150 | # E ============================================================================ |
| 2151 | # F ============================================================================ |
| 2152 | # G ============================================================================ |
| 2153 | # H ============================================================================ |
| 2154 | # I ============================================================================ |
| 2155 | def is_cmp(self,points): |
| 2156 | r""" |
| 2157 | Checks if the points in the list `points` are convenient multiple |
| 2158 | points of `V= \{ x\in CC^d : H(x) = 0\}`, where `H=self._H`. |
| 2159 | |
| 2160 | INPUT: |
| 2161 | |
| 2162 | - ``points`` - An individual or list of dictionaries with keys |
| 2163 | `self._variables` and values in some superfield of |
| 2164 | `self._R.base_ring()`. |
| 2165 | |
| 2166 | OUTPUT: |
| 2167 | |
| 2168 | A list of tuples `(p,verdict,comment)`, one for each point |
| 2169 | `p` in `points`, where `verdict` is True if `p` is a convenient |
| 2170 | multiple point and False otherwise, and where `comment` is a string |
| 2171 | comment relating to `verdict`, such as 'not a transverse intersection'. |
| 2172 | |
| 2173 | EXAMPLES:: |
| 2174 | |
| 2175 | sage: R.<x,y,z>= PolynomialRing(QQ) |
| 2176 | sage: G= 16 |
| 2177 | sage: H= (1-x*(1+y))*(1-z*x^2*(1+2*y)) |
| 2178 | sage: F= QuasiRationalExpression(G,H) |
| 2179 | sage: points= [{x:1/2,y:1,z:4/3},{x:1/2,y:1,z:2}] |
| 2180 | sage: F.is_cmp(points) |
| 2181 | [({y: 1, z: 4/3, x: 1/2}, True, 'all good'), ({y: 1, z: 2, x: 1/2}, False, 'not a singular point')] |
| 2182 | |
| 2183 | NOTES: |
| 2184 | |
| 2185 | A point `c` of `V` is a __convenient multiple point__ if `V` is locally |
| 2186 | a union of complex manifolds that intersect transversely at `c`; |
| 2187 | see [RaWi2011]_. |
| 2188 | |
| 2189 | REFERENCES: |
| 2190 | |
| 2191 | .. [RaWi2011] Alexander Raichev and Mark C. Wilson, "Asymptotics of |
| 2192 | coefficients of multivariate generating functions: improvements |
| 2193 | for smooth points", To appear. |
| 2194 | |
| 2195 | AUTHORS: |
| 2196 | |
| 2197 | - Alex Raichev (2011-04-18) |
| 2198 | """ |
| 2199 | H= self._H |
| 2200 | d= self._d |
| 2201 | Hs= [SR(h[Integer(0) ]) for h in self._Hfac] # irreducible factors of H |
| 2202 | X= self._variables |
| 2203 | J= [tuple([diff(h,x) for x in X]) for h in Hs] |
| 2204 | verdicts= [] |
| 2205 | if not isinstance(points,list): |
| 2206 | points= [points] |
| 2207 | for p in points: |
| 2208 | # Ensure variables in points lie in SR. |
| 2209 | pp= {} |
| 2210 | for x in p.keys(): |
| 2211 | pp[SR(x)]= p[x] |
| 2212 | # Test 1: Is p a zero of all factors of H? |
| 2213 | if [h.subs(pp) for h in Hs] != [Integer(0) for h in Hs]: |
| 2214 | # Failed test 1. Move on to next point. |
| 2215 | verdicts.append((p,False,'not a singular point')) |
| 2216 | continue |
| 2217 | # Test 2: Are the factors of H smooth and |
| 2218 | # do theyintersect transversely at p? |
| 2219 | J= [tuple([f.subs(pp) for f in dh]) for dh in J] |
| 2220 | l= len(J) |
| 2221 | if Set(J).cardinality() < l: |
| 2222 | # Fail. Move on to next point. |
| 2223 | verdicts.append((p,False,'not a transverse intersection')) |
| 2224 | continue |
| 2225 | temp= True |
| 2226 | for S in list(Set(J).subsets())[Integer(1) :]: # Subsets of size >= 1 |
| 2227 | k= len(list(S)) |
| 2228 | M = matrix(list(S)) |
| 2229 | if M.rank() != min(k,d): |
| 2230 | # Fail. |
| 2231 | temp= False |
| 2232 | verdicts.append((p,False,'not a transvere intersection')) |
| 2233 | break |
| 2234 | if not temp: |
| 2235 | # Move on to next point |
| 2236 | continue |
| 2237 | # Test 3: Is p convenient? |
| 2238 | Jlog= matrix([self._log_grad(h,X,pp) for h in Hs]) |
| 2239 | if [Integer(0) in f for f in Jlog.columns()] ==\ |
| 2240 | [True for f in Jlog.columns()]: |
| 2241 | # Fail. Move on to next point. |
| 2242 | verdict[p] = (False,'multiple point but not convenient') |
| 2243 | continue |
| 2244 | verdicts.append((p,True,'all good')) |
| 2245 | return verdicts |
| 2246 | # J ============================================================================ |
| 2247 | # K ============================================================================ |
| 2248 | # L ============================================================================ |
| 2249 | def _log_grad(self,f,X,c): |
| 2250 | r""" |
| 2251 | This function returns the logarithmic gradient of `f` with respect to the |
| 2252 | variables of `X` evalutated at `c`. |
| 2253 | |
| 2254 | Does not use `self`. |
| 2255 | |
| 2256 | INPUT: |
| 2257 | |
| 2258 | - ``f`` - An expression in the variables of `X`. |
| 2259 | - ``X`` - A list of variables. |
| 2260 | - ``c`` - A dictionary with keys `X` and values in a field `K`. |
| 2261 | |
| 2262 | OUTPUT: |
| 2263 | |
| 2264 | \code{[c[x] * diff(f,x).subs(c) for x in X]}. |
| 2265 | |
| 2266 | EXAMPLES:: |
| 2267 | |
| 2268 | sage: R.<x> = PolynomialRing(QQ) |
| 2269 | sage: G,H = 1,1 |
| 2270 | sage: F= QuasiRationalExpression(G,H) |
| 2271 | sage: X= var('x,y,z') |
| 2272 | sage: f= x*y*z^2 |
| 2273 | sage: c= {x:1,y:2,z:3} |
| 2274 | sage: f.gradient() |
| 2275 | (y*z^2, x*z^2, 2*x*y*z) |
| 2276 | sage: F._log_grad(f,X,c) |
| 2277 | (18, 18, 36) |
| 2278 | |
| 2279 | :: |
| 2280 | |
| 2281 | sage: R.<x,y,z>= PolynomialRing(QQ) |
| 2282 | sage: f= x*y*z^2 |
| 2283 | sage: c= {x:1,y:2,z:3} |
| 2284 | sage: F._log_grad(f,R.gens(),c) |
| 2285 | (18, 18, 36) |
| 2286 | |
| 2287 | AUTHORS: |
| 2288 | |
| 2289 | - Alex Raichev (2009-03-06) |
| 2290 | """ |
| 2291 | return tuple([SR(c[x] * diff(f,x).subs(c)).simplify() for x in X]) |
| 2292 | # M ============================================================================ |
| 2293 | def _make_var(self,L): |
| 2294 | r""" |
| 2295 | This function converts the list `L` to a string (without commas) and returns |
| 2296 | the string as a variable. |
| 2297 | For internal use by the function _diff_op() |
| 2298 | |
| 2299 | Does not use `self`. |
| 2300 | |
| 2301 | INPUT: |
| 2302 | |
| 2303 | - ``L`` - A list. |
| 2304 | |
| 2305 | OUTPUT: |
| 2306 | |
| 2307 | A variable whose name is the concatenation of the variable names in `L`. |
| 2308 | |
| 2309 | EXAMPLES:: |
| 2310 | |
| 2311 | sage: R.<x> = PolynomialRing(QQ) |
| 2312 | sage: G,H = 1,1 |
| 2313 | sage: F= QuasiRationalExpression(G,H) |
| 2314 | sage: L= list(var('x,y,hello')) |
| 2315 | sage: v= F._make_var(L) |
| 2316 | sage: print v, type(v) |
| 2317 | xyhello <type 'sage.symbolic.expression.Expression'> |
| 2318 | |
| 2319 | AUTHORS: |
| 2320 | |
| 2321 | - Alex Raichev (2010-01-21) |
| 2322 | """ |
| 2323 | return var(''.join([str(v) for v in L])) |
| 2324 | # N ============================================================================ |
| 2325 | def _new_var_name(self,name,V): |
| 2326 | r""" |
| 2327 | This function returns the first string in the sequence `name`, `name+name`, |
| 2328 | `name+name+name`,... that does not appear in the list `V`. |
| 2329 | It is for internal use by the function _asymptotics_main_multiple(). |
| 2330 | |
| 2331 | Does not use `self`. |
| 2332 | |
| 2333 | INPUT: |
| 2334 | |
| 2335 | - ``name`` - A string. |
| 2336 | - ``V`` - A list of variables. |
| 2337 | |
| 2338 | OUTPUT: |
| 2339 | |
| 2340 | The first string in the sequence `name`, `name+name`, |
| 2341 | `name+name+name`,... that does not appear in the list \code{str(V)}. |
| 2342 | |
| 2343 | EXAMPLES:: |
| 2344 | |
| 2345 | sage: R.<x> = PolynomialRing(QQ) |
| 2346 | sage: G,H = 1,1 |
| 2347 | sage: F= QuasiRationalExpression(G,H) |
| 2348 | sage: X= var('x,xx,y,z') |
| 2349 | sage: F._new_var_name('x',X) |
| 2350 | 'xxx' |
| 2351 | |
| 2352 | AUTHORS: |
| 2353 | |
| 2354 | - Alex Raichev (2008-10-01) |
| 2355 | """ |
| 2356 | newname= name |
| 2357 | while newname in str(V): |
| 2358 | newname= newname +name |
| 2359 | return newname |
| 2360 | # O ============================================================================ |
| 2361 | # P ============================================================================ |
| 2362 | # Q ============================================================================ |
| 2363 | # R ============================================================================ |
| 2364 | def relative_error(self,approx,alpha,interval,exp_scale=Integer(1) ): |
| 2365 | r""" |
| 2366 | Returns the relative error between the values of the Maclaurin |
| 2367 | coefficients of `self` with multi-indices `m alpha` for `m` in |
| 2368 | `interval` and the values of the functions in `approx`. |
| 2369 | |
| 2370 | INPUT: |
| 2371 | |
| 2372 | - ``approx`` - An individual or list of symbolic expressions in |
| 2373 | one variable. |
| 2374 | - ``alpha`` - A list of positive integers. |
| 2375 | - ``interval`` - A list of positive integers. |
| 2376 | - ``exp_scale`` - (optional) A number. Default: 1. |
| 2377 | |
| 2378 | OUTPUT: |
| 2379 | |
| 2380 | A list whose entries are of the form |
| 2381 | `[m\alpha,a_m,b_m,err_m]` for `m \in interval`. |
| 2382 | Here `m\alpha` is a tuple; `a_m` is the `m alpha` (multi-index) |
| 2383 | coefficient of the Maclaurin series for `F` divided by `exp_scale^m`; |
| 2384 | `b_m` is a list of the values of the functions in `approx` evaluated at |
| 2385 | `m` and divided by `exp_scale^m`; `err_m` is the list of relative errors |
| 2386 | `(a_m-f)/a_m` for `f` in `b_m`. |
| 2387 | All outputs are decimal approximations. |
| 2388 | |
| 2389 | EXAMPLES:: |
| 2390 | |
| 2391 | sage: R.<x,y>= PolynomialRing(QQ) |
| 2392 | sage: G=1 |
| 2393 | sage: H= 1-x-y-x*y |
| 2394 | sage: F= QuasiRationalExpression(G,H) |
| 2395 | sage: alpha= [1,1] |
| 2396 | sage: var('n') |
| 2397 | n |
| 2398 | sage: f= (0.573/sqrt(n))*5.83^n |
| 2399 | sage: es= 5.83 |
| 2400 | sage: F.relative_error(f,alpha,[1,2,4,8],es) |
| 2401 | Calculating errors table in the form |
| 2402 | exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... |
| 2403 | [[(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]]] |
| 2404 | sage: g= (0.573/sqrt(n) - 0.0674/n^(3/2))*5.83^n |
| 2405 | sage: F.relative_error([f,g],alpha,[1,2,4,8],es) |
| 2406 | Calculating errors table in the form |
| 2407 | exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... |
| 2408 | [[(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]]] |
| 2409 | |
| 2410 | AUTHORS: |
| 2411 | |
| 2412 | - Alex Raichev (2009-05-18, 2011-04-18) |
| 2413 | """ |
| 2414 | |
| 2415 | if not isinstance(approx,list): |
| 2416 | approx= [approx] |
| 2417 | av= approx[Integer(0) ].variables()[Integer(0) ] |
| 2418 | |
| 2419 | print "Calculating errors table in the form" |
| 2420 | print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..." |
| 2421 | |
| 2422 | # Get Maclaurin coefficients of self and scale them. |
| 2423 | # Then compute errors. |
| 2424 | n= interval[-Integer(1) ] |
| 2425 | mc= self.maclaurin_coefficients(alpha,n) |
| 2426 | mca={} |
| 2427 | stats=[] |
| 2428 | for m in interval: |
| 2429 | beta= tuple([m*a for a in alpha]) |
| 2430 | mc[beta]= mc[beta]/exp_scale**m |
| 2431 | mca[beta]= [f.subs({av:m})/exp_scale**m for f in approx] |
| 2432 | stats_row= [beta, mc[beta].n(), [a.n() for a in mca[beta]]] |
| 2433 | if mc[beta]==Integer(0) : |
| 2434 | stats_row.extend([None for a in mca[beta]]) |
| 2435 | else: |
| 2436 | stats_row.append([((mc[beta]-a)/mc[beta]).n() for a in mca[beta]]) |
| 2437 | stats.append(stats_row) |
| 2438 | return stats |
| 2439 | # S ============================================================================ |
| 2440 | def singular_points(self): |
| 2441 | r""" |
| 2442 | This function returns a Groebner basis ideal whose variety is the |
| 2443 | set of singular points of the algebraic variety |
| 2444 | `V= \{x\in\CC^d : H(x)=0\}`, where `H=sef._H`. |
| 2445 | |
| 2446 | INPUT: |
| 2447 | |
| 2448 | OUTPUT: |
| 2449 | |
| 2450 | A Groebner basis ideal whose variety is the set of singular points of |
| 2451 | the algebraic variety `V= \{x\in\CC^d : H(x)=0\}`. |
| 2452 | |
| 2453 | EXAMPLES:: |
| 2454 | |
| 2455 | sage: R.<x,y,z>= PolynomialRing(QQ) |
| 2456 | sage: G= 1 |
| 2457 | sage: H= (4-2*x-y-z)*(4-x-2*y-z) |
| 2458 | sage: F= QuasiRationalExpression(G,H) |
| 2459 | sage: F.singular_points() |
| 2460 | Ideal (x + 1/3*z - 4/3, y + 1/3*z - 4/3) of Multivariate Polynomial Ring in x, y, z over Rational Field |
| 2461 | |
| 2462 | AUTHORS: |
| 2463 | |
| 2464 | - Alex Raichev (2008-10-01, 2008-11-20, 2010-12-03, 2011-04-18) |
| 2465 | """ |
| 2466 | H= self._H |
| 2467 | R= self._R |
| 2468 | f= R.ideal(H).radical().gens()[Integer(0) ] # Compute the reduction of H. |
| 2469 | J= R.ideal([f] + f.gradient()) |
| 2470 | return R.ideal(J.groebner_basis()) |
| 2471 | #------------------------------------------------------------------------------- |
| 2472 | def smooth_critical(self,alpha): |
| 2473 | r""" |
| 2474 | This function returns a Groebner basis ideal whose variety is the set |
| 2475 | of smooth critical points of the algebraic variety |
| 2476 | `V= \{x\in\CC^d : H(x)=0\} for the direction `\alpha` where `H=self._H`. |
| 2477 | |
| 2478 | INPUT: |
| 2479 | |
| 2480 | - ``alpha`` - A `d`-tuple of positive integers and/or symbolic entries. |
| 2481 | |
| 2482 | OUTPUT: |
| 2483 | |
| 2484 | A Groebner basis ideal of smooth critical points of `V` for `\alpha`. |
| 2485 | |
| 2486 | EXAMPLES:: |
| 2487 | |
| 2488 | sage: R.<x,y> = PolynomialRing(QQ) |
| 2489 | sage: G=1 |
| 2490 | sage: H = (1-x-y-x*y)^2 |
| 2491 | sage: F= QuasiRationalExpression(G,H) |
| 2492 | sage: var('a1,a2') |
| 2493 | (a1, a2) |
| 2494 | sage: F.smooth_critical([a1,a2]) |
| 2495 | 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 |
| 2496 | |
| 2497 | NOTES: |
| 2498 | |
| 2499 | A point `c` of `V` is a __smooth critical point for `alpha`__ |
| 2500 | if the gradient of `f` at `c` is not identically zero and `\alpha` is in |
| 2501 | the span of the logarithmic gradient vector |
| 2502 | `(c[0] \partial_1 f(c)),\ldots,c[d-1] \partial_d f(c))`; see [RaWi2008a]_. |
| 2503 | |
| 2504 | REFERENCES: |
| 2505 | |
| 2506 | .. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of |
| 2507 | coefficients of multivariate generating functions: improvements |
| 2508 | for smooth points", Electronic Journal of Combinatorics, Vol. 15 |
| 2509 | (2008), R89. |
| 2510 | |
| 2511 | AUTHORS: |
| 2512 | |
| 2513 | - Alex Raichev (2008-10-01, 2008-11-20, 2009-03-09, 2010-12-02, 2011-04-18) |
| 2514 | """ |
| 2515 | H= self._H |
| 2516 | R= H.parent() |
| 2517 | B= R.base_ring() |
| 2518 | d= R.ngens() |
| 2519 | vars= R.gens() |
| 2520 | f= R.ideal(H).radical().gens()[Integer(0) ] # Compute the reduction of H. |
| 2521 | |
| 2522 | # Expand B by the variables of alpha if there are any. |
| 2523 | indets= [] |
| 2524 | indets_ind= [] |
| 2525 | for a in alpha: |
| 2526 | if not ((a in ZZ) and (a>Integer(0) )): |
| 2527 | try: |
| 2528 | CC(a) |
| 2529 | except: |
| 2530 | indets.append(var(a)) |
| 2531 | indets_ind.append(alpha.index(a)) |
| 2532 | else: |
| 2533 | print "The components of", alpha, \ |
| 2534 | "must be positive integers or symbolic variables." |
| 2535 | return |
| 2536 | indets= list(Set(indets)) # Delete duplicates in indets. |
| 2537 | if indets != []: |
| 2538 | BB= FractionField(PolynomialRing(B,tuple(indets))) |
| 2539 | S= R.change_ring(BB) |
| 2540 | vars= S.gens() |
| 2541 | # Coerce alpha into BB. |
| 2542 | for i in range(len(alpha)): |
| 2543 | alpha[i] = BB(alpha[i]) |
| 2544 | else: |
| 2545 | S= R |
| 2546 | |
| 2547 | # Find smooth, critical points for alpha. |
| 2548 | f= S(f) |
| 2549 | J= S.ideal([f] +[ alpha[d-Integer(1) ]*vars[i]*diff(f,vars[i]) \ |
| 2550 | -alpha[i]*vars[d-Integer(1) ]*diff(f,vars[d-Integer(1) ]) for i in range(d-Integer(1) )]) |
| 2551 | return S.ideal(J.groebner_basis()) |
| 2552 | #------------------------------------------------------------------------------- |
| 2553 | def _subs_all(self,f,sub,simplify=False): |
| 2554 | r""" |
| 2555 | This function returns the items of `f` substituted by the dictionaries of |
| 2556 | `sub` in order of their appearance in `sub`. |
| 2557 | |
| 2558 | Does not use `self`. |
| 2559 | |
| 2560 | INPUT: |
| 2561 | |
| 2562 | - ``f`` - An individual or list of symbolic expressions or dictionaries |
| 2563 | - ``sub`` - An individual or list of dictionaries. |
| 2564 | - ``simplify`` - Boolean (default: False). |
| 2565 | |
| 2566 | OUTPUT: |
| 2567 | |
| 2568 | The items of `f` substituted by the dictionaries of `sub` in order of |
| 2569 | their appearance in `sub`. The subs() command is used. |
| 2570 | If simplify is True, then simplify() is used after substitution. |
| 2571 | |
| 2572 | EXAMPLES:: |
| 2573 | |
| 2574 | sage: R.<x> = PolynomialRing(QQ) |
| 2575 | sage: G,H = 1,1 |
| 2576 | sage: F= QuasiRationalExpression(G,H) |
| 2577 | sage: var('x,y,z') |
| 2578 | (x, y, z) |
| 2579 | sage: a= {x:1} |
| 2580 | sage: b= {y:2} |
| 2581 | sage: c= {z:3} |
| 2582 | sage: F._subs_all(x+y+z,a) |
| 2583 | y + z + 1 |
| 2584 | sage: F._subs_all(x+y+z,[c,a]) |
| 2585 | y + 4 |
| 2586 | sage: F._subs_all([x+y+z,y^2],b) |
| 2587 | [x + z + 2, 4] |
| 2588 | sage: F._subs_all([x+y+z,y^2],[b,c]) |
| 2589 | [x + 5, 4] |
| 2590 | |
| 2591 | :: |
| 2592 | |
| 2593 | sage: var('x,y') |
| 2594 | (x, y) |
| 2595 | sage: a= {'foo':x^2+y^2, 'bar':x-y} |
| 2596 | sage: b= {x:1,y:2} |
| 2597 | sage: F._subs_all(a,b) |
| 2598 | {'foo': 5, 'bar': -1} |
| 2599 | |
| 2600 | AUTHORS: |
| 2601 | |
| 2602 | - Alex Raichev (2009-05-05) |
| 2603 | """ |
| 2604 | singleton= False |
| 2605 | if not isinstance(f,list): |
| 2606 | f= [f] |
| 2607 | singleton= True |
| 2608 | if not isinstance(sub,list): |
| 2609 | sub= [sub] |
| 2610 | g= [] |
| 2611 | for ff in f: |
| 2612 | for D in sub: |
| 2613 | if isinstance(ff,dict): |
| 2614 | ff= dict( [(k,ff[k].subs(D)) for k in ff.keys()] ) |
| 2615 | else: |
| 2616 | ff= ff.subs(D) |
| 2617 | g.append(ff) |
| 2618 | if singleton and simplify: |
| 2619 | if isinstance(g[Integer(0) ],dict): |
| 2620 | return g[Integer(0) ] |
| 2621 | else: |
| 2622 | return g[Integer(0) ].simplify() |
| 2623 | elif singleton and not simplify: |
| 2624 | return g[Integer(0) ] |
| 2625 | elif not singleton and simplify: |
| 2626 | G= [] |
| 2627 | for gg in g: |
| 2628 | if isinstance(gg,dict): |
| 2629 | G.append(gg) |
| 2630 | else: |
| 2631 | G.append(gg.simplify()) |
| 2632 | return G |
| 2633 | else: |
| 2634 | return g |
| 2635 | # T ============================================================================ |
| 2636 | def maclaurin_coefficients(self,alpha,n): |
| 2637 | r""" |
| 2638 | Returns the Maclaurin coefficients of self that have multi-indices |
| 2639 | `alpha`, `2*alpha`,...,`n*alpha`. |
| 2640 | |
| 2641 | INPUT: |
| 2642 | |
| 2643 | - ``n`` - positive integer |
| 2644 | - ``alpha`` - tuple of positive integers representing a multi-index. |
| 2645 | |
| 2646 | OUTPUT: |
| 2647 | |
| 2648 | A dictionary of the form (beta, beta Maclaurin coefficient of self). |
| 2649 | |
| 2650 | AUTHORS: |
| 2651 | |
| 2652 | - Alex Raichev (2011-04-08) |
| 2653 | """ |
| 2654 | # Do all computations in the Symbolic Ring. |
| 2655 | F= SR(self._G/self._H) |
| 2656 | v= F.variables() |
| 2657 | p= {} |
| 2658 | for x in v: |
| 2659 | p[x]=Integer(0) |
| 2660 | d= len(v) |
| 2661 | coeffs={} |
| 2662 | # Initialize sequence of variables to differentiate with. |
| 2663 | s=[] |
| 2664 | for i in range(d): |
| 2665 | s.extend([v[i] for j in range(alpha[i])]) |
| 2666 | F_deriv= diff(F,s) |
| 2667 | coeffs[tuple(alpha)]= F_deriv.subs(p)/ mul([factorial(a) for a in alpha]) |
| 2668 | old_beta= alpha |
| 2669 | for k in (ellipsis_range(Integer(2) ,Ellipsis,n)): |
| 2670 | # Update variable sequence to differentiate with. |
| 2671 | beta= [k*a for a in alpha] |
| 2672 | delta= [beta[i]-old_beta[i] for i in range(d)] |
| 2673 | s= [] |
| 2674 | for i in range(d): |
| 2675 | s.extend([v[i] for j in range(delta[i])]) |
| 2676 | F_deriv= diff(F_deriv,s) |
| 2677 | coeffs[tuple(beta)]= F_deriv.subs(p)/ mul([factorial(b) for b in beta]) |
| 2678 | old_beta= copy(beta) |
| 2679 | return coeffs |
| 2680 | # U ============================================================================ |
| 2681 | # V ============================================================================ |
| 2682 | def variables(self): |
| 2683 | r""" |
| 2684 | Returns the tuple of variables of `self`. |
| 2685 | |
| 2686 | EXAMPLES:: |
| 2687 | |
| 2688 | sage: R.<x,y>= PolynomialRing(QQ) |
| 2689 | sage: G= exp(x) |
| 2690 | sage: H= 1-y |
| 2691 | sage: F= QuasiRationalExpression(G,H) |
| 2692 | sage: F.variables() |
| 2693 | (x, y) |
| 2694 | |
| 2695 | sage: R.<x,y>= PolynomialRing(QQ,order='invlex') |
| 2696 | sage: G= exp(x) |
| 2697 | sage: H= 1-y |
| 2698 | sage: F= QuasiRationalExpression(G,H) |
| 2699 | sage: F.variables() |
| 2700 | (y, x) |
| 2701 | |
| 2702 | AUTHORS: |
| 2703 | |
| 2704 | - Alex Raichev (2011-04-01) |
| 2705 | """ |
| 2706 | return self._variables |
| 2707 | # W ============================================================================ |
| 2708 | # X ============================================================================ |
| 2709 | # Y ============================================================================ |
| 2710 | # Z ============================================================================ |