test_loggamma.py 2.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172
  1. from __future__ import division, print_function, absolute_import
  2. import numpy as np
  3. from numpy.testing import assert_allclose, assert_
  4. from scipy.special._testutils import FuncData
  5. from scipy.special import gamma, gammaln, loggamma
  6. def test_identities1():
  7. # test the identity exp(loggamma(z)) = gamma(z)
  8. x = np.array([-99.5, -9.5, -0.5, 0.5, 9.5, 99.5])
  9. y = x.copy()
  10. x, y = np.meshgrid(x, y)
  11. z = (x + 1J*y).flatten()
  12. dataset = np.vstack((z, gamma(z))).T
  13. def f(z):
  14. return np.exp(loggamma(z))
  15. FuncData(f, dataset, 0, 1, rtol=1e-14, atol=1e-14).check()
  16. def test_identities2():
  17. # test the identity loggamma(z + 1) = log(z) + loggamma(z)
  18. x = np.array([-99.5, -9.5, -0.5, 0.5, 9.5, 99.5])
  19. y = x.copy()
  20. x, y = np.meshgrid(x, y)
  21. z = (x + 1J*y).flatten()
  22. dataset = np.vstack((z, np.log(z) + loggamma(z))).T
  23. def f(z):
  24. return loggamma(z + 1)
  25. FuncData(f, dataset, 0, 1, rtol=1e-14, atol=1e-14).check()
  26. def test_complex_dispatch_realpart():
  27. # Test that the real parts of loggamma and gammaln agree on the
  28. # real axis.
  29. x = np.r_[-np.logspace(10, -10), np.logspace(-10, 10)] + 0.5
  30. dataset = np.vstack((x, gammaln(x))).T
  31. def f(z):
  32. z = np.array(z, dtype='complex128')
  33. return loggamma(z).real
  34. FuncData(f, dataset, 0, 1, rtol=1e-14, atol=1e-14).check()
  35. def test_real_dispatch():
  36. x = np.logspace(-10, 10) + 0.5
  37. dataset = np.vstack((x, gammaln(x))).T
  38. FuncData(loggamma, dataset, 0, 1, rtol=1e-14, atol=1e-14).check()
  39. assert_(loggamma(0) == np.inf)
  40. assert_(np.isnan(loggamma(-1)))
  41. def test_gh_6536():
  42. z = loggamma(complex(-3.4, +0.0))
  43. zbar = loggamma(complex(-3.4, -0.0))
  44. assert_allclose(z, zbar.conjugate(), rtol=1e-15, atol=0)
  45. def test_branch_cut():
  46. # Make sure negative zero is treated correctly
  47. x = -np.logspace(300, -30, 100)
  48. z = np.asarray([complex(x0, 0.0) for x0 in x])
  49. zbar = np.asarray([complex(x0, -0.0) for x0 in x])
  50. assert_allclose(z, zbar.conjugate(), rtol=1e-15, atol=0)