test_loggamma.py 1.9 KB

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