test_hyperexpand.py 40 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060
  1. from sympy.core.random import randrange
  2. from sympy.simplify.hyperexpand import (ShiftA, ShiftB, UnShiftA, UnShiftB,
  3. MeijerShiftA, MeijerShiftB, MeijerShiftC, MeijerShiftD,
  4. MeijerUnShiftA, MeijerUnShiftB, MeijerUnShiftC,
  5. MeijerUnShiftD,
  6. ReduceOrder, reduce_order, apply_operators,
  7. devise_plan, make_derivative_operator, Formula,
  8. hyperexpand, Hyper_Function, G_Function,
  9. reduce_order_meijer,
  10. build_hypergeometric_formula)
  11. from sympy.concrete.summations import Sum
  12. from sympy.core.containers import Tuple
  13. from sympy.core.expr import Expr
  14. from sympy.core.numbers import I
  15. from sympy.core.singleton import S
  16. from sympy.core.symbol import symbols
  17. from sympy.functions.combinatorial.factorials import binomial
  18. from sympy.functions.elementary.piecewise import Piecewise
  19. from sympy.functions.special.hyper import (hyper, meijerg)
  20. from sympy.abc import z, a, b, c
  21. from sympy.testing.pytest import XFAIL, raises, slow, ON_CI, skip
  22. from sympy.core.random import verify_numerically as tn
  23. from sympy.core.numbers import (Rational, pi)
  24. from sympy.functions.elementary.exponential import (exp, exp_polar, log)
  25. from sympy.functions.elementary.hyperbolic import atanh
  26. from sympy.functions.elementary.miscellaneous import sqrt
  27. from sympy.functions.elementary.trigonometric import (asin, cos, sin)
  28. from sympy.functions.special.bessel import besseli
  29. from sympy.functions.special.error_functions import erf
  30. from sympy.functions.special.gamma_functions import (gamma, lowergamma)
  31. def test_branch_bug():
  32. assert hyperexpand(hyper((Rational(-1, 3), S.Half), (Rational(2, 3), Rational(3, 2)), -z)) == \
  33. -z**S('1/3')*lowergamma(exp_polar(I*pi)/3, z)/5 \
  34. + sqrt(pi)*erf(sqrt(z))/(5*sqrt(z))
  35. assert hyperexpand(meijerg([Rational(7, 6), 1], [], [Rational(2, 3)], [Rational(1, 6), 0], z)) == \
  36. 2*z**S('2/3')*(2*sqrt(pi)*erf(sqrt(z))/sqrt(z) - 2*lowergamma(
  37. Rational(2, 3), z)/z**S('2/3'))*gamma(Rational(2, 3))/gamma(Rational(5, 3))
  38. def test_hyperexpand():
  39. # Luke, Y. L. (1969), The Special Functions and Their Approximations,
  40. # Volume 1, section 6.2
  41. assert hyperexpand(hyper([], [], z)) == exp(z)
  42. assert hyperexpand(hyper([1, 1], [2], -z)*z) == log(1 + z)
  43. assert hyperexpand(hyper([], [S.Half], -z**2/4)) == cos(z)
  44. assert hyperexpand(z*hyper([], [S('3/2')], -z**2/4)) == sin(z)
  45. assert hyperexpand(hyper([S('1/2'), S('1/2')], [S('3/2')], z**2)*z) \
  46. == asin(z)
  47. assert isinstance(Sum(binomial(2, z)*z**2, (z, 0, a)).doit(), Expr)
  48. def can_do(ap, bq, numerical=True, div=1, lowerplane=False):
  49. r = hyperexpand(hyper(ap, bq, z))
  50. if r.has(hyper):
  51. return False
  52. if not numerical:
  53. return True
  54. repl = {}
  55. randsyms = r.free_symbols - {z}
  56. while randsyms:
  57. # Only randomly generated parameters are checked.
  58. for n, ai in enumerate(randsyms):
  59. repl[ai] = randcplx(n)/div
  60. if not any(b.is_Integer and b <= 0 for b in Tuple(*bq).subs(repl)):
  61. break
  62. [a, b, c, d] = [2, -1, 3, 1]
  63. if lowerplane:
  64. [a, b, c, d] = [2, -2, 3, -1]
  65. return tn(
  66. hyper(ap, bq, z).subs(repl),
  67. r.replace(exp_polar, exp).subs(repl),
  68. z, a=a, b=b, c=c, d=d)
  69. def test_roach():
  70. # Kelly B. Roach. Meijer G Function Representations.
  71. # Section "Gallery"
  72. assert can_do([S.Half], [Rational(9, 2)])
  73. assert can_do([], [1, Rational(5, 2), 4])
  74. assert can_do([Rational(-1, 2), 1, 2], [3, 4])
  75. assert can_do([Rational(1, 3)], [Rational(-2, 3), Rational(-1, 2), S.Half, 1])
  76. assert can_do([Rational(-3, 2), Rational(-1, 2)], [Rational(-5, 2), 1])
  77. assert can_do([Rational(-3, 2), ], [Rational(-1, 2), S.Half]) # shine-integral
  78. assert can_do([Rational(-3, 2), Rational(-1, 2)], [2]) # elliptic integrals
  79. @XFAIL
  80. def test_roach_fail():
  81. assert can_do([Rational(-1, 2), 1], [Rational(1, 4), S.Half, Rational(3, 4)]) # PFDD
  82. assert can_do([Rational(3, 2)], [Rational(5, 2), 5]) # struve function
  83. assert can_do([Rational(-1, 2), S.Half, 1], [Rational(3, 2), Rational(5, 2)]) # polylog, pfdd
  84. assert can_do([1, 2, 3], [S.Half, 4]) # XXX ?
  85. assert can_do([S.Half], [Rational(-1, 3), Rational(-1, 2), Rational(-2, 3)]) # PFDD ?
  86. # For the long table tests, see end of file
  87. def test_polynomial():
  88. from sympy.core.numbers import oo
  89. assert hyperexpand(hyper([], [-1], z)) is oo
  90. assert hyperexpand(hyper([-2], [-1], z)) is oo
  91. assert hyperexpand(hyper([0, 0], [-1], z)) == 1
  92. assert can_do([-5, -2, randcplx(), randcplx()], [-10, randcplx()])
  93. assert hyperexpand(hyper((-1, 1), (-2,), z)) == 1 + z/2
  94. def test_hyperexpand_bases():
  95. assert hyperexpand(hyper([2], [a], z)) == \
  96. a + z**(-a + 1)*(-a**2 + 3*a + z*(a - 1) - 2)*exp(z)* \
  97. lowergamma(a - 1, z) - 1
  98. # TODO [a+1, aRational(-1, 2)], [2*a]
  99. assert hyperexpand(hyper([1, 2], [3], z)) == -2/z - 2*log(-z + 1)/z**2
  100. assert hyperexpand(hyper([S.Half, 2], [Rational(3, 2)], z)) == \
  101. -1/(2*z - 2) + atanh(sqrt(z))/sqrt(z)/2
  102. assert hyperexpand(hyper([S.Half, S.Half], [Rational(5, 2)], z)) == \
  103. (-3*z + 3)/4/(z*sqrt(-z + 1)) \
  104. + (6*z - 3)*asin(sqrt(z))/(4*z**Rational(3, 2))
  105. assert hyperexpand(hyper([1, 2], [Rational(3, 2)], z)) == -1/(2*z - 2) \
  106. - asin(sqrt(z))/(sqrt(z)*(2*z - 2)*sqrt(-z + 1))
  107. assert hyperexpand(hyper([Rational(-1, 2) - 1, 1, 2], [S.Half, 3], z)) == \
  108. sqrt(z)*(z*Rational(6, 7) - Rational(6, 5))*atanh(sqrt(z)) \
  109. + (-30*z**2 + 32*z - 6)/35/z - 6*log(-z + 1)/(35*z**2)
  110. assert hyperexpand(hyper([1 + S.Half, 1, 1], [2, 2], z)) == \
  111. -4*log(sqrt(-z + 1)/2 + S.Half)/z
  112. # TODO hyperexpand(hyper([a], [2*a + 1], z))
  113. # TODO [S.Half, a], [Rational(3, 2), a+1]
  114. assert hyperexpand(hyper([2], [b, 1], z)) == \
  115. z**(-b/2 + S.Half)*besseli(b - 1, 2*sqrt(z))*gamma(b) \
  116. + z**(-b/2 + 1)*besseli(b, 2*sqrt(z))*gamma(b)
  117. # TODO [a], [a - S.Half, 2*a]
  118. def test_hyperexpand_parametric():
  119. assert hyperexpand(hyper([a, S.Half + a], [S.Half], z)) \
  120. == (1 + sqrt(z))**(-2*a)/2 + (1 - sqrt(z))**(-2*a)/2
  121. assert hyperexpand(hyper([a, Rational(-1, 2) + a], [2*a], z)) \
  122. == 2**(2*a - 1)*((-z + 1)**S.Half + 1)**(-2*a + 1)
  123. def test_shifted_sum():
  124. from sympy.simplify.simplify import simplify
  125. assert simplify(hyperexpand(z**4*hyper([2], [3, S('3/2')], -z**2))) \
  126. == z*sin(2*z) + (-z**2 + S.Half)*cos(2*z) - S.Half
  127. def _randrat():
  128. """ Steer clear of integers. """
  129. return S(randrange(25) + 10)/50
  130. def randcplx(offset=-1):
  131. """ Polys is not good with real coefficients. """
  132. return _randrat() + I*_randrat() + I*(1 + offset)
  133. @slow
  134. def test_formulae():
  135. from sympy.simplify.hyperexpand import FormulaCollection
  136. formulae = FormulaCollection().formulae
  137. for formula in formulae:
  138. h = formula.func(formula.z)
  139. rep = {}
  140. for n, sym in enumerate(formula.symbols):
  141. rep[sym] = randcplx(n)
  142. # NOTE hyperexpand returns truly branched functions. We know we are
  143. # on the main sheet, but numerical evaluation can still go wrong
  144. # (e.g. if exp_polar cannot be evalf'd).
  145. # Just replace all exp_polar by exp, this usually works.
  146. # first test if the closed-form is actually correct
  147. h = h.subs(rep)
  148. closed_form = formula.closed_form.subs(rep).rewrite('nonrepsmall')
  149. z = formula.z
  150. assert tn(h, closed_form.replace(exp_polar, exp), z)
  151. # now test the computed matrix
  152. cl = (formula.C * formula.B)[0].subs(rep).rewrite('nonrepsmall')
  153. assert tn(closed_form.replace(
  154. exp_polar, exp), cl.replace(exp_polar, exp), z)
  155. deriv1 = z*formula.B.applyfunc(lambda t: t.rewrite(
  156. 'nonrepsmall')).diff(z)
  157. deriv2 = formula.M * formula.B
  158. for d1, d2 in zip(deriv1, deriv2):
  159. assert tn(d1.subs(rep).replace(exp_polar, exp),
  160. d2.subs(rep).rewrite('nonrepsmall').replace(exp_polar, exp), z)
  161. def test_meijerg_formulae():
  162. from sympy.simplify.hyperexpand import MeijerFormulaCollection
  163. formulae = MeijerFormulaCollection().formulae
  164. for sig in formulae:
  165. for formula in formulae[sig]:
  166. g = meijerg(formula.func.an, formula.func.ap,
  167. formula.func.bm, formula.func.bq,
  168. formula.z)
  169. rep = {}
  170. for sym in formula.symbols:
  171. rep[sym] = randcplx()
  172. # first test if the closed-form is actually correct
  173. g = g.subs(rep)
  174. closed_form = formula.closed_form.subs(rep)
  175. z = formula.z
  176. assert tn(g, closed_form, z)
  177. # now test the computed matrix
  178. cl = (formula.C * formula.B)[0].subs(rep)
  179. assert tn(closed_form, cl, z)
  180. deriv1 = z*formula.B.diff(z)
  181. deriv2 = formula.M * formula.B
  182. for d1, d2 in zip(deriv1, deriv2):
  183. assert tn(d1.subs(rep), d2.subs(rep), z)
  184. def op(f):
  185. return z*f.diff(z)
  186. def test_plan():
  187. assert devise_plan(Hyper_Function([0], ()),
  188. Hyper_Function([0], ()), z) == []
  189. with raises(ValueError):
  190. devise_plan(Hyper_Function([1], ()), Hyper_Function((), ()), z)
  191. with raises(ValueError):
  192. devise_plan(Hyper_Function([2], [1]), Hyper_Function([2], [2]), z)
  193. with raises(ValueError):
  194. devise_plan(Hyper_Function([2], []), Hyper_Function([S("1/2")], []), z)
  195. # We cannot use pi/(10000 + n) because polys is insanely slow.
  196. a1, a2, b1 = (randcplx(n) for n in range(3))
  197. b1 += 2*I
  198. h = hyper([a1, a2], [b1], z)
  199. h2 = hyper((a1 + 1, a2), [b1], z)
  200. assert tn(apply_operators(h,
  201. devise_plan(Hyper_Function((a1 + 1, a2), [b1]),
  202. Hyper_Function((a1, a2), [b1]), z), op),
  203. h2, z)
  204. h2 = hyper((a1 + 1, a2 - 1), [b1], z)
  205. assert tn(apply_operators(h,
  206. devise_plan(Hyper_Function((a1 + 1, a2 - 1), [b1]),
  207. Hyper_Function((a1, a2), [b1]), z), op),
  208. h2, z)
  209. def test_plan_derivatives():
  210. a1, a2, a3 = 1, 2, S('1/2')
  211. b1, b2 = 3, S('5/2')
  212. h = Hyper_Function((a1, a2, a3), (b1, b2))
  213. h2 = Hyper_Function((a1 + 1, a2 + 1, a3 + 2), (b1 + 1, b2 + 1))
  214. ops = devise_plan(h2, h, z)
  215. f = Formula(h, z, h(z), [])
  216. deriv = make_derivative_operator(f.M, z)
  217. assert tn((apply_operators(f.C, ops, deriv)*f.B)[0], h2(z), z)
  218. h2 = Hyper_Function((a1, a2 - 1, a3 - 2), (b1 - 1, b2 - 1))
  219. ops = devise_plan(h2, h, z)
  220. assert tn((apply_operators(f.C, ops, deriv)*f.B)[0], h2(z), z)
  221. def test_reduction_operators():
  222. a1, a2, b1 = (randcplx(n) for n in range(3))
  223. h = hyper([a1], [b1], z)
  224. assert ReduceOrder(2, 0) is None
  225. assert ReduceOrder(2, -1) is None
  226. assert ReduceOrder(1, S('1/2')) is None
  227. h2 = hyper((a1, a2), (b1, a2), z)
  228. assert tn(ReduceOrder(a2, a2).apply(h, op), h2, z)
  229. h2 = hyper((a1, a2 + 1), (b1, a2), z)
  230. assert tn(ReduceOrder(a2 + 1, a2).apply(h, op), h2, z)
  231. h2 = hyper((a2 + 4, a1), (b1, a2), z)
  232. assert tn(ReduceOrder(a2 + 4, a2).apply(h, op), h2, z)
  233. # test several step order reduction
  234. ap = (a2 + 4, a1, b1 + 1)
  235. bq = (a2, b1, b1)
  236. func, ops = reduce_order(Hyper_Function(ap, bq))
  237. assert func.ap == (a1,)
  238. assert func.bq == (b1,)
  239. assert tn(apply_operators(h, ops, op), hyper(ap, bq, z), z)
  240. def test_shift_operators():
  241. a1, a2, b1, b2, b3 = (randcplx(n) for n in range(5))
  242. h = hyper((a1, a2), (b1, b2, b3), z)
  243. raises(ValueError, lambda: ShiftA(0))
  244. raises(ValueError, lambda: ShiftB(1))
  245. assert tn(ShiftA(a1).apply(h, op), hyper((a1 + 1, a2), (b1, b2, b3), z), z)
  246. assert tn(ShiftA(a2).apply(h, op), hyper((a1, a2 + 1), (b1, b2, b3), z), z)
  247. assert tn(ShiftB(b1).apply(h, op), hyper((a1, a2), (b1 - 1, b2, b3), z), z)
  248. assert tn(ShiftB(b2).apply(h, op), hyper((a1, a2), (b1, b2 - 1, b3), z), z)
  249. assert tn(ShiftB(b3).apply(h, op), hyper((a1, a2), (b1, b2, b3 - 1), z), z)
  250. def test_ushift_operators():
  251. a1, a2, b1, b2, b3 = (randcplx(n) for n in range(5))
  252. h = hyper((a1, a2), (b1, b2, b3), z)
  253. raises(ValueError, lambda: UnShiftA((1,), (), 0, z))
  254. raises(ValueError, lambda: UnShiftB((), (-1,), 0, z))
  255. raises(ValueError, lambda: UnShiftA((1,), (0, -1, 1), 0, z))
  256. raises(ValueError, lambda: UnShiftB((0, 1), (1,), 0, z))
  257. s = UnShiftA((a1, a2), (b1, b2, b3), 0, z)
  258. assert tn(s.apply(h, op), hyper((a1 - 1, a2), (b1, b2, b3), z), z)
  259. s = UnShiftA((a1, a2), (b1, b2, b3), 1, z)
  260. assert tn(s.apply(h, op), hyper((a1, a2 - 1), (b1, b2, b3), z), z)
  261. s = UnShiftB((a1, a2), (b1, b2, b3), 0, z)
  262. assert tn(s.apply(h, op), hyper((a1, a2), (b1 + 1, b2, b3), z), z)
  263. s = UnShiftB((a1, a2), (b1, b2, b3), 1, z)
  264. assert tn(s.apply(h, op), hyper((a1, a2), (b1, b2 + 1, b3), z), z)
  265. s = UnShiftB((a1, a2), (b1, b2, b3), 2, z)
  266. assert tn(s.apply(h, op), hyper((a1, a2), (b1, b2, b3 + 1), z), z)
  267. def can_do_meijer(a1, a2, b1, b2, numeric=True):
  268. """
  269. This helper function tries to hyperexpand() the meijer g-function
  270. corresponding to the parameters a1, a2, b1, b2.
  271. It returns False if this expansion still contains g-functions.
  272. If numeric is True, it also tests the so-obtained formula numerically
  273. (at random values) and returns False if the test fails.
  274. Else it returns True.
  275. """
  276. from sympy.core.function import expand
  277. from sympy.functions.elementary.complexes import unpolarify
  278. r = hyperexpand(meijerg(a1, a2, b1, b2, z))
  279. if r.has(meijerg):
  280. return False
  281. # NOTE hyperexpand() returns a truly branched function, whereas numerical
  282. # evaluation only works on the main branch. Since we are evaluating on
  283. # the main branch, this should not be a problem, but expressions like
  284. # exp_polar(I*pi/2*x)**a are evaluated incorrectly. We thus have to get
  285. # rid of them. The expand heuristically does this...
  286. r = unpolarify(expand(r, force=True, power_base=True, power_exp=False,
  287. mul=False, log=False, multinomial=False, basic=False))
  288. if not numeric:
  289. return True
  290. repl = {}
  291. for n, ai in enumerate(meijerg(a1, a2, b1, b2, z).free_symbols - {z}):
  292. repl[ai] = randcplx(n)
  293. return tn(meijerg(a1, a2, b1, b2, z).subs(repl), r.subs(repl), z)
  294. @slow
  295. def test_meijerg_expand():
  296. from sympy.simplify.gammasimp import gammasimp
  297. from sympy.simplify.simplify import simplify
  298. # from mpmath docs
  299. assert hyperexpand(meijerg([[], []], [[0], []], -z)) == exp(z)
  300. assert hyperexpand(meijerg([[1, 1], []], [[1], [0]], z)) == \
  301. log(z + 1)
  302. assert hyperexpand(meijerg([[1, 1], []], [[1], [1]], z)) == \
  303. z/(z + 1)
  304. assert hyperexpand(meijerg([[], []], [[S.Half], [0]], (z/2)**2)) \
  305. == sin(z)/sqrt(pi)
  306. assert hyperexpand(meijerg([[], []], [[0], [S.Half]], (z/2)**2)) \
  307. == cos(z)/sqrt(pi)
  308. assert can_do_meijer([], [a], [a - 1, a - S.Half], [])
  309. assert can_do_meijer([], [], [a/2], [-a/2], False) # branches...
  310. assert can_do_meijer([a], [b], [a], [b, a - 1])
  311. # wikipedia
  312. assert hyperexpand(meijerg([1], [], [], [0], z)) == \
  313. Piecewise((0, abs(z) < 1), (1, abs(1/z) < 1),
  314. (meijerg([1], [], [], [0], z), True))
  315. assert hyperexpand(meijerg([], [1], [0], [], z)) == \
  316. Piecewise((1, abs(z) < 1), (0, abs(1/z) < 1),
  317. (meijerg([], [1], [0], [], z), True))
  318. # The Special Functions and their Approximations
  319. assert can_do_meijer([], [], [a + b/2], [a, a - b/2, a + S.Half])
  320. assert can_do_meijer(
  321. [], [], [a], [b], False) # branches only agree for small z
  322. assert can_do_meijer([], [S.Half], [a], [-a])
  323. assert can_do_meijer([], [], [a, b], [])
  324. assert can_do_meijer([], [], [a, b], [])
  325. assert can_do_meijer([], [], [a, a + S.Half], [b, b + S.Half])
  326. assert can_do_meijer([], [], [a, -a], [0, S.Half], False) # dito
  327. assert can_do_meijer([], [], [a, a + S.Half, b, b + S.Half], [])
  328. assert can_do_meijer([S.Half], [], [0], [a, -a])
  329. assert can_do_meijer([S.Half], [], [a], [0, -a], False) # dito
  330. assert can_do_meijer([], [a - S.Half], [a, b], [a - S.Half], False)
  331. assert can_do_meijer([], [a + S.Half], [a + b, a - b, a], [], False)
  332. assert can_do_meijer([a + S.Half], [], [b, 2*a - b, a], [], False)
  333. # This for example is actually zero.
  334. assert can_do_meijer([], [], [], [a, b])
  335. # Testing a bug:
  336. assert hyperexpand(meijerg([0, 2], [], [], [-1, 1], z)) == \
  337. Piecewise((0, abs(z) < 1),
  338. (z*(1 - 1/z**2)/2, abs(1/z) < 1),
  339. (meijerg([0, 2], [], [], [-1, 1], z), True))
  340. # Test that the simplest possible answer is returned:
  341. assert gammasimp(simplify(hyperexpand(
  342. meijerg([1], [1 - a], [-a/2, -a/2 + S.Half], [], 1/z)))) == \
  343. -2*sqrt(pi)*(sqrt(z + 1) + 1)**a/a
  344. # Test that hyper is returned
  345. assert hyperexpand(meijerg([1], [], [a], [0, 0], z)) == hyper(
  346. (a,), (a + 1, a + 1), z*exp_polar(I*pi))*z**a*gamma(a)/gamma(a + 1)**2
  347. # Test place option
  348. f = meijerg(((0, 1), ()), ((S.Half,), (0,)), z**2)
  349. assert hyperexpand(f) == sqrt(pi)/sqrt(1 + z**(-2))
  350. assert hyperexpand(f, place=0) == sqrt(pi)*z/sqrt(z**2 + 1)
  351. def test_meijerg_lookup():
  352. from sympy.functions.special.error_functions import (Ci, Si)
  353. from sympy.functions.special.gamma_functions import uppergamma
  354. assert hyperexpand(meijerg([a], [], [b, a], [], z)) == \
  355. z**b*exp(z)*gamma(-a + b + 1)*uppergamma(a - b, z)
  356. assert hyperexpand(meijerg([0], [], [0, 0], [], z)) == \
  357. exp(z)*uppergamma(0, z)
  358. assert can_do_meijer([a], [], [b, a + 1], [])
  359. assert can_do_meijer([a], [], [b + 2, a], [])
  360. assert can_do_meijer([a], [], [b - 2, a], [])
  361. assert hyperexpand(meijerg([a], [], [a, a, a - S.Half], [], z)) == \
  362. -sqrt(pi)*z**(a - S.Half)*(2*cos(2*sqrt(z))*(Si(2*sqrt(z)) - pi/2)
  363. - 2*sin(2*sqrt(z))*Ci(2*sqrt(z))) == \
  364. hyperexpand(meijerg([a], [], [a, a - S.Half, a], [], z)) == \
  365. hyperexpand(meijerg([a], [], [a - S.Half, a, a], [], z))
  366. assert can_do_meijer([a - 1], [], [a + 2, a - Rational(3, 2), a + 1], [])
  367. @XFAIL
  368. def test_meijerg_expand_fail():
  369. # These basically test hyper([], [1/2 - a, 1/2 + 1, 1/2], z),
  370. # which is *very* messy. But since the meijer g actually yields a
  371. # sum of bessel functions, things can sometimes be simplified a lot and
  372. # are then put into tables...
  373. assert can_do_meijer([], [], [a + S.Half], [a, a - b/2, a + b/2])
  374. assert can_do_meijer([], [], [0, S.Half], [a, -a])
  375. assert can_do_meijer([], [], [3*a - S.Half, a, -a - S.Half], [a - S.Half])
  376. assert can_do_meijer([], [], [0, a - S.Half, -a - S.Half], [S.Half])
  377. assert can_do_meijer([], [], [a, b + S.Half, b], [2*b - a])
  378. assert can_do_meijer([], [], [a, b + S.Half, b, 2*b - a])
  379. assert can_do_meijer([S.Half], [], [-a, a], [0])
  380. @slow
  381. def test_meijerg():
  382. # carefully set up the parameters.
  383. # NOTE: this used to fail sometimes. I believe it is fixed, but if you
  384. # hit an inexplicable test failure here, please let me know the seed.
  385. a1, a2 = (randcplx(n) - 5*I - n*I for n in range(2))
  386. b1, b2 = (randcplx(n) + 5*I + n*I for n in range(2))
  387. b3, b4, b5, a3, a4, a5 = (randcplx() for n in range(6))
  388. g = meijerg([a1], [a3, a4], [b1], [b3, b4], z)
  389. assert ReduceOrder.meijer_minus(3, 4) is None
  390. assert ReduceOrder.meijer_plus(4, 3) is None
  391. g2 = meijerg([a1, a2], [a3, a4], [b1], [b3, b4, a2], z)
  392. assert tn(ReduceOrder.meijer_plus(a2, a2).apply(g, op), g2, z)
  393. g2 = meijerg([a1, a2], [a3, a4], [b1], [b3, b4, a2 + 1], z)
  394. assert tn(ReduceOrder.meijer_plus(a2, a2 + 1).apply(g, op), g2, z)
  395. g2 = meijerg([a1, a2 - 1], [a3, a4], [b1], [b3, b4, a2 + 2], z)
  396. assert tn(ReduceOrder.meijer_plus(a2 - 1, a2 + 2).apply(g, op), g2, z)
  397. g2 = meijerg([a1], [a3, a4, b2 - 1], [b1, b2 + 2], [b3, b4], z)
  398. assert tn(ReduceOrder.meijer_minus(
  399. b2 + 2, b2 - 1).apply(g, op), g2, z, tol=1e-6)
  400. # test several-step reduction
  401. an = [a1, a2]
  402. bq = [b3, b4, a2 + 1]
  403. ap = [a3, a4, b2 - 1]
  404. bm = [b1, b2 + 1]
  405. niq, ops = reduce_order_meijer(G_Function(an, ap, bm, bq))
  406. assert niq.an == (a1,)
  407. assert set(niq.ap) == {a3, a4}
  408. assert niq.bm == (b1,)
  409. assert set(niq.bq) == {b3, b4}
  410. assert tn(apply_operators(g, ops, op), meijerg(an, ap, bm, bq, z), z)
  411. def test_meijerg_shift_operators():
  412. # carefully set up the parameters. XXX this still fails sometimes
  413. a1, a2, a3, a4, a5, b1, b2, b3, b4, b5 = (randcplx(n) for n in range(10))
  414. g = meijerg([a1], [a3, a4], [b1], [b3, b4], z)
  415. assert tn(MeijerShiftA(b1).apply(g, op),
  416. meijerg([a1], [a3, a4], [b1 + 1], [b3, b4], z), z)
  417. assert tn(MeijerShiftB(a1).apply(g, op),
  418. meijerg([a1 - 1], [a3, a4], [b1], [b3, b4], z), z)
  419. assert tn(MeijerShiftC(b3).apply(g, op),
  420. meijerg([a1], [a3, a4], [b1], [b3 + 1, b4], z), z)
  421. assert tn(MeijerShiftD(a3).apply(g, op),
  422. meijerg([a1], [a3 - 1, a4], [b1], [b3, b4], z), z)
  423. s = MeijerUnShiftA([a1], [a3, a4], [b1], [b3, b4], 0, z)
  424. assert tn(
  425. s.apply(g, op), meijerg([a1], [a3, a4], [b1 - 1], [b3, b4], z), z)
  426. s = MeijerUnShiftC([a1], [a3, a4], [b1], [b3, b4], 0, z)
  427. assert tn(
  428. s.apply(g, op), meijerg([a1], [a3, a4], [b1], [b3 - 1, b4], z), z)
  429. s = MeijerUnShiftB([a1], [a3, a4], [b1], [b3, b4], 0, z)
  430. assert tn(
  431. s.apply(g, op), meijerg([a1 + 1], [a3, a4], [b1], [b3, b4], z), z)
  432. s = MeijerUnShiftD([a1], [a3, a4], [b1], [b3, b4], 0, z)
  433. assert tn(
  434. s.apply(g, op), meijerg([a1], [a3 + 1, a4], [b1], [b3, b4], z), z)
  435. @slow
  436. def test_meijerg_confluence():
  437. def t(m, a, b):
  438. from sympy.core.sympify import sympify
  439. a, b = sympify([a, b])
  440. m_ = m
  441. m = hyperexpand(m)
  442. if not m == Piecewise((a, abs(z) < 1), (b, abs(1/z) < 1), (m_, True)):
  443. return False
  444. if not (m.args[0].args[0] == a and m.args[1].args[0] == b):
  445. return False
  446. z0 = randcplx()/10
  447. if abs(m.subs(z, z0).n() - a.subs(z, z0).n()).n() > 1e-10:
  448. return False
  449. if abs(m.subs(z, 1/z0).n() - b.subs(z, 1/z0).n()).n() > 1e-10:
  450. return False
  451. return True
  452. assert t(meijerg([], [1, 1], [0, 0], [], z), -log(z), 0)
  453. assert t(meijerg(
  454. [], [3, 1], [0, 0], [], z), -z**2/4 + z - log(z)/2 - Rational(3, 4), 0)
  455. assert t(meijerg([], [3, 1], [-1, 0], [], z),
  456. z**2/12 - z/2 + log(z)/2 + Rational(1, 4) + 1/(6*z), 0)
  457. assert t(meijerg([], [1, 1, 1, 1], [0, 0, 0, 0], [], z), -log(z)**3/6, 0)
  458. assert t(meijerg([1, 1], [], [], [0, 0], z), 0, -log(1/z))
  459. assert t(meijerg([1, 1], [2, 2], [1, 1], [0, 0], z),
  460. -z*log(z) + 2*z, -log(1/z) + 2)
  461. assert t(meijerg([S.Half], [1, 1], [0, 0], [Rational(3, 2)], z), log(z)/2 - 1, 0)
  462. def u(an, ap, bm, bq):
  463. m = meijerg(an, ap, bm, bq, z)
  464. m2 = hyperexpand(m, allow_hyper=True)
  465. if m2.has(meijerg) and not (m2.is_Piecewise and len(m2.args) == 3):
  466. return False
  467. return tn(m, m2, z)
  468. assert u([], [1], [0, 0], [])
  469. assert u([1, 1], [], [], [0])
  470. assert u([1, 1], [2, 2, 5], [1, 1, 6], [0, 0])
  471. assert u([1, 1], [2, 2, 5], [1, 1, 6], [0])
  472. def test_meijerg_with_Floats():
  473. # see issue #10681
  474. from sympy.polys.domains.realfield import RR
  475. f = meijerg(((3.0, 1), ()), ((Rational(3, 2),), (0,)), z)
  476. a = -2.3632718012073
  477. g = a*z**Rational(3, 2)*hyper((-0.5, Rational(3, 2)), (Rational(5, 2),), z*exp_polar(I*pi))
  478. assert RR.almosteq((hyperexpand(f)/g).n(), 1.0, 1e-12)
  479. def test_lerchphi():
  480. from sympy.functions.special.zeta_functions import (lerchphi, polylog)
  481. from sympy.simplify.gammasimp import gammasimp
  482. assert hyperexpand(hyper([1, a], [a + 1], z)/a) == lerchphi(z, 1, a)
  483. assert hyperexpand(
  484. hyper([1, a, a], [a + 1, a + 1], z)/a**2) == lerchphi(z, 2, a)
  485. assert hyperexpand(hyper([1, a, a, a], [a + 1, a + 1, a + 1], z)/a**3) == \
  486. lerchphi(z, 3, a)
  487. assert hyperexpand(hyper([1] + [a]*10, [a + 1]*10, z)/a**10) == \
  488. lerchphi(z, 10, a)
  489. assert gammasimp(hyperexpand(meijerg([0, 1 - a], [], [0],
  490. [-a], exp_polar(-I*pi)*z))) == lerchphi(z, 1, a)
  491. assert gammasimp(hyperexpand(meijerg([0, 1 - a, 1 - a], [], [0],
  492. [-a, -a], exp_polar(-I*pi)*z))) == lerchphi(z, 2, a)
  493. assert gammasimp(hyperexpand(meijerg([0, 1 - a, 1 - a, 1 - a], [], [0],
  494. [-a, -a, -a], exp_polar(-I*pi)*z))) == lerchphi(z, 3, a)
  495. assert hyperexpand(z*hyper([1, 1], [2], z)) == -log(1 + -z)
  496. assert hyperexpand(z*hyper([1, 1, 1], [2, 2], z)) == polylog(2, z)
  497. assert hyperexpand(z*hyper([1, 1, 1, 1], [2, 2, 2], z)) == polylog(3, z)
  498. assert hyperexpand(hyper([1, a, 1 + S.Half], [a + 1, S.Half], z)) == \
  499. -2*a/(z - 1) + (-2*a**2 + a)*lerchphi(z, 1, a)
  500. # Now numerical tests. These make sure reductions etc are carried out
  501. # correctly
  502. # a rational function (polylog at negative integer order)
  503. assert can_do([2, 2, 2], [1, 1])
  504. # NOTE these contain log(1-x) etc ... better make sure we have |z| < 1
  505. # reduction of order for polylog
  506. assert can_do([1, 1, 1, b + 5], [2, 2, b], div=10)
  507. # reduction of order for lerchphi
  508. # XXX lerchphi in mpmath is flaky
  509. assert can_do(
  510. [1, a, a, a, b + 5], [a + 1, a + 1, a + 1, b], numerical=False)
  511. # test a bug
  512. from sympy.functions.elementary.complexes import Abs
  513. assert hyperexpand(hyper([S.Half, S.Half, S.Half, 1],
  514. [Rational(3, 2), Rational(3, 2), Rational(3, 2)], Rational(1, 4))) == \
  515. Abs(-polylog(3, exp_polar(I*pi)/2) + polylog(3, S.Half))
  516. def test_partial_simp():
  517. # First test that hypergeometric function formulae work.
  518. a, b, c, d, e = (randcplx() for _ in range(5))
  519. for func in [Hyper_Function([a, b, c], [d, e]),
  520. Hyper_Function([], [a, b, c, d, e])]:
  521. f = build_hypergeometric_formula(func)
  522. z = f.z
  523. assert f.closed_form == func(z)
  524. deriv1 = f.B.diff(z)*z
  525. deriv2 = f.M*f.B
  526. for func1, func2 in zip(deriv1, deriv2):
  527. assert tn(func1, func2, z)
  528. # Now test that formulae are partially simplified.
  529. a, b, z = symbols('a b z')
  530. assert hyperexpand(hyper([3, a], [1, b], z)) == \
  531. (-a*b/2 + a*z/2 + 2*a)*hyper([a + 1], [b], z) \
  532. + (a*b/2 - 2*a + 1)*hyper([a], [b], z)
  533. assert tn(
  534. hyperexpand(hyper([3, d], [1, e], z)), hyper([3, d], [1, e], z), z)
  535. assert hyperexpand(hyper([3], [1, a, b], z)) == \
  536. hyper((), (a, b), z) \
  537. + z*hyper((), (a + 1, b), z)/(2*a) \
  538. - z*(b - 4)*hyper((), (a + 1, b + 1), z)/(2*a*b)
  539. assert tn(
  540. hyperexpand(hyper([3], [1, d, e], z)), hyper([3], [1, d, e], z), z)
  541. def test_hyperexpand_special():
  542. assert hyperexpand(hyper([a, b], [c], 1)) == \
  543. gamma(c)*gamma(c - a - b)/gamma(c - a)/gamma(c - b)
  544. assert hyperexpand(hyper([a, b], [1 + a - b], -1)) == \
  545. gamma(1 + a/2)*gamma(1 + a - b)/gamma(1 + a)/gamma(1 + a/2 - b)
  546. assert hyperexpand(hyper([a, b], [1 + b - a], -1)) == \
  547. gamma(1 + b/2)*gamma(1 + b - a)/gamma(1 + b)/gamma(1 + b/2 - a)
  548. assert hyperexpand(meijerg([1 - z - a/2], [1 - z + a/2], [b/2], [-b/2], 1)) == \
  549. gamma(1 - 2*z)*gamma(z + a/2 + b/2)/gamma(1 - z + a/2 - b/2) \
  550. /gamma(1 - z - a/2 + b/2)/gamma(1 - z + a/2 + b/2)
  551. assert hyperexpand(hyper([a], [b], 0)) == 1
  552. assert hyper([a], [b], 0) != 0
  553. def test_Mod1_behavior():
  554. from sympy.core.symbol import Symbol
  555. from sympy.simplify.simplify import simplify
  556. n = Symbol('n', integer=True)
  557. # Note: this should not hang.
  558. assert simplify(hyperexpand(meijerg([1], [], [n + 1], [0], z))) == \
  559. lowergamma(n + 1, z)
  560. @slow
  561. def test_prudnikov_misc():
  562. assert can_do([1, (3 + I)/2, (3 - I)/2], [Rational(3, 2), 2])
  563. assert can_do([S.Half, a - 1], [Rational(3, 2), a + 1], lowerplane=True)
  564. assert can_do([], [b + 1])
  565. assert can_do([a], [a - 1, b + 1])
  566. assert can_do([a], [a - S.Half, 2*a])
  567. assert can_do([a], [a - S.Half, 2*a + 1])
  568. assert can_do([a], [a - S.Half, 2*a - 1])
  569. assert can_do([a], [a + S.Half, 2*a])
  570. assert can_do([a], [a + S.Half, 2*a + 1])
  571. assert can_do([a], [a + S.Half, 2*a - 1])
  572. assert can_do([S.Half], [b, 2 - b])
  573. assert can_do([S.Half], [b, 3 - b])
  574. assert can_do([1], [2, b])
  575. assert can_do([a, a + S.Half], [2*a, b, 2*a - b + 1])
  576. assert can_do([a, a + S.Half], [S.Half, 2*a, 2*a + S.Half])
  577. assert can_do([a], [a + 1], lowerplane=True) # lowergamma
  578. def test_prudnikov_1():
  579. # A. P. Prudnikov, Yu. A. Brychkov and O. I. Marichev (1990).
  580. # Integrals and Series: More Special Functions, Vol. 3,.
  581. # Gordon and Breach Science Publisher
  582. # 7.3.1
  583. assert can_do([a, -a], [S.Half])
  584. assert can_do([a, 1 - a], [S.Half])
  585. assert can_do([a, 1 - a], [Rational(3, 2)])
  586. assert can_do([a, 2 - a], [S.Half])
  587. assert can_do([a, 2 - a], [Rational(3, 2)])
  588. assert can_do([a, 2 - a], [Rational(3, 2)])
  589. assert can_do([a, a + S.Half], [2*a - 1])
  590. assert can_do([a, a + S.Half], [2*a])
  591. assert can_do([a, a + S.Half], [2*a + 1])
  592. assert can_do([a, a + S.Half], [S.Half])
  593. assert can_do([a, a + S.Half], [Rational(3, 2)])
  594. assert can_do([a, a/2 + 1], [a/2])
  595. assert can_do([1, b], [2])
  596. assert can_do([1, b], [b + 1], numerical=False) # Lerch Phi
  597. # NOTE: branches are complicated for |z| > 1
  598. assert can_do([a], [2*a])
  599. assert can_do([a], [2*a + 1])
  600. assert can_do([a], [2*a - 1])
  601. @slow
  602. def test_prudnikov_2():
  603. h = S.Half
  604. assert can_do([-h, -h], [h])
  605. assert can_do([-h, h], [3*h])
  606. assert can_do([-h, h], [5*h])
  607. assert can_do([-h, h], [7*h])
  608. assert can_do([-h, 1], [h])
  609. for p in [-h, h]:
  610. for n in [-h, h, 1, 3*h, 2, 5*h, 3, 7*h, 4]:
  611. for m in [-h, h, 3*h, 5*h, 7*h]:
  612. assert can_do([p, n], [m])
  613. for n in [1, 2, 3, 4]:
  614. for m in [1, 2, 3, 4]:
  615. assert can_do([p, n], [m])
  616. @slow
  617. def test_prudnikov_3():
  618. if ON_CI:
  619. # See https://github.com/sympy/sympy/pull/12795
  620. skip("Too slow for CI.")
  621. h = S.Half
  622. assert can_do([Rational(1, 4), Rational(3, 4)], [h])
  623. assert can_do([Rational(1, 4), Rational(3, 4)], [3*h])
  624. assert can_do([Rational(1, 3), Rational(2, 3)], [3*h])
  625. assert can_do([Rational(3, 4), Rational(5, 4)], [h])
  626. assert can_do([Rational(3, 4), Rational(5, 4)], [3*h])
  627. for p in [1, 2, 3, 4]:
  628. for n in [-h, h, 1, 3*h, 2, 5*h, 3, 7*h, 4, 9*h]:
  629. for m in [1, 3*h, 2, 5*h, 3, 7*h, 4]:
  630. assert can_do([p, m], [n])
  631. @slow
  632. def test_prudnikov_4():
  633. h = S.Half
  634. for p in [3*h, 5*h, 7*h]:
  635. for n in [-h, h, 3*h, 5*h, 7*h]:
  636. for m in [3*h, 2, 5*h, 3, 7*h, 4]:
  637. assert can_do([p, m], [n])
  638. for n in [1, 2, 3, 4]:
  639. for m in [2, 3, 4]:
  640. assert can_do([p, m], [n])
  641. @slow
  642. def test_prudnikov_5():
  643. h = S.Half
  644. for p in [1, 2, 3]:
  645. for q in range(p, 4):
  646. for r in [1, 2, 3]:
  647. for s in range(r, 4):
  648. assert can_do([-h, p, q], [r, s])
  649. for p in [h, 1, 3*h, 2, 5*h, 3]:
  650. for q in [h, 3*h, 5*h]:
  651. for r in [h, 3*h, 5*h]:
  652. for s in [h, 3*h, 5*h]:
  653. if s <= q and s <= r:
  654. assert can_do([-h, p, q], [r, s])
  655. for p in [h, 1, 3*h, 2, 5*h, 3]:
  656. for q in [1, 2, 3]:
  657. for r in [h, 3*h, 5*h]:
  658. for s in [1, 2, 3]:
  659. assert can_do([-h, p, q], [r, s])
  660. @slow
  661. def test_prudnikov_6():
  662. h = S.Half
  663. for m in [3*h, 5*h]:
  664. for n in [1, 2, 3]:
  665. for q in [h, 1, 2]:
  666. for p in [1, 2, 3]:
  667. assert can_do([h, q, p], [m, n])
  668. for q in [1, 2, 3]:
  669. for p in [3*h, 5*h]:
  670. assert can_do([h, q, p], [m, n])
  671. for q in [1, 2]:
  672. for p in [1, 2, 3]:
  673. for m in [1, 2, 3]:
  674. for n in [1, 2, 3]:
  675. assert can_do([h, q, p], [m, n])
  676. assert can_do([h, h, 5*h], [3*h, 3*h])
  677. assert can_do([h, 1, 5*h], [3*h, 3*h])
  678. assert can_do([h, 2, 2], [1, 3])
  679. # pages 435 to 457 contain more PFDD and stuff like this
  680. @slow
  681. def test_prudnikov_7():
  682. assert can_do([3], [6])
  683. h = S.Half
  684. for n in [h, 3*h, 5*h, 7*h]:
  685. assert can_do([-h], [n])
  686. for m in [-h, h, 1, 3*h, 2, 5*h, 3, 7*h, 4]: # HERE
  687. for n in [-h, h, 3*h, 5*h, 7*h, 1, 2, 3, 4]:
  688. assert can_do([m], [n])
  689. @slow
  690. def test_prudnikov_8():
  691. h = S.Half
  692. # 7.12.2
  693. for ai in [1, 2, 3]:
  694. for bi in [1, 2, 3]:
  695. for ci in range(1, ai + 1):
  696. for di in [h, 1, 3*h, 2, 5*h, 3]:
  697. assert can_do([ai, bi], [ci, di])
  698. for bi in [3*h, 5*h]:
  699. for ci in [h, 1, 3*h, 2, 5*h, 3]:
  700. for di in [1, 2, 3]:
  701. assert can_do([ai, bi], [ci, di])
  702. for ai in [-h, h, 3*h, 5*h]:
  703. for bi in [1, 2, 3]:
  704. for ci in [h, 1, 3*h, 2, 5*h, 3]:
  705. for di in [1, 2, 3]:
  706. assert can_do([ai, bi], [ci, di])
  707. for bi in [h, 3*h, 5*h]:
  708. for ci in [h, 3*h, 5*h, 3]:
  709. for di in [h, 1, 3*h, 2, 5*h, 3]:
  710. if ci <= bi:
  711. assert can_do([ai, bi], [ci, di])
  712. def test_prudnikov_9():
  713. # 7.13.1 [we have a general formula ... so this is a bit pointless]
  714. for i in range(9):
  715. assert can_do([], [(S(i) + 1)/2])
  716. for i in range(5):
  717. assert can_do([], [-(2*S(i) + 1)/2])
  718. @slow
  719. def test_prudnikov_10():
  720. # 7.14.2
  721. h = S.Half
  722. for p in [-h, h, 1, 3*h, 2, 5*h, 3, 7*h, 4]:
  723. for m in [1, 2, 3, 4]:
  724. for n in range(m, 5):
  725. assert can_do([p], [m, n])
  726. for p in [1, 2, 3, 4]:
  727. for n in [h, 3*h, 5*h, 7*h]:
  728. for m in [1, 2, 3, 4]:
  729. assert can_do([p], [n, m])
  730. for p in [3*h, 5*h, 7*h]:
  731. for m in [h, 1, 2, 5*h, 3, 7*h, 4]:
  732. assert can_do([p], [h, m])
  733. assert can_do([p], [3*h, m])
  734. for m in [h, 1, 2, 5*h, 3, 7*h, 4]:
  735. assert can_do([7*h], [5*h, m])
  736. assert can_do([Rational(-1, 2)], [S.Half, S.Half]) # shine-integral shi
  737. def test_prudnikov_11():
  738. # 7.15
  739. assert can_do([a, a + S.Half], [2*a, b, 2*a - b])
  740. assert can_do([a, a + S.Half], [Rational(3, 2), 2*a, 2*a - S.Half])
  741. assert can_do([Rational(1, 4), Rational(3, 4)], [S.Half, S.Half, 1])
  742. assert can_do([Rational(5, 4), Rational(3, 4)], [Rational(3, 2), S.Half, 2])
  743. assert can_do([Rational(5, 4), Rational(3, 4)], [Rational(3, 2), Rational(3, 2), 1])
  744. assert can_do([Rational(5, 4), Rational(7, 4)], [Rational(3, 2), Rational(5, 2), 2])
  745. assert can_do([1, 1], [Rational(3, 2), 2, 2]) # cosh-integral chi
  746. def test_prudnikov_12():
  747. # 7.16
  748. assert can_do(
  749. [], [a, a + S.Half, 2*a], False) # branches only agree for some z!
  750. assert can_do([], [a, a + S.Half, 2*a + 1], False) # dito
  751. assert can_do([], [S.Half, a, a + S.Half])
  752. assert can_do([], [Rational(3, 2), a, a + S.Half])
  753. assert can_do([], [Rational(1, 4), S.Half, Rational(3, 4)])
  754. assert can_do([], [S.Half, S.Half, 1])
  755. assert can_do([], [S.Half, Rational(3, 2), 1])
  756. assert can_do([], [Rational(3, 4), Rational(3, 2), Rational(5, 4)])
  757. assert can_do([], [1, 1, Rational(3, 2)])
  758. assert can_do([], [1, 2, Rational(3, 2)])
  759. assert can_do([], [1, Rational(3, 2), Rational(3, 2)])
  760. assert can_do([], [Rational(5, 4), Rational(3, 2), Rational(7, 4)])
  761. assert can_do([], [2, Rational(3, 2), Rational(3, 2)])
  762. @slow
  763. def test_prudnikov_2F1():
  764. h = S.Half
  765. # Elliptic integrals
  766. for p in [-h, h]:
  767. for m in [h, 3*h, 5*h, 7*h]:
  768. for n in [1, 2, 3, 4]:
  769. assert can_do([p, m], [n])
  770. @XFAIL
  771. def test_prudnikov_fail_2F1():
  772. assert can_do([a, b], [b + 1]) # incomplete beta function
  773. assert can_do([-1, b], [c]) # Poly. also -2, -3 etc
  774. # TODO polys
  775. # Legendre functions:
  776. assert can_do([a, b], [a + b + S.Half])
  777. assert can_do([a, b], [a + b - S.Half])
  778. assert can_do([a, b], [a + b + Rational(3, 2)])
  779. assert can_do([a, b], [(a + b + 1)/2])
  780. assert can_do([a, b], [(a + b)/2 + 1])
  781. assert can_do([a, b], [a - b + 1])
  782. assert can_do([a, b], [a - b + 2])
  783. assert can_do([a, b], [2*b])
  784. assert can_do([a, b], [S.Half])
  785. assert can_do([a, b], [Rational(3, 2)])
  786. assert can_do([a, 1 - a], [c])
  787. assert can_do([a, 2 - a], [c])
  788. assert can_do([a, 3 - a], [c])
  789. assert can_do([a, a + S.Half], [c])
  790. assert can_do([1, b], [c])
  791. assert can_do([1, b], [Rational(3, 2)])
  792. assert can_do([Rational(1, 4), Rational(3, 4)], [1])
  793. # PFDD
  794. o = S.One
  795. assert can_do([o/8, 1], [o/8*9])
  796. assert can_do([o/6, 1], [o/6*7])
  797. assert can_do([o/6, 1], [o/6*13])
  798. assert can_do([o/5, 1], [o/5*6])
  799. assert can_do([o/5, 1], [o/5*11])
  800. assert can_do([o/4, 1], [o/4*5])
  801. assert can_do([o/4, 1], [o/4*9])
  802. assert can_do([o/3, 1], [o/3*4])
  803. assert can_do([o/3, 1], [o/3*7])
  804. assert can_do([o/8*3, 1], [o/8*11])
  805. assert can_do([o/5*2, 1], [o/5*7])
  806. assert can_do([o/5*2, 1], [o/5*12])
  807. assert can_do([o/5*3, 1], [o/5*8])
  808. assert can_do([o/5*3, 1], [o/5*13])
  809. assert can_do([o/8*5, 1], [o/8*13])
  810. assert can_do([o/4*3, 1], [o/4*7])
  811. assert can_do([o/4*3, 1], [o/4*11])
  812. assert can_do([o/3*2, 1], [o/3*5])
  813. assert can_do([o/3*2, 1], [o/3*8])
  814. assert can_do([o/5*4, 1], [o/5*9])
  815. assert can_do([o/5*4, 1], [o/5*14])
  816. assert can_do([o/6*5, 1], [o/6*11])
  817. assert can_do([o/6*5, 1], [o/6*17])
  818. assert can_do([o/8*7, 1], [o/8*15])
  819. @XFAIL
  820. def test_prudnikov_fail_3F2():
  821. assert can_do([a, a + Rational(1, 3), a + Rational(2, 3)], [Rational(1, 3), Rational(2, 3)])
  822. assert can_do([a, a + Rational(1, 3), a + Rational(2, 3)], [Rational(2, 3), Rational(4, 3)])
  823. assert can_do([a, a + Rational(1, 3), a + Rational(2, 3)], [Rational(4, 3), Rational(5, 3)])
  824. # page 421
  825. assert can_do([a, a + Rational(1, 3), a + Rational(2, 3)], [a*Rational(3, 2), (3*a + 1)/2])
  826. # pages 422 ...
  827. assert can_do([Rational(-1, 2), S.Half, S.Half], [1, 1]) # elliptic integrals
  828. assert can_do([Rational(-1, 2), S.Half, 1], [Rational(3, 2), Rational(3, 2)])
  829. # TODO LOTS more
  830. # PFDD
  831. assert can_do([Rational(1, 8), Rational(3, 8), 1], [Rational(9, 8), Rational(11, 8)])
  832. assert can_do([Rational(1, 8), Rational(5, 8), 1], [Rational(9, 8), Rational(13, 8)])
  833. assert can_do([Rational(1, 8), Rational(7, 8), 1], [Rational(9, 8), Rational(15, 8)])
  834. assert can_do([Rational(1, 6), Rational(1, 3), 1], [Rational(7, 6), Rational(4, 3)])
  835. assert can_do([Rational(1, 6), Rational(2, 3), 1], [Rational(7, 6), Rational(5, 3)])
  836. assert can_do([Rational(1, 6), Rational(2, 3), 1], [Rational(5, 3), Rational(13, 6)])
  837. assert can_do([S.Half, 1, 1], [Rational(1, 4), Rational(3, 4)])
  838. # LOTS more
  839. @XFAIL
  840. def test_prudnikov_fail_other():
  841. # 7.11.2
  842. # 7.12.1
  843. assert can_do([1, a], [b, 1 - 2*a + b]) # ???
  844. # 7.14.2
  845. assert can_do([Rational(-1, 2)], [S.Half, 1]) # struve
  846. assert can_do([1], [S.Half, S.Half]) # struve
  847. assert can_do([Rational(1, 4)], [S.Half, Rational(5, 4)]) # PFDD
  848. assert can_do([Rational(3, 4)], [Rational(3, 2), Rational(7, 4)]) # PFDD
  849. assert can_do([1], [Rational(1, 4), Rational(3, 4)]) # PFDD
  850. assert can_do([1], [Rational(3, 4), Rational(5, 4)]) # PFDD
  851. assert can_do([1], [Rational(5, 4), Rational(7, 4)]) # PFDD
  852. # TODO LOTS more
  853. # 7.15.2
  854. assert can_do([S.Half, 1], [Rational(3, 4), Rational(5, 4), Rational(3, 2)]) # PFDD
  855. assert can_do([S.Half, 1], [Rational(7, 4), Rational(5, 4), Rational(3, 2)]) # PFDD
  856. # 7.16.1
  857. assert can_do([], [Rational(1, 3), S(2/3)]) # PFDD
  858. assert can_do([], [Rational(2, 3), S(4/3)]) # PFDD
  859. assert can_do([], [Rational(5, 3), S(4/3)]) # PFDD
  860. # XXX this does not *evaluate* right??
  861. assert can_do([], [a, a + S.Half, 2*a - 1])
  862. def test_bug():
  863. h = hyper([-1, 1], [z], -1)
  864. assert hyperexpand(h) == (z + 1)/z
  865. def test_omgissue_203():
  866. h = hyper((-5, -3, -4), (-6, -6), 1)
  867. assert hyperexpand(h) == Rational(1, 30)
  868. h = hyper((-6, -7, -5), (-6, -6), 1)
  869. assert hyperexpand(h) == Rational(-1, 6)