1""" 2Expand Hypergeometric (and Meijer G) functions into named 3special functions. 4 5The algorithm for doing this uses a collection of lookup tables of 6hypergeometric functions, and various of their properties, to expand 7many hypergeometric functions in terms of special functions. 8 9It is based on the following paper: 10 Kelly B. Roach. Meijer G Function Representations. 11 In: Proceedings of the 1997 International Symposium on Symbolic and 12 Algebraic Computation, pages 205-211, New York, 1997. ACM. 13 14It is described in great(er) detail in the Sphinx documentation. 15""" 16# SUMMARY OF EXTENSIONS FOR MEIJER G FUNCTIONS 17# 18# o z**rho G(ap, bq; z) = G(ap + rho, bq + rho; z) 19# 20# o denote z*d/dz by D 21# 22# o It is helpful to keep in mind that ap and bq play essentially symmetric 23# roles: G(1/z) has slightly altered parameters, with ap and bq interchanged. 24# 25# o There are four shift operators: 26# A_J = b_J - D, J = 1, ..., n 27# B_J = 1 - a_j + D, J = 1, ..., m 28# C_J = -b_J + D, J = m+1, ..., q 29# D_J = a_J - 1 - D, J = n+1, ..., p 30# 31# A_J, C_J increment b_J 32# B_J, D_J decrement a_J 33# 34# o The corresponding four inverse-shift operators are defined if there 35# is no cancellation. Thus e.g. an index a_J (upper or lower) can be 36# incremented if a_J != b_i for i = 1, ..., q. 37# 38# o Order reduction: if b_j - a_i is a non-negative integer, where 39# j <= m and i > n, the corresponding quotient of gamma functions reduces 40# to a polynomial. Hence the G function can be expressed using a G-function 41# of lower order. 42# Similarly if j > m and i <= n. 43# 44# Secondly, there are paired index theorems [Adamchik, The evaluation of 45# integrals of Bessel functions via G-function identities]. Suppose there 46# are three parameters a, b, c, where a is an a_i, i <= n, b is a b_j, 47# j <= m and c is a denominator parameter (i.e. a_i, i > n or b_j, j > m). 48# Suppose further all three differ by integers. 49# Then the order can be reduced. 50# TODO work this out in detail. 51# 52# o An index quadruple is called suitable if its order cannot be reduced. 53# If there exists a sequence of shift operators transforming one index 54# quadruple into another, we say one is reachable from the other. 55# 56# o Deciding if one index quadruple is reachable from another is tricky. For 57# this reason, we use hand-built routines to match and instantiate formulas. 58# 59from collections import defaultdict 60from itertools import product 61 62from sympy import SYMPY_DEBUG 63from sympy.core import (S, Dummy, symbols, sympify, Tuple, expand, I, pi, Mul, 64 EulerGamma, oo, zoo, expand_func, Add, nan, Expr, Rational) 65from sympy.core.compatibility import default_sort_key, reduce 66from sympy.core.mod import Mod 67from sympy.functions import (exp, sqrt, root, log, lowergamma, cos, 68 besseli, gamma, uppergamma, expint, erf, sin, besselj, Ei, Ci, Si, Shi, 69 sinh, cosh, Chi, fresnels, fresnelc, polar_lift, exp_polar, floor, ceiling, 70 rf, factorial, lerchphi, Piecewise, re, elliptic_k, elliptic_e) 71from sympy.functions.elementary.complexes import polarify, unpolarify 72from sympy.functions.special.hyper import (hyper, HyperRep_atanh, 73 HyperRep_power1, HyperRep_power2, HyperRep_log1, HyperRep_asin1, 74 HyperRep_asin2, HyperRep_sqrts1, HyperRep_sqrts2, HyperRep_log2, 75 HyperRep_cosasin, HyperRep_sinasin, meijerg) 76from sympy.polys import poly, Poly 77from sympy.series import residue 78from sympy.simplify import simplify # type: ignore 79from sympy.simplify.powsimp import powdenest 80from sympy.utilities.iterables import sift 81 82# function to define "buckets" 83def _mod1(x): 84 # TODO see if this can work as Mod(x, 1); this will require 85 # different handling of the "buckets" since these need to 86 # be sorted and that fails when there is a mixture of 87 # integers and expressions with parameters. With the current 88 # Mod behavior, Mod(k, 1) == Mod(1, 1) == 0 if k is an integer. 89 # Although the sorting can be done with Basic.compare, this may 90 # still require different handling of the sorted buckets. 91 if x.is_Number: 92 return Mod(x, 1) 93 c, x = x.as_coeff_Add() 94 return Mod(c, 1) + x 95 96 97# leave add formulae at the top for easy reference 98def add_formulae(formulae): 99 """ Create our knowledge base. """ 100 from sympy.matrices import Matrix 101 102 a, b, c, z = symbols('a b c, z', cls=Dummy) 103 104 def add(ap, bq, res): 105 func = Hyper_Function(ap, bq) 106 formulae.append(Formula(func, z, res, (a, b, c))) 107 108 def addb(ap, bq, B, C, M): 109 func = Hyper_Function(ap, bq) 110 formulae.append(Formula(func, z, None, (a, b, c), B, C, M)) 111 112 # Luke, Y. L. (1969), The Special Functions and Their Approximations, 113 # Volume 1, section 6.2 114 115 # 0F0 116 add((), (), exp(z)) 117 118 # 1F0 119 add((a, ), (), HyperRep_power1(-a, z)) 120 121 # 2F1 122 addb((a, a - S.Half), (2*a, ), 123 Matrix([HyperRep_power2(a, z), 124 HyperRep_power2(a + S.Half, z)/2]), 125 Matrix([[1, 0]]), 126 Matrix([[(a - S.Half)*z/(1 - z), (S.Half - a)*z/(1 - z)], 127 [a/(1 - z), a*(z - 2)/(1 - z)]])) 128 addb((1, 1), (2, ), 129 Matrix([HyperRep_log1(z), 1]), Matrix([[-1/z, 0]]), 130 Matrix([[0, z/(z - 1)], [0, 0]])) 131 addb((S.Half, 1), (S('3/2'), ), 132 Matrix([HyperRep_atanh(z), 1]), 133 Matrix([[1, 0]]), 134 Matrix([[Rational(-1, 2), 1/(1 - z)/2], [0, 0]])) 135 addb((S.Half, S.Half), (S('3/2'), ), 136 Matrix([HyperRep_asin1(z), HyperRep_power1(Rational(-1, 2), z)]), 137 Matrix([[1, 0]]), 138 Matrix([[Rational(-1, 2), S.Half], [0, z/(1 - z)/2]])) 139 addb((a, S.Half + a), (S.Half, ), 140 Matrix([HyperRep_sqrts1(-a, z), -HyperRep_sqrts2(-a - S.Half, z)]), 141 Matrix([[1, 0]]), 142 Matrix([[0, -a], 143 [z*(-2*a - 1)/2/(1 - z), S.Half - z*(-2*a - 1)/(1 - z)]])) 144 145 # A. P. Prudnikov, Yu. A. Brychkov and O. I. Marichev (1990). 146 # Integrals and Series: More Special Functions, Vol. 3,. 147 # Gordon and Breach Science Publisher 148 addb([a, -a], [S.Half], 149 Matrix([HyperRep_cosasin(a, z), HyperRep_sinasin(a, z)]), 150 Matrix([[1, 0]]), 151 Matrix([[0, -a], [a*z/(1 - z), 1/(1 - z)/2]])) 152 addb([1, 1], [3*S.Half], 153 Matrix([HyperRep_asin2(z), 1]), Matrix([[1, 0]]), 154 Matrix([[(z - S.Half)/(1 - z), 1/(1 - z)/2], [0, 0]])) 155 156 # Complete elliptic integrals K(z) and E(z), both a 2F1 function 157 addb([S.Half, S.Half], [S.One], 158 Matrix([elliptic_k(z), elliptic_e(z)]), 159 Matrix([[2/pi, 0]]), 160 Matrix([[Rational(-1, 2), -1/(2*z-2)], 161 [Rational(-1, 2), S.Half]])) 162 addb([Rational(-1, 2), S.Half], [S.One], 163 Matrix([elliptic_k(z), elliptic_e(z)]), 164 Matrix([[0, 2/pi]]), 165 Matrix([[Rational(-1, 2), -1/(2*z-2)], 166 [Rational(-1, 2), S.Half]])) 167 168 # 3F2 169 addb([Rational(-1, 2), 1, 1], [S.Half, 2], 170 Matrix([z*HyperRep_atanh(z), HyperRep_log1(z), 1]), 171 Matrix([[Rational(-2, 3), -S.One/(3*z), Rational(2, 3)]]), 172 Matrix([[S.Half, 0, z/(1 - z)/2], 173 [0, 0, z/(z - 1)], 174 [0, 0, 0]])) 175 # actually the formula for 3/2 is much nicer ... 176 addb([Rational(-1, 2), 1, 1], [2, 2], 177 Matrix([HyperRep_power1(S.Half, z), HyperRep_log2(z), 1]), 178 Matrix([[Rational(4, 9) - 16/(9*z), 4/(3*z), 16/(9*z)]]), 179 Matrix([[z/2/(z - 1), 0, 0], [1/(2*(z - 1)), 0, S.Half], [0, 0, 0]])) 180 181 # 1F1 182 addb([1], [b], Matrix([z**(1 - b) * exp(z) * lowergamma(b - 1, z), 1]), 183 Matrix([[b - 1, 0]]), Matrix([[1 - b + z, 1], [0, 0]])) 184 addb([a], [2*a], 185 Matrix([z**(S.Half - a)*exp(z/2)*besseli(a - S.Half, z/2) 186 * gamma(a + S.Half)/4**(S.Half - a), 187 z**(S.Half - a)*exp(z/2)*besseli(a + S.Half, z/2) 188 * gamma(a + S.Half)/4**(S.Half - a)]), 189 Matrix([[1, 0]]), 190 Matrix([[z/2, z/2], [z/2, (z/2 - 2*a)]])) 191 mz = polar_lift(-1)*z 192 addb([a], [a + 1], 193 Matrix([mz**(-a)*a*lowergamma(a, mz), a*exp(z)]), 194 Matrix([[1, 0]]), 195 Matrix([[-a, 1], [0, z]])) 196 # This one is redundant. 197 add([Rational(-1, 2)], [S.Half], exp(z) - sqrt(pi*z)*(-I)*erf(I*sqrt(z))) 198 199 # Added to get nice results for Laplace transform of Fresnel functions 200 # http://functions.wolfram.com/07.22.03.6437.01 201 # Basic rule 202 #add([1], [Rational(3, 4), Rational(5, 4)], 203 # sqrt(pi) * (cos(2*sqrt(polar_lift(-1)*z))*fresnelc(2*root(polar_lift(-1)*z,4)/sqrt(pi)) + 204 # sin(2*sqrt(polar_lift(-1)*z))*fresnels(2*root(polar_lift(-1)*z,4)/sqrt(pi))) 205 # / (2*root(polar_lift(-1)*z,4))) 206 # Manually tuned rule 207 addb([1], [Rational(3, 4), Rational(5, 4)], 208 Matrix([ sqrt(pi)*(I*sinh(2*sqrt(z))*fresnels(2*root(z, 4)*exp(I*pi/4)/sqrt(pi)) 209 + cosh(2*sqrt(z))*fresnelc(2*root(z, 4)*exp(I*pi/4)/sqrt(pi))) 210 * exp(-I*pi/4)/(2*root(z, 4)), 211 sqrt(pi)*root(z, 4)*(sinh(2*sqrt(z))*fresnelc(2*root(z, 4)*exp(I*pi/4)/sqrt(pi)) 212 + I*cosh(2*sqrt(z))*fresnels(2*root(z, 4)*exp(I*pi/4)/sqrt(pi))) 213 *exp(-I*pi/4)/2, 214 1 ]), 215 Matrix([[1, 0, 0]]), 216 Matrix([[Rational(-1, 4), 1, Rational(1, 4)], 217 [ z, Rational(1, 4), 0], 218 [ 0, 0, 0]])) 219 220 # 2F2 221 addb([S.Half, a], [Rational(3, 2), a + 1], 222 Matrix([a/(2*a - 1)*(-I)*sqrt(pi/z)*erf(I*sqrt(z)), 223 a/(2*a - 1)*(polar_lift(-1)*z)**(-a)* 224 lowergamma(a, polar_lift(-1)*z), 225 a/(2*a - 1)*exp(z)]), 226 Matrix([[1, -1, 0]]), 227 Matrix([[Rational(-1, 2), 0, 1], [0, -a, 1], [0, 0, z]])) 228 # We make a "basis" of four functions instead of three, and give EulerGamma 229 # an extra slot (it could just be a coefficient to 1). The advantage is 230 # that this way Polys will not see multivariate polynomials (it treats 231 # EulerGamma as an indeterminate), which is *way* faster. 232 addb([1, 1], [2, 2], 233 Matrix([Ei(z) - log(z), exp(z), 1, EulerGamma]), 234 Matrix([[1/z, 0, 0, -1/z]]), 235 Matrix([[0, 1, -1, 0], [0, z, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]])) 236 237 # 0F1 238 add((), (S.Half, ), cosh(2*sqrt(z))) 239 addb([], [b], 240 Matrix([gamma(b)*z**((1 - b)/2)*besseli(b - 1, 2*sqrt(z)), 241 gamma(b)*z**(1 - b/2)*besseli(b, 2*sqrt(z))]), 242 Matrix([[1, 0]]), Matrix([[0, 1], [z, (1 - b)]])) 243 244 # 0F3 245 x = 4*z**Rational(1, 4) 246 247 def fp(a, z): 248 return besseli(a, x) + besselj(a, x) 249 250 def fm(a, z): 251 return besseli(a, x) - besselj(a, x) 252 253 # TODO branching 254 addb([], [S.Half, a, a + S.Half], 255 Matrix([fp(2*a - 1, z), fm(2*a, z)*z**Rational(1, 4), 256 fm(2*a - 1, z)*sqrt(z), fp(2*a, z)*z**Rational(3, 4)]) 257 * 2**(-2*a)*gamma(2*a)*z**((1 - 2*a)/4), 258 Matrix([[1, 0, 0, 0]]), 259 Matrix([[0, 1, 0, 0], 260 [0, S.Half - a, 1, 0], 261 [0, 0, S.Half, 1], 262 [z, 0, 0, 1 - a]])) 263 x = 2*(4*z)**Rational(1, 4)*exp_polar(I*pi/4) 264 addb([], [a, a + S.Half, 2*a], 265 (2*sqrt(polar_lift(-1)*z))**(1 - 2*a)*gamma(2*a)**2 * 266 Matrix([besselj(2*a - 1, x)*besseli(2*a - 1, x), 267 x*(besseli(2*a, x)*besselj(2*a - 1, x) 268 - besseli(2*a - 1, x)*besselj(2*a, x)), 269 x**2*besseli(2*a, x)*besselj(2*a, x), 270 x**3*(besseli(2*a, x)*besselj(2*a - 1, x) 271 + besseli(2*a - 1, x)*besselj(2*a, x))]), 272 Matrix([[1, 0, 0, 0]]), 273 Matrix([[0, Rational(1, 4), 0, 0], 274 [0, (1 - 2*a)/2, Rational(-1, 2), 0], 275 [0, 0, 1 - 2*a, Rational(1, 4)], 276 [-32*z, 0, 0, 1 - a]])) 277 278 # 1F2 279 addb([a], [a - S.Half, 2*a], 280 Matrix([z**(S.Half - a)*besseli(a - S.Half, sqrt(z))**2, 281 z**(1 - a)*besseli(a - S.Half, sqrt(z)) 282 *besseli(a - Rational(3, 2), sqrt(z)), 283 z**(Rational(3, 2) - a)*besseli(a - Rational(3, 2), sqrt(z))**2]), 284 Matrix([[-gamma(a + S.Half)**2/4**(S.Half - a), 285 2*gamma(a - S.Half)*gamma(a + S.Half)/4**(1 - a), 286 0]]), 287 Matrix([[1 - 2*a, 1, 0], [z/2, S.Half - a, S.Half], [0, z, 0]])) 288 addb([S.Half], [b, 2 - b], 289 pi*(1 - b)/sin(pi*b)* 290 Matrix([besseli(1 - b, sqrt(z))*besseli(b - 1, sqrt(z)), 291 sqrt(z)*(besseli(-b, sqrt(z))*besseli(b - 1, sqrt(z)) 292 + besseli(1 - b, sqrt(z))*besseli(b, sqrt(z))), 293 besseli(-b, sqrt(z))*besseli(b, sqrt(z))]), 294 Matrix([[1, 0, 0]]), 295 Matrix([[b - 1, S.Half, 0], 296 [z, 0, z], 297 [0, S.Half, -b]])) 298 addb([S.Half], [Rational(3, 2), Rational(3, 2)], 299 Matrix([Shi(2*sqrt(z))/2/sqrt(z), sinh(2*sqrt(z))/2/sqrt(z), 300 cosh(2*sqrt(z))]), 301 Matrix([[1, 0, 0]]), 302 Matrix([[Rational(-1, 2), S.Half, 0], [0, Rational(-1, 2), S.Half], [0, 2*z, 0]])) 303 304 # FresnelS 305 # Basic rule 306 #add([Rational(3, 4)], [Rational(3, 2),Rational(7, 4)], 6*fresnels( exp(pi*I/4)*root(z,4)*2/sqrt(pi) ) / ( pi * (exp(pi*I/4)*root(z,4)*2/sqrt(pi))**3 ) ) 307 # Manually tuned rule 308 addb([Rational(3, 4)], [Rational(3, 2), Rational(7, 4)], 309 Matrix( 310 [ fresnels( 311 exp( 312 pi*I/4)*root( 313 z, 4)*2/sqrt( 314 pi) ) / ( 315 pi * (exp(pi*I/4)*root(z, 4)*2/sqrt(pi))**3 ), 316 sinh(2*sqrt(z))/sqrt(z), 317 cosh(2*sqrt(z)) ]), 318 Matrix([[6, 0, 0]]), 319 Matrix([[Rational(-3, 4), Rational(1, 16), 0], 320 [ 0, Rational(-1, 2), 1], 321 [ 0, z, 0]])) 322 323 # FresnelC 324 # Basic rule 325 #add([Rational(1, 4)], [S.Half,Rational(5, 4)], fresnelc( exp(pi*I/4)*root(z,4)*2/sqrt(pi) ) / ( exp(pi*I/4)*root(z,4)*2/sqrt(pi) ) ) 326 # Manually tuned rule 327 addb([Rational(1, 4)], [S.Half, Rational(5, 4)], 328 Matrix( 329 [ sqrt( 330 pi)*exp( 331 -I*pi/4)*fresnelc( 332 2*root(z, 4)*exp(I*pi/4)/sqrt(pi))/(2*root(z, 4)), 333 cosh(2*sqrt(z)), 334 sinh(2*sqrt(z))*sqrt(z) ]), 335 Matrix([[1, 0, 0]]), 336 Matrix([[Rational(-1, 4), Rational(1, 4), 0 ], 337 [ 0, 0, 1 ], 338 [ 0, z, S.Half]])) 339 340 # 2F3 341 # XXX with this five-parameter formula is pretty slow with the current 342 # Formula.find_instantiations (creates 2!*3!*3**(2+3) ~ 3000 343 # instantiations ... But it's not too bad. 344 addb([a, a + S.Half], [2*a, b, 2*a - b + 1], 345 gamma(b)*gamma(2*a - b + 1) * (sqrt(z)/2)**(1 - 2*a) * 346 Matrix([besseli(b - 1, sqrt(z))*besseli(2*a - b, sqrt(z)), 347 sqrt(z)*besseli(b, sqrt(z))*besseli(2*a - b, sqrt(z)), 348 sqrt(z)*besseli(b - 1, sqrt(z))*besseli(2*a - b + 1, sqrt(z)), 349 besseli(b, sqrt(z))*besseli(2*a - b + 1, sqrt(z))]), 350 Matrix([[1, 0, 0, 0]]), 351 Matrix([[0, S.Half, S.Half, 0], 352 [z/2, 1 - b, 0, z/2], 353 [z/2, 0, b - 2*a, z/2], 354 [0, S.Half, S.Half, -2*a]])) 355 # (C/f above comment about eulergamma in the basis). 356 addb([1, 1], [2, 2, Rational(3, 2)], 357 Matrix([Chi(2*sqrt(z)) - log(2*sqrt(z)), 358 cosh(2*sqrt(z)), sqrt(z)*sinh(2*sqrt(z)), 1, EulerGamma]), 359 Matrix([[1/z, 0, 0, 0, -1/z]]), 360 Matrix([[0, S.Half, 0, Rational(-1, 2), 0], 361 [0, 0, 1, 0, 0], 362 [0, z, S.Half, 0, 0], 363 [0, 0, 0, 0, 0], 364 [0, 0, 0, 0, 0]])) 365 366 # 3F3 367 # This is rule: http://functions.wolfram.com/07.31.03.0134.01 368 # Initial reason to add it was a nice solution for 369 # integrate(erf(a*z)/z**2, z) and same for erfc and erfi. 370 # Basic rule 371 # add([1, 1, a], [2, 2, a+1], (a/(z*(a-1)**2)) * 372 # (1 - (-z)**(1-a) * (gamma(a) - uppergamma(a,-z)) 373 # - (a-1) * (EulerGamma + uppergamma(0,-z) + log(-z)) 374 # - exp(z))) 375 # Manually tuned rule 376 addb([1, 1, a], [2, 2, a+1], 377 Matrix([a*(log(-z) + expint(1, -z) + EulerGamma)/(z*(a**2 - 2*a + 1)), 378 a*(-z)**(-a)*(gamma(a) - uppergamma(a, -z))/(a - 1)**2, 379 a*exp(z)/(a**2 - 2*a + 1), 380 a/(z*(a**2 - 2*a + 1))]), 381 Matrix([[1-a, 1, -1/z, 1]]), 382 Matrix([[-1,0,-1/z,1], 383 [0,-a,1,0], 384 [0,0,z,0], 385 [0,0,0,-1]])) 386 387 388def add_meijerg_formulae(formulae): 389 from sympy.matrices import Matrix 390 391 a, b, c, z = list(map(Dummy, 'abcz')) 392 rho = Dummy('rho') 393 394 def add(an, ap, bm, bq, B, C, M, matcher): 395 formulae.append(MeijerFormula(an, ap, bm, bq, z, [a, b, c, rho], 396 B, C, M, matcher)) 397 398 def detect_uppergamma(func): 399 x = func.an[0] 400 y, z = func.bm 401 swapped = False 402 if not _mod1((x - y).simplify()): 403 swapped = True 404 (y, z) = (z, y) 405 if _mod1((x - z).simplify()) or x - z > 0: 406 return None 407 l = [y, x] 408 if swapped: 409 l = [x, y] 410 return {rho: y, a: x - y}, G_Function([x], [], l, []) 411 412 add([a + rho], [], [rho, a + rho], [], 413 Matrix([gamma(1 - a)*z**rho*exp(z)*uppergamma(a, z), 414 gamma(1 - a)*z**(a + rho)]), 415 Matrix([[1, 0]]), 416 Matrix([[rho + z, -1], [0, a + rho]]), 417 detect_uppergamma) 418 419 def detect_3113(func): 420 """http://functions.wolfram.com/07.34.03.0984.01""" 421 x = func.an[0] 422 u, v, w = func.bm 423 if _mod1((u - v).simplify()) == 0: 424 if _mod1((v - w).simplify()) == 0: 425 return 426 sig = (S.Half, S.Half, S.Zero) 427 x1, x2, y = u, v, w 428 else: 429 if _mod1((x - u).simplify()) == 0: 430 sig = (S.Half, S.Zero, S.Half) 431 x1, y, x2 = u, v, w 432 else: 433 sig = (S.Zero, S.Half, S.Half) 434 y, x1, x2 = u, v, w 435 436 if (_mod1((x - x1).simplify()) != 0 or 437 _mod1((x - x2).simplify()) != 0 or 438 _mod1((x - y).simplify()) != S.Half or 439 x - x1 > 0 or x - x2 > 0): 440 return 441 442 return {a: x}, G_Function([x], [], [x - S.Half + t for t in sig], []) 443 444 s = sin(2*sqrt(z)) 445 c_ = cos(2*sqrt(z)) 446 S_ = Si(2*sqrt(z)) - pi/2 447 C = Ci(2*sqrt(z)) 448 add([a], [], [a, a, a - S.Half], [], 449 Matrix([sqrt(pi)*z**(a - S.Half)*(c_*S_ - s*C), 450 sqrt(pi)*z**a*(s*S_ + c_*C), 451 sqrt(pi)*z**a]), 452 Matrix([[-2, 0, 0]]), 453 Matrix([[a - S.Half, -1, 0], [z, a, S.Half], [0, 0, a]]), 454 detect_3113) 455 456 457def make_simp(z): 458 """ Create a function that simplifies rational functions in ``z``. """ 459 460 def simp(expr): 461 """ Efficiently simplify the rational function ``expr``. """ 462 numer, denom = expr.as_numer_denom() 463 numer = numer.expand() 464 # denom = denom.expand() # is this needed? 465 c, numer, denom = poly(numer, z).cancel(poly(denom, z)) 466 return c * numer.as_expr() / denom.as_expr() 467 468 return simp 469 470 471def debug(*args): 472 if SYMPY_DEBUG: 473 for a in args: 474 print(a, end="") 475 print() 476 477 478class Hyper_Function(Expr): 479 """ A generalized hypergeometric function. """ 480 481 def __new__(cls, ap, bq): 482 obj = super().__new__(cls) 483 obj.ap = Tuple(*list(map(expand, ap))) 484 obj.bq = Tuple(*list(map(expand, bq))) 485 return obj 486 487 @property 488 def args(self): 489 return (self.ap, self.bq) 490 491 @property 492 def sizes(self): 493 return (len(self.ap), len(self.bq)) 494 495 @property 496 def gamma(self): 497 """ 498 Number of upper parameters that are negative integers 499 500 This is a transformation invariant. 501 """ 502 return sum(bool(x.is_integer and x.is_negative) for x in self.ap) 503 504 def _hashable_content(self): 505 return super()._hashable_content() + (self.ap, 506 self.bq) 507 508 def __call__(self, arg): 509 return hyper(self.ap, self.bq, arg) 510 511 def build_invariants(self): 512 """ 513 Compute the invariant vector. 514 515 Explanation 516 =========== 517 518 The invariant vector is: 519 (gamma, ((s1, n1), ..., (sk, nk)), ((t1, m1), ..., (tr, mr))) 520 where gamma is the number of integer a < 0, 521 s1 < ... < sk 522 nl is the number of parameters a_i congruent to sl mod 1 523 t1 < ... < tr 524 ml is the number of parameters b_i congruent to tl mod 1 525 526 If the index pair contains parameters, then this is not truly an 527 invariant, since the parameters cannot be sorted uniquely mod1. 528 529 Examples 530 ======== 531 532 >>> from sympy.simplify.hyperexpand import Hyper_Function 533 >>> from sympy import S 534 >>> ap = (S.Half, S.One/3, S(-1)/2, -2) 535 >>> bq = (1, 2) 536 537 Here gamma = 1, 538 k = 3, s1 = 0, s2 = 1/3, s3 = 1/2 539 n1 = 1, n2 = 1, n2 = 2 540 r = 1, t1 = 0 541 m1 = 2: 542 543 >>> Hyper_Function(ap, bq).build_invariants() 544 (1, ((0, 1), (1/3, 1), (1/2, 2)), ((0, 2),)) 545 """ 546 abuckets, bbuckets = sift(self.ap, _mod1), sift(self.bq, _mod1) 547 548 def tr(bucket): 549 bucket = list(bucket.items()) 550 if not any(isinstance(x[0], Mod) for x in bucket): 551 bucket.sort(key=lambda x: default_sort_key(x[0])) 552 bucket = tuple([(mod, len(values)) for mod, values in bucket if 553 values]) 554 return bucket 555 556 return (self.gamma, tr(abuckets), tr(bbuckets)) 557 558 def difficulty(self, func): 559 """ Estimate how many steps it takes to reach ``func`` from self. 560 Return -1 if impossible. """ 561 if self.gamma != func.gamma: 562 return -1 563 oabuckets, obbuckets, abuckets, bbuckets = [sift(params, _mod1) for 564 params in (self.ap, self.bq, func.ap, func.bq)] 565 566 diff = 0 567 for bucket, obucket in [(abuckets, oabuckets), (bbuckets, obbuckets)]: 568 for mod in set(list(bucket.keys()) + list(obucket.keys())): 569 if (not mod in bucket) or (not mod in obucket) \ 570 or len(bucket[mod]) != len(obucket[mod]): 571 return -1 572 l1 = list(bucket[mod]) 573 l2 = list(obucket[mod]) 574 l1.sort() 575 l2.sort() 576 for i, j in zip(l1, l2): 577 diff += abs(i - j) 578 579 return diff 580 581 def _is_suitable_origin(self): 582 """ 583 Decide if ``self`` is a suitable origin. 584 585 Explanation 586 =========== 587 588 A function is a suitable origin iff: 589 * none of the ai equals bj + n, with n a non-negative integer 590 * none of the ai is zero 591 * none of the bj is a non-positive integer 592 593 Note that this gives meaningful results only when none of the indices 594 are symbolic. 595 596 """ 597 for a in self.ap: 598 for b in self.bq: 599 if (a - b).is_integer and (a - b).is_negative is False: 600 return False 601 for a in self.ap: 602 if a == 0: 603 return False 604 for b in self.bq: 605 if b.is_integer and b.is_nonpositive: 606 return False 607 return True 608 609 610class G_Function(Expr): 611 """ A Meijer G-function. """ 612 613 def __new__(cls, an, ap, bm, bq): 614 obj = super().__new__(cls) 615 obj.an = Tuple(*list(map(expand, an))) 616 obj.ap = Tuple(*list(map(expand, ap))) 617 obj.bm = Tuple(*list(map(expand, bm))) 618 obj.bq = Tuple(*list(map(expand, bq))) 619 return obj 620 621 @property 622 def args(self): 623 return (self.an, self.ap, self.bm, self.bq) 624 625 def _hashable_content(self): 626 return super()._hashable_content() + self.args 627 628 def __call__(self, z): 629 return meijerg(self.an, self.ap, self.bm, self.bq, z) 630 631 def compute_buckets(self): 632 """ 633 Compute buckets for the fours sets of parameters. 634 635 Explanation 636 =========== 637 638 We guarantee that any two equal Mod objects returned are actually the 639 same, and that the buckets are sorted by real part (an and bq 640 descendending, bm and ap ascending). 641 642 Examples 643 ======== 644 645 >>> from sympy.simplify.hyperexpand import G_Function 646 >>> from sympy.abc import y 647 >>> from sympy import S 648 649 >>> a, b = [1, 3, 2, S(3)/2], [1 + y, y, 2, y + 3] 650 >>> G_Function(a, b, [2], [y]).compute_buckets() 651 ({0: [3, 2, 1], 1/2: [3/2]}, 652 {0: [2], y: [y, y + 1, y + 3]}, {0: [2]}, {y: [y]}) 653 654 """ 655 dicts = pan, pap, pbm, pbq = [defaultdict(list) for i in range(4)] 656 for dic, lis in zip(dicts, (self.an, self.ap, self.bm, self.bq)): 657 for x in lis: 658 dic[_mod1(x)].append(x) 659 660 for dic, flip in zip(dicts, (True, False, False, True)): 661 for m, items in dic.items(): 662 x0 = items[0] 663 items.sort(key=lambda x: x - x0, reverse=flip) 664 dic[m] = items 665 666 return tuple([dict(w) for w in dicts]) 667 668 @property 669 def signature(self): 670 return (len(self.an), len(self.ap), len(self.bm), len(self.bq)) 671 672 673# Dummy variable. 674_x = Dummy('x') 675 676class Formula: 677 """ 678 This class represents hypergeometric formulae. 679 680 Explanation 681 =========== 682 683 Its data members are: 684 - z, the argument 685 - closed_form, the closed form expression 686 - symbols, the free symbols (parameters) in the formula 687 - func, the function 688 - B, C, M (see _compute_basis) 689 690 Examples 691 ======== 692 693 >>> from sympy.abc import a, b, z 694 >>> from sympy.simplify.hyperexpand import Formula, Hyper_Function 695 >>> func = Hyper_Function((a/2, a/3 + b, (1+a)/2), (a, b, (a+b)/7)) 696 >>> f = Formula(func, z, None, [a, b]) 697 698 """ 699 700 def _compute_basis(self, closed_form): 701 """ 702 Compute a set of functions B=(f1, ..., fn), a nxn matrix M 703 and a 1xn matrix C such that: 704 closed_form = C B 705 z d/dz B = M B. 706 """ 707 from sympy.matrices import Matrix, eye, zeros 708 709 afactors = [_x + a for a in self.func.ap] 710 bfactors = [_x + b - 1 for b in self.func.bq] 711 expr = _x*Mul(*bfactors) - self.z*Mul(*afactors) 712 poly = Poly(expr, _x) 713 714 n = poly.degree() - 1 715 b = [closed_form] 716 for _ in range(n): 717 b.append(self.z*b[-1].diff(self.z)) 718 719 self.B = Matrix(b) 720 self.C = Matrix([[1] + [0]*n]) 721 722 m = eye(n) 723 m = m.col_insert(0, zeros(n, 1)) 724 l = poly.all_coeffs()[1:] 725 l.reverse() 726 self.M = m.row_insert(n, -Matrix([l])/poly.all_coeffs()[0]) 727 728 def __init__(self, func, z, res, symbols, B=None, C=None, M=None): 729 z = sympify(z) 730 res = sympify(res) 731 symbols = [x for x in sympify(symbols) if func.has(x)] 732 733 self.z = z 734 self.symbols = symbols 735 self.B = B 736 self.C = C 737 self.M = M 738 self.func = func 739 740 # TODO with symbolic parameters, it could be advantageous 741 # (for prettier answers) to compute a basis only *after* 742 # instantiation 743 if res is not None: 744 self._compute_basis(res) 745 746 @property 747 def closed_form(self): 748 return reduce(lambda s,m: s+m[0]*m[1], zip(self.C, self.B), S.Zero) 749 750 def find_instantiations(self, func): 751 """ 752 Find substitutions of the free symbols that match ``func``. 753 754 Return the substitution dictionaries as a list. Note that the returned 755 instantiations need not actually match, or be valid! 756 757 """ 758 from sympy.solvers import solve 759 ap = func.ap 760 bq = func.bq 761 if len(ap) != len(self.func.ap) or len(bq) != len(self.func.bq): 762 raise TypeError('Cannot instantiate other number of parameters') 763 symbol_values = [] 764 for a in self.symbols: 765 if a in self.func.ap.args: 766 symbol_values.append(ap) 767 elif a in self.func.bq.args: 768 symbol_values.append(bq) 769 else: 770 raise ValueError("At least one of the parameters of the " 771 "formula must be equal to %s" % (a,)) 772 base_repl = [dict(list(zip(self.symbols, values))) 773 for values in product(*symbol_values)] 774 abuckets, bbuckets = [sift(params, _mod1) for params in [ap, bq]] 775 a_inv, b_inv = [{a: len(vals) for a, vals in bucket.items()} 776 for bucket in [abuckets, bbuckets]] 777 critical_values = [[0] for _ in self.symbols] 778 result = [] 779 _n = Dummy() 780 for repl in base_repl: 781 symb_a, symb_b = [sift(params, lambda x: _mod1(x.xreplace(repl))) 782 for params in [self.func.ap, self.func.bq]] 783 for bucket, obucket in [(abuckets, symb_a), (bbuckets, symb_b)]: 784 for mod in set(list(bucket.keys()) + list(obucket.keys())): 785 if (not mod in bucket) or (not mod in obucket) \ 786 or len(bucket[mod]) != len(obucket[mod]): 787 break 788 for a, vals in zip(self.symbols, critical_values): 789 if repl[a].free_symbols: 790 continue 791 exprs = [expr for expr in obucket[mod] if expr.has(a)] 792 repl0 = repl.copy() 793 repl0[a] += _n 794 for expr in exprs: 795 for target in bucket[mod]: 796 n0, = solve(expr.xreplace(repl0) - target, _n) 797 if n0.free_symbols: 798 raise ValueError("Value should not be true") 799 vals.append(n0) 800 else: 801 values = [] 802 for a, vals in zip(self.symbols, critical_values): 803 a0 = repl[a] 804 min_ = floor(min(vals)) 805 max_ = ceiling(max(vals)) 806 values.append([a0 + n for n in range(min_, max_ + 1)]) 807 result.extend(dict(list(zip(self.symbols, l))) for l in product(*values)) 808 return result 809 810 811 812 813class FormulaCollection: 814 """ A collection of formulae to use as origins. """ 815 816 def __init__(self): 817 """ Doing this globally at module init time is a pain ... """ 818 self.symbolic_formulae = {} 819 self.concrete_formulae = {} 820 self.formulae = [] 821 822 add_formulae(self.formulae) 823 824 # Now process the formulae into a helpful form. 825 # These dicts are indexed by (p, q). 826 827 for f in self.formulae: 828 sizes = f.func.sizes 829 if len(f.symbols) > 0: 830 self.symbolic_formulae.setdefault(sizes, []).append(f) 831 else: 832 inv = f.func.build_invariants() 833 self.concrete_formulae.setdefault(sizes, {})[inv] = f 834 835 def lookup_origin(self, func): 836 """ 837 Given the suitable target ``func``, try to find an origin in our 838 knowledge base. 839 840 Examples 841 ======== 842 843 >>> from sympy.simplify.hyperexpand import (FormulaCollection, 844 ... Hyper_Function) 845 >>> f = FormulaCollection() 846 >>> f.lookup_origin(Hyper_Function((), ())).closed_form 847 exp(_z) 848 >>> f.lookup_origin(Hyper_Function([1], ())).closed_form 849 HyperRep_power1(-1, _z) 850 851 >>> from sympy import S 852 >>> i = Hyper_Function([S('1/4'), S('3/4 + 4')], [S.Half]) 853 >>> f.lookup_origin(i).closed_form 854 HyperRep_sqrts1(-1/4, _z) 855 """ 856 inv = func.build_invariants() 857 sizes = func.sizes 858 if sizes in self.concrete_formulae and \ 859 inv in self.concrete_formulae[sizes]: 860 return self.concrete_formulae[sizes][inv] 861 862 # We don't have a concrete formula. Try to instantiate. 863 if not sizes in self.symbolic_formulae: 864 return None # Too bad... 865 866 possible = [] 867 for f in self.symbolic_formulae[sizes]: 868 repls = f.find_instantiations(func) 869 for repl in repls: 870 func2 = f.func.xreplace(repl) 871 if not func2._is_suitable_origin(): 872 continue 873 diff = func2.difficulty(func) 874 if diff == -1: 875 continue 876 possible.append((diff, repl, f, func2)) 877 878 # find the nearest origin 879 possible.sort(key=lambda x: x[0]) 880 for _, repl, f, func2 in possible: 881 f2 = Formula(func2, f.z, None, [], f.B.subs(repl), 882 f.C.subs(repl), f.M.subs(repl)) 883 if not any(e.has(S.NaN, oo, -oo, zoo) for e in [f2.B, f2.M, f2.C]): 884 return f2 885 886 return None 887 888 889class MeijerFormula: 890 """ 891 This class represents a Meijer G-function formula. 892 893 Its data members are: 894 - z, the argument 895 - symbols, the free symbols (parameters) in the formula 896 - func, the function 897 - B, C, M (c/f ordinary Formula) 898 """ 899 900 def __init__(self, an, ap, bm, bq, z, symbols, B, C, M, matcher): 901 an, ap, bm, bq = [Tuple(*list(map(expand, w))) for w in [an, ap, bm, bq]] 902 self.func = G_Function(an, ap, bm, bq) 903 self.z = z 904 self.symbols = symbols 905 self._matcher = matcher 906 self.B = B 907 self.C = C 908 self.M = M 909 910 @property 911 def closed_form(self): 912 return reduce(lambda s,m: s+m[0]*m[1], zip(self.C, self.B), S.Zero) 913 914 def try_instantiate(self, func): 915 """ 916 Try to instantiate the current formula to (almost) match func. 917 This uses the _matcher passed on init. 918 """ 919 if func.signature != self.func.signature: 920 return None 921 res = self._matcher(func) 922 if res is not None: 923 subs, newfunc = res 924 return MeijerFormula(newfunc.an, newfunc.ap, newfunc.bm, newfunc.bq, 925 self.z, [], 926 self.B.subs(subs), self.C.subs(subs), 927 self.M.subs(subs), None) 928 929 930class MeijerFormulaCollection: 931 """ 932 This class holds a collection of meijer g formulae. 933 """ 934 935 def __init__(self): 936 formulae = [] 937 add_meijerg_formulae(formulae) 938 self.formulae = defaultdict(list) 939 for formula in formulae: 940 self.formulae[formula.func.signature].append(formula) 941 self.formulae = dict(self.formulae) 942 943 def lookup_origin(self, func): 944 """ Try to find a formula that matches func. """ 945 if not func.signature in self.formulae: 946 return None 947 for formula in self.formulae[func.signature]: 948 res = formula.try_instantiate(func) 949 if res is not None: 950 return res 951 952 953class Operator: 954 """ 955 Base class for operators to be applied to our functions. 956 957 Explanation 958 =========== 959 960 These operators are differential operators. They are by convention 961 expressed in the variable D = z*d/dz (although this base class does 962 not actually care). 963 Note that when the operator is applied to an object, we typically do 964 *not* blindly differentiate but instead use a different representation 965 of the z*d/dz operator (see make_derivative_operator). 966 967 To subclass from this, define a __init__ method that initializes a 968 self._poly variable. This variable stores a polynomial. By convention 969 the generator is z*d/dz, and acts to the right of all coefficients. 970 971 Thus this poly 972 x**2 + 2*z*x + 1 973 represents the differential operator 974 (z*d/dz)**2 + 2*z**2*d/dz. 975 976 This class is used only in the implementation of the hypergeometric 977 function expansion algorithm. 978 """ 979 980 def apply(self, obj, op): 981 """ 982 Apply ``self`` to the object ``obj``, where the generator is ``op``. 983 984 Examples 985 ======== 986 987 >>> from sympy.simplify.hyperexpand import Operator 988 >>> from sympy.polys.polytools import Poly 989 >>> from sympy.abc import x, y, z 990 >>> op = Operator() 991 >>> op._poly = Poly(x**2 + z*x + y, x) 992 >>> op.apply(z**7, lambda f: f.diff(z)) 993 y*z**7 + 7*z**7 + 42*z**5 994 """ 995 coeffs = self._poly.all_coeffs() 996 coeffs.reverse() 997 diffs = [obj] 998 for c in coeffs[1:]: 999 diffs.append(op(diffs[-1])) 1000 r = coeffs[0]*diffs[0] 1001 for c, d in zip(coeffs[1:], diffs[1:]): 1002 r += c*d 1003 return r 1004 1005 1006class MultOperator(Operator): 1007 """ Simply multiply by a "constant" """ 1008 1009 def __init__(self, p): 1010 self._poly = Poly(p, _x) 1011 1012 1013class ShiftA(Operator): 1014 """ Increment an upper index. """ 1015 1016 def __init__(self, ai): 1017 ai = sympify(ai) 1018 if ai == 0: 1019 raise ValueError('Cannot increment zero upper index.') 1020 self._poly = Poly(_x/ai + 1, _x) 1021 1022 def __str__(self): 1023 return '<Increment upper %s.>' % (1/self._poly.all_coeffs()[0]) 1024 1025 1026class ShiftB(Operator): 1027 """ Decrement a lower index. """ 1028 1029 def __init__(self, bi): 1030 bi = sympify(bi) 1031 if bi == 1: 1032 raise ValueError('Cannot decrement unit lower index.') 1033 self._poly = Poly(_x/(bi - 1) + 1, _x) 1034 1035 def __str__(self): 1036 return '<Decrement lower %s.>' % (1/self._poly.all_coeffs()[0] + 1) 1037 1038 1039class UnShiftA(Operator): 1040 """ Decrement an upper index. """ 1041 1042 def __init__(self, ap, bq, i, z): 1043 """ Note: i counts from zero! """ 1044 ap, bq, i = list(map(sympify, [ap, bq, i])) 1045 1046 self._ap = ap 1047 self._bq = bq 1048 self._i = i 1049 1050 ap = list(ap) 1051 bq = list(bq) 1052 ai = ap.pop(i) - 1 1053 1054 if ai == 0: 1055 raise ValueError('Cannot decrement unit upper index.') 1056 1057 m = Poly(z*ai, _x) 1058 for a in ap: 1059 m *= Poly(_x + a, _x) 1060 1061 A = Dummy('A') 1062 n = D = Poly(ai*A - ai, A) 1063 for b in bq: 1064 n *= D + (b - 1).as_poly(A) 1065 1066 b0 = -n.nth(0) 1067 if b0 == 0: 1068 raise ValueError('Cannot decrement upper index: ' 1069 'cancels with lower') 1070 1071 n = Poly(Poly(n.all_coeffs()[:-1], A).as_expr().subs(A, _x/ai + 1), _x) 1072 1073 self._poly = Poly((n - m)/b0, _x) 1074 1075 def __str__(self): 1076 return '<Decrement upper index #%s of %s, %s.>' % (self._i, 1077 self._ap, self._bq) 1078 1079 1080class UnShiftB(Operator): 1081 """ Increment a lower index. """ 1082 1083 def __init__(self, ap, bq, i, z): 1084 """ Note: i counts from zero! """ 1085 ap, bq, i = list(map(sympify, [ap, bq, i])) 1086 1087 self._ap = ap 1088 self._bq = bq 1089 self._i = i 1090 1091 ap = list(ap) 1092 bq = list(bq) 1093 bi = bq.pop(i) + 1 1094 1095 if bi == 0: 1096 raise ValueError('Cannot increment -1 lower index.') 1097 1098 m = Poly(_x*(bi - 1), _x) 1099 for b in bq: 1100 m *= Poly(_x + b - 1, _x) 1101 1102 B = Dummy('B') 1103 D = Poly((bi - 1)*B - bi + 1, B) 1104 n = Poly(z, B) 1105 for a in ap: 1106 n *= (D + a.as_poly(B)) 1107 1108 b0 = n.nth(0) 1109 if b0 == 0: 1110 raise ValueError('Cannot increment index: cancels with upper') 1111 1112 n = Poly(Poly(n.all_coeffs()[:-1], B).as_expr().subs( 1113 B, _x/(bi - 1) + 1), _x) 1114 1115 self._poly = Poly((m - n)/b0, _x) 1116 1117 def __str__(self): 1118 return '<Increment lower index #%s of %s, %s.>' % (self._i, 1119 self._ap, self._bq) 1120 1121 1122class MeijerShiftA(Operator): 1123 """ Increment an upper b index. """ 1124 1125 def __init__(self, bi): 1126 bi = sympify(bi) 1127 self._poly = Poly(bi - _x, _x) 1128 1129 def __str__(self): 1130 return '<Increment upper b=%s.>' % (self._poly.all_coeffs()[1]) 1131 1132 1133class MeijerShiftB(Operator): 1134 """ Decrement an upper a index. """ 1135 1136 def __init__(self, bi): 1137 bi = sympify(bi) 1138 self._poly = Poly(1 - bi + _x, _x) 1139 1140 def __str__(self): 1141 return '<Decrement upper a=%s.>' % (1 - self._poly.all_coeffs()[1]) 1142 1143 1144class MeijerShiftC(Operator): 1145 """ Increment a lower b index. """ 1146 1147 def __init__(self, bi): 1148 bi = sympify(bi) 1149 self._poly = Poly(-bi + _x, _x) 1150 1151 def __str__(self): 1152 return '<Increment lower b=%s.>' % (-self._poly.all_coeffs()[1]) 1153 1154 1155class MeijerShiftD(Operator): 1156 """ Decrement a lower a index. """ 1157 1158 def __init__(self, bi): 1159 bi = sympify(bi) 1160 self._poly = Poly(bi - 1 - _x, _x) 1161 1162 def __str__(self): 1163 return '<Decrement lower a=%s.>' % (self._poly.all_coeffs()[1] + 1) 1164 1165 1166class MeijerUnShiftA(Operator): 1167 """ Decrement an upper b index. """ 1168 1169 def __init__(self, an, ap, bm, bq, i, z): 1170 """ Note: i counts from zero! """ 1171 an, ap, bm, bq, i = list(map(sympify, [an, ap, bm, bq, i])) 1172 1173 self._an = an 1174 self._ap = ap 1175 self._bm = bm 1176 self._bq = bq 1177 self._i = i 1178 1179 an = list(an) 1180 ap = list(ap) 1181 bm = list(bm) 1182 bq = list(bq) 1183 bi = bm.pop(i) - 1 1184 1185 m = Poly(1, _x) 1186 for b in bm: 1187 m *= Poly(b - _x, _x) 1188 for b in bq: 1189 m *= Poly(_x - b, _x) 1190 1191 A = Dummy('A') 1192 D = Poly(bi - A, A) 1193 n = Poly(z, A) 1194 for a in an: 1195 n *= (D + 1 - a) 1196 for a in ap: 1197 n *= (-D + a - 1) 1198 1199 b0 = n.nth(0) 1200 if b0 == 0: 1201 raise ValueError('Cannot decrement upper b index (cancels)') 1202 1203 n = Poly(Poly(n.all_coeffs()[:-1], A).as_expr().subs(A, bi - _x), _x) 1204 1205 self._poly = Poly((m - n)/b0, _x) 1206 1207 def __str__(self): 1208 return '<Decrement upper b index #%s of %s, %s, %s, %s.>' % (self._i, 1209 self._an, self._ap, self._bm, self._bq) 1210 1211 1212class MeijerUnShiftB(Operator): 1213 """ Increment an upper a index. """ 1214 1215 def __init__(self, an, ap, bm, bq, i, z): 1216 """ Note: i counts from zero! """ 1217 an, ap, bm, bq, i = list(map(sympify, [an, ap, bm, bq, i])) 1218 1219 self._an = an 1220 self._ap = ap 1221 self._bm = bm 1222 self._bq = bq 1223 self._i = i 1224 1225 an = list(an) 1226 ap = list(ap) 1227 bm = list(bm) 1228 bq = list(bq) 1229 ai = an.pop(i) + 1 1230 1231 m = Poly(z, _x) 1232 for a in an: 1233 m *= Poly(1 - a + _x, _x) 1234 for a in ap: 1235 m *= Poly(a - 1 - _x, _x) 1236 1237 B = Dummy('B') 1238 D = Poly(B + ai - 1, B) 1239 n = Poly(1, B) 1240 for b in bm: 1241 n *= (-D + b) 1242 for b in bq: 1243 n *= (D - b) 1244 1245 b0 = n.nth(0) 1246 if b0 == 0: 1247 raise ValueError('Cannot increment upper a index (cancels)') 1248 1249 n = Poly(Poly(n.all_coeffs()[:-1], B).as_expr().subs( 1250 B, 1 - ai + _x), _x) 1251 1252 self._poly = Poly((m - n)/b0, _x) 1253 1254 def __str__(self): 1255 return '<Increment upper a index #%s of %s, %s, %s, %s.>' % (self._i, 1256 self._an, self._ap, self._bm, self._bq) 1257 1258 1259class MeijerUnShiftC(Operator): 1260 """ Decrement a lower b index. """ 1261 # XXX this is "essentially" the same as MeijerUnShiftA. This "essentially" 1262 # can be made rigorous using the functional equation G(1/z) = G'(z), 1263 # where G' denotes a G function of slightly altered parameters. 1264 # However, sorting out the details seems harder than just coding it 1265 # again. 1266 1267 def __init__(self, an, ap, bm, bq, i, z): 1268 """ Note: i counts from zero! """ 1269 an, ap, bm, bq, i = list(map(sympify, [an, ap, bm, bq, i])) 1270 1271 self._an = an 1272 self._ap = ap 1273 self._bm = bm 1274 self._bq = bq 1275 self._i = i 1276 1277 an = list(an) 1278 ap = list(ap) 1279 bm = list(bm) 1280 bq = list(bq) 1281 bi = bq.pop(i) - 1 1282 1283 m = Poly(1, _x) 1284 for b in bm: 1285 m *= Poly(b - _x, _x) 1286 for b in bq: 1287 m *= Poly(_x - b, _x) 1288 1289 C = Dummy('C') 1290 D = Poly(bi + C, C) 1291 n = Poly(z, C) 1292 for a in an: 1293 n *= (D + 1 - a) 1294 for a in ap: 1295 n *= (-D + a - 1) 1296 1297 b0 = n.nth(0) 1298 if b0 == 0: 1299 raise ValueError('Cannot decrement lower b index (cancels)') 1300 1301 n = Poly(Poly(n.all_coeffs()[:-1], C).as_expr().subs(C, _x - bi), _x) 1302 1303 self._poly = Poly((m - n)/b0, _x) 1304 1305 def __str__(self): 1306 return '<Decrement lower b index #%s of %s, %s, %s, %s.>' % (self._i, 1307 self._an, self._ap, self._bm, self._bq) 1308 1309 1310class MeijerUnShiftD(Operator): 1311 """ Increment a lower a index. """ 1312 # XXX This is essentially the same as MeijerUnShiftA. 1313 # See comment at MeijerUnShiftC. 1314 1315 def __init__(self, an, ap, bm, bq, i, z): 1316 """ Note: i counts from zero! """ 1317 an, ap, bm, bq, i = list(map(sympify, [an, ap, bm, bq, i])) 1318 1319 self._an = an 1320 self._ap = ap 1321 self._bm = bm 1322 self._bq = bq 1323 self._i = i 1324 1325 an = list(an) 1326 ap = list(ap) 1327 bm = list(bm) 1328 bq = list(bq) 1329 ai = ap.pop(i) + 1 1330 1331 m = Poly(z, _x) 1332 for a in an: 1333 m *= Poly(1 - a + _x, _x) 1334 for a in ap: 1335 m *= Poly(a - 1 - _x, _x) 1336 1337 B = Dummy('B') # - this is the shift operator `D_I` 1338 D = Poly(ai - 1 - B, B) 1339 n = Poly(1, B) 1340 for b in bm: 1341 n *= (-D + b) 1342 for b in bq: 1343 n *= (D - b) 1344 1345 b0 = n.nth(0) 1346 if b0 == 0: 1347 raise ValueError('Cannot increment lower a index (cancels)') 1348 1349 n = Poly(Poly(n.all_coeffs()[:-1], B).as_expr().subs( 1350 B, ai - 1 - _x), _x) 1351 1352 self._poly = Poly((m - n)/b0, _x) 1353 1354 def __str__(self): 1355 return '<Increment lower a index #%s of %s, %s, %s, %s.>' % (self._i, 1356 self._an, self._ap, self._bm, self._bq) 1357 1358 1359class ReduceOrder(Operator): 1360 """ Reduce Order by cancelling an upper and a lower index. """ 1361 1362 def __new__(cls, ai, bj): 1363 """ For convenience if reduction is not possible, return None. """ 1364 ai = sympify(ai) 1365 bj = sympify(bj) 1366 n = ai - bj 1367 if not n.is_Integer or n < 0: 1368 return None 1369 if bj.is_integer and bj.is_nonpositive: 1370 return None 1371 1372 expr = Operator.__new__(cls) 1373 1374 p = S.One 1375 for k in range(n): 1376 p *= (_x + bj + k)/(bj + k) 1377 1378 expr._poly = Poly(p, _x) 1379 expr._a = ai 1380 expr._b = bj 1381 1382 return expr 1383 1384 @classmethod 1385 def _meijer(cls, b, a, sign): 1386 """ Cancel b + sign*s and a + sign*s 1387 This is for meijer G functions. """ 1388 b = sympify(b) 1389 a = sympify(a) 1390 n = b - a 1391 if n.is_negative or not n.is_Integer: 1392 return None 1393 1394 expr = Operator.__new__(cls) 1395 1396 p = S.One 1397 for k in range(n): 1398 p *= (sign*_x + a + k) 1399 1400 expr._poly = Poly(p, _x) 1401 if sign == -1: 1402 expr._a = b 1403 expr._b = a 1404 else: 1405 expr._b = Add(1, a - 1, evaluate=False) 1406 expr._a = Add(1, b - 1, evaluate=False) 1407 1408 return expr 1409 1410 @classmethod 1411 def meijer_minus(cls, b, a): 1412 return cls._meijer(b, a, -1) 1413 1414 @classmethod 1415 def meijer_plus(cls, a, b): 1416 return cls._meijer(1 - a, 1 - b, 1) 1417 1418 def __str__(self): 1419 return '<Reduce order by cancelling upper %s with lower %s.>' % \ 1420 (self._a, self._b) 1421 1422 1423def _reduce_order(ap, bq, gen, key): 1424 """ Order reduction algorithm used in Hypergeometric and Meijer G """ 1425 ap = list(ap) 1426 bq = list(bq) 1427 1428 ap.sort(key=key) 1429 bq.sort(key=key) 1430 1431 nap = [] 1432 # we will edit bq in place 1433 operators = [] 1434 for a in ap: 1435 op = None 1436 for i in range(len(bq)): 1437 op = gen(a, bq[i]) 1438 if op is not None: 1439 bq.pop(i) 1440 break 1441 if op is None: 1442 nap.append(a) 1443 else: 1444 operators.append(op) 1445 1446 return nap, bq, operators 1447 1448 1449def reduce_order(func): 1450 """ 1451 Given the hypergeometric function ``func``, find a sequence of operators to 1452 reduces order as much as possible. 1453 1454 Explanation 1455 =========== 1456 1457 Return (newfunc, [operators]), where applying the operators to the 1458 hypergeometric function newfunc yields func. 1459 1460 Examples 1461 ======== 1462 1463 >>> from sympy.simplify.hyperexpand import reduce_order, Hyper_Function 1464 >>> reduce_order(Hyper_Function((1, 2), (3, 4))) 1465 (Hyper_Function((1, 2), (3, 4)), []) 1466 >>> reduce_order(Hyper_Function((1,), (1,))) 1467 (Hyper_Function((), ()), [<Reduce order by cancelling upper 1 with lower 1.>]) 1468 >>> reduce_order(Hyper_Function((2, 4), (3, 3))) 1469 (Hyper_Function((2,), (3,)), [<Reduce order by cancelling 1470 upper 4 with lower 3.>]) 1471 """ 1472 nap, nbq, operators = _reduce_order(func.ap, func.bq, ReduceOrder, default_sort_key) 1473 1474 return Hyper_Function(Tuple(*nap), Tuple(*nbq)), operators 1475 1476 1477def reduce_order_meijer(func): 1478 """ 1479 Given the Meijer G function parameters, ``func``, find a sequence of 1480 operators that reduces order as much as possible. 1481 1482 Return newfunc, [operators]. 1483 1484 Examples 1485 ======== 1486 1487 >>> from sympy.simplify.hyperexpand import (reduce_order_meijer, 1488 ... G_Function) 1489 >>> reduce_order_meijer(G_Function([3, 4], [5, 6], [3, 4], [1, 2]))[0] 1490 G_Function((4, 3), (5, 6), (3, 4), (2, 1)) 1491 >>> reduce_order_meijer(G_Function([3, 4], [5, 6], [3, 4], [1, 8]))[0] 1492 G_Function((3,), (5, 6), (3, 4), (1,)) 1493 >>> reduce_order_meijer(G_Function([3, 4], [5, 6], [7, 5], [1, 5]))[0] 1494 G_Function((3,), (), (), (1,)) 1495 >>> reduce_order_meijer(G_Function([3, 4], [5, 6], [7, 5], [5, 3]))[0] 1496 G_Function((), (), (), ()) 1497 """ 1498 1499 nan, nbq, ops1 = _reduce_order(func.an, func.bq, ReduceOrder.meijer_plus, 1500 lambda x: default_sort_key(-x)) 1501 nbm, nap, ops2 = _reduce_order(func.bm, func.ap, ReduceOrder.meijer_minus, 1502 default_sort_key) 1503 1504 return G_Function(nan, nap, nbm, nbq), ops1 + ops2 1505 1506 1507def make_derivative_operator(M, z): 1508 """ Create a derivative operator, to be passed to Operator.apply. """ 1509 def doit(C): 1510 r = z*C.diff(z) + C*M 1511 r = r.applyfunc(make_simp(z)) 1512 return r 1513 return doit 1514 1515 1516def apply_operators(obj, ops, op): 1517 """ 1518 Apply the list of operators ``ops`` to object ``obj``, substituting 1519 ``op`` for the generator. 1520 """ 1521 res = obj 1522 for o in reversed(ops): 1523 res = o.apply(res, op) 1524 return res 1525 1526 1527def devise_plan(target, origin, z): 1528 """ 1529 Devise a plan (consisting of shift and un-shift operators) to be applied 1530 to the hypergeometric function ``target`` to yield ``origin``. 1531 Returns a list of operators. 1532 1533 Examples 1534 ======== 1535 1536 >>> from sympy.simplify.hyperexpand import devise_plan, Hyper_Function 1537 >>> from sympy.abc import z 1538 1539 Nothing to do: 1540 1541 >>> devise_plan(Hyper_Function((1, 2), ()), Hyper_Function((1, 2), ()), z) 1542 [] 1543 >>> devise_plan(Hyper_Function((), (1, 2)), Hyper_Function((), (1, 2)), z) 1544 [] 1545 1546 Very simple plans: 1547 1548 >>> devise_plan(Hyper_Function((2,), ()), Hyper_Function((1,), ()), z) 1549 [<Increment upper 1.>] 1550 >>> devise_plan(Hyper_Function((), (2,)), Hyper_Function((), (1,)), z) 1551 [<Increment lower index #0 of [], [1].>] 1552 1553 Several buckets: 1554 1555 >>> from sympy import S 1556 >>> devise_plan(Hyper_Function((1, S.Half), ()), 1557 ... Hyper_Function((2, S('3/2')), ()), z) #doctest: +NORMALIZE_WHITESPACE 1558 [<Decrement upper index #0 of [3/2, 1], [].>, 1559 <Decrement upper index #0 of [2, 3/2], [].>] 1560 1561 A slightly more complicated plan: 1562 1563 >>> devise_plan(Hyper_Function((1, 3), ()), Hyper_Function((2, 2), ()), z) 1564 [<Increment upper 2.>, <Decrement upper index #0 of [2, 2], [].>] 1565 1566 Another more complicated plan: (note that the ap have to be shifted first!) 1567 1568 >>> devise_plan(Hyper_Function((1, -1), (2,)), Hyper_Function((3, -2), (4,)), z) 1569 [<Decrement lower 3.>, <Decrement lower 4.>, 1570 <Decrement upper index #1 of [-1, 2], [4].>, 1571 <Decrement upper index #1 of [-1, 3], [4].>, <Increment upper -2.>] 1572 """ 1573 abuckets, bbuckets, nabuckets, nbbuckets = [sift(params, _mod1) for 1574 params in (target.ap, target.bq, origin.ap, origin.bq)] 1575 1576 if len(list(abuckets.keys())) != len(list(nabuckets.keys())) or \ 1577 len(list(bbuckets.keys())) != len(list(nbbuckets.keys())): 1578 raise ValueError('%s not reachable from %s' % (target, origin)) 1579 1580 ops = [] 1581 1582 def do_shifts(fro, to, inc, dec): 1583 ops = [] 1584 for i in range(len(fro)): 1585 if to[i] - fro[i] > 0: 1586 sh = inc 1587 ch = 1 1588 else: 1589 sh = dec 1590 ch = -1 1591 1592 while to[i] != fro[i]: 1593 ops += [sh(fro, i)] 1594 fro[i] += ch 1595 1596 return ops 1597 1598 def do_shifts_a(nal, nbk, al, aother, bother): 1599 """ Shift us from (nal, nbk) to (al, nbk). """ 1600 return do_shifts(nal, al, lambda p, i: ShiftA(p[i]), 1601 lambda p, i: UnShiftA(p + aother, nbk + bother, i, z)) 1602 1603 def do_shifts_b(nal, nbk, bk, aother, bother): 1604 """ Shift us from (nal, nbk) to (nal, bk). """ 1605 return do_shifts(nbk, bk, 1606 lambda p, i: UnShiftB(nal + aother, p + bother, i, z), 1607 lambda p, i: ShiftB(p[i])) 1608 1609 for r in sorted(list(abuckets.keys()) + list(bbuckets.keys()), key=default_sort_key): 1610 al = () 1611 nal = () 1612 bk = () 1613 nbk = () 1614 if r in abuckets: 1615 al = abuckets[r] 1616 nal = nabuckets[r] 1617 if r in bbuckets: 1618 bk = bbuckets[r] 1619 nbk = nbbuckets[r] 1620 if len(al) != len(nal) or len(bk) != len(nbk): 1621 raise ValueError('%s not reachable from %s' % (target, origin)) 1622 1623 al, nal, bk, nbk = [sorted(list(w), key=default_sort_key) 1624 for w in [al, nal, bk, nbk]] 1625 1626 def others(dic, key): 1627 l = [] 1628 for k, value in dic.items(): 1629 if k != key: 1630 l += list(dic[k]) 1631 return l 1632 aother = others(nabuckets, r) 1633 bother = others(nbbuckets, r) 1634 1635 if len(al) == 0: 1636 # there can be no complications, just shift the bs as we please 1637 ops += do_shifts_b([], nbk, bk, aother, bother) 1638 elif len(bk) == 0: 1639 # there can be no complications, just shift the as as we please 1640 ops += do_shifts_a(nal, [], al, aother, bother) 1641 else: 1642 namax = nal[-1] 1643 amax = al[-1] 1644 1645 if nbk[0] - namax <= 0 or bk[0] - amax <= 0: 1646 raise ValueError('Non-suitable parameters.') 1647 1648 if namax - amax > 0: 1649 # we are going to shift down - first do the as, then the bs 1650 ops += do_shifts_a(nal, nbk, al, aother, bother) 1651 ops += do_shifts_b(al, nbk, bk, aother, bother) 1652 else: 1653 # we are going to shift up - first do the bs, then the as 1654 ops += do_shifts_b(nal, nbk, bk, aother, bother) 1655 ops += do_shifts_a(nal, bk, al, aother, bother) 1656 1657 nabuckets[r] = al 1658 nbbuckets[r] = bk 1659 1660 ops.reverse() 1661 return ops 1662 1663 1664def try_shifted_sum(func, z): 1665 """ Try to recognise a hypergeometric sum that starts from k > 0. """ 1666 abuckets, bbuckets = sift(func.ap, _mod1), sift(func.bq, _mod1) 1667 if len(abuckets[S.Zero]) != 1: 1668 return None 1669 r = abuckets[S.Zero][0] 1670 if r <= 0: 1671 return None 1672 if not S.Zero in bbuckets: 1673 return None 1674 l = list(bbuckets[S.Zero]) 1675 l.sort() 1676 k = l[0] 1677 if k <= 0: 1678 return None 1679 1680 nap = list(func.ap) 1681 nap.remove(r) 1682 nbq = list(func.bq) 1683 nbq.remove(k) 1684 k -= 1 1685 nap = [x - k for x in nap] 1686 nbq = [x - k for x in nbq] 1687 1688 ops = [] 1689 for n in range(r - 1): 1690 ops.append(ShiftA(n + 1)) 1691 ops.reverse() 1692 1693 fac = factorial(k)/z**k 1694 for a in nap: 1695 fac /= rf(a, k) 1696 for b in nbq: 1697 fac *= rf(b, k) 1698 1699 ops += [MultOperator(fac)] 1700 1701 p = 0 1702 for n in range(k): 1703 m = z**n/factorial(n) 1704 for a in nap: 1705 m *= rf(a, n) 1706 for b in nbq: 1707 m /= rf(b, n) 1708 p += m 1709 1710 return Hyper_Function(nap, nbq), ops, -p 1711 1712 1713def try_polynomial(func, z): 1714 """ Recognise polynomial cases. Returns None if not such a case. 1715 Requires order to be fully reduced. """ 1716 abuckets, bbuckets = sift(func.ap, _mod1), sift(func.bq, _mod1) 1717 a0 = abuckets[S.Zero] 1718 b0 = bbuckets[S.Zero] 1719 a0.sort() 1720 b0.sort() 1721 al0 = [x for x in a0 if x <= 0] 1722 bl0 = [x for x in b0 if x <= 0] 1723 1724 if bl0 and all(a < bl0[-1] for a in al0): 1725 return oo 1726 if not al0: 1727 return None 1728 1729 a = al0[-1] 1730 fac = 1 1731 res = S.One 1732 for n in Tuple(*list(range(-a))): 1733 fac *= z 1734 fac /= n + 1 1735 for a in func.ap: 1736 fac *= a + n 1737 for b in func.bq: 1738 fac /= b + n 1739 res += fac 1740 return res 1741 1742 1743def try_lerchphi(func): 1744 """ 1745 Try to find an expression for Hyper_Function ``func`` in terms of Lerch 1746 Transcendents. 1747 1748 Return None if no such expression can be found. 1749 """ 1750 # This is actually quite simple, and is described in Roach's paper, 1751 # section 18. 1752 # We don't need to implement the reduction to polylog here, this 1753 # is handled by expand_func. 1754 from sympy.matrices import Matrix, zeros 1755 from sympy.polys import apart 1756 1757 # First we need to figure out if the summation coefficient is a rational 1758 # function of the summation index, and construct that rational function. 1759 abuckets, bbuckets = sift(func.ap, _mod1), sift(func.bq, _mod1) 1760 1761 paired = {} 1762 for key, value in abuckets.items(): 1763 if key != 0 and not key in bbuckets: 1764 return None 1765 bvalue = bbuckets[key] 1766 paired[key] = (list(value), list(bvalue)) 1767 bbuckets.pop(key, None) 1768 if bbuckets != {}: 1769 return None 1770 if not S.Zero in abuckets: 1771 return None 1772 aints, bints = paired[S.Zero] 1773 # Account for the additional n! in denominator 1774 paired[S.Zero] = (aints, bints + [1]) 1775 1776 t = Dummy('t') 1777 numer = S.One 1778 denom = S.One 1779 for key, (avalue, bvalue) in paired.items(): 1780 if len(avalue) != len(bvalue): 1781 return None 1782 # Note that since order has been reduced fully, all the b are 1783 # bigger than all the a they differ from by an integer. In particular 1784 # if there are any negative b left, this function is not well-defined. 1785 for a, b in zip(avalue, bvalue): 1786 if (a - b).is_positive: 1787 k = a - b 1788 numer *= rf(b + t, k) 1789 denom *= rf(b, k) 1790 else: 1791 k = b - a 1792 numer *= rf(a, k) 1793 denom *= rf(a + t, k) 1794 1795 # Now do a partial fraction decomposition. 1796 # We assemble two structures: a list monomials of pairs (a, b) representing 1797 # a*t**b (b a non-negative integer), and a dict terms, where 1798 # terms[a] = [(b, c)] means that there is a term b/(t-a)**c. 1799 part = apart(numer/denom, t) 1800 args = Add.make_args(part) 1801 monomials = [] 1802 terms = {} 1803 for arg in args: 1804 numer, denom = arg.as_numer_denom() 1805 if not denom.has(t): 1806 p = Poly(numer, t) 1807 if not p.is_monomial: 1808 raise TypeError("p should be monomial") 1809 ((b, ), a) = p.LT() 1810 monomials += [(a/denom, b)] 1811 continue 1812 if numer.has(t): 1813 raise NotImplementedError('Need partial fraction decomposition' 1814 ' with linear denominators') 1815 indep, [dep] = denom.as_coeff_mul(t) 1816 n = 1 1817 if dep.is_Pow: 1818 n = dep.exp 1819 dep = dep.base 1820 if dep == t: 1821 a == 0 1822 elif dep.is_Add: 1823 a, tmp = dep.as_independent(t) 1824 b = 1 1825 if tmp != t: 1826 b, _ = tmp.as_independent(t) 1827 if dep != b*t + a: 1828 raise NotImplementedError('unrecognised form %s' % dep) 1829 a /= b 1830 indep *= b**n 1831 else: 1832 raise NotImplementedError('unrecognised form of partial fraction') 1833 terms.setdefault(a, []).append((numer/indep, n)) 1834 1835 # Now that we have this information, assemble our formula. All the 1836 # monomials yield rational functions and go into one basis element. 1837 # The terms[a] are related by differentiation. If the largest exponent is 1838 # n, we need lerchphi(z, k, a) for k = 1, 2, ..., n. 1839 # deriv maps a basis to its derivative, expressed as a C(z)-linear 1840 # combination of other basis elements. 1841 deriv = {} 1842 coeffs = {} 1843 z = Dummy('z') 1844 monomials.sort(key=lambda x: x[1]) 1845 mon = {0: 1/(1 - z)} 1846 if monomials: 1847 for k in range(monomials[-1][1]): 1848 mon[k + 1] = z*mon[k].diff(z) 1849 for a, n in monomials: 1850 coeffs.setdefault(S.One, []).append(a*mon[n]) 1851 for a, l in terms.items(): 1852 for c, k in l: 1853 coeffs.setdefault(lerchphi(z, k, a), []).append(c) 1854 l.sort(key=lambda x: x[1]) 1855 for k in range(2, l[-1][1] + 1): 1856 deriv[lerchphi(z, k, a)] = [(-a, lerchphi(z, k, a)), 1857 (1, lerchphi(z, k - 1, a))] 1858 deriv[lerchphi(z, 1, a)] = [(-a, lerchphi(z, 1, a)), 1859 (1/(1 - z), S.One)] 1860 trans = {} 1861 for n, b in enumerate([S.One] + list(deriv.keys())): 1862 trans[b] = n 1863 basis = [expand_func(b) for (b, _) in sorted(list(trans.items()), 1864 key=lambda x:x[1])] 1865 B = Matrix(basis) 1866 C = Matrix([[0]*len(B)]) 1867 for b, c in coeffs.items(): 1868 C[trans[b]] = Add(*c) 1869 M = zeros(len(B)) 1870 for b, l in deriv.items(): 1871 for c, b2 in l: 1872 M[trans[b], trans[b2]] = c 1873 return Formula(func, z, None, [], B, C, M) 1874 1875 1876def build_hypergeometric_formula(func): 1877 """ 1878 Create a formula object representing the hypergeometric function ``func``. 1879 1880 """ 1881 # We know that no `ap` are negative integers, otherwise "detect poly" 1882 # would have kicked in. However, `ap` could be empty. In this case we can 1883 # use a different basis. 1884 # I'm not aware of a basis that works in all cases. 1885 from sympy import zeros, Matrix, eye 1886 z = Dummy('z') 1887 if func.ap: 1888 afactors = [_x + a for a in func.ap] 1889 bfactors = [_x + b - 1 for b in func.bq] 1890 expr = _x*Mul(*bfactors) - z*Mul(*afactors) 1891 poly = Poly(expr, _x) 1892 n = poly.degree() 1893 basis = [] 1894 M = zeros(n) 1895 for k in range(n): 1896 a = func.ap[0] + k 1897 basis += [hyper([a] + list(func.ap[1:]), func.bq, z)] 1898 if k < n - 1: 1899 M[k, k] = -a 1900 M[k, k + 1] = a 1901 B = Matrix(basis) 1902 C = Matrix([[1] + [0]*(n - 1)]) 1903 derivs = [eye(n)] 1904 for k in range(n): 1905 derivs.append(M*derivs[k]) 1906 l = poly.all_coeffs() 1907 l.reverse() 1908 res = [0]*n 1909 for k, c in enumerate(l): 1910 for r, d in enumerate(C*derivs[k]): 1911 res[r] += c*d 1912 for k, c in enumerate(res): 1913 M[n - 1, k] = -c/derivs[n - 1][0, n - 1]/poly.all_coeffs()[0] 1914 return Formula(func, z, None, [], B, C, M) 1915 else: 1916 # Since there are no `ap`, none of the `bq` can be non-positive 1917 # integers. 1918 basis = [] 1919 bq = list(func.bq[:]) 1920 for i in range(len(bq)): 1921 basis += [hyper([], bq, z)] 1922 bq[i] += 1 1923 basis += [hyper([], bq, z)] 1924 B = Matrix(basis) 1925 n = len(B) 1926 C = Matrix([[1] + [0]*(n - 1)]) 1927 M = zeros(n) 1928 M[0, n - 1] = z/Mul(*func.bq) 1929 for k in range(1, n): 1930 M[k, k - 1] = func.bq[k - 1] 1931 M[k, k] = -func.bq[k - 1] 1932 return Formula(func, z, None, [], B, C, M) 1933 1934 1935def hyperexpand_special(ap, bq, z): 1936 """ 1937 Try to find a closed-form expression for hyper(ap, bq, z), where ``z`` 1938 is supposed to be a "special" value, e.g. 1. 1939 1940 This function tries various of the classical summation formulae 1941 (Gauss, Saalschuetz, etc). 1942 """ 1943 # This code is very ad-hoc. There are many clever algorithms 1944 # (notably Zeilberger's) related to this problem. 1945 # For now we just want a few simple cases to work. 1946 p, q = len(ap), len(bq) 1947 z_ = z 1948 z = unpolarify(z) 1949 if z == 0: 1950 return S.One 1951 if p == 2 and q == 1: 1952 # 2F1 1953 a, b, c = ap + bq 1954 if z == 1: 1955 # Gauss 1956 return gamma(c - a - b)*gamma(c)/gamma(c - a)/gamma(c - b) 1957 if z == -1 and simplify(b - a + c) == 1: 1958 b, a = a, b 1959 if z == -1 and simplify(a - b + c) == 1: 1960 # Kummer 1961 if b.is_integer and b.is_negative: 1962 return 2*cos(pi*b/2)*gamma(-b)*gamma(b - a + 1) \ 1963 /gamma(-b/2)/gamma(b/2 - a + 1) 1964 else: 1965 return gamma(b/2 + 1)*gamma(b - a + 1) \ 1966 /gamma(b + 1)/gamma(b/2 - a + 1) 1967 # TODO tons of more formulae 1968 # investigate what algorithms exist 1969 return hyper(ap, bq, z_) 1970 1971_collection = None 1972 1973 1974def _hyperexpand(func, z, ops0=[], z0=Dummy('z0'), premult=1, prem=0, 1975 rewrite='default'): 1976 """ 1977 Try to find an expression for the hypergeometric function ``func``. 1978 1979 Explanation 1980 =========== 1981 1982 The result is expressed in terms of a dummy variable ``z0``. Then it 1983 is multiplied by ``premult``. Then ``ops0`` is applied. 1984 ``premult`` must be a*z**prem for some a independent of ``z``. 1985 """ 1986 1987 if z.is_zero: 1988 return S.One 1989 1990 z = polarify(z, subs=False) 1991 if rewrite == 'default': 1992 rewrite = 'nonrepsmall' 1993 1994 def carryout_plan(f, ops): 1995 C = apply_operators(f.C.subs(f.z, z0), ops, 1996 make_derivative_operator(f.M.subs(f.z, z0), z0)) 1997 from sympy import eye 1998 C = apply_operators(C, ops0, 1999 make_derivative_operator(f.M.subs(f.z, z0) 2000 + prem*eye(f.M.shape[0]), z0)) 2001 2002 if premult == 1: 2003 C = C.applyfunc(make_simp(z0)) 2004 r = reduce(lambda s,m: s+m[0]*m[1], zip(C, f.B.subs(f.z, z0)), S.Zero)*premult 2005 res = r.subs(z0, z) 2006 if rewrite: 2007 res = res.rewrite(rewrite) 2008 return res 2009 2010 # TODO 2011 # The following would be possible: 2012 # *) PFD Duplication (see Kelly Roach's paper) 2013 # *) In a similar spirit, try_lerchphi() can be generalised considerably. 2014 2015 global _collection 2016 if _collection is None: 2017 _collection = FormulaCollection() 2018 2019 debug('Trying to expand hypergeometric function ', func) 2020 2021 # First reduce order as much as possible. 2022 func, ops = reduce_order(func) 2023 if ops: 2024 debug(' Reduced order to ', func) 2025 else: 2026 debug(' Could not reduce order.') 2027 2028 # Now try polynomial cases 2029 res = try_polynomial(func, z0) 2030 if res is not None: 2031 debug(' Recognised polynomial.') 2032 p = apply_operators(res, ops, lambda f: z0*f.diff(z0)) 2033 p = apply_operators(p*premult, ops0, lambda f: z0*f.diff(z0)) 2034 return unpolarify(simplify(p).subs(z0, z)) 2035 2036 # Try to recognise a shifted sum. 2037 p = S.Zero 2038 res = try_shifted_sum(func, z0) 2039 if res is not None: 2040 func, nops, p = res 2041 debug(' Recognised shifted sum, reduced order to ', func) 2042 ops += nops 2043 2044 # apply the plan for poly 2045 p = apply_operators(p, ops, lambda f: z0*f.diff(z0)) 2046 p = apply_operators(p*premult, ops0, lambda f: z0*f.diff(z0)) 2047 p = simplify(p).subs(z0, z) 2048 2049 # Try special expansions early. 2050 if unpolarify(z) in [1, -1] and (len(func.ap), len(func.bq)) == (2, 1): 2051 f = build_hypergeometric_formula(func) 2052 r = carryout_plan(f, ops).replace(hyper, hyperexpand_special) 2053 if not r.has(hyper): 2054 return r + p 2055 2056 # Try to find a formula in our collection 2057 formula = _collection.lookup_origin(func) 2058 2059 # Now try a lerch phi formula 2060 if formula is None: 2061 formula = try_lerchphi(func) 2062 2063 if formula is None: 2064 debug(' Could not find an origin. ', 2065 'Will return answer in terms of ' 2066 'simpler hypergeometric functions.') 2067 formula = build_hypergeometric_formula(func) 2068 2069 debug(' Found an origin: ', formula.closed_form, ' ', formula.func) 2070 2071 # We need to find the operators that convert formula into func. 2072 ops += devise_plan(func, formula.func, z0) 2073 2074 # Now carry out the plan. 2075 r = carryout_plan(formula, ops) + p 2076 2077 return powdenest(r, polar=True).replace(hyper, hyperexpand_special) 2078 2079 2080def devise_plan_meijer(fro, to, z): 2081 """ 2082 Find operators to convert G-function ``fro`` into G-function ``to``. 2083 2084 Explanation 2085 =========== 2086 2087 It is assumed that ``fro`` and ``to`` have the same signatures, and that in fact 2088 any corresponding pair of parameters differs by integers, and a direct path 2089 is possible. I.e. if there are parameters a1 b1 c1 and a2 b2 c2 it is 2090 assumed that a1 can be shifted to a2, etc. The only thing this routine 2091 determines is the order of shifts to apply, nothing clever will be tried. 2092 It is also assumed that ``fro`` is suitable. 2093 2094 Examples 2095 ======== 2096 2097 >>> from sympy.simplify.hyperexpand import (devise_plan_meijer, 2098 ... G_Function) 2099 >>> from sympy.abc import z 2100 2101 Empty plan: 2102 2103 >>> devise_plan_meijer(G_Function([1], [2], [3], [4]), 2104 ... G_Function([1], [2], [3], [4]), z) 2105 [] 2106 2107 Very simple plans: 2108 2109 >>> devise_plan_meijer(G_Function([0], [], [], []), 2110 ... G_Function([1], [], [], []), z) 2111 [<Increment upper a index #0 of [0], [], [], [].>] 2112 >>> devise_plan_meijer(G_Function([0], [], [], []), 2113 ... G_Function([-1], [], [], []), z) 2114 [<Decrement upper a=0.>] 2115 >>> devise_plan_meijer(G_Function([], [1], [], []), 2116 ... G_Function([], [2], [], []), z) 2117 [<Increment lower a index #0 of [], [1], [], [].>] 2118 2119 Slightly more complicated plans: 2120 2121 >>> devise_plan_meijer(G_Function([0], [], [], []), 2122 ... G_Function([2], [], [], []), z) 2123 [<Increment upper a index #0 of [1], [], [], [].>, 2124 <Increment upper a index #0 of [0], [], [], [].>] 2125 >>> devise_plan_meijer(G_Function([0], [], [0], []), 2126 ... G_Function([-1], [], [1], []), z) 2127 [<Increment upper b=0.>, <Decrement upper a=0.>] 2128 2129 Order matters: 2130 2131 >>> devise_plan_meijer(G_Function([0], [], [0], []), 2132 ... G_Function([1], [], [1], []), z) 2133 [<Increment upper a index #0 of [0], [], [1], [].>, <Increment upper b=0.>] 2134 """ 2135 # TODO for now, we use the following simple heuristic: inverse-shift 2136 # when possible, shift otherwise. Give up if we cannot make progress. 2137 2138 def try_shift(f, t, shifter, diff, counter): 2139 """ Try to apply ``shifter`` in order to bring some element in ``f`` 2140 nearer to its counterpart in ``to``. ``diff`` is +/- 1 and 2141 determines the effect of ``shifter``. Counter is a list of elements 2142 blocking the shift. 2143 2144 Return an operator if change was possible, else None. 2145 """ 2146 for idx, (a, b) in enumerate(zip(f, t)): 2147 if ( 2148 (a - b).is_integer and (b - a)/diff > 0 and 2149 all(a != x for x in counter)): 2150 sh = shifter(idx) 2151 f[idx] += diff 2152 return sh 2153 fan = list(fro.an) 2154 fap = list(fro.ap) 2155 fbm = list(fro.bm) 2156 fbq = list(fro.bq) 2157 ops = [] 2158 change = True 2159 while change: 2160 change = False 2161 op = try_shift(fan, to.an, 2162 lambda i: MeijerUnShiftB(fan, fap, fbm, fbq, i, z), 2163 1, fbm + fbq) 2164 if op is not None: 2165 ops += [op] 2166 change = True 2167 continue 2168 op = try_shift(fap, to.ap, 2169 lambda i: MeijerUnShiftD(fan, fap, fbm, fbq, i, z), 2170 1, fbm + fbq) 2171 if op is not None: 2172 ops += [op] 2173 change = True 2174 continue 2175 op = try_shift(fbm, to.bm, 2176 lambda i: MeijerUnShiftA(fan, fap, fbm, fbq, i, z), 2177 -1, fan + fap) 2178 if op is not None: 2179 ops += [op] 2180 change = True 2181 continue 2182 op = try_shift(fbq, to.bq, 2183 lambda i: MeijerUnShiftC(fan, fap, fbm, fbq, i, z), 2184 -1, fan + fap) 2185 if op is not None: 2186 ops += [op] 2187 change = True 2188 continue 2189 op = try_shift(fan, to.an, lambda i: MeijerShiftB(fan[i]), -1, []) 2190 if op is not None: 2191 ops += [op] 2192 change = True 2193 continue 2194 op = try_shift(fap, to.ap, lambda i: MeijerShiftD(fap[i]), -1, []) 2195 if op is not None: 2196 ops += [op] 2197 change = True 2198 continue 2199 op = try_shift(fbm, to.bm, lambda i: MeijerShiftA(fbm[i]), 1, []) 2200 if op is not None: 2201 ops += [op] 2202 change = True 2203 continue 2204 op = try_shift(fbq, to.bq, lambda i: MeijerShiftC(fbq[i]), 1, []) 2205 if op is not None: 2206 ops += [op] 2207 change = True 2208 continue 2209 if fan != list(to.an) or fap != list(to.ap) or fbm != list(to.bm) or \ 2210 fbq != list(to.bq): 2211 raise NotImplementedError('Could not devise plan.') 2212 ops.reverse() 2213 return ops 2214 2215_meijercollection = None 2216 2217 2218def _meijergexpand(func, z0, allow_hyper=False, rewrite='default', 2219 place=None): 2220 """ 2221 Try to find an expression for the Meijer G function specified 2222 by the G_Function ``func``. If ``allow_hyper`` is True, then returning 2223 an expression in terms of hypergeometric functions is allowed. 2224 2225 Currently this just does Slater's theorem. 2226 If expansions exist both at zero and at infinity, ``place`` 2227 can be set to ``0`` or ``zoo`` for the preferred choice. 2228 """ 2229 global _meijercollection 2230 if _meijercollection is None: 2231 _meijercollection = MeijerFormulaCollection() 2232 if rewrite == 'default': 2233 rewrite = None 2234 2235 func0 = func 2236 debug('Try to expand Meijer G function corresponding to ', func) 2237 2238 # We will play games with analytic continuation - rather use a fresh symbol 2239 z = Dummy('z') 2240 2241 func, ops = reduce_order_meijer(func) 2242 if ops: 2243 debug(' Reduced order to ', func) 2244 else: 2245 debug(' Could not reduce order.') 2246 2247 # Try to find a direct formula 2248 f = _meijercollection.lookup_origin(func) 2249 if f is not None: 2250 debug(' Found a Meijer G formula: ', f.func) 2251 ops += devise_plan_meijer(f.func, func, z) 2252 2253 # Now carry out the plan. 2254 C = apply_operators(f.C.subs(f.z, z), ops, 2255 make_derivative_operator(f.M.subs(f.z, z), z)) 2256 2257 C = C.applyfunc(make_simp(z)) 2258 r = C*f.B.subs(f.z, z) 2259 r = r[0].subs(z, z0) 2260 return powdenest(r, polar=True) 2261 2262 debug(" Could not find a direct formula. Trying Slater's theorem.") 2263 2264 # TODO the following would be possible: 2265 # *) Paired Index Theorems 2266 # *) PFD Duplication 2267 # (See Kelly Roach's paper for details on either.) 2268 # 2269 # TODO Also, we tend to create combinations of gamma functions that can be 2270 # simplified. 2271 2272 def can_do(pbm, pap): 2273 """ Test if slater applies. """ 2274 for i in pbm: 2275 if len(pbm[i]) > 1: 2276 l = 0 2277 if i in pap: 2278 l = len(pap[i]) 2279 if l + 1 < len(pbm[i]): 2280 return False 2281 return True 2282 2283 def do_slater(an, bm, ap, bq, z, zfinal): 2284 # zfinal is the value that will eventually be substituted for z. 2285 # We pass it to _hyperexpand to improve performance. 2286 func = G_Function(an, bm, ap, bq) 2287 _, pbm, pap, _ = func.compute_buckets() 2288 if not can_do(pbm, pap): 2289 return S.Zero, False 2290 2291 cond = len(an) + len(ap) < len(bm) + len(bq) 2292 if len(an) + len(ap) == len(bm) + len(bq): 2293 cond = abs(z) < 1 2294 if cond is False: 2295 return S.Zero, False 2296 2297 res = S.Zero 2298 for m in pbm: 2299 if len(pbm[m]) == 1: 2300 bh = pbm[m][0] 2301 fac = 1 2302 bo = list(bm) 2303 bo.remove(bh) 2304 for bj in bo: 2305 fac *= gamma(bj - bh) 2306 for aj in an: 2307 fac *= gamma(1 + bh - aj) 2308 for bj in bq: 2309 fac /= gamma(1 + bh - bj) 2310 for aj in ap: 2311 fac /= gamma(aj - bh) 2312 nap = [1 + bh - a for a in list(an) + list(ap)] 2313 nbq = [1 + bh - b for b in list(bo) + list(bq)] 2314 2315 k = polar_lift(S.NegativeOne**(len(ap) - len(bm))) 2316 harg = k*zfinal 2317 # NOTE even though k "is" +-1, this has to be t/k instead of 2318 # t*k ... we are using polar numbers for consistency! 2319 premult = (t/k)**bh 2320 hyp = _hyperexpand(Hyper_Function(nap, nbq), harg, ops, 2321 t, premult, bh, rewrite=None) 2322 res += fac * hyp 2323 else: 2324 b_ = pbm[m][0] 2325 ki = [bi - b_ for bi in pbm[m][1:]] 2326 u = len(ki) 2327 li = [ai - b_ for ai in pap[m][:u + 1]] 2328 bo = list(bm) 2329 for b in pbm[m]: 2330 bo.remove(b) 2331 ao = list(ap) 2332 for a in pap[m][:u]: 2333 ao.remove(a) 2334 lu = li[-1] 2335 di = [l - k for (l, k) in zip(li, ki)] 2336 2337 # We first work out the integrand: 2338 s = Dummy('s') 2339 integrand = z**s 2340 for b in bm: 2341 if not Mod(b, 1) and b.is_Number: 2342 b = int(round(b)) 2343 integrand *= gamma(b - s) 2344 for a in an: 2345 integrand *= gamma(1 - a + s) 2346 for b in bq: 2347 integrand /= gamma(1 - b + s) 2348 for a in ap: 2349 integrand /= gamma(a - s) 2350 2351 # Now sum the finitely many residues: 2352 # XXX This speeds up some cases - is it a good idea? 2353 integrand = expand_func(integrand) 2354 for r in range(int(round(lu))): 2355 resid = residue(integrand, s, b_ + r) 2356 resid = apply_operators(resid, ops, lambda f: z*f.diff(z)) 2357 res -= resid 2358 2359 # Now the hypergeometric term. 2360 au = b_ + lu 2361 k = polar_lift(S.NegativeOne**(len(ao) + len(bo) + 1)) 2362 harg = k*zfinal 2363 premult = (t/k)**au 2364 nap = [1 + au - a for a in list(an) + list(ap)] + [1] 2365 nbq = [1 + au - b for b in list(bm) + list(bq)] 2366 2367 hyp = _hyperexpand(Hyper_Function(nap, nbq), harg, ops, 2368 t, premult, au, rewrite=None) 2369 2370 C = S.NegativeOne**(lu)/factorial(lu) 2371 for i in range(u): 2372 C *= S.NegativeOne**di[i]/rf(lu - li[i] + 1, di[i]) 2373 for a in an: 2374 C *= gamma(1 - a + au) 2375 for b in bo: 2376 C *= gamma(b - au) 2377 for a in ao: 2378 C /= gamma(a - au) 2379 for b in bq: 2380 C /= gamma(1 - b + au) 2381 2382 res += C*hyp 2383 2384 return res, cond 2385 2386 t = Dummy('t') 2387 slater1, cond1 = do_slater(func.an, func.bm, func.ap, func.bq, z, z0) 2388 2389 def tr(l): 2390 return [1 - x for x in l] 2391 2392 for op in ops: 2393 op._poly = Poly(op._poly.subs({z: 1/t, _x: -_x}), _x) 2394 slater2, cond2 = do_slater(tr(func.bm), tr(func.an), tr(func.bq), tr(func.ap), 2395 t, 1/z0) 2396 2397 slater1 = powdenest(slater1.subs(z, z0), polar=True) 2398 slater2 = powdenest(slater2.subs(t, 1/z0), polar=True) 2399 if not isinstance(cond2, bool): 2400 cond2 = cond2.subs(t, 1/z) 2401 2402 m = func(z) 2403 if m.delta > 0 or \ 2404 (m.delta == 0 and len(m.ap) == len(m.bq) and 2405 (re(m.nu) < -1) is not False and polar_lift(z0) == polar_lift(1)): 2406 # The condition delta > 0 means that the convergence region is 2407 # connected. Any expression we find can be continued analytically 2408 # to the entire convergence region. 2409 # The conditions delta==0, p==q, re(nu) < -1 imply that G is continuous 2410 # on the positive reals, so the values at z=1 agree. 2411 if cond1 is not False: 2412 cond1 = True 2413 if cond2 is not False: 2414 cond2 = True 2415 2416 if cond1 is True: 2417 slater1 = slater1.rewrite(rewrite or 'nonrep') 2418 else: 2419 slater1 = slater1.rewrite(rewrite or 'nonrepsmall') 2420 if cond2 is True: 2421 slater2 = slater2.rewrite(rewrite or 'nonrep') 2422 else: 2423 slater2 = slater2.rewrite(rewrite or 'nonrepsmall') 2424 2425 if cond1 is not False and cond2 is not False: 2426 # If one condition is False, there is no choice. 2427 if place == 0: 2428 cond2 = False 2429 if place == zoo: 2430 cond1 = False 2431 2432 if not isinstance(cond1, bool): 2433 cond1 = cond1.subs(z, z0) 2434 if not isinstance(cond2, bool): 2435 cond2 = cond2.subs(z, z0) 2436 2437 def weight(expr, cond): 2438 if cond is True: 2439 c0 = 0 2440 elif cond is False: 2441 c0 = 1 2442 else: 2443 c0 = 2 2444 if expr.has(oo, zoo, -oo, nan): 2445 # XXX this actually should not happen, but consider 2446 # S('meijerg(((0, -1/2, 0, -1/2, 1/2), ()), ((0,), 2447 # (-1/2, -1/2, -1/2, -1)), exp_polar(I*pi))/4') 2448 c0 = 3 2449 return (c0, expr.count(hyper), expr.count_ops()) 2450 2451 w1 = weight(slater1, cond1) 2452 w2 = weight(slater2, cond2) 2453 if min(w1, w2) <= (0, 1, oo): 2454 if w1 < w2: 2455 return slater1 2456 else: 2457 return slater2 2458 if max(w1[0], w2[0]) <= 1 and max(w1[1], w2[1]) <= 1: 2459 return Piecewise((slater1, cond1), (slater2, cond2), (func0(z0), True)) 2460 2461 # We couldn't find an expression without hypergeometric functions. 2462 # TODO it would be helpful to give conditions under which the integral 2463 # is known to diverge. 2464 r = Piecewise((slater1, cond1), (slater2, cond2), (func0(z0), True)) 2465 if r.has(hyper) and not allow_hyper: 2466 debug(' Could express using hypergeometric functions, ' 2467 'but not allowed.') 2468 if not r.has(hyper) or allow_hyper: 2469 return r 2470 2471 return func0(z0) 2472 2473 2474def hyperexpand(f, allow_hyper=False, rewrite='default', place=None): 2475 """ 2476 Expand hypergeometric functions. If allow_hyper is True, allow partial 2477 simplification (that is a result different from input, 2478 but still containing hypergeometric functions). 2479 2480 If a G-function has expansions both at zero and at infinity, 2481 ``place`` can be set to ``0`` or ``zoo`` to indicate the 2482 preferred choice. 2483 2484 Examples 2485 ======== 2486 2487 >>> from sympy.simplify.hyperexpand import hyperexpand 2488 >>> from sympy.functions import hyper 2489 >>> from sympy.abc import z 2490 >>> hyperexpand(hyper([], [], z)) 2491 exp(z) 2492 2493 Non-hyperegeometric parts of the expression and hypergeometric expressions 2494 that are not recognised are left unchanged: 2495 2496 >>> hyperexpand(1 + hyper([1, 1, 1], [], z)) 2497 hyper((1, 1, 1), (), z) + 1 2498 """ 2499 f = sympify(f) 2500 2501 def do_replace(ap, bq, z): 2502 r = _hyperexpand(Hyper_Function(ap, bq), z, rewrite=rewrite) 2503 if r is None: 2504 return hyper(ap, bq, z) 2505 else: 2506 return r 2507 2508 def do_meijer(ap, bq, z): 2509 r = _meijergexpand(G_Function(ap[0], ap[1], bq[0], bq[1]), z, 2510 allow_hyper, rewrite=rewrite, place=place) 2511 if not r.has(nan, zoo, oo, -oo): 2512 return r 2513 return f.replace(hyper, do_replace).replace(meijerg, do_meijer) 2514