test_hyperexpand.py 40 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063
  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, tooslow
  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. def test_prudnikov_3():
  617. h = S.Half
  618. assert can_do([Rational(1, 4), Rational(3, 4)], [h])
  619. assert can_do([Rational(1, 4), Rational(3, 4)], [3*h])
  620. assert can_do([Rational(1, 3), Rational(2, 3)], [3*h])
  621. assert can_do([Rational(3, 4), Rational(5, 4)], [h])
  622. assert can_do([Rational(3, 4), Rational(5, 4)], [3*h])
  623. @tooslow
  624. def test_prudnikov_3_slow():
  625. # XXX: This is marked as tooslow and hence skipped in CI. None of the
  626. # individual cases below fails or hangs. Some cases are slow and the loops
  627. # below generate 280 different cases. Is it really necessary to test all
  628. # 280 cases here?
  629. h = S.Half
  630. for p in [1, 2, 3, 4]:
  631. for n in [-h, h, 1, 3*h, 2, 5*h, 3, 7*h, 4, 9*h]:
  632. for m in [1, 3*h, 2, 5*h, 3, 7*h, 4]:
  633. assert can_do([p, m], [n])
  634. @slow
  635. def test_prudnikov_4():
  636. h = S.Half
  637. for p in [3*h, 5*h, 7*h]:
  638. for n in [-h, h, 3*h, 5*h, 7*h]:
  639. for m in [3*h, 2, 5*h, 3, 7*h, 4]:
  640. assert can_do([p, m], [n])
  641. for n in [1, 2, 3, 4]:
  642. for m in [2, 3, 4]:
  643. assert can_do([p, m], [n])
  644. @slow
  645. def test_prudnikov_5():
  646. h = S.Half
  647. for p in [1, 2, 3]:
  648. for q in range(p, 4):
  649. for r in [1, 2, 3]:
  650. for s in range(r, 4):
  651. assert can_do([-h, p, q], [r, s])
  652. for p in [h, 1, 3*h, 2, 5*h, 3]:
  653. for q in [h, 3*h, 5*h]:
  654. for r in [h, 3*h, 5*h]:
  655. for s in [h, 3*h, 5*h]:
  656. if s <= q and s <= r:
  657. assert can_do([-h, p, q], [r, s])
  658. for p in [h, 1, 3*h, 2, 5*h, 3]:
  659. for q in [1, 2, 3]:
  660. for r in [h, 3*h, 5*h]:
  661. for s in [1, 2, 3]:
  662. assert can_do([-h, p, q], [r, s])
  663. @slow
  664. def test_prudnikov_6():
  665. h = S.Half
  666. for m in [3*h, 5*h]:
  667. for n in [1, 2, 3]:
  668. for q in [h, 1, 2]:
  669. for p in [1, 2, 3]:
  670. assert can_do([h, q, p], [m, n])
  671. for q in [1, 2, 3]:
  672. for p in [3*h, 5*h]:
  673. assert can_do([h, q, p], [m, n])
  674. for q in [1, 2]:
  675. for p in [1, 2, 3]:
  676. for m in [1, 2, 3]:
  677. for n in [1, 2, 3]:
  678. assert can_do([h, q, p], [m, n])
  679. assert can_do([h, h, 5*h], [3*h, 3*h])
  680. assert can_do([h, 1, 5*h], [3*h, 3*h])
  681. assert can_do([h, 2, 2], [1, 3])
  682. # pages 435 to 457 contain more PFDD and stuff like this
  683. @slow
  684. def test_prudnikov_7():
  685. assert can_do([3], [6])
  686. h = S.Half
  687. for n in [h, 3*h, 5*h, 7*h]:
  688. assert can_do([-h], [n])
  689. for m in [-h, h, 1, 3*h, 2, 5*h, 3, 7*h, 4]: # HERE
  690. for n in [-h, h, 3*h, 5*h, 7*h, 1, 2, 3, 4]:
  691. assert can_do([m], [n])
  692. @slow
  693. def test_prudnikov_8():
  694. h = S.Half
  695. # 7.12.2
  696. for ai in [1, 2, 3]:
  697. for bi in [1, 2, 3]:
  698. for ci in range(1, ai + 1):
  699. for di in [h, 1, 3*h, 2, 5*h, 3]:
  700. assert can_do([ai, bi], [ci, di])
  701. for bi in [3*h, 5*h]:
  702. for ci in [h, 1, 3*h, 2, 5*h, 3]:
  703. for di in [1, 2, 3]:
  704. assert can_do([ai, bi], [ci, di])
  705. for ai in [-h, h, 3*h, 5*h]:
  706. for bi in [1, 2, 3]:
  707. for ci in [h, 1, 3*h, 2, 5*h, 3]:
  708. for di in [1, 2, 3]:
  709. assert can_do([ai, bi], [ci, di])
  710. for bi in [h, 3*h, 5*h]:
  711. for ci in [h, 3*h, 5*h, 3]:
  712. for di in [h, 1, 3*h, 2, 5*h, 3]:
  713. if ci <= bi:
  714. assert can_do([ai, bi], [ci, di])
  715. def test_prudnikov_9():
  716. # 7.13.1 [we have a general formula ... so this is a bit pointless]
  717. for i in range(9):
  718. assert can_do([], [(S(i) + 1)/2])
  719. for i in range(5):
  720. assert can_do([], [-(2*S(i) + 1)/2])
  721. @slow
  722. def test_prudnikov_10():
  723. # 7.14.2
  724. h = S.Half
  725. for p in [-h, h, 1, 3*h, 2, 5*h, 3, 7*h, 4]:
  726. for m in [1, 2, 3, 4]:
  727. for n in range(m, 5):
  728. assert can_do([p], [m, n])
  729. for p in [1, 2, 3, 4]:
  730. for n in [h, 3*h, 5*h, 7*h]:
  731. for m in [1, 2, 3, 4]:
  732. assert can_do([p], [n, m])
  733. for p in [3*h, 5*h, 7*h]:
  734. for m in [h, 1, 2, 5*h, 3, 7*h, 4]:
  735. assert can_do([p], [h, m])
  736. assert can_do([p], [3*h, m])
  737. for m in [h, 1, 2, 5*h, 3, 7*h, 4]:
  738. assert can_do([7*h], [5*h, m])
  739. assert can_do([Rational(-1, 2)], [S.Half, S.Half]) # shine-integral shi
  740. def test_prudnikov_11():
  741. # 7.15
  742. assert can_do([a, a + S.Half], [2*a, b, 2*a - b])
  743. assert can_do([a, a + S.Half], [Rational(3, 2), 2*a, 2*a - S.Half])
  744. assert can_do([Rational(1, 4), Rational(3, 4)], [S.Half, S.Half, 1])
  745. assert can_do([Rational(5, 4), Rational(3, 4)], [Rational(3, 2), S.Half, 2])
  746. assert can_do([Rational(5, 4), Rational(3, 4)], [Rational(3, 2), Rational(3, 2), 1])
  747. assert can_do([Rational(5, 4), Rational(7, 4)], [Rational(3, 2), Rational(5, 2), 2])
  748. assert can_do([1, 1], [Rational(3, 2), 2, 2]) # cosh-integral chi
  749. def test_prudnikov_12():
  750. # 7.16
  751. assert can_do(
  752. [], [a, a + S.Half, 2*a], False) # branches only agree for some z!
  753. assert can_do([], [a, a + S.Half, 2*a + 1], False) # dito
  754. assert can_do([], [S.Half, a, a + S.Half])
  755. assert can_do([], [Rational(3, 2), a, a + S.Half])
  756. assert can_do([], [Rational(1, 4), S.Half, Rational(3, 4)])
  757. assert can_do([], [S.Half, S.Half, 1])
  758. assert can_do([], [S.Half, Rational(3, 2), 1])
  759. assert can_do([], [Rational(3, 4), Rational(3, 2), Rational(5, 4)])
  760. assert can_do([], [1, 1, Rational(3, 2)])
  761. assert can_do([], [1, 2, Rational(3, 2)])
  762. assert can_do([], [1, Rational(3, 2), Rational(3, 2)])
  763. assert can_do([], [Rational(5, 4), Rational(3, 2), Rational(7, 4)])
  764. assert can_do([], [2, Rational(3, 2), Rational(3, 2)])
  765. @slow
  766. def test_prudnikov_2F1():
  767. h = S.Half
  768. # Elliptic integrals
  769. for p in [-h, h]:
  770. for m in [h, 3*h, 5*h, 7*h]:
  771. for n in [1, 2, 3, 4]:
  772. assert can_do([p, m], [n])
  773. @XFAIL
  774. def test_prudnikov_fail_2F1():
  775. assert can_do([a, b], [b + 1]) # incomplete beta function
  776. assert can_do([-1, b], [c]) # Poly. also -2, -3 etc
  777. # TODO polys
  778. # Legendre functions:
  779. assert can_do([a, b], [a + b + S.Half])
  780. assert can_do([a, b], [a + b - S.Half])
  781. assert can_do([a, b], [a + b + Rational(3, 2)])
  782. assert can_do([a, b], [(a + b + 1)/2])
  783. assert can_do([a, b], [(a + b)/2 + 1])
  784. assert can_do([a, b], [a - b + 1])
  785. assert can_do([a, b], [a - b + 2])
  786. assert can_do([a, b], [2*b])
  787. assert can_do([a, b], [S.Half])
  788. assert can_do([a, b], [Rational(3, 2)])
  789. assert can_do([a, 1 - a], [c])
  790. assert can_do([a, 2 - a], [c])
  791. assert can_do([a, 3 - a], [c])
  792. assert can_do([a, a + S.Half], [c])
  793. assert can_do([1, b], [c])
  794. assert can_do([1, b], [Rational(3, 2)])
  795. assert can_do([Rational(1, 4), Rational(3, 4)], [1])
  796. # PFDD
  797. o = S.One
  798. assert can_do([o/8, 1], [o/8*9])
  799. assert can_do([o/6, 1], [o/6*7])
  800. assert can_do([o/6, 1], [o/6*13])
  801. assert can_do([o/5, 1], [o/5*6])
  802. assert can_do([o/5, 1], [o/5*11])
  803. assert can_do([o/4, 1], [o/4*5])
  804. assert can_do([o/4, 1], [o/4*9])
  805. assert can_do([o/3, 1], [o/3*4])
  806. assert can_do([o/3, 1], [o/3*7])
  807. assert can_do([o/8*3, 1], [o/8*11])
  808. assert can_do([o/5*2, 1], [o/5*7])
  809. assert can_do([o/5*2, 1], [o/5*12])
  810. assert can_do([o/5*3, 1], [o/5*8])
  811. assert can_do([o/5*3, 1], [o/5*13])
  812. assert can_do([o/8*5, 1], [o/8*13])
  813. assert can_do([o/4*3, 1], [o/4*7])
  814. assert can_do([o/4*3, 1], [o/4*11])
  815. assert can_do([o/3*2, 1], [o/3*5])
  816. assert can_do([o/3*2, 1], [o/3*8])
  817. assert can_do([o/5*4, 1], [o/5*9])
  818. assert can_do([o/5*4, 1], [o/5*14])
  819. assert can_do([o/6*5, 1], [o/6*11])
  820. assert can_do([o/6*5, 1], [o/6*17])
  821. assert can_do([o/8*7, 1], [o/8*15])
  822. @XFAIL
  823. def test_prudnikov_fail_3F2():
  824. assert can_do([a, a + Rational(1, 3), a + Rational(2, 3)], [Rational(1, 3), Rational(2, 3)])
  825. assert can_do([a, a + Rational(1, 3), a + Rational(2, 3)], [Rational(2, 3), Rational(4, 3)])
  826. assert can_do([a, a + Rational(1, 3), a + Rational(2, 3)], [Rational(4, 3), Rational(5, 3)])
  827. # page 421
  828. assert can_do([a, a + Rational(1, 3), a + Rational(2, 3)], [a*Rational(3, 2), (3*a + 1)/2])
  829. # pages 422 ...
  830. assert can_do([Rational(-1, 2), S.Half, S.Half], [1, 1]) # elliptic integrals
  831. assert can_do([Rational(-1, 2), S.Half, 1], [Rational(3, 2), Rational(3, 2)])
  832. # TODO LOTS more
  833. # PFDD
  834. assert can_do([Rational(1, 8), Rational(3, 8), 1], [Rational(9, 8), Rational(11, 8)])
  835. assert can_do([Rational(1, 8), Rational(5, 8), 1], [Rational(9, 8), Rational(13, 8)])
  836. assert can_do([Rational(1, 8), Rational(7, 8), 1], [Rational(9, 8), Rational(15, 8)])
  837. assert can_do([Rational(1, 6), Rational(1, 3), 1], [Rational(7, 6), Rational(4, 3)])
  838. assert can_do([Rational(1, 6), Rational(2, 3), 1], [Rational(7, 6), Rational(5, 3)])
  839. assert can_do([Rational(1, 6), Rational(2, 3), 1], [Rational(5, 3), Rational(13, 6)])
  840. assert can_do([S.Half, 1, 1], [Rational(1, 4), Rational(3, 4)])
  841. # LOTS more
  842. @XFAIL
  843. def test_prudnikov_fail_other():
  844. # 7.11.2
  845. # 7.12.1
  846. assert can_do([1, a], [b, 1 - 2*a + b]) # ???
  847. # 7.14.2
  848. assert can_do([Rational(-1, 2)], [S.Half, 1]) # struve
  849. assert can_do([1], [S.Half, S.Half]) # struve
  850. assert can_do([Rational(1, 4)], [S.Half, Rational(5, 4)]) # PFDD
  851. assert can_do([Rational(3, 4)], [Rational(3, 2), Rational(7, 4)]) # PFDD
  852. assert can_do([1], [Rational(1, 4), Rational(3, 4)]) # PFDD
  853. assert can_do([1], [Rational(3, 4), Rational(5, 4)]) # PFDD
  854. assert can_do([1], [Rational(5, 4), Rational(7, 4)]) # PFDD
  855. # TODO LOTS more
  856. # 7.15.2
  857. assert can_do([S.Half, 1], [Rational(3, 4), Rational(5, 4), Rational(3, 2)]) # PFDD
  858. assert can_do([S.Half, 1], [Rational(7, 4), Rational(5, 4), Rational(3, 2)]) # PFDD
  859. # 7.16.1
  860. assert can_do([], [Rational(1, 3), S(2/3)]) # PFDD
  861. assert can_do([], [Rational(2, 3), S(4/3)]) # PFDD
  862. assert can_do([], [Rational(5, 3), S(4/3)]) # PFDD
  863. # XXX this does not *evaluate* right??
  864. assert can_do([], [a, a + S.Half, 2*a - 1])
  865. def test_bug():
  866. h = hyper([-1, 1], [z], -1)
  867. assert hyperexpand(h) == (z + 1)/z
  868. def test_omgissue_203():
  869. h = hyper((-5, -3, -4), (-6, -6), 1)
  870. assert hyperexpand(h) == Rational(1, 30)
  871. h = hyper((-6, -7, -5), (-6, -6), 1)
  872. assert hyperexpand(h) == Rational(-1, 6)