Coverage for /usr/lib/python3/dist-packages/sympy/functions/special/hyper.py: 36%

1"""Hypergeometric and Meijer G-functions""" 

2from functools import reduce 

3 

4from sympy.core import S, ilcm, Mod 

5from sympy.core.add import Add 

6from sympy.core.expr import Expr 

7from sympy.core.function import Function, Derivative, ArgumentIndexError 

8 

9from sympy.core.containers import Tuple 

10from sympy.core.mul import Mul 

11from sympy.core.numbers import I, pi, oo, zoo 

12from sympy.core.relational import Ne 

13from sympy.core.sorting import default_sort_key 

14from sympy.core.symbol import Dummy 

15 

16from sympy.functions import (sqrt, exp, log, sin, cos, asin, atan, 

17 sinh, cosh, asinh, acosh, atanh, acoth) 

18from sympy.functions import factorial, RisingFactorial 

19from sympy.functions.elementary.complexes import Abs, re, unpolarify 

20from sympy.functions.elementary.exponential import exp_polar 

21from sympy.functions.elementary.integers import ceiling 

22from sympy.functions.elementary.piecewise import Piecewise 

23from sympy.logic.boolalg import (And, Or) 

24 

25class TupleArg(Tuple): 

26 def limit(self, x, xlim, dir='+'): 

27 """ Compute limit x->xlim. 

28 """ 

29 from sympy.series.limits import limit 

30 return TupleArg(*[limit(f, x, xlim, dir) for f in self.args]) 

31 

32 

33# TODO should __new__ accept **options? 

34# TODO should constructors should check if parameters are sensible? 

35 

36 

37def _prep_tuple(v): 

38 """ 

39 Turn an iterable argument *v* into a tuple and unpolarify, since both 

40 hypergeometric and meijer g-functions are unbranched in their parameters. 

41 

42 Examples 

43 ======== 

44 

45 >>> from sympy.functions.special.hyper import _prep_tuple 

46 >>> _prep_tuple([1, 2, 3]) 

47 (1, 2, 3) 

48 >>> _prep_tuple((4, 5)) 

49 (4, 5) 

50 >>> _prep_tuple((7, 8, 9)) 

51 (7, 8, 9) 

52 

53 """ 

54 return TupleArg(*[unpolarify(x) for x in v]) 

55 

56 

57class TupleParametersBase(Function): 

58 """ Base class that takes care of differentiation, when some of 

59 the arguments are actually tuples. """ 

60 # This is not deduced automatically since there are Tuples as arguments. 

61 is_commutative = True 

62 

63 def _eval_derivative(self, s): 

64 try: 

65 res = 0 

66 if self.args[0].has(s) or self.args[1].has(s): 

67 for i, p in enumerate(self._diffargs): 

68 m = self._diffargs[i].diff(s) 

69 if m != 0: 

70 res += self.fdiff((1, i))*m 

71 return res + self.fdiff(3)*self.args[2].diff(s) 

72 except (ArgumentIndexError, NotImplementedError): 

73 return Derivative(self, s) 

74 

75 

76class hyper(TupleParametersBase): 

77 r""" 

78 The generalized hypergeometric function is defined by a series where 

79 the ratios of successive terms are a rational function of the summation 

80 index. When convergent, it is continued analytically to the largest 

81 possible domain. 

82 

83 Explanation 

84 =========== 

85 

86 The hypergeometric function depends on two vectors of parameters, called 

87 the numerator parameters $a_p$, and the denominator parameters 

88 $b_q$. It also has an argument $z$. The series definition is 

89 

90 .. math :: 

91 {}_pF_q\left(\begin{matrix} a_1, \cdots, a_p \\ b_1, \cdots, b_q \end{matrix} 

92 \middle| z \right) 

93 = \sum_{n=0}^\infty \frac{(a_1)_n \cdots (a_p)_n}{(b_1)_n \cdots (b_q)_n} 

94 \frac{z^n}{n!}, 

95 

96 where $(a)_n = (a)(a+1)\cdots(a+n-1)$ denotes the rising factorial. 

97 

98 If one of the $b_q$ is a non-positive integer then the series is 

99 undefined unless one of the $a_p$ is a larger (i.e., smaller in 

100 magnitude) non-positive integer. If none of the $b_q$ is a 

101 non-positive integer and one of the $a_p$ is a non-positive 

102 integer, then the series reduces to a polynomial. To simplify the 

103 following discussion, we assume that none of the $a_p$ or 

104 $b_q$ is a non-positive integer. For more details, see the 

105 references. 

106 

107 The series converges for all $z$ if $p \le q$, and thus 

108 defines an entire single-valued function in this case. If $p = 

109 q+1$ the series converges for $|z| < 1$, and can be continued 

110 analytically into a half-plane. If $p > q+1$ the series is 

111 divergent for all $z$. 

112 

113 Please note the hypergeometric function constructor currently does *not* 

114 check if the parameters actually yield a well-defined function. 

115 

116 Examples 

117 ======== 

118 

119 The parameters $a_p$ and $b_q$ can be passed as arbitrary 

120 iterables, for example: 

121 

122 >>> from sympy import hyper 

123 >>> from sympy.abc import x, n, a 

124 >>> hyper((1, 2, 3), [3, 4], x) 

125 hyper((1, 2, 3), (3, 4), x) 

126 

127 There is also pretty printing (it looks better using Unicode): 

128 

129 >>> from sympy import pprint 

130 >>> pprint(hyper((1, 2, 3), [3, 4], x), use_unicode=False) 

131 _ 

132 |_ /1, 2, 3 | \ 

133 | | | x| 

134 3 2 \ 3, 4 | / 

135 

136 The parameters must always be iterables, even if they are vectors of 

137 length one or zero: 

138 

139 >>> hyper((1, ), [], x) 

140 hyper((1,), (), x) 

141 

142 But of course they may be variables (but if they depend on $x$ then you 

143 should not expect much implemented functionality): 

144 

145 >>> hyper((n, a), (n**2,), x) 

146 hyper((n, a), (n**2,), x) 

147 

148 The hypergeometric function generalizes many named special functions. 

149 The function ``hyperexpand()`` tries to express a hypergeometric function 

150 using named special functions. For example: 

151 

152 >>> from sympy import hyperexpand 

153 >>> hyperexpand(hyper([], [], x)) 

154 exp(x) 

155 

156 You can also use ``expand_func()``: 

157 

158 >>> from sympy import expand_func 

159 >>> expand_func(x*hyper([1, 1], [2], -x)) 

160 log(x + 1) 

161 

162 More examples: 

163 

164 >>> from sympy import S 

165 >>> hyperexpand(hyper([], [S(1)/2], -x**2/4)) 

166 cos(x) 

167 >>> hyperexpand(x*hyper([S(1)/2, S(1)/2], [S(3)/2], x**2)) 

168 asin(x) 

169 

170 We can also sometimes ``hyperexpand()`` parametric functions: 

171 

172 >>> from sympy.abc import a 

173 >>> hyperexpand(hyper([-a], [], x)) 

174 (1 - x)**a 

175 

176 See Also 

177 ======== 

178 

179 sympy.simplify.hyperexpand 

180 gamma 

181 meijerg 

182 

183 References 

184 ========== 

185 

186 .. [1] Luke, Y. L. (1969), The Special Functions and Their Approximations, 

187 Volume 1 

188 .. [2] https://en.wikipedia.org/wiki/Generalized_hypergeometric_function 

189 

190 """ 

191 

192 

193 def __new__(cls, ap, bq, z, **kwargs): 

194 # TODO should we check convergence conditions? 

195 return Function.__new__(cls, _prep_tuple(ap), _prep_tuple(bq), z, **kwargs) 

196 

197 @classmethod 

198 def eval(cls, ap, bq, z): 

199 if len(ap) <= len(bq) or (len(ap) == len(bq) + 1 and (Abs(z) <= 1) == True): 

200 nz = unpolarify(z) 

201 if z != nz: 

202 return hyper(ap, bq, nz) 

203 

204 def fdiff(self, argindex=3): 

205 if argindex != 3: 

206 raise ArgumentIndexError(self, argindex) 

207 nap = Tuple(*[a + 1 for a in self.ap]) 

208 nbq = Tuple(*[b + 1 for b in self.bq]) 

209 fac = Mul(*self.ap)/Mul(*self.bq) 

210 return fac*hyper(nap, nbq, self.argument) 

211 

212 def _eval_expand_func(self, **hints): 

213 from sympy.functions.special.gamma_functions import gamma 

214 from sympy.simplify.hyperexpand import hyperexpand 

215 if len(self.ap) == 2 and len(self.bq) == 1 and self.argument == 1: 

216 a, b = self.ap 

217 c = self.bq[0] 

218 return gamma(c)*gamma(c - a - b)/gamma(c - a)/gamma(c - b) 

219 return hyperexpand(self) 

220 

221 def _eval_rewrite_as_Sum(self, ap, bq, z, **kwargs): 

222 from sympy.concrete.summations import Sum 

223 n = Dummy("n", integer=True) 

224 rfap = [RisingFactorial(a, n) for a in ap] 

225 rfbq = [RisingFactorial(b, n) for b in bq] 

226 coeff = Mul(*rfap) / Mul(*rfbq) 

227 return Piecewise((Sum(coeff * z**n / factorial(n), (n, 0, oo)), 

228 self.convergence_statement), (self, True)) 

229 

230 def _eval_as_leading_term(self, x, logx=None, cdir=0): 

231 arg = self.args[2] 

232 x0 = arg.subs(x, 0) 

233 if x0 is S.NaN: 

234 x0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+') 

235 

236 if x0 is S.Zero: 

237 return S.One 

238 return super()._eval_as_leading_term(x, logx=logx, cdir=cdir) 

239 

240 def _eval_nseries(self, x, n, logx, cdir=0): 

241 

242 from sympy.series.order import Order 

243 

244 arg = self.args[2] 

245 x0 = arg.limit(x, 0) 

246 ap = self.args[0] 

247 bq = self.args[1] 

248 

249 if x0 != 0: 

250 return super()._eval_nseries(x, n, logx) 

251 

252 terms = [] 

253 

254 for i in range(n): 

255 num = Mul(*[RisingFactorial(a, i) for a in ap]) 

256 den = Mul(*[RisingFactorial(b, i) for b in bq]) 

257 terms.append(((num/den) * (arg**i)) / factorial(i)) 

258 

259 return (Add(*terms) + Order(x**n,x)) 

260 

261 @property 

262 def argument(self): 

263 """ Argument of the hypergeometric function. """ 

264 return self.args[2] 

265 

266 @property 

267 def ap(self): 

268 """ Numerator parameters of the hypergeometric function. """ 

269 return Tuple(*self.args[0]) 

270 

271 @property 

272 def bq(self): 

273 """ Denominator parameters of the hypergeometric function. """ 

274 return Tuple(*self.args[1]) 

275 

276 @property 

277 def _diffargs(self): 

278 return self.ap + self.bq 

279 

280 @property 

281 def eta(self): 

282 """ A quantity related to the convergence of the series. """ 

283 return sum(self.ap) - sum(self.bq) 

284 

285 @property 

286 def radius_of_convergence(self): 

287 """ 

288 Compute the radius of convergence of the defining series. 

289 

290 Explanation 

291 =========== 

292 

293 Note that even if this is not ``oo``, the function may still be 

294 evaluated outside of the radius of convergence by analytic 

295 continuation. But if this is zero, then the function is not actually 

296 defined anywhere else. 

297 

298 Examples 

299 ======== 

300 

301 >>> from sympy import hyper 

302 >>> from sympy.abc import z 

303 >>> hyper((1, 2), [3], z).radius_of_convergence 

304 1 

305 >>> hyper((1, 2, 3), [4], z).radius_of_convergence 

306 0 

307 >>> hyper((1, 2), (3, 4), z).radius_of_convergence 

308 oo 

309 

310 """ 

311 if any(a.is_integer and (a <= 0) == True for a in self.ap + self.bq): 

312 aints = [a for a in self.ap if a.is_Integer and (a <= 0) == True] 

313 bints = [a for a in self.bq if a.is_Integer and (a <= 0) == True] 

314 if len(aints) < len(bints): 

315 return S.Zero 

316 popped = False 

317 for b in bints: 

318 cancelled = False 

319 while aints: 

320 a = aints.pop() 

321 if a >= b: 

322 cancelled = True 

323 break 

324 popped = True 

325 if not cancelled: 

326 return S.Zero 

327 if aints or popped: 

328 # There are still non-positive numerator parameters. 

329 # This is a polynomial. 

330 return oo 

331 if len(self.ap) == len(self.bq) + 1: 

332 return S.One 

333 elif len(self.ap) <= len(self.bq): 

334 return oo 

335 else: 

336 return S.Zero 

337 

338 @property 

339 def convergence_statement(self): 

340 """ Return a condition on z under which the series converges. """ 

341 R = self.radius_of_convergence 

342 if R == 0: 

343 return False 

344 if R == oo: 

345 return True 

346 # The special functions and their approximations, page 44 

347 e = self.eta 

348 z = self.argument 

349 c1 = And(re(e) < 0, abs(z) <= 1) 

350 c2 = And(0 <= re(e), re(e) < 1, abs(z) <= 1, Ne(z, 1)) 

351 c3 = And(re(e) >= 1, abs(z) < 1) 

352 return Or(c1, c2, c3) 

353 

354 def _eval_simplify(self, **kwargs): 

355 from sympy.simplify.hyperexpand import hyperexpand 

356 return hyperexpand(self) 

357 

358 

359class meijerg(TupleParametersBase): 

360 r""" 

361 The Meijer G-function is defined by a Mellin-Barnes type integral that 

362 resembles an inverse Mellin transform. It generalizes the hypergeometric 

363 functions. 

364 

365 Explanation 

366 =========== 

367 

368 The Meijer G-function depends on four sets of parameters. There are 

369 "*numerator parameters*" 

370 $a_1, \ldots, a_n$ and $a_{n+1}, \ldots, a_p$, and there are 

371 "*denominator parameters*" 

372 $b_1, \ldots, b_m$ and $b_{m+1}, \ldots, b_q$. 

373 Confusingly, it is traditionally denoted as follows (note the position 

374 of $m$, $n$, $p$, $q$, and how they relate to the lengths of the four 

375 parameter vectors): 

376 

377 .. math :: 

378 G_{p,q}^{m,n} \left(\begin{matrix}a_1, \cdots, a_n & a_{n+1}, \cdots, a_p \\ 

379 b_1, \cdots, b_m & b_{m+1}, \cdots, b_q 

380 \end{matrix} \middle| z \right). 

381 

382 However, in SymPy the four parameter vectors are always available 

383 separately (see examples), so that there is no need to keep track of the 

384 decorating sub- and super-scripts on the G symbol. 

385 

386 The G function is defined as the following integral: 

387 

388 .. math :: 

389 \frac{1}{2 \pi i} \int_L \frac{\prod_{j=1}^m \Gamma(b_j - s) 

390 \prod_{j=1}^n \Gamma(1 - a_j + s)}{\prod_{j=m+1}^q \Gamma(1- b_j +s) 

391 \prod_{j=n+1}^p \Gamma(a_j - s)} z^s \mathrm{d}s, 

392 

393 where $\Gamma(z)$ is the gamma function. There are three possible 

394 contours which we will not describe in detail here (see the references). 

395 If the integral converges along more than one of them, the definitions 

396 agree. The contours all separate the poles of $\Gamma(1-a_j+s)$ 

397 from the poles of $\Gamma(b_k-s)$, so in particular the G function 

398 is undefined if $a_j - b_k \in \mathbb{Z}_{>0}$ for some 

399 $j \le n$ and $k \le m$. 

400 

401 The conditions under which one of the contours yields a convergent integral 

402 are complicated and we do not state them here, see the references. 

403 

404 Please note currently the Meijer G-function constructor does *not* check any 

405 convergence conditions. 

406 

407 Examples 

408 ======== 

409 

410 You can pass the parameters either as four separate vectors: 

411 

412 >>> from sympy import meijerg, Tuple, pprint 

413 >>> from sympy.abc import x, a 

414 >>> pprint(meijerg((1, 2), (a, 4), (5,), [], x), use_unicode=False) 

415 __1, 2 /1, 2 a, 4 | \ 

416 /__ | | x| 

417 \_|4, 1 \ 5 | / 

418 

419 Or as two nested vectors: 

420 

421 >>> pprint(meijerg([(1, 2), (3, 4)], ([5], Tuple()), x), use_unicode=False) 

422 __1, 2 /1, 2 3, 4 | \ 

423 /__ | | x| 

424 \_|4, 1 \ 5 | / 

425 

426 As with the hypergeometric function, the parameters may be passed as 

427 arbitrary iterables. Vectors of length zero and one also have to be 

428 passed as iterables. The parameters need not be constants, but if they 

429 depend on the argument then not much implemented functionality should be 

430 expected. 

431 

432 All the subvectors of parameters are available: 

433 

434 >>> from sympy import pprint 

435 >>> g = meijerg([1], [2], [3], [4], x) 

436 >>> pprint(g, use_unicode=False) 

437 __1, 1 /1 2 | \ 

438 /__ | | x| 

439 \_|2, 2 \3 4 | / 

440 >>> g.an 

441 (1,) 

442 >>> g.ap 

443 (1, 2) 

444 >>> g.aother 

445 (2,) 

446 >>> g.bm 

447 (3,) 

448 >>> g.bq 

449 (3, 4) 

450 >>> g.bother 

451 (4,) 

452 

453 The Meijer G-function generalizes the hypergeometric functions. 

454 In some cases it can be expressed in terms of hypergeometric functions, 

455 using Slater's theorem. For example: 

456 

457 >>> from sympy import hyperexpand 

458 >>> from sympy.abc import a, b, c 

459 >>> hyperexpand(meijerg([a], [], [c], [b], x), allow_hyper=True) 

460 x**c*gamma(-a + c + 1)*hyper((-a + c + 1,), 

461 (-b + c + 1,), -x)/gamma(-b + c + 1) 

462 

463 Thus the Meijer G-function also subsumes many named functions as special 

464 cases. You can use ``expand_func()`` or ``hyperexpand()`` to (try to) 

465 rewrite a Meijer G-function in terms of named special functions. For 

466 example: 

467 

468 >>> from sympy import expand_func, S 

469 >>> expand_func(meijerg([[],[]], [[0],[]], -x)) 

470 exp(x) 

471 >>> hyperexpand(meijerg([[],[]], [[S(1)/2],[0]], (x/2)**2)) 

472 sin(x)/sqrt(pi) 

473 

474 See Also 

475 ======== 

476 

477 hyper 

478 sympy.simplify.hyperexpand 

479 

480 References 

481 ========== 

482 

483 .. [1] Luke, Y. L. (1969), The Special Functions and Their Approximations, 

484 Volume 1 

485 .. [2] https://en.wikipedia.org/wiki/Meijer_G-function 

486 

487 """ 

488 

489 

490 def __new__(cls, *args, **kwargs): 

491 if len(args) == 5: 

492 args = [(args[0], args[1]), (args[2], args[3]), args[4]] 

493 if len(args) != 3: 

494 raise TypeError("args must be either as, as', bs, bs', z or " 

495 "as, bs, z") 

496 

497 def tr(p): 

498 if len(p) != 2: 

499 raise TypeError("wrong argument") 

500 return TupleArg(_prep_tuple(p[0]), _prep_tuple(p[1])) 

501 

502 arg0, arg1 = tr(args[0]), tr(args[1]) 

503 if Tuple(arg0, arg1).has(oo, zoo, -oo): 

504 raise ValueError("G-function parameters must be finite") 

505 if any((a - b).is_Integer and a - b > 0 

506 for a in arg0[0] for b in arg1[0]): 

507 raise ValueError("no parameter a1, ..., an may differ from " 

508 "any b1, ..., bm by a positive integer") 

509 

510 # TODO should we check convergence conditions? 

511 return Function.__new__(cls, arg0, arg1, args[2], **kwargs) 

512 

513 def fdiff(self, argindex=3): 

514 if argindex != 3: 

515 return self._diff_wrt_parameter(argindex[1]) 

516 if len(self.an) >= 1: 

517 a = list(self.an) 

518 a[0] -= 1 

519 G = meijerg(a, self.aother, self.bm, self.bother, self.argument) 

520 return 1/self.argument * ((self.an[0] - 1)*self + G) 

521 elif len(self.bm) >= 1: 

522 b = list(self.bm) 

523 b[0] += 1 

524 G = meijerg(self.an, self.aother, b, self.bother, self.argument) 

525 return 1/self.argument * (self.bm[0]*self - G) 

526 else: 

527 return S.Zero 

528 

529 def _diff_wrt_parameter(self, idx): 

530 # Differentiation wrt a parameter can only be done in very special 

531 # cases. In particular, if we want to differentiate with respect to 

532 # `a`, all other gamma factors have to reduce to rational functions. 

533 # 

534 # Let MT denote mellin transform. Suppose T(-s) is the gamma factor 

535 # appearing in the definition of G. Then 

536 # 

537 # MT(log(z)G(z)) = d/ds T(s) = d/da T(s) + ... 

538 # 

539 # Thus d/da G(z) = log(z)G(z) - ... 

540 # The ... can be evaluated as a G function under the above conditions, 

541 # the formula being most easily derived by using 

542 # 

543 # d Gamma(s + n) Gamma(s + n) / 1 1 1 \ 

544 # -- ------------ = ------------ | - + ---- + ... + --------- | 

545 # ds Gamma(s) Gamma(s) \ s s + 1 s + n - 1 / 

546 # 

547 # which follows from the difference equation of the digamma function. 

548 # (There is a similar equation for -n instead of +n). 

549 

550 # We first figure out how to pair the parameters. 

551 an = list(self.an) 

552 ap = list(self.aother) 

553 bm = list(self.bm) 

554 bq = list(self.bother) 

555 if idx < len(an): 

556 an.pop(idx) 

557 else: 

558 idx -= len(an) 

559 if idx < len(ap): 

560 ap.pop(idx) 

561 else: 

562 idx -= len(ap) 

563 if idx < len(bm): 

564 bm.pop(idx) 

565 else: 

566 bq.pop(idx - len(bm)) 

567 pairs1 = [] 

568 pairs2 = [] 

569 for l1, l2, pairs in [(an, bq, pairs1), (ap, bm, pairs2)]: 

570 while l1: 

571 x = l1.pop() 

572 found = None 

573 for i, y in enumerate(l2): 

574 if not Mod((x - y).simplify(), 1): 

575 found = i 

576 break 

577 if found is None: 

578 raise NotImplementedError('Derivative not expressible ' 

579 'as G-function?') 

580 y = l2[i] 

581 l2.pop(i) 

582 pairs.append((x, y)) 

583 

584 # Now build the result. 

585 res = log(self.argument)*self 

586 

587 for a, b in pairs1: 

588 sign = 1 

589 n = a - b 

590 base = b 

591 if n < 0: 

592 sign = -1 

593 n = b - a 

594 base = a 

595 for k in range(n): 

596 res -= sign*meijerg(self.an + (base + k + 1,), self.aother, 

597 self.bm, self.bother + (base + k + 0,), 

598 self.argument) 

599 

600 for a, b in pairs2: 

601 sign = 1 

602 n = b - a 

603 base = a 

604 if n < 0: 

605 sign = -1 

606 n = a - b 

607 base = b 

608 for k in range(n): 

609 res -= sign*meijerg(self.an, self.aother + (base + k + 1,), 

610 self.bm + (base + k + 0,), self.bother, 

611 self.argument) 

612 

613 return res 

614 

615 def get_period(self): 

616 """ 

617 Return a number $P$ such that $G(x*exp(I*P)) == G(x)$. 

618 

619 Examples 

620 ======== 

621 

622 >>> from sympy import meijerg, pi, S 

623 >>> from sympy.abc import z 

624 

625 >>> meijerg([1], [], [], [], z).get_period() 

626 2*pi 

627 >>> meijerg([pi], [], [], [], z).get_period() 

628 oo 

629 >>> meijerg([1, 2], [], [], [], z).get_period() 

630 oo 

631 >>> meijerg([1,1], [2], [1, S(1)/2, S(1)/3], [1], z).get_period() 

632 12*pi 

633 

634 """ 

635 # This follows from slater's theorem. 

636 def compute(l): 

637 # first check that no two differ by an integer 

638 for i, b in enumerate(l): 

639 if not b.is_Rational: 

640 return oo 

641 for j in range(i + 1, len(l)): 

642 if not Mod((b - l[j]).simplify(), 1): 

643 return oo 

644 return reduce(ilcm, (x.q for x in l), 1) 

645 beta = compute(self.bm) 

646 alpha = compute(self.an) 

647 p, q = len(self.ap), len(self.bq) 

648 if p == q: 

649 if oo in (alpha, beta): 

650 return oo 

651 return 2*pi*ilcm(alpha, beta) 

652 elif p < q: 

653 return 2*pi*beta 

654 else: 

655 return 2*pi*alpha 

656 

657 def _eval_expand_func(self, **hints): 

658 from sympy.simplify.hyperexpand import hyperexpand 

659 return hyperexpand(self) 

660 

661 def _eval_evalf(self, prec): 

662 # The default code is insufficient for polar arguments. 

663 # mpmath provides an optional argument "r", which evaluates 

664 # G(z**(1/r)). I am not sure what its intended use is, but we hijack it 

665 # here in the following way: to evaluate at a number z of |argument| 

666 # less than (say) n*pi, we put r=1/n, compute z' = root(z, n) 

667 # (carefully so as not to loose the branch information), and evaluate 

668 # G(z'**(1/r)) = G(z'**n) = G(z). 

669 import mpmath 

670 znum = self.argument._eval_evalf(prec) 

671 if znum.has(exp_polar): 

672 znum, branch = znum.as_coeff_mul(exp_polar) 

673 if len(branch) != 1: 

674 return 

675 branch = branch[0].args[0]/I 

676 else: 

677 branch = S.Zero 

678 n = ceiling(abs(branch/pi)) + 1 

679 znum = znum**(S.One/n)*exp(I*branch / n) 

680 

681 # Convert all args to mpf or mpc 

682 try: 

683 [z, r, ap, bq] = [arg._to_mpmath(prec) 

684 for arg in [znum, 1/n, self.args[0], self.args[1]]] 

685 except ValueError: 

686 return 

687 

688 with mpmath.workprec(prec): 

689 v = mpmath.meijerg(ap, bq, z, r) 

690 

691 return Expr._from_mpmath(v, prec) 

692 

693 def _eval_as_leading_term(self, x, logx=None, cdir=0): 

694 from sympy.simplify.hyperexpand import hyperexpand 

695 return hyperexpand(self).as_leading_term(x, logx=logx, cdir=cdir) 

696 

697 def integrand(self, s): 

698 """ Get the defining integrand D(s). """ 

699 from sympy.functions.special.gamma_functions import gamma 

700 return self.argument**s \ 

701 * Mul(*(gamma(b - s) for b in self.bm)) \ 

702 * Mul(*(gamma(1 - a + s) for a in self.an)) \ 

703 / Mul(*(gamma(1 - b + s) for b in self.bother)) \ 

704 / Mul(*(gamma(a - s) for a in self.aother)) 

705 

706 @property 

707 def argument(self): 

708 """ Argument of the Meijer G-function. """ 

709 return self.args[2] 

710 

711 @property 

712 def an(self): 

713 """ First set of numerator parameters. """ 

714 return Tuple(*self.args[0][0]) 

715 

716 @property 

717 def ap(self): 

718 """ Combined numerator parameters. """ 

719 return Tuple(*(self.args[0][0] + self.args[0][1])) 

720 

721 @property 

722 def aother(self): 

723 """ Second set of numerator parameters. """ 

724 return Tuple(*self.args[0][1]) 

725 

726 @property 

727 def bm(self): 

728 """ First set of denominator parameters. """ 

729 return Tuple(*self.args[1][0]) 

730 

731 @property 

732 def bq(self): 

733 """ Combined denominator parameters. """ 

734 return Tuple(*(self.args[1][0] + self.args[1][1])) 

735 

736 @property 

737 def bother(self): 

738 """ Second set of denominator parameters. """ 

739 return Tuple(*self.args[1][1]) 

740 

741 @property 

742 def _diffargs(self): 

743 return self.ap + self.bq 

744 

745 @property 

746 def nu(self): 

747 """ A quantity related to the convergence region of the integral, 

748 c.f. references. """ 

749 return sum(self.bq) - sum(self.ap) 

750 

751 @property 

752 def delta(self): 

753 """ A quantity related to the convergence region of the integral, 

754 c.f. references. """ 

755 return len(self.bm) + len(self.an) - S(len(self.ap) + len(self.bq))/2 

756 

757 @property 

758 def is_number(self): 

759 """ Returns true if expression has numeric data only. """ 

760 return not self.free_symbols 

761 

762 

763class HyperRep(Function): 

764 """ 

765 A base class for "hyper representation functions". 

766 

767 This is used exclusively in ``hyperexpand()``, but fits more logically here. 

768 

769 pFq is branched at 1 if p == q+1. For use with slater-expansion, we want 

770 define an "analytic continuation" to all polar numbers, which is 

771 continuous on circles and on the ray t*exp_polar(I*pi). Moreover, we want 

772 a "nice" expression for the various cases. 

773 

774 This base class contains the core logic, concrete derived classes only 

775 supply the actual functions. 

776 

777 """ 

778 

779 

780 @classmethod 

781 def eval(cls, *args): 

782 newargs = tuple(map(unpolarify, args[:-1])) + args[-1:] 

783 if args != newargs: 

784 return cls(*newargs) 

785 

786 @classmethod 

787 def _expr_small(cls, x): 

788 """ An expression for F(x) which holds for |x| < 1. """ 

789 raise NotImplementedError 

790 

791 @classmethod 

792 def _expr_small_minus(cls, x): 

793 """ An expression for F(-x) which holds for |x| < 1. """ 

794 raise NotImplementedError 

795 

796 @classmethod 

797 def _expr_big(cls, x, n): 

798 """ An expression for F(exp_polar(2*I*pi*n)*x), |x| > 1. """ 

799 raise NotImplementedError 

800 

801 @classmethod 

802 def _expr_big_minus(cls, x, n): 

803 """ An expression for F(exp_polar(2*I*pi*n + pi*I)*x), |x| > 1. """ 

804 raise NotImplementedError 

805 

806 def _eval_rewrite_as_nonrep(self, *args, **kwargs): 

807 x, n = self.args[-1].extract_branch_factor(allow_half=True) 

808 minus = False 

809 newargs = self.args[:-1] + (x,) 

810 if not n.is_Integer: 

811 minus = True 

812 n -= S.Half 

813 newerargs = newargs + (n,) 

814 if minus: 

815 small = self._expr_small_minus(*newargs) 

816 big = self._expr_big_minus(*newerargs) 

817 else: 

818 small = self._expr_small(*newargs) 

819 big = self._expr_big(*newerargs) 

820 

821 if big == small: 

822 return small 

823 return Piecewise((big, abs(x) > 1), (small, True)) 

824 

825 def _eval_rewrite_as_nonrepsmall(self, *args, **kwargs): 

826 x, n = self.args[-1].extract_branch_factor(allow_half=True) 

827 args = self.args[:-1] + (x,) 

828 if not n.is_Integer: 

829 return self._expr_small_minus(*args) 

830 return self._expr_small(*args) 

831 

832 

833class HyperRep_power1(HyperRep): 

834 """ Return a representative for hyper([-a], [], z) == (1 - z)**a. """ 

835 

836 @classmethod 

837 def _expr_small(cls, a, x): 

838 return (1 - x)**a 

839 

840 @classmethod 

841 def _expr_small_minus(cls, a, x): 

842 return (1 + x)**a 

843 

844 @classmethod 

845 def _expr_big(cls, a, x, n): 

846 if a.is_integer: 

847 return cls._expr_small(a, x) 

848 return (x - 1)**a*exp((2*n - 1)*pi*I*a) 

849 

850 @classmethod 

851 def _expr_big_minus(cls, a, x, n):