test_data.py 28 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617
  1. import os
  2. import numpy as np
  3. from numpy.testing import suppress_warnings
  4. import pytest
  5. from scipy.special import (
  6. lpn, lpmn, lpmv, lqn, lqmn, sph_harm, eval_legendre, eval_hermite,
  7. eval_laguerre, eval_genlaguerre, binom, cbrt, expm1, log1p, zeta,
  8. jn, jv, jvp, yn, yv, yvp, iv, ivp, kn, kv, kvp,
  9. gamma, gammaln, gammainc, gammaincc, gammaincinv, gammainccinv, digamma,
  10. beta, betainc, betaincinv, poch,
  11. ellipe, ellipeinc, ellipk, ellipkm1, ellipkinc, ellipj,
  12. elliprc, elliprd, elliprf, elliprg, elliprj,
  13. erf, erfc, erfinv, erfcinv, exp1, expi, expn,
  14. bdtrik, btdtr, btdtri, btdtria, btdtrib, chndtr, gdtr, gdtrc, gdtrix, gdtrib,
  15. nbdtrik, pdtrik, owens_t,
  16. mathieu_a, mathieu_b, mathieu_cem, mathieu_sem, mathieu_modcem1,
  17. mathieu_modsem1, mathieu_modcem2, mathieu_modsem2,
  18. ellip_harm, ellip_harm_2, spherical_jn, spherical_yn, wright_bessel
  19. )
  20. from scipy.integrate import IntegrationWarning
  21. from scipy.special._testutils import FuncData
  22. DATASETS_BOOST = np.load(os.path.join(os.path.dirname(__file__),
  23. "data", "boost.npz"))
  24. DATASETS_GSL = np.load(os.path.join(os.path.dirname(__file__),
  25. "data", "gsl.npz"))
  26. DATASETS_LOCAL = np.load(os.path.join(os.path.dirname(__file__),
  27. "data", "local.npz"))
  28. def data(func, dataname, *a, **kw):
  29. kw.setdefault('dataname', dataname)
  30. return FuncData(func, DATASETS_BOOST[dataname], *a, **kw)
  31. def data_gsl(func, dataname, *a, **kw):
  32. kw.setdefault('dataname', dataname)
  33. return FuncData(func, DATASETS_GSL[dataname], *a, **kw)
  34. def data_local(func, dataname, *a, **kw):
  35. kw.setdefault('dataname', dataname)
  36. return FuncData(func, DATASETS_LOCAL[dataname], *a, **kw)
  37. def ellipk_(k):
  38. return ellipk(k*k)
  39. def ellipkinc_(f, k):
  40. return ellipkinc(f, k*k)
  41. def ellipe_(k):
  42. return ellipe(k*k)
  43. def ellipeinc_(f, k):
  44. return ellipeinc(f, k*k)
  45. def ellipj_(k):
  46. return ellipj(k*k)
  47. def zeta_(x):
  48. return zeta(x, 1.)
  49. def assoc_legendre_p_boost_(nu, mu, x):
  50. # the boost test data is for integer orders only
  51. return lpmv(mu, nu.astype(int), x)
  52. def legendre_p_via_assoc_(nu, x):
  53. return lpmv(0, nu, x)
  54. def lpn_(n, x):
  55. return lpn(n.astype('l'), x)[0][-1]
  56. def lqn_(n, x):
  57. return lqn(n.astype('l'), x)[0][-1]
  58. def legendre_p_via_lpmn(n, x):
  59. return lpmn(0, n, x)[0][0,-1]
  60. def legendre_q_via_lqmn(n, x):
  61. return lqmn(0, n, x)[0][0,-1]
  62. def mathieu_ce_rad(m, q, x):
  63. return mathieu_cem(m, q, x*180/np.pi)[0]
  64. def mathieu_se_rad(m, q, x):
  65. return mathieu_sem(m, q, x*180/np.pi)[0]
  66. def mathieu_mc1_scaled(m, q, x):
  67. # GSL follows a different normalization.
  68. # We follow Abramowitz & Stegun, they apparently something else.
  69. return mathieu_modcem1(m, q, x)[0] * np.sqrt(np.pi/2)
  70. def mathieu_ms1_scaled(m, q, x):
  71. return mathieu_modsem1(m, q, x)[0] * np.sqrt(np.pi/2)
  72. def mathieu_mc2_scaled(m, q, x):
  73. return mathieu_modcem2(m, q, x)[0] * np.sqrt(np.pi/2)
  74. def mathieu_ms2_scaled(m, q, x):
  75. return mathieu_modsem2(m, q, x)[0] * np.sqrt(np.pi/2)
  76. def eval_legendre_ld(n, x):
  77. return eval_legendre(n.astype('l'), x)
  78. def eval_legendre_dd(n, x):
  79. return eval_legendre(n.astype('d'), x)
  80. def eval_hermite_ld(n, x):
  81. return eval_hermite(n.astype('l'), x)
  82. def eval_laguerre_ld(n, x):
  83. return eval_laguerre(n.astype('l'), x)
  84. def eval_laguerre_dd(n, x):
  85. return eval_laguerre(n.astype('d'), x)
  86. def eval_genlaguerre_ldd(n, a, x):
  87. return eval_genlaguerre(n.astype('l'), a, x)
  88. def eval_genlaguerre_ddd(n, a, x):
  89. return eval_genlaguerre(n.astype('d'), a, x)
  90. def bdtrik_comp(y, n, p):
  91. return bdtrik(1-y, n, p)
  92. def btdtri_comp(a, b, p):
  93. return btdtri(a, b, 1-p)
  94. def btdtria_comp(p, b, x):
  95. return btdtria(1-p, b, x)
  96. def btdtrib_comp(a, p, x):
  97. return btdtrib(a, 1-p, x)
  98. def gdtr_(p, x):
  99. return gdtr(1.0, p, x)
  100. def gdtrc_(p, x):
  101. return gdtrc(1.0, p, x)
  102. def gdtrix_(b, p):
  103. return gdtrix(1.0, b, p)
  104. def gdtrix_comp(b, p):
  105. return gdtrix(1.0, b, 1-p)
  106. def gdtrib_(p, x):
  107. return gdtrib(1.0, p, x)
  108. def gdtrib_comp(p, x):
  109. return gdtrib(1.0, 1-p, x)
  110. def nbdtrik_comp(y, n, p):
  111. return nbdtrik(1-y, n, p)
  112. def pdtrik_comp(p, m):
  113. return pdtrik(1-p, m)
  114. def poch_(z, m):
  115. return 1.0 / poch(z, m)
  116. def poch_minus(z, m):
  117. return 1.0 / poch(z, -m)
  118. def spherical_jn_(n, x):
  119. return spherical_jn(n.astype('l'), x)
  120. def spherical_yn_(n, x):
  121. return spherical_yn(n.astype('l'), x)
  122. def sph_harm_(m, n, theta, phi):
  123. y = sph_harm(m, n, theta, phi)
  124. return (y.real, y.imag)
  125. def cexpm1(x, y):
  126. z = expm1(x + 1j*y)
  127. return z.real, z.imag
  128. def clog1p(x, y):
  129. z = log1p(x + 1j*y)
  130. return z.real, z.imag
  131. BOOST_TESTS = [
  132. data(assoc_legendre_p_boost_, 'assoc_legendre_p_ipp-assoc_legendre_p', (0,1,2), 3, rtol=1e-11),
  133. data(legendre_p_via_assoc_, 'legendre_p_ipp-legendre_p', (0,1), 2, rtol=1e-11),
  134. data(legendre_p_via_assoc_, 'legendre_p_large_ipp-legendre_p_large', (0,1), 2, rtol=9.6e-14),
  135. data(legendre_p_via_lpmn, 'legendre_p_ipp-legendre_p', (0,1), 2, rtol=5e-14, vectorized=False),
  136. data(legendre_p_via_lpmn, 'legendre_p_large_ipp-legendre_p_large', (0,1), 2, rtol=9.6e-14, vectorized=False),
  137. data(lpn_, 'legendre_p_ipp-legendre_p', (0,1), 2, rtol=5e-14, vectorized=False),
  138. data(lpn_, 'legendre_p_large_ipp-legendre_p_large', (0,1), 2, rtol=3e-13, vectorized=False),
  139. data(eval_legendre_ld, 'legendre_p_ipp-legendre_p', (0,1), 2, rtol=6e-14),
  140. data(eval_legendre_ld, 'legendre_p_large_ipp-legendre_p_large', (0,1), 2, rtol=2e-13),
  141. data(eval_legendre_dd, 'legendre_p_ipp-legendre_p', (0,1), 2, rtol=2e-14),
  142. data(eval_legendre_dd, 'legendre_p_large_ipp-legendre_p_large', (0,1), 2, rtol=2e-13),
  143. data(lqn_, 'legendre_p_ipp-legendre_p', (0,1), 3, rtol=2e-14, vectorized=False),
  144. data(lqn_, 'legendre_p_large_ipp-legendre_p_large', (0,1), 3, rtol=2e-12, vectorized=False),
  145. data(legendre_q_via_lqmn, 'legendre_p_ipp-legendre_p', (0,1), 3, rtol=2e-14, vectorized=False),
  146. data(legendre_q_via_lqmn, 'legendre_p_large_ipp-legendre_p_large', (0,1), 3, rtol=2e-12, vectorized=False),
  147. data(beta, 'beta_exp_data_ipp-beta_exp_data', (0,1), 2, rtol=1e-13),
  148. data(beta, 'beta_exp_data_ipp-beta_exp_data', (0,1), 2, rtol=1e-13),
  149. data(beta, 'beta_med_data_ipp-beta_med_data', (0,1), 2, rtol=5e-13),
  150. data(betainc, 'ibeta_small_data_ipp-ibeta_small_data', (0,1,2), 5, rtol=6e-15),
  151. data(betainc, 'ibeta_data_ipp-ibeta_data', (0,1,2), 5, rtol=5e-13),
  152. data(betainc, 'ibeta_int_data_ipp-ibeta_int_data', (0,1,2), 5, rtol=2e-14),
  153. data(betainc, 'ibeta_large_data_ipp-ibeta_large_data', (0,1,2), 5, rtol=4e-10),
  154. data(betaincinv, 'ibeta_inv_data_ipp-ibeta_inv_data', (0,1,2), 3, rtol=1e-5),
  155. data(btdtr, 'ibeta_small_data_ipp-ibeta_small_data', (0,1,2), 5, rtol=6e-15),
  156. data(btdtr, 'ibeta_data_ipp-ibeta_data', (0,1,2), 5, rtol=4e-13),
  157. data(btdtr, 'ibeta_int_data_ipp-ibeta_int_data', (0,1,2), 5, rtol=2e-14),
  158. data(btdtr, 'ibeta_large_data_ipp-ibeta_large_data', (0,1,2), 5, rtol=4e-10),
  159. data(btdtri, 'ibeta_inv_data_ipp-ibeta_inv_data', (0,1,2), 3, rtol=1e-5),
  160. data(btdtri_comp, 'ibeta_inv_data_ipp-ibeta_inv_data', (0,1,2), 4, rtol=8e-7),
  161. data(btdtria, 'ibeta_inva_data_ipp-ibeta_inva_data', (2,0,1), 3, rtol=5e-9),
  162. data(btdtria_comp, 'ibeta_inva_data_ipp-ibeta_inva_data', (2,0,1), 4, rtol=5e-9),
  163. data(btdtrib, 'ibeta_inva_data_ipp-ibeta_inva_data', (0,2,1), 5, rtol=5e-9),
  164. data(btdtrib_comp, 'ibeta_inva_data_ipp-ibeta_inva_data', (0,2,1), 6, rtol=5e-9),
  165. data(binom, 'binomial_data_ipp-binomial_data', (0,1), 2, rtol=1e-13),
  166. data(binom, 'binomial_large_data_ipp-binomial_large_data', (0,1), 2, rtol=5e-13),
  167. data(bdtrik, 'binomial_quantile_ipp-binomial_quantile_data', (2,0,1), 3, rtol=5e-9),
  168. data(bdtrik_comp, 'binomial_quantile_ipp-binomial_quantile_data', (2,0,1), 4, rtol=5e-9),
  169. data(nbdtrik, 'negative_binomial_quantile_ipp-negative_binomial_quantile_data', (2,0,1), 3, rtol=4e-9),
  170. data(nbdtrik_comp, 'negative_binomial_quantile_ipp-negative_binomial_quantile_data', (2,0,1), 4, rtol=4e-9),
  171. data(pdtrik, 'poisson_quantile_ipp-poisson_quantile_data', (1,0), 2, rtol=3e-9),
  172. data(pdtrik_comp, 'poisson_quantile_ipp-poisson_quantile_data', (1,0), 3, rtol=4e-9),
  173. data(cbrt, 'cbrt_data_ipp-cbrt_data', 1, 0),
  174. data(digamma, 'digamma_data_ipp-digamma_data', 0, 1),
  175. data(digamma, 'digamma_data_ipp-digamma_data', 0j, 1),
  176. data(digamma, 'digamma_neg_data_ipp-digamma_neg_data', 0, 1, rtol=2e-13),
  177. data(digamma, 'digamma_neg_data_ipp-digamma_neg_data', 0j, 1, rtol=1e-13),
  178. data(digamma, 'digamma_root_data_ipp-digamma_root_data', 0, 1, rtol=1e-15),
  179. data(digamma, 'digamma_root_data_ipp-digamma_root_data', 0j, 1, rtol=1e-15),
  180. data(digamma, 'digamma_small_data_ipp-digamma_small_data', 0, 1, rtol=1e-15),
  181. data(digamma, 'digamma_small_data_ipp-digamma_small_data', 0j, 1, rtol=1e-14),
  182. data(ellipk_, 'ellint_k_data_ipp-ellint_k_data', 0, 1),
  183. data(ellipkinc_, 'ellint_f_data_ipp-ellint_f_data', (0,1), 2, rtol=1e-14),
  184. data(ellipe_, 'ellint_e_data_ipp-ellint_e_data', 0, 1),
  185. data(ellipeinc_, 'ellint_e2_data_ipp-ellint_e2_data', (0,1), 2, rtol=1e-14),
  186. data(erf, 'erf_data_ipp-erf_data', 0, 1),
  187. data(erf, 'erf_data_ipp-erf_data', 0j, 1, rtol=1e-13),
  188. data(erfc, 'erf_data_ipp-erf_data', 0, 2, rtol=6e-15),
  189. data(erf, 'erf_large_data_ipp-erf_large_data', 0, 1),
  190. data(erf, 'erf_large_data_ipp-erf_large_data', 0j, 1),
  191. data(erfc, 'erf_large_data_ipp-erf_large_data', 0, 2, rtol=4e-14),
  192. data(erf, 'erf_small_data_ipp-erf_small_data', 0, 1),
  193. data(erf, 'erf_small_data_ipp-erf_small_data', 0j, 1, rtol=1e-13),
  194. data(erfc, 'erf_small_data_ipp-erf_small_data', 0, 2),
  195. data(erfinv, 'erf_inv_data_ipp-erf_inv_data', 0, 1),
  196. data(erfcinv, 'erfc_inv_data_ipp-erfc_inv_data', 0, 1),
  197. data(erfcinv, 'erfc_inv_big_data_ipp-erfc_inv_big_data', 0, 1, param_filter=(lambda s: s > 0)),
  198. data(exp1, 'expint_1_data_ipp-expint_1_data', 1, 2, rtol=1e-13),
  199. data(exp1, 'expint_1_data_ipp-expint_1_data', 1j, 2, rtol=5e-9),
  200. data(expi, 'expinti_data_ipp-expinti_data', 0, 1, rtol=1e-13),
  201. data(expi, 'expinti_data_double_ipp-expinti_data_double', 0, 1, rtol=1e-13),
  202. data(expi, 'expinti_data_long_ipp-expinti_data_long', 0, 1),
  203. data(expn, 'expint_small_data_ipp-expint_small_data', (0,1), 2),
  204. data(expn, 'expint_data_ipp-expint_data', (0,1), 2, rtol=1e-14),
  205. data(gamma, 'test_gamma_data_ipp-near_0', 0, 1),
  206. data(gamma, 'test_gamma_data_ipp-near_1', 0, 1),
  207. data(gamma, 'test_gamma_data_ipp-near_2', 0, 1),
  208. data(gamma, 'test_gamma_data_ipp-near_m10', 0, 1),
  209. data(gamma, 'test_gamma_data_ipp-near_m55', 0, 1, rtol=7e-12),
  210. data(gamma, 'test_gamma_data_ipp-factorials', 0, 1, rtol=4e-14),
  211. data(gamma, 'test_gamma_data_ipp-near_0', 0j, 1, rtol=2e-9),
  212. data(gamma, 'test_gamma_data_ipp-near_1', 0j, 1, rtol=2e-9),
  213. data(gamma, 'test_gamma_data_ipp-near_2', 0j, 1, rtol=2e-9),
  214. data(gamma, 'test_gamma_data_ipp-near_m10', 0j, 1, rtol=2e-9),
  215. data(gamma, 'test_gamma_data_ipp-near_m55', 0j, 1, rtol=2e-9),
  216. data(gamma, 'test_gamma_data_ipp-factorials', 0j, 1, rtol=2e-13),
  217. data(gammaln, 'test_gamma_data_ipp-near_0', 0, 2, rtol=5e-11),
  218. data(gammaln, 'test_gamma_data_ipp-near_1', 0, 2, rtol=5e-11),
  219. data(gammaln, 'test_gamma_data_ipp-near_2', 0, 2, rtol=2e-10),
  220. data(gammaln, 'test_gamma_data_ipp-near_m10', 0, 2, rtol=5e-11),
  221. data(gammaln, 'test_gamma_data_ipp-near_m55', 0, 2, rtol=5e-11),
  222. data(gammaln, 'test_gamma_data_ipp-factorials', 0, 2),
  223. data(gammainc, 'igamma_small_data_ipp-igamma_small_data', (0,1), 5, rtol=5e-15),
  224. data(gammainc, 'igamma_med_data_ipp-igamma_med_data', (0,1), 5, rtol=2e-13),
  225. data(gammainc, 'igamma_int_data_ipp-igamma_int_data', (0,1), 5, rtol=2e-13),
  226. data(gammainc, 'igamma_big_data_ipp-igamma_big_data', (0,1), 5, rtol=1e-12),
  227. data(gdtr_, 'igamma_small_data_ipp-igamma_small_data', (0,1), 5, rtol=1e-13),
  228. data(gdtr_, 'igamma_med_data_ipp-igamma_med_data', (0,1), 5, rtol=2e-13),
  229. data(gdtr_, 'igamma_int_data_ipp-igamma_int_data', (0,1), 5, rtol=2e-13),
  230. data(gdtr_, 'igamma_big_data_ipp-igamma_big_data', (0,1), 5, rtol=2e-9),
  231. data(gammaincc, 'igamma_small_data_ipp-igamma_small_data', (0,1), 3, rtol=1e-13),
  232. data(gammaincc, 'igamma_med_data_ipp-igamma_med_data', (0,1), 3, rtol=2e-13),
  233. data(gammaincc, 'igamma_int_data_ipp-igamma_int_data', (0,1), 3, rtol=4e-14),
  234. data(gammaincc, 'igamma_big_data_ipp-igamma_big_data', (0,1), 3, rtol=1e-11),
  235. data(gdtrc_, 'igamma_small_data_ipp-igamma_small_data', (0,1), 3, rtol=1e-13),
  236. data(gdtrc_, 'igamma_med_data_ipp-igamma_med_data', (0,1), 3, rtol=2e-13),
  237. data(gdtrc_, 'igamma_int_data_ipp-igamma_int_data', (0,1), 3, rtol=4e-14),
  238. data(gdtrc_, 'igamma_big_data_ipp-igamma_big_data', (0,1), 3, rtol=1e-11),
  239. data(gdtrib_, 'igamma_inva_data_ipp-igamma_inva_data', (1,0), 2, rtol=5e-9),
  240. data(gdtrib_comp, 'igamma_inva_data_ipp-igamma_inva_data', (1,0), 3, rtol=5e-9),
  241. data(poch_, 'tgamma_delta_ratio_data_ipp-tgamma_delta_ratio_data', (0,1), 2, rtol=2e-13),
  242. data(poch_, 'tgamma_delta_ratio_int_ipp-tgamma_delta_ratio_int', (0,1), 2,),
  243. data(poch_, 'tgamma_delta_ratio_int2_ipp-tgamma_delta_ratio_int2', (0,1), 2,),
  244. data(poch_minus, 'tgamma_delta_ratio_data_ipp-tgamma_delta_ratio_data', (0,1), 3, rtol=2e-13),
  245. data(poch_minus, 'tgamma_delta_ratio_int_ipp-tgamma_delta_ratio_int', (0,1), 3),
  246. data(poch_minus, 'tgamma_delta_ratio_int2_ipp-tgamma_delta_ratio_int2', (0,1), 3),
  247. data(eval_hermite_ld, 'hermite_ipp-hermite', (0,1), 2, rtol=2e-14),
  248. data(eval_laguerre_ld, 'laguerre2_ipp-laguerre2', (0,1), 2, rtol=7e-12),
  249. data(eval_laguerre_dd, 'laguerre2_ipp-laguerre2', (0,1), 2, knownfailure='hyp2f1 insufficiently accurate.'),
  250. data(eval_genlaguerre_ldd, 'laguerre3_ipp-laguerre3', (0,1,2), 3, rtol=2e-13),
  251. data(eval_genlaguerre_ddd, 'laguerre3_ipp-laguerre3', (0,1,2), 3, knownfailure='hyp2f1 insufficiently accurate.'),
  252. data(log1p, 'log1p_expm1_data_ipp-log1p_expm1_data', 0, 1),
  253. data(expm1, 'log1p_expm1_data_ipp-log1p_expm1_data', 0, 2),
  254. data(iv, 'bessel_i_data_ipp-bessel_i_data', (0,1), 2, rtol=1e-12),
  255. data(iv, 'bessel_i_data_ipp-bessel_i_data', (0,1j), 2, rtol=2e-10, atol=1e-306),
  256. data(iv, 'bessel_i_int_data_ipp-bessel_i_int_data', (0,1), 2, rtol=1e-9),
  257. data(iv, 'bessel_i_int_data_ipp-bessel_i_int_data', (0,1j), 2, rtol=2e-10),
  258. data(ivp, 'bessel_i_prime_int_data_ipp-bessel_i_prime_int_data', (0,1), 2, rtol=1.2e-13),
  259. data(ivp, 'bessel_i_prime_int_data_ipp-bessel_i_prime_int_data', (0,1j), 2, rtol=1.2e-13, atol=1e-300),
  260. data(jn, 'bessel_j_int_data_ipp-bessel_j_int_data', (0,1), 2, rtol=1e-12),
  261. data(jn, 'bessel_j_int_data_ipp-bessel_j_int_data', (0,1j), 2, rtol=1e-12),
  262. data(jn, 'bessel_j_large_data_ipp-bessel_j_large_data', (0,1), 2, rtol=6e-11),
  263. data(jn, 'bessel_j_large_data_ipp-bessel_j_large_data', (0,1j), 2, rtol=6e-11),
  264. data(jv, 'bessel_j_int_data_ipp-bessel_j_int_data', (0,1), 2, rtol=1e-12),
  265. data(jv, 'bessel_j_int_data_ipp-bessel_j_int_data', (0,1j), 2, rtol=1e-12),
  266. data(jv, 'bessel_j_data_ipp-bessel_j_data', (0,1), 2, rtol=1e-12),
  267. data(jv, 'bessel_j_data_ipp-bessel_j_data', (0,1j), 2, rtol=1e-12),
  268. data(jvp, 'bessel_j_prime_int_data_ipp-bessel_j_prime_int_data', (0,1), 2, rtol=1e-13),
  269. data(jvp, 'bessel_j_prime_int_data_ipp-bessel_j_prime_int_data', (0,1j), 2, rtol=1e-13),
  270. data(jvp, 'bessel_j_prime_large_data_ipp-bessel_j_prime_large_data', (0,1), 2, rtol=1e-11),
  271. data(jvp, 'bessel_j_prime_large_data_ipp-bessel_j_prime_large_data', (0,1j), 2, rtol=1e-11),
  272. data(kn, 'bessel_k_int_data_ipp-bessel_k_int_data', (0,1), 2, rtol=1e-12),
  273. data(kv, 'bessel_k_int_data_ipp-bessel_k_int_data', (0,1), 2, rtol=1e-12),
  274. data(kv, 'bessel_k_int_data_ipp-bessel_k_int_data', (0,1j), 2, rtol=1e-12),
  275. data(kv, 'bessel_k_data_ipp-bessel_k_data', (0,1), 2, rtol=1e-12),
  276. data(kv, 'bessel_k_data_ipp-bessel_k_data', (0,1j), 2, rtol=1e-12),
  277. data(kvp, 'bessel_k_prime_int_data_ipp-bessel_k_prime_int_data', (0,1), 2, rtol=3e-14),
  278. data(kvp, 'bessel_k_prime_int_data_ipp-bessel_k_prime_int_data', (0,1j), 2, rtol=3e-14),
  279. data(kvp, 'bessel_k_prime_data_ipp-bessel_k_prime_data', (0,1), 2, rtol=7e-14),
  280. data(kvp, 'bessel_k_prime_data_ipp-bessel_k_prime_data', (0,1j), 2, rtol=7e-14),
  281. data(yn, 'bessel_y01_data_ipp-bessel_y01_data', (0,1), 2, rtol=1e-12),
  282. data(yn, 'bessel_yn_data_ipp-bessel_yn_data', (0,1), 2, rtol=1e-12),
  283. data(yv, 'bessel_yn_data_ipp-bessel_yn_data', (0,1), 2, rtol=1e-12),
  284. data(yv, 'bessel_yn_data_ipp-bessel_yn_data', (0,1j), 2, rtol=1e-12),
  285. data(yv, 'bessel_yv_data_ipp-bessel_yv_data', (0,1), 2, rtol=1e-10),
  286. data(yv, 'bessel_yv_data_ipp-bessel_yv_data', (0,1j), 2, rtol=1e-10),
  287. data(yvp, 'bessel_yv_prime_data_ipp-bessel_yv_prime_data', (0, 1), 2, rtol=4e-9),
  288. data(yvp, 'bessel_yv_prime_data_ipp-bessel_yv_prime_data', (0, 1j), 2, rtol=4e-9),
  289. data(zeta_, 'zeta_data_ipp-zeta_data', 0, 1, param_filter=(lambda s: s > 1)),
  290. data(zeta_, 'zeta_neg_data_ipp-zeta_neg_data', 0, 1, param_filter=(lambda s: s > 1)),
  291. data(zeta_, 'zeta_1_up_data_ipp-zeta_1_up_data', 0, 1, param_filter=(lambda s: s > 1)),
  292. data(zeta_, 'zeta_1_below_data_ipp-zeta_1_below_data', 0, 1, param_filter=(lambda s: s > 1)),
  293. data(gammaincinv, 'gamma_inv_small_data_ipp-gamma_inv_small_data', (0,1), 2, rtol=1e-11),
  294. data(gammaincinv, 'gamma_inv_data_ipp-gamma_inv_data', (0,1), 2, rtol=1e-14),
  295. data(gammaincinv, 'gamma_inv_big_data_ipp-gamma_inv_big_data', (0,1), 2, rtol=1e-11),
  296. data(gammainccinv, 'gamma_inv_small_data_ipp-gamma_inv_small_data', (0,1), 3, rtol=1e-12),
  297. data(gammainccinv, 'gamma_inv_data_ipp-gamma_inv_data', (0,1), 3, rtol=1e-14),
  298. data(gammainccinv, 'gamma_inv_big_data_ipp-gamma_inv_big_data', (0,1), 3, rtol=1e-14),
  299. data(gdtrix_, 'gamma_inv_small_data_ipp-gamma_inv_small_data', (0,1), 2, rtol=3e-13, knownfailure='gdtrix unflow some points'),
  300. data(gdtrix_, 'gamma_inv_data_ipp-gamma_inv_data', (0,1), 2, rtol=3e-15),
  301. data(gdtrix_, 'gamma_inv_big_data_ipp-gamma_inv_big_data', (0,1), 2),
  302. data(gdtrix_comp, 'gamma_inv_small_data_ipp-gamma_inv_small_data', (0,1), 2, knownfailure='gdtrix bad some points'),
  303. data(gdtrix_comp, 'gamma_inv_data_ipp-gamma_inv_data', (0,1), 3, rtol=6e-15),
  304. data(gdtrix_comp, 'gamma_inv_big_data_ipp-gamma_inv_big_data', (0,1), 3),
  305. data(chndtr, 'nccs_ipp-nccs', (2,0,1), 3, rtol=3e-5),
  306. data(chndtr, 'nccs_big_ipp-nccs_big', (2,0,1), 3, rtol=5e-4, knownfailure='chndtr inaccurate some points'),
  307. data(sph_harm_, 'spherical_harmonic_ipp-spherical_harmonic', (1,0,3,2), (4,5), rtol=5e-11,
  308. param_filter=(lambda p: np.ones(p.shape, '?'),
  309. lambda p: np.ones(p.shape, '?'),
  310. lambda p: np.logical_and(p < 2*np.pi, p >= 0),
  311. lambda p: np.logical_and(p < np.pi, p >= 0))),
  312. data(spherical_jn_, 'sph_bessel_data_ipp-sph_bessel_data', (0,1), 2, rtol=1e-13),
  313. data(spherical_yn_, 'sph_neumann_data_ipp-sph_neumann_data', (0,1), 2, rtol=8e-15),
  314. data(owens_t, 'owens_t_ipp-owens_t', (0, 1), 2, rtol=5e-14),
  315. data(owens_t, 'owens_t_large_data_ipp-owens_t_large_data', (0, 1), 2, rtol=8e-12),
  316. # -- test data exists in boost but is not used in scipy --
  317. # ibeta_derivative_data_ipp/ibeta_derivative_data.txt
  318. # ibeta_derivative_int_data_ipp/ibeta_derivative_int_data.txt
  319. # ibeta_derivative_large_data_ipp/ibeta_derivative_large_data.txt
  320. # ibeta_derivative_small_data_ipp/ibeta_derivative_small_data.txt
  321. # bessel_y01_prime_data_ipp/bessel_y01_prime_data.txt
  322. # bessel_yn_prime_data_ipp/bessel_yn_prime_data.txt
  323. # sph_bessel_prime_data_ipp/sph_bessel_prime_data.txt
  324. # sph_neumann_prime_data_ipp/sph_neumann_prime_data.txt
  325. # ellint_d2_data_ipp/ellint_d2_data.txt
  326. # ellint_d_data_ipp/ellint_d_data.txt
  327. # ellint_pi2_data_ipp/ellint_pi2_data.txt
  328. # ellint_pi3_data_ipp/ellint_pi3_data.txt
  329. # ellint_pi3_large_data_ipp/ellint_pi3_large_data.txt
  330. data(elliprc, 'ellint_rc_data_ipp-ellint_rc_data', (0, 1), 2,
  331. rtol=5e-16),
  332. data(elliprd, 'ellint_rd_data_ipp-ellint_rd_data', (0, 1, 2), 3,
  333. rtol=5e-16),
  334. data(elliprd, 'ellint_rd_0xy_ipp-ellint_rd_0xy', (0, 1, 2), 3,
  335. rtol=5e-16),
  336. data(elliprd, 'ellint_rd_0yy_ipp-ellint_rd_0yy', (0, 1, 2), 3,
  337. rtol=5e-16),
  338. data(elliprd, 'ellint_rd_xxx_ipp-ellint_rd_xxx', (0, 1, 2), 3,
  339. rtol=5e-16),
  340. # Some of the following rtol for elliprd may be larger than 5e-16 to
  341. # work around some hard cases in the Boost test where we get slightly
  342. # larger error than the ideal bound when the x (==y) input is close to
  343. # zero.
  344. # Also the accuracy on 32-bit buids with g++ may suffer from excess
  345. # loss of precision; see GCC bugzilla 323
  346. # https://gcc.gnu.org/bugzilla/show_bug.cgi?id=323
  347. data(elliprd, 'ellint_rd_xxz_ipp-ellint_rd_xxz', (0, 1, 2), 3,
  348. rtol=6.5e-16),
  349. data(elliprd, 'ellint_rd_xyy_ipp-ellint_rd_xyy', (0, 1, 2), 3,
  350. rtol=6e-16),
  351. data(elliprf, 'ellint_rf_data_ipp-ellint_rf_data', (0, 1, 2), 3,
  352. rtol=5e-16),
  353. data(elliprf, 'ellint_rf_xxx_ipp-ellint_rf_xxx', (0, 1, 2), 3,
  354. rtol=5e-16),
  355. data(elliprf, 'ellint_rf_xyy_ipp-ellint_rf_xyy', (0, 1, 2), 3,
  356. rtol=5e-16),
  357. data(elliprf, 'ellint_rf_xy0_ipp-ellint_rf_xy0', (0, 1, 2), 3,
  358. rtol=5e-16),
  359. data(elliprf, 'ellint_rf_0yy_ipp-ellint_rf_0yy', (0, 1, 2), 3,
  360. rtol=5e-16),
  361. # The accuracy of R_G is primarily limited by R_D that is used
  362. # internally. It is generally worse than R_D. Notice that we increased
  363. # the rtol for R_G here. The cases with duplicate arguments are
  364. # slightly less likely to be unbalanced (at least two arguments are
  365. # already balanced) so the error bound is slightly better. Again,
  366. # precision with g++ 32-bit is even worse.
  367. data(elliprg, 'ellint_rg_ipp-ellint_rg', (0, 1, 2), 3,
  368. rtol=8.0e-16),
  369. data(elliprg, 'ellint_rg_xxx_ipp-ellint_rg_xxx', (0, 1, 2), 3,
  370. rtol=6e-16),
  371. data(elliprg, 'ellint_rg_xyy_ipp-ellint_rg_xyy', (0, 1, 2), 3,
  372. rtol=7.5e-16),
  373. data(elliprg, 'ellint_rg_xy0_ipp-ellint_rg_xy0', (0, 1, 2), 3,
  374. rtol=5e-16),
  375. data(elliprg, 'ellint_rg_00x_ipp-ellint_rg_00x', (0, 1, 2), 3,
  376. rtol=5e-16),
  377. data(elliprj, 'ellint_rj_data_ipp-ellint_rj_data', (0, 1, 2, 3), 4,
  378. rtol=5e-16, atol=1e-25,
  379. param_filter=(lambda s: s <= 5e-26,)),
  380. # ellint_rc_data_ipp/ellint_rc_data.txt
  381. # ellint_rd_0xy_ipp/ellint_rd_0xy.txt
  382. # ellint_rd_0yy_ipp/ellint_rd_0yy.txt
  383. # ellint_rd_data_ipp/ellint_rd_data.txt
  384. # ellint_rd_xxx_ipp/ellint_rd_xxx.txt
  385. # ellint_rd_xxz_ipp/ellint_rd_xxz.txt
  386. # ellint_rd_xyy_ipp/ellint_rd_xyy.txt
  387. # ellint_rf_0yy_ipp/ellint_rf_0yy.txt
  388. # ellint_rf_data_ipp/ellint_rf_data.txt
  389. # ellint_rf_xxx_ipp/ellint_rf_xxx.txt
  390. # ellint_rf_xy0_ipp/ellint_rf_xy0.txt
  391. # ellint_rf_xyy_ipp/ellint_rf_xyy.txt
  392. # ellint_rg_00x_ipp/ellint_rg_00x.txt
  393. # ellint_rg_ipp/ellint_rg.txt
  394. # ellint_rg_xxx_ipp/ellint_rg_xxx.txt
  395. # ellint_rg_xy0_ipp/ellint_rg_xy0.txt
  396. # ellint_rg_xyy_ipp/ellint_rg_xyy.txt
  397. # ellint_rj_data_ipp/ellint_rj_data.txt
  398. # ellint_rj_e2_ipp/ellint_rj_e2.txt
  399. # ellint_rj_e3_ipp/ellint_rj_e3.txt
  400. # ellint_rj_e4_ipp/ellint_rj_e4.txt
  401. # ellint_rj_zp_ipp/ellint_rj_zp.txt
  402. # jacobi_elliptic_ipp/jacobi_elliptic.txt
  403. # jacobi_elliptic_small_ipp/jacobi_elliptic_small.txt
  404. # jacobi_large_phi_ipp/jacobi_large_phi.txt
  405. # jacobi_near_1_ipp/jacobi_near_1.txt
  406. # jacobi_zeta_big_phi_ipp/jacobi_zeta_big_phi.txt
  407. # jacobi_zeta_data_ipp/jacobi_zeta_data.txt
  408. # heuman_lambda_data_ipp/heuman_lambda_data.txt
  409. # hypergeometric_0F2_ipp/hypergeometric_0F2.txt
  410. # hypergeometric_1F1_big_ipp/hypergeometric_1F1_big.txt
  411. # hypergeometric_1F1_ipp/hypergeometric_1F1.txt
  412. # hypergeometric_1F1_small_random_ipp/hypergeometric_1F1_small_random.txt
  413. # hypergeometric_1F2_ipp/hypergeometric_1F2.txt
  414. # hypergeometric_1f1_large_regularized_ipp/hypergeometric_1f1_large_regularized.txt
  415. # hypergeometric_1f1_log_large_unsolved_ipp/hypergeometric_1f1_log_large_unsolved.txt
  416. # hypergeometric_2F0_half_ipp/hypergeometric_2F0_half.txt
  417. # hypergeometric_2F0_integer_a2_ipp/hypergeometric_2F0_integer_a2.txt
  418. # hypergeometric_2F0_ipp/hypergeometric_2F0.txt
  419. # hypergeometric_2F0_large_z_ipp/hypergeometric_2F0_large_z.txt
  420. # hypergeometric_2F1_ipp/hypergeometric_2F1.txt
  421. # hypergeometric_2F2_ipp/hypergeometric_2F2.txt
  422. # ncbeta_big_ipp/ncbeta_big.txt
  423. # nct_small_delta_ipp/nct_small_delta.txt
  424. # nct_asym_ipp/nct_asym.txt
  425. # ncbeta_ipp/ncbeta.txt
  426. # powm1_data_ipp/powm1_big_data.txt
  427. # powm1_sqrtp1m1_test_hpp/sqrtp1m1_data.txt
  428. # sinc_data_ipp/sinc_data.txt
  429. # test_gamma_data_ipp/gammap1m1_data.txt
  430. # tgamma_ratio_data_ipp/tgamma_ratio_data.txt
  431. # trig_data_ipp/trig_data.txt
  432. # trig_data2_ipp/trig_data2.txt
  433. ]
  434. @pytest.mark.parametrize('test', BOOST_TESTS, ids=repr)
  435. def test_boost(test):
  436. _test_factory(test)
  437. GSL_TESTS = [
  438. data_gsl(mathieu_a, 'mathieu_ab', (0, 1), 2, rtol=1e-13, atol=1e-13),
  439. data_gsl(mathieu_b, 'mathieu_ab', (0, 1), 3, rtol=1e-13, atol=1e-13),
  440. # Also the GSL output has limited accuracy...
  441. data_gsl(mathieu_ce_rad, 'mathieu_ce_se', (0, 1, 2), 3, rtol=1e-7, atol=1e-13),
  442. data_gsl(mathieu_se_rad, 'mathieu_ce_se', (0, 1, 2), 4, rtol=1e-7, atol=1e-13),
  443. data_gsl(mathieu_mc1_scaled, 'mathieu_mc_ms', (0, 1, 2), 3, rtol=1e-7, atol=1e-13),
  444. data_gsl(mathieu_ms1_scaled, 'mathieu_mc_ms', (0, 1, 2), 4, rtol=1e-7, atol=1e-13),
  445. data_gsl(mathieu_mc2_scaled, 'mathieu_mc_ms', (0, 1, 2), 5, rtol=1e-7, atol=1e-13),
  446. data_gsl(mathieu_ms2_scaled, 'mathieu_mc_ms', (0, 1, 2), 6, rtol=1e-7, atol=1e-13),
  447. ]
  448. @pytest.mark.parametrize('test', GSL_TESTS, ids=repr)
  449. def test_gsl(test):
  450. _test_factory(test)
  451. LOCAL_TESTS = [
  452. data_local(ellipkinc, 'ellipkinc_neg_m', (0, 1), 2),
  453. data_local(ellipkm1, 'ellipkm1', 0, 1),
  454. data_local(ellipeinc, 'ellipeinc_neg_m', (0, 1), 2),
  455. data_local(clog1p, 'log1p_expm1_complex', (0,1), (2,3), rtol=1e-14),
  456. data_local(cexpm1, 'log1p_expm1_complex', (0,1), (4,5), rtol=1e-14),
  457. data_local(gammainc, 'gammainc', (0, 1), 2, rtol=1e-12),
  458. data_local(gammaincc, 'gammaincc', (0, 1), 2, rtol=1e-11),
  459. data_local(ellip_harm_2, 'ellip',(0, 1, 2, 3, 4), 6, rtol=1e-10, atol=1e-13),
  460. data_local(ellip_harm, 'ellip',(0, 1, 2, 3, 4), 5, rtol=1e-10, atol=1e-13),
  461. data_local(wright_bessel, 'wright_bessel', (0, 1, 2), 3, rtol=1e-11),
  462. ]
  463. @pytest.mark.parametrize('test', LOCAL_TESTS, ids=repr)
  464. def test_local(test):
  465. _test_factory(test)
  466. def _test_factory(test, dtype=np.double):
  467. """Boost test"""
  468. with suppress_warnings() as sup:
  469. sup.filter(IntegrationWarning, "The occurrence of roundoff error is detected")
  470. with np.errstate(all='ignore'):
  471. test.check(dtype=dtype)