图片解析应用
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

389 lines
16 KiB

  1. from sympy.core.containers import Tuple
  2. from sympy.core.function import Derivative
  3. from sympy.core.numbers import (I, Rational, oo, pi)
  4. from sympy.core.singleton import S
  5. from sympy.core.symbol import symbols
  6. from sympy.functions.elementary.exponential import (exp, log)
  7. from sympy.functions.elementary.miscellaneous import sqrt
  8. from sympy.functions.elementary.trigonometric import cos
  9. from sympy.functions.special.gamma_functions import gamma
  10. from sympy.functions.special.hyper import (appellf1, hyper, meijerg)
  11. from sympy.series.order import O
  12. from sympy.abc import x, z, k
  13. from sympy.series.limits import limit
  14. from sympy.testing.pytest import raises, slow
  15. from sympy.core.random import (
  16. random_complex_number as randcplx,
  17. verify_numerically as tn,
  18. test_derivative_numerically as td)
  19. def test_TupleParametersBase():
  20. # test that our implementation of the chain rule works
  21. p = hyper((), (), z**2)
  22. assert p.diff(z) == p*2*z
  23. def test_hyper():
  24. raises(TypeError, lambda: hyper(1, 2, z))
  25. assert hyper((1, 2), (1,), z) == hyper(Tuple(1, 2), Tuple(1), z)
  26. h = hyper((1, 2), (3, 4, 5), z)
  27. assert h.ap == Tuple(1, 2)
  28. assert h.bq == Tuple(3, 4, 5)
  29. assert h.argument == z
  30. assert h.is_commutative is True
  31. # just a few checks to make sure that all arguments go where they should
  32. assert tn(hyper(Tuple(), Tuple(), z), exp(z), z)
  33. assert tn(z*hyper((1, 1), Tuple(2), -z), log(1 + z), z)
  34. # differentiation
  35. h = hyper(
  36. (randcplx(), randcplx(), randcplx()), (randcplx(), randcplx()), z)
  37. assert td(h, z)
  38. a1, a2, b1, b2, b3 = symbols('a1:3, b1:4')
  39. assert hyper((a1, a2), (b1, b2, b3), z).diff(z) == \
  40. a1*a2/(b1*b2*b3) * hyper((a1 + 1, a2 + 1), (b1 + 1, b2 + 1, b3 + 1), z)
  41. # differentiation wrt parameters is not supported
  42. assert hyper([z], [], z).diff(z) == Derivative(hyper([z], [], z), z)
  43. # hyper is unbranched wrt parameters
  44. from sympy.functions.elementary.complexes import polar_lift
  45. assert hyper([polar_lift(z)], [polar_lift(k)], polar_lift(x)) == \
  46. hyper([z], [k], polar_lift(x))
  47. # hyper does not automatically evaluate anyway, but the test is to make
  48. # sure that the evaluate keyword is accepted
  49. assert hyper((1, 2), (1,), z, evaluate=False).func is hyper
  50. def test_expand_func():
  51. # evaluation at 1 of Gauss' hypergeometric function:
  52. from sympy.abc import a, b, c
  53. from sympy.core.function import expand_func
  54. a1, b1, c1 = randcplx(), randcplx(), randcplx() + 5
  55. assert expand_func(hyper([a, b], [c], 1)) == \
  56. gamma(c)*gamma(-a - b + c)/(gamma(-a + c)*gamma(-b + c))
  57. assert abs(expand_func(hyper([a1, b1], [c1], 1)).n()
  58. - hyper([a1, b1], [c1], 1).n()) < 1e-10
  59. # hyperexpand wrapper for hyper:
  60. assert expand_func(hyper([], [], z)) == exp(z)
  61. assert expand_func(hyper([1, 2, 3], [], z)) == hyper([1, 2, 3], [], z)
  62. assert expand_func(meijerg([[1, 1], []], [[1], [0]], z)) == log(z + 1)
  63. assert expand_func(meijerg([[1, 1], []], [[], []], z)) == \
  64. meijerg([[1, 1], []], [[], []], z)
  65. def replace_dummy(expr, sym):
  66. from sympy.core.symbol import Dummy
  67. dum = expr.atoms(Dummy)
  68. if not dum:
  69. return expr
  70. assert len(dum) == 1
  71. return expr.xreplace({dum.pop(): sym})
  72. def test_hyper_rewrite_sum():
  73. from sympy.concrete.summations import Sum
  74. from sympy.core.symbol import Dummy
  75. from sympy.functions.combinatorial.factorials import (RisingFactorial, factorial)
  76. _k = Dummy("k")
  77. assert replace_dummy(hyper((1, 2), (1, 3), x).rewrite(Sum), _k) == \
  78. Sum(x**_k / factorial(_k) * RisingFactorial(2, _k) /
  79. RisingFactorial(3, _k), (_k, 0, oo))
  80. assert hyper((1, 2, 3), (-1, 3), z).rewrite(Sum) == \
  81. hyper((1, 2, 3), (-1, 3), z)
  82. def test_radius_of_convergence():
  83. assert hyper((1, 2), [3], z).radius_of_convergence == 1
  84. assert hyper((1, 2), [3, 4], z).radius_of_convergence is oo
  85. assert hyper((1, 2, 3), [4], z).radius_of_convergence == 0
  86. assert hyper((0, 1, 2), [4], z).radius_of_convergence is oo
  87. assert hyper((-1, 1, 2), [-4], z).radius_of_convergence == 0
  88. assert hyper((-1, -2, 2), [-1], z).radius_of_convergence is oo
  89. assert hyper((-1, 2), [-1, -2], z).radius_of_convergence == 0
  90. assert hyper([-1, 1, 3], [-2, 2], z).radius_of_convergence == 1
  91. assert hyper([-1, 1], [-2, 2], z).radius_of_convergence is oo
  92. assert hyper([-1, 1, 3], [-2], z).radius_of_convergence == 0
  93. assert hyper((-1, 2, 3, 4), [], z).radius_of_convergence is oo
  94. assert hyper([1, 1], [3], 1).convergence_statement == True
  95. assert hyper([1, 1], [2], 1).convergence_statement == False
  96. assert hyper([1, 1], [2], -1).convergence_statement == True
  97. assert hyper([1, 1], [1], -1).convergence_statement == False
  98. def test_meijer():
  99. raises(TypeError, lambda: meijerg(1, z))
  100. raises(TypeError, lambda: meijerg(((1,), (2,)), (3,), (4,), z))
  101. assert meijerg(((1, 2), (3,)), ((4,), (5,)), z) == \
  102. meijerg(Tuple(1, 2), Tuple(3), Tuple(4), Tuple(5), z)
  103. g = meijerg((1, 2), (3, 4, 5), (6, 7, 8, 9), (10, 11, 12, 13, 14), z)
  104. assert g.an == Tuple(1, 2)
  105. assert g.ap == Tuple(1, 2, 3, 4, 5)
  106. assert g.aother == Tuple(3, 4, 5)
  107. assert g.bm == Tuple(6, 7, 8, 9)
  108. assert g.bq == Tuple(6, 7, 8, 9, 10, 11, 12, 13, 14)
  109. assert g.bother == Tuple(10, 11, 12, 13, 14)
  110. assert g.argument == z
  111. assert g.nu == 75
  112. assert g.delta == -1
  113. assert g.is_commutative is True
  114. assert g.is_number is False
  115. #issue 13071
  116. assert meijerg([[],[]], [[S.Half],[0]], 1).is_number is True
  117. assert meijerg([1, 2], [3], [4], [5], z).delta == S.Half
  118. # just a few checks to make sure that all arguments go where they should
  119. assert tn(meijerg(Tuple(), Tuple(), Tuple(0), Tuple(), -z), exp(z), z)
  120. assert tn(sqrt(pi)*meijerg(Tuple(), Tuple(),
  121. Tuple(0), Tuple(S.Half), z**2/4), cos(z), z)
  122. assert tn(meijerg(Tuple(1, 1), Tuple(), Tuple(1), Tuple(0), z),
  123. log(1 + z), z)
  124. # test exceptions
  125. raises(ValueError, lambda: meijerg(((3, 1), (2,)), ((oo,), (2, 0)), x))
  126. raises(ValueError, lambda: meijerg(((3, 1), (2,)), ((1,), (2, 0)), x))
  127. # differentiation
  128. g = meijerg((randcplx(),), (randcplx() + 2*I,), Tuple(),
  129. (randcplx(), randcplx()), z)
  130. assert td(g, z)
  131. g = meijerg(Tuple(), (randcplx(),), Tuple(),
  132. (randcplx(), randcplx()), z)
  133. assert td(g, z)
  134. g = meijerg(Tuple(), Tuple(), Tuple(randcplx()),
  135. Tuple(randcplx(), randcplx()), z)
  136. assert td(g, z)
  137. a1, a2, b1, b2, c1, c2, d1, d2 = symbols('a1:3, b1:3, c1:3, d1:3')
  138. assert meijerg((a1, a2), (b1, b2), (c1, c2), (d1, d2), z).diff(z) == \
  139. (meijerg((a1 - 1, a2), (b1, b2), (c1, c2), (d1, d2), z)
  140. + (a1 - 1)*meijerg((a1, a2), (b1, b2), (c1, c2), (d1, d2), z))/z
  141. assert meijerg([z, z], [], [], [], z).diff(z) == \
  142. Derivative(meijerg([z, z], [], [], [], z), z)
  143. # meijerg is unbranched wrt parameters
  144. from sympy.functions.elementary.complexes import polar_lift as pl
  145. assert meijerg([pl(a1)], [pl(a2)], [pl(b1)], [pl(b2)], pl(z)) == \
  146. meijerg([a1], [a2], [b1], [b2], pl(z))
  147. # integrand
  148. from sympy.abc import a, b, c, d, s
  149. assert meijerg([a], [b], [c], [d], z).integrand(s) == \
  150. z**s*gamma(c - s)*gamma(-a + s + 1)/(gamma(b - s)*gamma(-d + s + 1))
  151. def test_meijerg_derivative():
  152. assert meijerg([], [1, 1], [0, 0, x], [], z).diff(x) == \
  153. log(z)*meijerg([], [1, 1], [0, 0, x], [], z) \
  154. + 2*meijerg([], [1, 1, 1], [0, 0, x, 0], [], z)
  155. y = randcplx()
  156. a = 5 # mpmath chokes with non-real numbers, and Mod1 with floats
  157. assert td(meijerg([x], [], [], [], y), x)
  158. assert td(meijerg([x**2], [], [], [], y), x)
  159. assert td(meijerg([], [x], [], [], y), x)
  160. assert td(meijerg([], [], [x], [], y), x)
  161. assert td(meijerg([], [], [], [x], y), x)
  162. assert td(meijerg([x], [a], [a + 1], [], y), x)
  163. assert td(meijerg([x], [a + 1], [a], [], y), x)
  164. assert td(meijerg([x, a], [], [], [a + 1], y), x)
  165. assert td(meijerg([x, a + 1], [], [], [a], y), x)
  166. b = Rational(3, 2)
  167. assert td(meijerg([a + 2], [b], [b - 3, x], [a], y), x)
  168. def test_meijerg_period():
  169. assert meijerg([], [1], [0], [], x).get_period() == 2*pi
  170. assert meijerg([1], [], [], [0], x).get_period() == 2*pi
  171. assert meijerg([], [], [0], [], x).get_period() == 2*pi # exp(x)
  172. assert meijerg(
  173. [], [], [0], [S.Half], x).get_period() == 2*pi # cos(sqrt(x))
  174. assert meijerg(
  175. [], [], [S.Half], [0], x).get_period() == 4*pi # sin(sqrt(x))
  176. assert meijerg([1, 1], [], [1], [0], x).get_period() is oo # log(1 + x)
  177. def test_hyper_unpolarify():
  178. from sympy.functions.elementary.exponential import exp_polar
  179. a = exp_polar(2*pi*I)*x
  180. b = x
  181. assert hyper([], [], a).argument == b
  182. assert hyper([0], [], a).argument == a
  183. assert hyper([0], [0], a).argument == b
  184. assert hyper([0, 1], [0], a).argument == a
  185. assert hyper([0, 1], [0], exp_polar(2*pi*I)).argument == 1
  186. @slow
  187. def test_hyperrep():
  188. from sympy.functions.special.hyper import (HyperRep, HyperRep_atanh,
  189. HyperRep_power1, HyperRep_power2, HyperRep_log1, HyperRep_asin1,
  190. HyperRep_asin2, HyperRep_sqrts1, HyperRep_sqrts2, HyperRep_log2,
  191. HyperRep_cosasin, HyperRep_sinasin)
  192. # First test the base class works.
  193. from sympy.functions.elementary.exponential import exp_polar
  194. from sympy.functions.elementary.piecewise import Piecewise
  195. a, b, c, d, z = symbols('a b c d z')
  196. class myrep(HyperRep):
  197. @classmethod
  198. def _expr_small(cls, x):
  199. return a
  200. @classmethod
  201. def _expr_small_minus(cls, x):
  202. return b
  203. @classmethod
  204. def _expr_big(cls, x, n):
  205. return c*n
  206. @classmethod
  207. def _expr_big_minus(cls, x, n):
  208. return d*n
  209. assert myrep(z).rewrite('nonrep') == Piecewise((0, abs(z) > 1), (a, True))
  210. assert myrep(exp_polar(I*pi)*z).rewrite('nonrep') == \
  211. Piecewise((0, abs(z) > 1), (b, True))
  212. assert myrep(exp_polar(2*I*pi)*z).rewrite('nonrep') == \
  213. Piecewise((c, abs(z) > 1), (a, True))
  214. assert myrep(exp_polar(3*I*pi)*z).rewrite('nonrep') == \
  215. Piecewise((d, abs(z) > 1), (b, True))
  216. assert myrep(exp_polar(4*I*pi)*z).rewrite('nonrep') == \
  217. Piecewise((2*c, abs(z) > 1), (a, True))
  218. assert myrep(exp_polar(5*I*pi)*z).rewrite('nonrep') == \
  219. Piecewise((2*d, abs(z) > 1), (b, True))
  220. assert myrep(z).rewrite('nonrepsmall') == a
  221. assert myrep(exp_polar(I*pi)*z).rewrite('nonrepsmall') == b
  222. def t(func, hyp, z):
  223. """ Test that func is a valid representation of hyp. """
  224. # First test that func agrees with hyp for small z
  225. if not tn(func.rewrite('nonrepsmall'), hyp, z,
  226. a=Rational(-1, 2), b=Rational(-1, 2), c=S.Half, d=S.Half):
  227. return False
  228. # Next check that the two small representations agree.
  229. if not tn(
  230. func.rewrite('nonrepsmall').subs(
  231. z, exp_polar(I*pi)*z).replace(exp_polar, exp),
  232. func.subs(z, exp_polar(I*pi)*z).rewrite('nonrepsmall'),
  233. z, a=Rational(-1, 2), b=Rational(-1, 2), c=S.Half, d=S.Half):
  234. return False
  235. # Next check continuity along exp_polar(I*pi)*t
  236. expr = func.subs(z, exp_polar(I*pi)*z).rewrite('nonrep')
  237. if abs(expr.subs(z, 1 + 1e-15).n() - expr.subs(z, 1 - 1e-15).n()) > 1e-10:
  238. return False
  239. # Finally check continuity of the big reps.
  240. def dosubs(func, a, b):
  241. rv = func.subs(z, exp_polar(a)*z).rewrite('nonrep')
  242. return rv.subs(z, exp_polar(b)*z).replace(exp_polar, exp)
  243. for n in [0, 1, 2, 3, 4, -1, -2, -3, -4]:
  244. expr1 = dosubs(func, 2*I*pi*n, I*pi/2)
  245. expr2 = dosubs(func, 2*I*pi*n + I*pi, -I*pi/2)
  246. if not tn(expr1, expr2, z):
  247. return False
  248. expr1 = dosubs(func, 2*I*pi*(n + 1), -I*pi/2)
  249. expr2 = dosubs(func, 2*I*pi*n + I*pi, I*pi/2)
  250. if not tn(expr1, expr2, z):
  251. return False
  252. return True
  253. # Now test the various representatives.
  254. a = Rational(1, 3)
  255. assert t(HyperRep_atanh(z), hyper([S.Half, 1], [Rational(3, 2)], z), z)
  256. assert t(HyperRep_power1(a, z), hyper([-a], [], z), z)
  257. assert t(HyperRep_power2(a, z), hyper([a, a - S.Half], [2*a], z), z)
  258. assert t(HyperRep_log1(z), -z*hyper([1, 1], [2], z), z)
  259. assert t(HyperRep_asin1(z), hyper([S.Half, S.Half], [Rational(3, 2)], z), z)
  260. assert t(HyperRep_asin2(z), hyper([1, 1], [Rational(3, 2)], z), z)
  261. assert t(HyperRep_sqrts1(a, z), hyper([-a, S.Half - a], [S.Half], z), z)
  262. assert t(HyperRep_sqrts2(a, z),
  263. -2*z/(2*a + 1)*hyper([-a - S.Half, -a], [S.Half], z).diff(z), z)
  264. assert t(HyperRep_log2(z), -z/4*hyper([Rational(3, 2), 1, 1], [2, 2], z), z)
  265. assert t(HyperRep_cosasin(a, z), hyper([-a, a], [S.Half], z), z)
  266. assert t(HyperRep_sinasin(a, z), 2*a*z*hyper([1 - a, 1 + a], [Rational(3, 2)], z), z)
  267. @slow
  268. def test_meijerg_eval():
  269. from sympy.functions.elementary.exponential import exp_polar
  270. from sympy.functions.special.bessel import besseli
  271. from sympy.abc import l
  272. a = randcplx()
  273. arg = x*exp_polar(k*pi*I)
  274. expr1 = pi*meijerg([[], [(a + 1)/2]], [[a/2], [-a/2, (a + 1)/2]], arg**2/4)
  275. expr2 = besseli(a, arg)
  276. # Test that the two expressions agree for all arguments.
  277. for x_ in [0.5, 1.5]:
  278. for k_ in [0.0, 0.1, 0.3, 0.5, 0.8, 1, 5.751, 15.3]:
  279. assert abs((expr1 - expr2).n(subs={x: x_, k: k_})) < 1e-10
  280. assert abs((expr1 - expr2).n(subs={x: x_, k: -k_})) < 1e-10
  281. # Test continuity independently
  282. eps = 1e-13
  283. expr2 = expr1.subs(k, l)
  284. for x_ in [0.5, 1.5]:
  285. for k_ in [0.5, Rational(1, 3), 0.25, 0.75, Rational(2, 3), 1.0, 1.5]:
  286. assert abs((expr1 - expr2).n(
  287. subs={x: x_, k: k_ + eps, l: k_ - eps})) < 1e-10
  288. assert abs((expr1 - expr2).n(
  289. subs={x: x_, k: -k_ + eps, l: -k_ - eps})) < 1e-10
  290. expr = (meijerg(((0.5,), ()), ((0.5, 0, 0.5), ()), exp_polar(-I*pi)/4)
  291. + meijerg(((0.5,), ()), ((0.5, 0, 0.5), ()), exp_polar(I*pi)/4)) \
  292. /(2*sqrt(pi))
  293. assert (expr - pi/exp(1)).n(chop=True) == 0
  294. def test_limits():
  295. k, x = symbols('k, x')
  296. assert hyper((1,), (Rational(4, 3), Rational(5, 3)), k**2).series(k) == \
  297. 1 + 9*k**2/20 + 81*k**4/1120 + O(k**6) # issue 6350
  298. assert limit(meijerg((), (), (1,), (0,), -x), x, 0) == \
  299. meijerg(((), ()), ((1,), (0,)), 0) # issue 6052
  300. # https://github.com/sympy/sympy/issues/11465
  301. assert limit(1/hyper((1, ), (1, ), x), x, 0) == 1
  302. def test_appellf1():
  303. a, b1, b2, c, x, y = symbols('a b1 b2 c x y')
  304. assert appellf1(a, b2, b1, c, y, x) == appellf1(a, b1, b2, c, x, y)
  305. assert appellf1(a, b1, b1, c, y, x) == appellf1(a, b1, b1, c, x, y)
  306. assert appellf1(a, b1, b2, c, S.Zero, S.Zero) is S.One
  307. f = appellf1(a, b1, b2, c, S.Zero, S.Zero, evaluate=False)
  308. assert f.func is appellf1
  309. assert f.doit() is S.One
  310. def test_derivative_appellf1():
  311. from sympy.core.function import diff
  312. a, b1, b2, c, x, y, z = symbols('a b1 b2 c x y z')
  313. assert diff(appellf1(a, b1, b2, c, x, y), x) == a*b1*appellf1(a + 1, b2, b1 + 1, c + 1, y, x)/c
  314. assert diff(appellf1(a, b1, b2, c, x, y), y) == a*b2*appellf1(a + 1, b1, b2 + 1, c + 1, x, y)/c
  315. assert diff(appellf1(a, b1, b2, c, x, y), z) == 0
  316. assert diff(appellf1(a, b1, b2, c, x, y), a) == Derivative(appellf1(a, b1, b2, c, x, y), a)
  317. def test_eval_nseries():
  318. a1, b1, a2, b2 = symbols('a1 b1 a2 b2')
  319. assert hyper((1,2), (1,2,3), x**2)._eval_nseries(x, 7, None) == 1 + x**2/3 + x**4/24 + x**6/360 + O(x**7)
  320. assert exp(x)._eval_nseries(x,7,None) == hyper((a1, b1), (a1, b1), x)._eval_nseries(x, 7, None)
  321. assert hyper((a1, a2), (b1, b2), x)._eval_nseries(z, 7, None) == hyper((a1, a2), (b1, b2), x) + O(z**7)