test_precompute_gammainc.py 4.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109
  1. import numpy as np # np is actually used, in the decorators below.
  2. import pytest
  3. from scipy.special._testutils import MissingModule, check_version
  4. from scipy.special._mptestutils import (
  5. Arg, IntArg, mp_assert_allclose, assert_mpmath_equal)
  6. from scipy.special._precompute.gammainc_asy import (
  7. compute_g, compute_alpha, compute_d)
  8. from scipy.special._precompute.gammainc_data import gammainc, gammaincc
  9. try:
  10. import sympy
  11. except ImportError:
  12. sympy = MissingModule('sympy')
  13. try:
  14. import mpmath as mp
  15. except ImportError:
  16. mp = MissingModule('mpmath')
  17. @check_version(mp, '0.19')
  18. def test_g():
  19. # Test data for the g_k. See DLMF 5.11.4.
  20. with mp.workdps(30):
  21. g = [mp.mpf(1), mp.mpf(1)/12, mp.mpf(1)/288,
  22. -mp.mpf(139)/51840, -mp.mpf(571)/2488320,
  23. mp.mpf(163879)/209018880, mp.mpf(5246819)/75246796800]
  24. mp_assert_allclose(compute_g(7), g)
  25. @pytest.mark.slow
  26. @check_version(mp, '0.19')
  27. @check_version(sympy, '0.7')
  28. @pytest.mark.xfail_on_32bit("rtol only 2e-11, see gh-6938")
  29. def test_alpha():
  30. # Test data for the alpha_k. See DLMF 8.12.14.
  31. with mp.workdps(30):
  32. alpha = [mp.mpf(0), mp.mpf(1), mp.mpf(1)/3, mp.mpf(1)/36,
  33. -mp.mpf(1)/270, mp.mpf(1)/4320, mp.mpf(1)/17010,
  34. -mp.mpf(139)/5443200, mp.mpf(1)/204120]
  35. mp_assert_allclose(compute_alpha(9), alpha)
  36. @pytest.mark.xslow
  37. @check_version(mp, '0.19')
  38. @check_version(sympy, '0.7')
  39. def test_d():
  40. # Compare the d_{k, n} to the results in appendix F of [1].
  41. #
  42. # Sources
  43. # -------
  44. # [1] DiDonato and Morris, Computation of the Incomplete Gamma
  45. # Function Ratios and their Inverse, ACM Transactions on
  46. # Mathematical Software, 1986.
  47. with mp.workdps(50):
  48. dataset = [(0, 0, -mp.mpf('0.333333333333333333333333333333')),
  49. (0, 12, mp.mpf('0.102618097842403080425739573227e-7')),
  50. (1, 0, -mp.mpf('0.185185185185185185185185185185e-2')),
  51. (1, 12, mp.mpf('0.119516285997781473243076536700e-7')),
  52. (2, 0, mp.mpf('0.413359788359788359788359788360e-2')),
  53. (2, 12, -mp.mpf('0.140925299108675210532930244154e-7')),
  54. (3, 0, mp.mpf('0.649434156378600823045267489712e-3')),
  55. (3, 12, -mp.mpf('0.191111684859736540606728140873e-7')),
  56. (4, 0, -mp.mpf('0.861888290916711698604702719929e-3')),
  57. (4, 12, mp.mpf('0.288658297427087836297341274604e-7')),
  58. (5, 0, -mp.mpf('0.336798553366358150308767592718e-3')),
  59. (5, 12, mp.mpf('0.482409670378941807563762631739e-7')),
  60. (6, 0, mp.mpf('0.531307936463992223165748542978e-3')),
  61. (6, 12, -mp.mpf('0.882860074633048352505085243179e-7')),
  62. (7, 0, mp.mpf('0.344367606892377671254279625109e-3')),
  63. (7, 12, -mp.mpf('0.175629733590604619378669693914e-6')),
  64. (8, 0, -mp.mpf('0.652623918595309418922034919727e-3')),
  65. (8, 12, mp.mpf('0.377358774161109793380344937299e-6')),
  66. (9, 0, -mp.mpf('0.596761290192746250124390067179e-3')),
  67. (9, 12, mp.mpf('0.870823417786464116761231237189e-6'))]
  68. d = compute_d(10, 13)
  69. res = [d[k][n] for k, n, std in dataset]
  70. std = [x[2] for x in dataset]
  71. mp_assert_allclose(res, std)
  72. @check_version(mp, '0.19')
  73. def test_gammainc():
  74. # Quick check that the gammainc in
  75. # special._precompute.gammainc_data agrees with mpmath's
  76. # gammainc.
  77. assert_mpmath_equal(gammainc,
  78. lambda a, x: mp.gammainc(a, b=x, regularized=True),
  79. [Arg(0, 100, inclusive_a=False), Arg(0, 100)],
  80. nan_ok=False, rtol=1e-17, n=50, dps=50)
  81. @pytest.mark.xslow
  82. @check_version(mp, '0.19')
  83. def test_gammaincc():
  84. # Check that the gammaincc in special._precompute.gammainc_data
  85. # agrees with mpmath's gammainc.
  86. assert_mpmath_equal(lambda a, x: gammaincc(a, x, dps=1000),
  87. lambda a, x: mp.gammainc(a, a=x, regularized=True),
  88. [Arg(20, 100), Arg(20, 100)],
  89. nan_ok=False, rtol=1e-17, n=50, dps=1000)
  90. # Test the fast integer path
  91. assert_mpmath_equal(gammaincc,
  92. lambda a, x: mp.gammainc(a, a=x, regularized=True),
  93. [IntArg(1, 100), Arg(0, 100)],
  94. nan_ok=False, rtol=1e-17, n=50, dps=50)