test_error_functions.py 30 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816
  1. from sympy.core.function import (diff, expand, expand_func)
  2. from sympy.core import EulerGamma
  3. from sympy.core.numbers import (E, Float, I, Rational, nan, oo, pi)
  4. from sympy.core.singleton import S
  5. from sympy.core.symbol import (Symbol, symbols)
  6. from sympy.functions.elementary.complexes import (conjugate, im, polar_lift, re)
  7. from sympy.functions.elementary.exponential import (exp, exp_polar, log)
  8. from sympy.functions.elementary.hyperbolic import (cosh, sinh)
  9. from sympy.functions.elementary.miscellaneous import sqrt
  10. from sympy.functions.elementary.trigonometric import (cos, sin, sinc)
  11. from sympy.functions.special.error_functions import (Chi, Ci, E1, Ei, Li, Shi, Si, erf, erf2, erf2inv, erfc, erfcinv, erfi, erfinv, expint, fresnelc, fresnels, li)
  12. from sympy.functions.special.gamma_functions import (gamma, uppergamma)
  13. from sympy.functions.special.hyper import (hyper, meijerg)
  14. from sympy.integrals.integrals import (Integral, integrate)
  15. from sympy.series.gruntz import gruntz
  16. from sympy.series.limits import limit
  17. from sympy.series.order import O
  18. from sympy.core.expr import unchanged
  19. from sympy.core.function import ArgumentIndexError
  20. from sympy.functions.special.error_functions import _erfs, _eis
  21. from sympy.testing.pytest import raises
  22. x, y, z = symbols('x,y,z')
  23. w = Symbol("w", real=True)
  24. n = Symbol("n", integer=True)
  25. def test_erf():
  26. assert erf(nan) is nan
  27. assert erf(oo) == 1
  28. assert erf(-oo) == -1
  29. assert erf(0) is S.Zero
  30. assert erf(I*oo) == oo*I
  31. assert erf(-I*oo) == -oo*I
  32. assert erf(-2) == -erf(2)
  33. assert erf(-x*y) == -erf(x*y)
  34. assert erf(-x - y) == -erf(x + y)
  35. assert erf(erfinv(x)) == x
  36. assert erf(erfcinv(x)) == 1 - x
  37. assert erf(erf2inv(0, x)) == x
  38. assert erf(erf2inv(0, x, evaluate=False)) == x # To cover code in erf
  39. assert erf(erf2inv(0, erf(erfcinv(1 - erf(erfinv(x)))))) == x
  40. assert erf(I).is_real is False
  41. assert erf(0, evaluate=False).is_real
  42. assert erf(0, evaluate=False).is_zero
  43. assert conjugate(erf(z)) == erf(conjugate(z))
  44. assert erf(x).as_leading_term(x) == 2*x/sqrt(pi)
  45. assert erf(x*y).as_leading_term(y) == 2*x*y/sqrt(pi)
  46. assert (erf(x*y)/erf(y)).as_leading_term(y) == x
  47. assert erf(1/x).as_leading_term(x) == S.One
  48. assert erf(z).rewrite('uppergamma') == sqrt(z**2)*(1 - erfc(sqrt(z**2)))/z
  49. assert erf(z).rewrite('erfc') == S.One - erfc(z)
  50. assert erf(z).rewrite('erfi') == -I*erfi(I*z)
  51. assert erf(z).rewrite('fresnels') == (1 + I)*(fresnelc(z*(1 - I)/sqrt(pi)) -
  52. I*fresnels(z*(1 - I)/sqrt(pi)))
  53. assert erf(z).rewrite('fresnelc') == (1 + I)*(fresnelc(z*(1 - I)/sqrt(pi)) -
  54. I*fresnels(z*(1 - I)/sqrt(pi)))
  55. assert erf(z).rewrite('hyper') == 2*z*hyper([S.Half], [3*S.Half], -z**2)/sqrt(pi)
  56. assert erf(z).rewrite('meijerg') == z*meijerg([S.Half], [], [0], [Rational(-1, 2)], z**2)/sqrt(pi)
  57. assert erf(z).rewrite('expint') == sqrt(z**2)/z - z*expint(S.Half, z**2)/sqrt(S.Pi)
  58. assert limit(exp(x)*exp(x**2)*(erf(x + 1/exp(x)) - erf(x)), x, oo) == \
  59. 2/sqrt(pi)
  60. assert limit((1 - erf(z))*exp(z**2)*z, z, oo) == 1/sqrt(pi)
  61. assert limit((1 - erf(x))*exp(x**2)*sqrt(pi)*x, x, oo) == 1
  62. assert limit(((1 - erf(x))*exp(x**2)*sqrt(pi)*x - 1)*2*x**2, x, oo) == -1
  63. assert limit(erf(x)/x, x, 0) == 2/sqrt(pi)
  64. assert limit(x**(-4) - sqrt(pi)*erf(x**2) / (2*x**6), x, 0) == S(1)/3
  65. assert erf(x).as_real_imag() == \
  66. (erf(re(x) - I*im(x))/2 + erf(re(x) + I*im(x))/2,
  67. -I*(-erf(re(x) - I*im(x)) + erf(re(x) + I*im(x)))/2)
  68. assert erf(x).as_real_imag(deep=False) == \
  69. (erf(re(x) - I*im(x))/2 + erf(re(x) + I*im(x))/2,
  70. -I*(-erf(re(x) - I*im(x)) + erf(re(x) + I*im(x)))/2)
  71. assert erf(w).as_real_imag() == (erf(w), 0)
  72. assert erf(w).as_real_imag(deep=False) == (erf(w), 0)
  73. # issue 13575
  74. assert erf(I).as_real_imag() == (0, -I*erf(I))
  75. raises(ArgumentIndexError, lambda: erf(x).fdiff(2))
  76. assert erf(x).inverse() == erfinv
  77. def test_erf_series():
  78. assert erf(x).series(x, 0, 7) == 2*x/sqrt(pi) - \
  79. 2*x**3/3/sqrt(pi) + x**5/5/sqrt(pi) + O(x**7)
  80. assert erf(x).series(x, oo) == \
  81. -exp(-x**2)*(3/(4*x**5) - 1/(2*x**3) + 1/x + O(x**(-6), (x, oo)))/sqrt(pi) + 1
  82. assert erf(x**2).series(x, oo, n=8) == \
  83. (-1/(2*x**6) + x**(-2) + O(x**(-8), (x, oo)))*exp(-x**4)/sqrt(pi)*-1 + 1
  84. assert erf(sqrt(x)).series(x, oo, n=3) == (sqrt(1/x) - (1/x)**(S(3)/2)/2\
  85. + 3*(1/x)**(S(5)/2)/4 + O(x**(-3), (x, oo)))*exp(-x)/sqrt(pi)*-1 + 1
  86. def test_erf_evalf():
  87. assert abs( erf(Float(2.0)) - 0.995322265 ) < 1E-8 # XXX
  88. def test__erfs():
  89. assert _erfs(z).diff(z) == -2/sqrt(S.Pi) + 2*z*_erfs(z)
  90. assert _erfs(1/z).series(z) == \
  91. z/sqrt(pi) - z**3/(2*sqrt(pi)) + 3*z**5/(4*sqrt(pi)) + O(z**6)
  92. assert expand(erf(z).rewrite('tractable').diff(z).rewrite('intractable')) \
  93. == erf(z).diff(z)
  94. assert _erfs(z).rewrite("intractable") == (-erf(z) + 1)*exp(z**2)
  95. raises(ArgumentIndexError, lambda: _erfs(z).fdiff(2))
  96. def test_erfc():
  97. assert erfc(nan) is nan
  98. assert erfc(oo) is S.Zero
  99. assert erfc(-oo) == 2
  100. assert erfc(0) == 1
  101. assert erfc(I*oo) == -oo*I
  102. assert erfc(-I*oo) == oo*I
  103. assert erfc(-x) == S(2) - erfc(x)
  104. assert erfc(erfcinv(x)) == x
  105. assert erfc(I).is_real is False
  106. assert erfc(0, evaluate=False).is_real
  107. assert erfc(0, evaluate=False).is_zero is False
  108. assert erfc(erfinv(x)) == 1 - x
  109. assert conjugate(erfc(z)) == erfc(conjugate(z))
  110. assert erfc(x).as_leading_term(x) is S.One
  111. assert erfc(1/x).as_leading_term(x) == S.Zero
  112. assert erfc(z).rewrite('erf') == 1 - erf(z)
  113. assert erfc(z).rewrite('erfi') == 1 + I*erfi(I*z)
  114. assert erfc(z).rewrite('fresnels') == 1 - (1 + I)*(fresnelc(z*(1 - I)/sqrt(pi)) -
  115. I*fresnels(z*(1 - I)/sqrt(pi)))
  116. assert erfc(z).rewrite('fresnelc') == 1 - (1 + I)*(fresnelc(z*(1 - I)/sqrt(pi)) -
  117. I*fresnels(z*(1 - I)/sqrt(pi)))
  118. assert erfc(z).rewrite('hyper') == 1 - 2*z*hyper([S.Half], [3*S.Half], -z**2)/sqrt(pi)
  119. assert erfc(z).rewrite('meijerg') == 1 - z*meijerg([S.Half], [], [0], [Rational(-1, 2)], z**2)/sqrt(pi)
  120. assert erfc(z).rewrite('uppergamma') == 1 - sqrt(z**2)*(1 - erfc(sqrt(z**2)))/z
  121. assert erfc(z).rewrite('expint') == S.One - sqrt(z**2)/z + z*expint(S.Half, z**2)/sqrt(S.Pi)
  122. assert erfc(z).rewrite('tractable') == _erfs(z)*exp(-z**2)
  123. assert expand_func(erf(x) + erfc(x)) is S.One
  124. assert erfc(x).as_real_imag() == \
  125. (erfc(re(x) - I*im(x))/2 + erfc(re(x) + I*im(x))/2,
  126. -I*(-erfc(re(x) - I*im(x)) + erfc(re(x) + I*im(x)))/2)
  127. assert erfc(x).as_real_imag(deep=False) == \
  128. (erfc(re(x) - I*im(x))/2 + erfc(re(x) + I*im(x))/2,
  129. -I*(-erfc(re(x) - I*im(x)) + erfc(re(x) + I*im(x)))/2)
  130. assert erfc(w).as_real_imag() == (erfc(w), 0)
  131. assert erfc(w).as_real_imag(deep=False) == (erfc(w), 0)
  132. raises(ArgumentIndexError, lambda: erfc(x).fdiff(2))
  133. assert erfc(x).inverse() == erfcinv
  134. def test_erfc_series():
  135. assert erfc(x).series(x, 0, 7) == 1 - 2*x/sqrt(pi) + \
  136. 2*x**3/3/sqrt(pi) - x**5/5/sqrt(pi) + O(x**7)
  137. assert erfc(x).series(x, oo) == \
  138. (3/(4*x**5) - 1/(2*x**3) + 1/x + O(x**(-6), (x, oo)))*exp(-x**2)/sqrt(pi)
  139. def test_erfc_evalf():
  140. assert abs( erfc(Float(2.0)) - 0.00467773 ) < 1E-8 # XXX
  141. def test_erfi():
  142. assert erfi(nan) is nan
  143. assert erfi(oo) is S.Infinity
  144. assert erfi(-oo) is S.NegativeInfinity
  145. assert erfi(0) is S.Zero
  146. assert erfi(I*oo) == I
  147. assert erfi(-I*oo) == -I
  148. assert erfi(-x) == -erfi(x)
  149. assert erfi(I*erfinv(x)) == I*x
  150. assert erfi(I*erfcinv(x)) == I*(1 - x)
  151. assert erfi(I*erf2inv(0, x)) == I*x
  152. assert erfi(I*erf2inv(0, x, evaluate=False)) == I*x # To cover code in erfi
  153. assert erfi(I).is_real is False
  154. assert erfi(0, evaluate=False).is_real
  155. assert erfi(0, evaluate=False).is_zero
  156. assert conjugate(erfi(z)) == erfi(conjugate(z))
  157. assert erfi(x).as_leading_term(x) == 2*x/sqrt(pi)
  158. assert erfi(x*y).as_leading_term(y) == 2*x*y/sqrt(pi)
  159. assert (erfi(x*y)/erfi(y)).as_leading_term(y) == x
  160. assert erfi(1/x).as_leading_term(x) == erfi(1/x)
  161. assert erfi(z).rewrite('erf') == -I*erf(I*z)
  162. assert erfi(z).rewrite('erfc') == I*erfc(I*z) - I
  163. assert erfi(z).rewrite('fresnels') == (1 - I)*(fresnelc(z*(1 + I)/sqrt(pi)) -
  164. I*fresnels(z*(1 + I)/sqrt(pi)))
  165. assert erfi(z).rewrite('fresnelc') == (1 - I)*(fresnelc(z*(1 + I)/sqrt(pi)) -
  166. I*fresnels(z*(1 + I)/sqrt(pi)))
  167. assert erfi(z).rewrite('hyper') == 2*z*hyper([S.Half], [3*S.Half], z**2)/sqrt(pi)
  168. assert erfi(z).rewrite('meijerg') == z*meijerg([S.Half], [], [0], [Rational(-1, 2)], -z**2)/sqrt(pi)
  169. assert erfi(z).rewrite('uppergamma') == (sqrt(-z**2)/z*(uppergamma(S.Half,
  170. -z**2)/sqrt(S.Pi) - S.One))
  171. assert erfi(z).rewrite('expint') == sqrt(-z**2)/z - z*expint(S.Half, -z**2)/sqrt(S.Pi)
  172. assert erfi(z).rewrite('tractable') == -I*(-_erfs(I*z)*exp(z**2) + 1)
  173. assert expand_func(erfi(I*z)) == I*erf(z)
  174. assert erfi(x).as_real_imag() == \
  175. (erfi(re(x) - I*im(x))/2 + erfi(re(x) + I*im(x))/2,
  176. -I*(-erfi(re(x) - I*im(x)) + erfi(re(x) + I*im(x)))/2)
  177. assert erfi(x).as_real_imag(deep=False) == \
  178. (erfi(re(x) - I*im(x))/2 + erfi(re(x) + I*im(x))/2,
  179. -I*(-erfi(re(x) - I*im(x)) + erfi(re(x) + I*im(x)))/2)
  180. assert erfi(w).as_real_imag() == (erfi(w), 0)
  181. assert erfi(w).as_real_imag(deep=False) == (erfi(w), 0)
  182. raises(ArgumentIndexError, lambda: erfi(x).fdiff(2))
  183. def test_erfi_series():
  184. assert erfi(x).series(x, 0, 7) == 2*x/sqrt(pi) + \
  185. 2*x**3/3/sqrt(pi) + x**5/5/sqrt(pi) + O(x**7)
  186. assert erfi(x).series(x, oo) == \
  187. (3/(4*x**5) + 1/(2*x**3) + 1/x + O(x**(-6), (x, oo)))*exp(x**2)/sqrt(pi) - I
  188. def test_erfi_evalf():
  189. assert abs( erfi(Float(2.0)) - 18.5648024145756 ) < 1E-13 # XXX
  190. def test_erf2():
  191. assert erf2(0, 0) is S.Zero
  192. assert erf2(x, x) is S.Zero
  193. assert erf2(nan, 0) is nan
  194. assert erf2(-oo, y) == erf(y) + 1
  195. assert erf2( oo, y) == erf(y) - 1
  196. assert erf2( x, oo) == 1 - erf(x)
  197. assert erf2( x,-oo) == -1 - erf(x)
  198. assert erf2(x, erf2inv(x, y)) == y
  199. assert erf2(-x, -y) == -erf2(x,y)
  200. assert erf2(-x, y) == erf(y) + erf(x)
  201. assert erf2( x, -y) == -erf(y) - erf(x)
  202. assert erf2(x, y).rewrite('fresnels') == erf(y).rewrite(fresnels)-erf(x).rewrite(fresnels)
  203. assert erf2(x, y).rewrite('fresnelc') == erf(y).rewrite(fresnelc)-erf(x).rewrite(fresnelc)
  204. assert erf2(x, y).rewrite('hyper') == erf(y).rewrite(hyper)-erf(x).rewrite(hyper)
  205. assert erf2(x, y).rewrite('meijerg') == erf(y).rewrite(meijerg)-erf(x).rewrite(meijerg)
  206. assert erf2(x, y).rewrite('uppergamma') == erf(y).rewrite(uppergamma) - erf(x).rewrite(uppergamma)
  207. assert erf2(x, y).rewrite('expint') == erf(y).rewrite(expint)-erf(x).rewrite(expint)
  208. assert erf2(I, 0).is_real is False
  209. assert erf2(0, 0, evaluate=False).is_real
  210. assert erf2(0, 0, evaluate=False).is_zero
  211. assert erf2(x, x, evaluate=False).is_zero
  212. assert erf2(x, y).is_zero is None
  213. assert expand_func(erf(x) + erf2(x, y)) == erf(y)
  214. assert conjugate(erf2(x, y)) == erf2(conjugate(x), conjugate(y))
  215. assert erf2(x, y).rewrite('erf') == erf(y) - erf(x)
  216. assert erf2(x, y).rewrite('erfc') == erfc(x) - erfc(y)
  217. assert erf2(x, y).rewrite('erfi') == I*(erfi(I*x) - erfi(I*y))
  218. assert erf2(x, y).diff(x) == erf2(x, y).fdiff(1)
  219. assert erf2(x, y).diff(y) == erf2(x, y).fdiff(2)
  220. assert erf2(x, y).diff(x) == -2*exp(-x**2)/sqrt(pi)
  221. assert erf2(x, y).diff(y) == 2*exp(-y**2)/sqrt(pi)
  222. raises(ArgumentIndexError, lambda: erf2(x, y).fdiff(3))
  223. assert erf2(x, y).is_extended_real is None
  224. xr, yr = symbols('xr yr', extended_real=True)
  225. assert erf2(xr, yr).is_extended_real is True
  226. def test_erfinv():
  227. assert erfinv(0) is S.Zero
  228. assert erfinv(1) is S.Infinity
  229. assert erfinv(nan) is S.NaN
  230. assert erfinv(-1) is S.NegativeInfinity
  231. assert erfinv(erf(w)) == w
  232. assert erfinv(erf(-w)) == -w
  233. assert erfinv(x).diff() == sqrt(pi)*exp(erfinv(x)**2)/2
  234. raises(ArgumentIndexError, lambda: erfinv(x).fdiff(2))
  235. assert erfinv(z).rewrite('erfcinv') == erfcinv(1-z)
  236. assert erfinv(z).inverse() == erf
  237. def test_erfinv_evalf():
  238. assert abs( erfinv(Float(0.2)) - 0.179143454621292 ) < 1E-13
  239. def test_erfcinv():
  240. assert erfcinv(1) is S.Zero
  241. assert erfcinv(0) is S.Infinity
  242. assert erfcinv(nan) is S.NaN
  243. assert erfcinv(x).diff() == -sqrt(pi)*exp(erfcinv(x)**2)/2
  244. raises(ArgumentIndexError, lambda: erfcinv(x).fdiff(2))
  245. assert erfcinv(z).rewrite('erfinv') == erfinv(1-z)
  246. assert erfcinv(z).inverse() == erfc
  247. def test_erf2inv():
  248. assert erf2inv(0, 0) is S.Zero
  249. assert erf2inv(0, 1) is S.Infinity
  250. assert erf2inv(1, 0) is S.One
  251. assert erf2inv(0, y) == erfinv(y)
  252. assert erf2inv(oo, y) == erfcinv(-y)
  253. assert erf2inv(x, 0) == x
  254. assert erf2inv(x, oo) == erfinv(x)
  255. assert erf2inv(nan, 0) is nan
  256. assert erf2inv(0, nan) is nan
  257. assert erf2inv(x, y).diff(x) == exp(-x**2 + erf2inv(x, y)**2)
  258. assert erf2inv(x, y).diff(y) == sqrt(pi)*exp(erf2inv(x, y)**2)/2
  259. raises(ArgumentIndexError, lambda: erf2inv(x, y).fdiff(3))
  260. # NOTE we multiply by exp_polar(I*pi) and need this to be on the principal
  261. # branch, hence take x in the lower half plane (d=0).
  262. def mytn(expr1, expr2, expr3, x, d=0):
  263. from sympy.core.random import verify_numerically, random_complex_number
  264. subs = {}
  265. for a in expr1.free_symbols:
  266. if a != x:
  267. subs[a] = random_complex_number()
  268. return expr2 == expr3 and verify_numerically(expr1.subs(subs),
  269. expr2.subs(subs), x, d=d)
  270. def mytd(expr1, expr2, x):
  271. from sympy.core.random import test_derivative_numerically, \
  272. random_complex_number
  273. subs = {}
  274. for a in expr1.free_symbols:
  275. if a != x:
  276. subs[a] = random_complex_number()
  277. return expr1.diff(x) == expr2 and test_derivative_numerically(expr1.subs(subs), x)
  278. def tn_branch(func, s=None):
  279. from sympy.core.random import uniform
  280. def fn(x):
  281. if s is None:
  282. return func(x)
  283. return func(s, x)
  284. c = uniform(1, 5)
  285. expr = fn(c*exp_polar(I*pi)) - fn(c*exp_polar(-I*pi))
  286. eps = 1e-15
  287. expr2 = fn(-c + eps*I) - fn(-c - eps*I)
  288. return abs(expr.n() - expr2.n()).n() < 1e-10
  289. def test_ei():
  290. assert Ei(0) is S.NegativeInfinity
  291. assert Ei(oo) is S.Infinity
  292. assert Ei(-oo) is S.Zero
  293. assert tn_branch(Ei)
  294. assert mytd(Ei(x), exp(x)/x, x)
  295. assert mytn(Ei(x), Ei(x).rewrite(uppergamma),
  296. -uppergamma(0, x*polar_lift(-1)) - I*pi, x)
  297. assert mytn(Ei(x), Ei(x).rewrite(expint),
  298. -expint(1, x*polar_lift(-1)) - I*pi, x)
  299. assert Ei(x).rewrite(expint).rewrite(Ei) == Ei(x)
  300. assert Ei(x*exp_polar(2*I*pi)) == Ei(x) + 2*I*pi
  301. assert Ei(x*exp_polar(-2*I*pi)) == Ei(x) - 2*I*pi
  302. assert mytn(Ei(x), Ei(x).rewrite(Shi), Chi(x) + Shi(x), x)
  303. assert mytn(Ei(x*polar_lift(I)), Ei(x*polar_lift(I)).rewrite(Si),
  304. Ci(x) + I*Si(x) + I*pi/2, x)
  305. assert Ei(log(x)).rewrite(li) == li(x)
  306. assert Ei(2*log(x)).rewrite(li) == li(x**2)
  307. assert gruntz(Ei(x+exp(-x))*exp(-x)*x, x, oo) == 1
  308. assert Ei(x).series(x) == EulerGamma + log(x) + x + x**2/4 + \
  309. x**3/18 + x**4/96 + x**5/600 + O(x**6)
  310. assert Ei(x).series(x, 1, 3) == Ei(1) + E*(x - 1) + O((x - 1)**3, (x, 1))
  311. assert Ei(x).series(x, oo) == \
  312. (120/x**5 + 24/x**4 + 6/x**3 + 2/x**2 + 1/x + 1 + O(x**(-6), (x, oo)))*exp(x)/x
  313. assert str(Ei(cos(2)).evalf(n=10)) == '-0.6760647401'
  314. raises(ArgumentIndexError, lambda: Ei(x).fdiff(2))
  315. def test_expint():
  316. assert mytn(expint(x, y), expint(x, y).rewrite(uppergamma),
  317. y**(x - 1)*uppergamma(1 - x, y), x)
  318. assert mytd(
  319. expint(x, y), -y**(x - 1)*meijerg([], [1, 1], [0, 0, 1 - x], [], y), x)
  320. assert mytd(expint(x, y), -expint(x - 1, y), y)
  321. assert mytn(expint(1, x), expint(1, x).rewrite(Ei),
  322. -Ei(x*polar_lift(-1)) + I*pi, x)
  323. assert expint(-4, x) == exp(-x)/x + 4*exp(-x)/x**2 + 12*exp(-x)/x**3 \
  324. + 24*exp(-x)/x**4 + 24*exp(-x)/x**5
  325. assert expint(Rational(-3, 2), x) == \
  326. exp(-x)/x + 3*exp(-x)/(2*x**2) + 3*sqrt(pi)*erfc(sqrt(x))/(4*x**S('5/2'))
  327. assert tn_branch(expint, 1)
  328. assert tn_branch(expint, 2)
  329. assert tn_branch(expint, 3)
  330. assert tn_branch(expint, 1.7)
  331. assert tn_branch(expint, pi)
  332. assert expint(y, x*exp_polar(2*I*pi)) == \
  333. x**(y - 1)*(exp(2*I*pi*y) - 1)*gamma(-y + 1) + expint(y, x)
  334. assert expint(y, x*exp_polar(-2*I*pi)) == \
  335. x**(y - 1)*(exp(-2*I*pi*y) - 1)*gamma(-y + 1) + expint(y, x)
  336. assert expint(2, x*exp_polar(2*I*pi)) == 2*I*pi*x + expint(2, x)
  337. assert expint(2, x*exp_polar(-2*I*pi)) == -2*I*pi*x + expint(2, x)
  338. assert expint(1, x).rewrite(Ei).rewrite(expint) == expint(1, x)
  339. assert expint(x, y).rewrite(Ei) == expint(x, y)
  340. assert expint(x, y).rewrite(Ci) == expint(x, y)
  341. assert mytn(E1(x), E1(x).rewrite(Shi), Shi(x) - Chi(x), x)
  342. assert mytn(E1(polar_lift(I)*x), E1(polar_lift(I)*x).rewrite(Si),
  343. -Ci(x) + I*Si(x) - I*pi/2, x)
  344. assert mytn(expint(2, x), expint(2, x).rewrite(Ei).rewrite(expint),
  345. -x*E1(x) + exp(-x), x)
  346. assert mytn(expint(3, x), expint(3, x).rewrite(Ei).rewrite(expint),
  347. x**2*E1(x)/2 + (1 - x)*exp(-x)/2, x)
  348. assert expint(Rational(3, 2), z).nseries(z) == \
  349. 2 + 2*z - z**2/3 + z**3/15 - z**4/84 + z**5/540 - \
  350. 2*sqrt(pi)*sqrt(z) + O(z**6)
  351. assert E1(z).series(z) == -EulerGamma - log(z) + z - \
  352. z**2/4 + z**3/18 - z**4/96 + z**5/600 + O(z**6)
  353. assert expint(4, z).series(z) == Rational(1, 3) - z/2 + z**2/2 + \
  354. z**3*(log(z)/6 - Rational(11, 36) + EulerGamma/6 - I*pi/6) - z**4/24 + \
  355. z**5/240 + O(z**6)
  356. assert expint(n, x).series(x, oo, n=3) == \
  357. (n*(n + 1)/x**2 - n/x + 1 + O(x**(-3), (x, oo)))*exp(-x)/x
  358. assert expint(z, y).series(z, 0, 2) == exp(-y)/y - z*meijerg(((), (1, 1)),
  359. ((0, 0, 1), ()), y)/y + O(z**2)
  360. raises(ArgumentIndexError, lambda: expint(x, y).fdiff(3))
  361. neg = Symbol('neg', negative=True)
  362. assert Ei(neg).rewrite(Si) == Shi(neg) + Chi(neg) - I*pi
  363. def test__eis():
  364. assert _eis(z).diff(z) == -_eis(z) + 1/z
  365. assert _eis(1/z).series(z) == \
  366. z + z**2 + 2*z**3 + 6*z**4 + 24*z**5 + O(z**6)
  367. assert Ei(z).rewrite('tractable') == exp(z)*_eis(z)
  368. assert li(z).rewrite('tractable') == z*_eis(log(z))
  369. assert _eis(z).rewrite('intractable') == exp(-z)*Ei(z)
  370. assert expand(li(z).rewrite('tractable').diff(z).rewrite('intractable')) \
  371. == li(z).diff(z)
  372. assert expand(Ei(z).rewrite('tractable').diff(z).rewrite('intractable')) \
  373. == Ei(z).diff(z)
  374. assert _eis(z).series(z, n=3) == EulerGamma + log(z) + z*(-log(z) - \
  375. EulerGamma + 1) + z**2*(log(z)/2 - Rational(3, 4) + EulerGamma/2)\
  376. + O(z**3*log(z))
  377. raises(ArgumentIndexError, lambda: _eis(z).fdiff(2))
  378. def tn_arg(func):
  379. def test(arg, e1, e2):
  380. from sympy.core.random import uniform
  381. v = uniform(1, 5)
  382. v1 = func(arg*x).subs(x, v).n()
  383. v2 = func(e1*v + e2*1e-15).n()
  384. return abs(v1 - v2).n() < 1e-10
  385. return test(exp_polar(I*pi/2), I, 1) and \
  386. test(exp_polar(-I*pi/2), -I, 1) and \
  387. test(exp_polar(I*pi), -1, I) and \
  388. test(exp_polar(-I*pi), -1, -I)
  389. def test_li():
  390. z = Symbol("z")
  391. zr = Symbol("z", real=True)
  392. zp = Symbol("z", positive=True)
  393. zn = Symbol("z", negative=True)
  394. assert li(0) is S.Zero
  395. assert li(1) is -oo
  396. assert li(oo) is oo
  397. assert isinstance(li(z), li)
  398. assert unchanged(li, -zp)
  399. assert unchanged(li, zn)
  400. assert diff(li(z), z) == 1/log(z)
  401. assert conjugate(li(z)) == li(conjugate(z))
  402. assert conjugate(li(-zr)) == li(-zr)
  403. assert unchanged(conjugate, li(-zp))
  404. assert unchanged(conjugate, li(zn))
  405. assert li(z).rewrite(Li) == Li(z) + li(2)
  406. assert li(z).rewrite(Ei) == Ei(log(z))
  407. assert li(z).rewrite(uppergamma) == (-log(1/log(z))/2 - log(-log(z)) +
  408. log(log(z))/2 - expint(1, -log(z)))
  409. assert li(z).rewrite(Si) == (-log(I*log(z)) - log(1/log(z))/2 +
  410. log(log(z))/2 + Ci(I*log(z)) + Shi(log(z)))
  411. assert li(z).rewrite(Ci) == (-log(I*log(z)) - log(1/log(z))/2 +
  412. log(log(z))/2 + Ci(I*log(z)) + Shi(log(z)))
  413. assert li(z).rewrite(Shi) == (-log(1/log(z))/2 + log(log(z))/2 +
  414. Chi(log(z)) - Shi(log(z)))
  415. assert li(z).rewrite(Chi) == (-log(1/log(z))/2 + log(log(z))/2 +
  416. Chi(log(z)) - Shi(log(z)))
  417. assert li(z).rewrite(hyper) ==(log(z)*hyper((1, 1), (2, 2), log(z)) -
  418. log(1/log(z))/2 + log(log(z))/2 + EulerGamma)
  419. assert li(z).rewrite(meijerg) == (-log(1/log(z))/2 - log(-log(z)) + log(log(z))/2 -
  420. meijerg(((), (1,)), ((0, 0), ()), -log(z)))
  421. assert gruntz(1/li(z), z, oo) is S.Zero
  422. assert li(z).series(z) == log(z)**5/600 + log(z)**4/96 + log(z)**3/18 + log(z)**2/4 + \
  423. log(z) + log(log(z)) + EulerGamma
  424. raises(ArgumentIndexError, lambda: li(z).fdiff(2))
  425. def test_Li():
  426. assert Li(2) is S.Zero
  427. assert Li(oo) is oo
  428. assert isinstance(Li(z), Li)
  429. assert diff(Li(z), z) == 1/log(z)
  430. assert gruntz(1/Li(z), z, oo) is S.Zero
  431. assert Li(z).rewrite(li) == li(z) - li(2)
  432. assert Li(z).series(z) == \
  433. log(z)**5/600 + log(z)**4/96 + log(z)**3/18 + log(z)**2/4 + log(z) + log(log(z)) - li(2) + EulerGamma
  434. raises(ArgumentIndexError, lambda: Li(z).fdiff(2))
  435. def test_si():
  436. assert Si(I*x) == I*Shi(x)
  437. assert Shi(I*x) == I*Si(x)
  438. assert Si(-I*x) == -I*Shi(x)
  439. assert Shi(-I*x) == -I*Si(x)
  440. assert Si(-x) == -Si(x)
  441. assert Shi(-x) == -Shi(x)
  442. assert Si(exp_polar(2*pi*I)*x) == Si(x)
  443. assert Si(exp_polar(-2*pi*I)*x) == Si(x)
  444. assert Shi(exp_polar(2*pi*I)*x) == Shi(x)
  445. assert Shi(exp_polar(-2*pi*I)*x) == Shi(x)
  446. assert Si(oo) == pi/2
  447. assert Si(-oo) == -pi/2
  448. assert Shi(oo) is oo
  449. assert Shi(-oo) is -oo
  450. assert mytd(Si(x), sin(x)/x, x)
  451. assert mytd(Shi(x), sinh(x)/x, x)
  452. assert mytn(Si(x), Si(x).rewrite(Ei),
  453. -I*(-Ei(x*exp_polar(-I*pi/2))/2
  454. + Ei(x*exp_polar(I*pi/2))/2 - I*pi) + pi/2, x)
  455. assert mytn(Si(x), Si(x).rewrite(expint),
  456. -I*(-expint(1, x*exp_polar(-I*pi/2))/2 +
  457. expint(1, x*exp_polar(I*pi/2))/2) + pi/2, x)
  458. assert mytn(Shi(x), Shi(x).rewrite(Ei),
  459. Ei(x)/2 - Ei(x*exp_polar(I*pi))/2 + I*pi/2, x)
  460. assert mytn(Shi(x), Shi(x).rewrite(expint),
  461. expint(1, x)/2 - expint(1, x*exp_polar(I*pi))/2 - I*pi/2, x)
  462. assert tn_arg(Si)
  463. assert tn_arg(Shi)
  464. assert Si(x).nseries(x, n=8) == \
  465. x - x**3/18 + x**5/600 - x**7/35280 + O(x**9)
  466. assert Shi(x).nseries(x, n=8) == \
  467. x + x**3/18 + x**5/600 + x**7/35280 + O(x**9)
  468. assert Si(sin(x)).nseries(x, n=5) == x - 2*x**3/9 + 17*x**5/450 + O(x**6)
  469. assert Si(x).nseries(x, 1, n=3) == \
  470. Si(1) + (x - 1)*sin(1) + (x - 1)**2*(-sin(1)/2 + cos(1)/2) + O((x - 1)**3, (x, 1))
  471. assert Si(x).series(x, oo) == pi/2 - (- 6/x**3 + 1/x \
  472. + O(x**(-7), (x, oo)))*sin(x)/x - (24/x**4 - 2/x**2 + 1 \
  473. + O(x**(-7), (x, oo)))*cos(x)/x
  474. t = Symbol('t', Dummy=True)
  475. assert Si(x).rewrite(sinc) == Integral(sinc(t), (t, 0, x))
  476. assert limit(Shi(x), x, S.Infinity) == S.Infinity
  477. assert limit(Shi(x), x, S.NegativeInfinity) == S.NegativeInfinity
  478. def test_ci():
  479. m1 = exp_polar(I*pi)
  480. m1_ = exp_polar(-I*pi)
  481. pI = exp_polar(I*pi/2)
  482. mI = exp_polar(-I*pi/2)
  483. assert Ci(m1*x) == Ci(x) + I*pi
  484. assert Ci(m1_*x) == Ci(x) - I*pi
  485. assert Ci(pI*x) == Chi(x) + I*pi/2
  486. assert Ci(mI*x) == Chi(x) - I*pi/2
  487. assert Chi(m1*x) == Chi(x) + I*pi
  488. assert Chi(m1_*x) == Chi(x) - I*pi
  489. assert Chi(pI*x) == Ci(x) + I*pi/2
  490. assert Chi(mI*x) == Ci(x) - I*pi/2
  491. assert Ci(exp_polar(2*I*pi)*x) == Ci(x) + 2*I*pi
  492. assert Chi(exp_polar(-2*I*pi)*x) == Chi(x) - 2*I*pi
  493. assert Chi(exp_polar(2*I*pi)*x) == Chi(x) + 2*I*pi
  494. assert Ci(exp_polar(-2*I*pi)*x) == Ci(x) - 2*I*pi
  495. assert Ci(oo) is S.Zero
  496. assert Ci(-oo) == I*pi
  497. assert Chi(oo) is oo
  498. assert Chi(-oo) is oo
  499. assert mytd(Ci(x), cos(x)/x, x)
  500. assert mytd(Chi(x), cosh(x)/x, x)
  501. assert mytn(Ci(x), Ci(x).rewrite(Ei),
  502. Ei(x*exp_polar(-I*pi/2))/2 + Ei(x*exp_polar(I*pi/2))/2, x)
  503. assert mytn(Chi(x), Chi(x).rewrite(Ei),
  504. Ei(x)/2 + Ei(x*exp_polar(I*pi))/2 - I*pi/2, x)
  505. assert tn_arg(Ci)
  506. assert tn_arg(Chi)
  507. assert Ci(x).nseries(x, n=4) == \
  508. EulerGamma + log(x) - x**2/4 + x**4/96 + O(x**5)
  509. assert Chi(x).nseries(x, n=4) == \
  510. EulerGamma + log(x) + x**2/4 + x**4/96 + O(x**5)
  511. assert Ci(x).series(x, oo) == -cos(x)*(-6/x**3 + 1/x \
  512. + O(x**(-7), (x, oo)))/x + (24/x**4 - 2/x**2 + 1 \
  513. + O(x**(-7), (x, oo)))*sin(x)/x
  514. assert limit(log(x) - Ci(2*x), x, 0) == -log(2) - EulerGamma
  515. assert Ci(x).rewrite(uppergamma) == -expint(1, x*exp_polar(-I*pi/2))/2 -\
  516. expint(1, x*exp_polar(I*pi/2))/2
  517. assert Ci(x).rewrite(expint) == -expint(1, x*exp_polar(-I*pi/2))/2 -\
  518. expint(1, x*exp_polar(I*pi/2))/2
  519. raises(ArgumentIndexError, lambda: Ci(x).fdiff(2))
  520. def test_fresnel():
  521. assert fresnels(0) is S.Zero
  522. assert fresnels(oo) is S.Half
  523. assert fresnels(-oo) == Rational(-1, 2)
  524. assert fresnels(I*oo) == -I*S.Half
  525. assert unchanged(fresnels, z)
  526. assert fresnels(-z) == -fresnels(z)
  527. assert fresnels(I*z) == -I*fresnels(z)
  528. assert fresnels(-I*z) == I*fresnels(z)
  529. assert conjugate(fresnels(z)) == fresnels(conjugate(z))
  530. assert fresnels(z).diff(z) == sin(pi*z**2/2)
  531. assert fresnels(z).rewrite(erf) == (S.One + I)/4 * (
  532. erf((S.One + I)/2*sqrt(pi)*z) - I*erf((S.One - I)/2*sqrt(pi)*z))
  533. assert fresnels(z).rewrite(hyper) == \
  534. pi*z**3/6 * hyper([Rational(3, 4)], [Rational(3, 2), Rational(7, 4)], -pi**2*z**4/16)
  535. assert fresnels(z).series(z, n=15) == \
  536. pi*z**3/6 - pi**3*z**7/336 + pi**5*z**11/42240 + O(z**15)
  537. assert fresnels(w).is_extended_real is True
  538. assert fresnels(w).is_finite is True
  539. assert fresnels(z).is_extended_real is None
  540. assert fresnels(z).is_finite is None
  541. assert fresnels(z).as_real_imag() == (fresnels(re(z) - I*im(z))/2 +
  542. fresnels(re(z) + I*im(z))/2,
  543. -I*(-fresnels(re(z) - I*im(z)) + fresnels(re(z) + I*im(z)))/2)
  544. assert fresnels(z).as_real_imag(deep=False) == (fresnels(re(z) - I*im(z))/2 +
  545. fresnels(re(z) + I*im(z))/2,
  546. -I*(-fresnels(re(z) - I*im(z)) + fresnels(re(z) + I*im(z)))/2)
  547. assert fresnels(w).as_real_imag() == (fresnels(w), 0)
  548. assert fresnels(w).as_real_imag(deep=True) == (fresnels(w), 0)
  549. assert fresnels(2 + 3*I).as_real_imag() == (
  550. fresnels(2 + 3*I)/2 + fresnels(2 - 3*I)/2,
  551. -I*(fresnels(2 + 3*I) - fresnels(2 - 3*I))/2
  552. )
  553. assert expand_func(integrate(fresnels(z), z)) == \
  554. z*fresnels(z) + cos(pi*z**2/2)/pi
  555. assert fresnels(z).rewrite(meijerg) == sqrt(2)*pi*z**Rational(9, 4) * \
  556. meijerg(((), (1,)), ((Rational(3, 4),),
  557. (Rational(1, 4), 0)), -pi**2*z**4/16)/(2*(-z)**Rational(3, 4)*(z**2)**Rational(3, 4))
  558. assert fresnelc(0) is S.Zero
  559. assert fresnelc(oo) == S.Half
  560. assert fresnelc(-oo) == Rational(-1, 2)
  561. assert fresnelc(I*oo) == I*S.Half
  562. assert unchanged(fresnelc, z)
  563. assert fresnelc(-z) == -fresnelc(z)
  564. assert fresnelc(I*z) == I*fresnelc(z)
  565. assert fresnelc(-I*z) == -I*fresnelc(z)
  566. assert conjugate(fresnelc(z)) == fresnelc(conjugate(z))
  567. assert fresnelc(z).diff(z) == cos(pi*z**2/2)
  568. assert fresnelc(z).rewrite(erf) == (S.One - I)/4 * (
  569. erf((S.One + I)/2*sqrt(pi)*z) + I*erf((S.One - I)/2*sqrt(pi)*z))
  570. assert fresnelc(z).rewrite(hyper) == \
  571. z * hyper([Rational(1, 4)], [S.Half, Rational(5, 4)], -pi**2*z**4/16)
  572. assert fresnelc(w).is_extended_real is True
  573. assert fresnelc(z).as_real_imag() == \
  574. (fresnelc(re(z) - I*im(z))/2 + fresnelc(re(z) + I*im(z))/2,
  575. -I*(-fresnelc(re(z) - I*im(z)) + fresnelc(re(z) + I*im(z)))/2)
  576. assert fresnelc(z).as_real_imag(deep=False) == \
  577. (fresnelc(re(z) - I*im(z))/2 + fresnelc(re(z) + I*im(z))/2,
  578. -I*(-fresnelc(re(z) - I*im(z)) + fresnelc(re(z) + I*im(z)))/2)
  579. assert fresnelc(2 + 3*I).as_real_imag() == (
  580. fresnelc(2 - 3*I)/2 + fresnelc(2 + 3*I)/2,
  581. -I*(fresnelc(2 + 3*I) - fresnelc(2 - 3*I))/2
  582. )
  583. assert expand_func(integrate(fresnelc(z), z)) == \
  584. z*fresnelc(z) - sin(pi*z**2/2)/pi
  585. assert fresnelc(z).rewrite(meijerg) == sqrt(2)*pi*z**Rational(3, 4) * \
  586. meijerg(((), (1,)), ((Rational(1, 4),),
  587. (Rational(3, 4), 0)), -pi**2*z**4/16)/(2*(-z)**Rational(1, 4)*(z**2)**Rational(1, 4))
  588. from sympy.core.random import verify_numerically
  589. verify_numerically(re(fresnels(z)), fresnels(z).as_real_imag()[0], z)
  590. verify_numerically(im(fresnels(z)), fresnels(z).as_real_imag()[1], z)
  591. verify_numerically(fresnels(z), fresnels(z).rewrite(hyper), z)
  592. verify_numerically(fresnels(z), fresnels(z).rewrite(meijerg), z)
  593. verify_numerically(re(fresnelc(z)), fresnelc(z).as_real_imag()[0], z)
  594. verify_numerically(im(fresnelc(z)), fresnelc(z).as_real_imag()[1], z)
  595. verify_numerically(fresnelc(z), fresnelc(z).rewrite(hyper), z)
  596. verify_numerically(fresnelc(z), fresnelc(z).rewrite(meijerg), z)
  597. raises(ArgumentIndexError, lambda: fresnels(z).fdiff(2))
  598. raises(ArgumentIndexError, lambda: fresnelc(z).fdiff(2))
  599. assert fresnels(x).taylor_term(-1, x) is S.Zero
  600. assert fresnelc(x).taylor_term(-1, x) is S.Zero
  601. assert fresnelc(x).taylor_term(1, x) == -pi**2*x**5/40
  602. def test_fresnel_series():
  603. assert fresnelc(z).series(z, n=15) == \
  604. z - pi**2*z**5/40 + pi**4*z**9/3456 - pi**6*z**13/599040 + O(z**15)
  605. # issues 6510, 10102
  606. fs = (S.Half - sin(pi*z**2/2)/(pi**2*z**3)
  607. + (-1/(pi*z) + 3/(pi**3*z**5))*cos(pi*z**2/2))
  608. fc = (S.Half - cos(pi*z**2/2)/(pi**2*z**3)
  609. + (1/(pi*z) - 3/(pi**3*z**5))*sin(pi*z**2/2))
  610. assert fresnels(z).series(z, oo) == fs + O(z**(-6), (z, oo))
  611. assert fresnelc(z).series(z, oo) == fc + O(z**(-6), (z, oo))
  612. assert (fresnels(z).series(z, -oo) + fs.subs(z, -z)).expand().is_Order
  613. assert (fresnelc(z).series(z, -oo) + fc.subs(z, -z)).expand().is_Order
  614. assert (fresnels(1/z).series(z) - fs.subs(z, 1/z)).expand().is_Order
  615. assert (fresnelc(1/z).series(z) - fc.subs(z, 1/z)).expand().is_Order
  616. assert ((2*fresnels(3*z)).series(z, oo) - 2*fs.subs(z, 3*z)).expand().is_Order
  617. assert ((3*fresnelc(2*z)).series(z, oo) - 3*fc.subs(z, 2*z)).expand().is_Order