diff --git a/Doc/library/math.rst b/Doc/library/math.rst index fdd3a1fbd7b..e903c5f9022 100644 --- a/Doc/library/math.rst +++ b/Doc/library/math.rst @@ -288,6 +288,14 @@ Special functions .. versionadded:: 3.2 +.. function:: lgamma(x) + + Return the natural logarithm of the absolute value of the Gamma + function at *x*. + + .. versionadded:: 2.7 + + Constants --------- diff --git a/Lib/test/math_testcases.txt b/Lib/test/math_testcases.txt index 764c0e02cbe..2b1fc694700 100644 --- a/Lib/test/math_testcases.txt +++ b/Lib/test/math_testcases.txt @@ -47,6 +47,111 @@ -- MPFR homepage at http://www.mpfr.org for more information about the -- MPFR project. +--------------------------------------------------------- +-- lgamma: log of absolute value of the gamma function -- +--------------------------------------------------------- + +-- special values +lgam0000 lgamma 0.0 -> inf divide-by-zero +lgam0001 lgamma -0.0 -> inf divide-by-zero +lgam0002 lgamma inf -> inf +lgam0003 lgamma -inf -> inf +lgam0004 lgamma nan -> nan + +-- negative integers +lgam0010 lgamma -1 -> inf divide-by-zero +lgam0011 lgamma -2 -> inf divide-by-zero +lgam0012 lgamma -1e16 -> inf divide-by-zero +lgam0013 lgamma -1e300 -> inf divide-by-zero +lgam0014 lgamma -1.79e308 -> inf divide-by-zero + +-- small positive integers give factorials +lgam0020 lgamma 1 -> 0.0 +lgam0021 lgamma 2 -> 0.0 +lgam0022 lgamma 3 -> 0.69314718055994529 +lgam0023 lgamma 4 -> 1.791759469228055 +lgam0024 lgamma 5 -> 3.1780538303479458 +lgam0025 lgamma 6 -> 4.7874917427820458 + +-- half integers +lgam0030 lgamma 0.5 -> 0.57236494292470008 +lgam0031 lgamma 1.5 -> -0.12078223763524522 +lgam0032 lgamma 2.5 -> 0.28468287047291918 +lgam0033 lgamma 3.5 -> 1.2009736023470743 +lgam0034 lgamma -0.5 -> 1.2655121234846454 +lgam0035 lgamma -1.5 -> 0.86004701537648098 +lgam0036 lgamma -2.5 -> -0.056243716497674054 +lgam0037 lgamma -3.5 -> -1.309006684993042 + +-- values near 0 +lgam0040 lgamma 0.1 -> 2.252712651734206 +lgam0041 lgamma 0.01 -> 4.5994798780420219 +lgam0042 lgamma 1e-8 -> 18.420680738180209 +lgam0043 lgamma 1e-16 -> 36.841361487904734 +lgam0044 lgamma 1e-30 -> 69.077552789821368 +lgam0045 lgamma 1e-160 -> 368.41361487904732 +lgam0046 lgamma 1e-308 -> 709.19620864216608 +lgam0047 lgamma 5.6e-309 -> 709.77602713741896 +lgam0048 lgamma 5.5e-309 -> 709.79404564292167 +lgam0049 lgamma 1e-309 -> 711.49879373516012 +lgam0050 lgamma 1e-323 -> 743.74692474082133 +lgam0051 lgamma 5e-324 -> 744.44007192138122 +lgam0060 lgamma -0.1 -> 2.3689613327287886 +lgam0061 lgamma -0.01 -> 4.6110249927528013 +lgam0062 lgamma -1e-8 -> 18.420680749724522 +lgam0063 lgamma -1e-16 -> 36.841361487904734 +lgam0064 lgamma -1e-30 -> 69.077552789821368 +lgam0065 lgamma -1e-160 -> 368.41361487904732 +lgam0066 lgamma -1e-308 -> 709.19620864216608 +lgam0067 lgamma -5.6e-309 -> 709.77602713741896 +lgam0068 lgamma -5.5e-309 -> 709.79404564292167 +lgam0069 lgamma -1e-309 -> 711.49879373516012 +lgam0070 lgamma -1e-323 -> 743.74692474082133 +lgam0071 lgamma -5e-324 -> 744.44007192138122 + +-- values near negative integers +lgam0080 lgamma -0.99999999999999989 -> 36.736800569677101 +lgam0081 lgamma -1.0000000000000002 -> 36.043653389117154 +lgam0082 lgamma -1.9999999999999998 -> 35.350506208557213 +lgam0083 lgamma -2.0000000000000004 -> 34.657359027997266 +lgam0084 lgamma -100.00000000000001 -> -331.85460524980607 +lgam0085 lgamma -99.999999999999986 -> -331.85460524980596 + +-- large inputs +lgam0100 lgamma 170 -> 701.43726380873704 +lgam0101 lgamma 171 -> 706.57306224578736 +lgam0102 lgamma 171.624 -> 709.78077443669895 +lgam0103 lgamma 171.625 -> 709.78591682948365 +lgam0104 lgamma 172 -> 711.71472580228999 +lgam0105 lgamma 2000 -> 13198.923448054265 +lgam0106 lgamma 2.55998332785163e305 -> 1.7976931348623099e+308 +lgam0107 lgamma 2.55998332785164e305 -> inf overflow +lgam0108 lgamma 1.7e308 -> inf overflow + +-- inputs for which gamma(x) is tiny +lgam0120 lgamma -100.5 -> -364.90096830942736 +lgam0121 lgamma -160.5 -> -656.88005261126432 +lgam0122 lgamma -170.5 -> -707.99843314507882 +lgam0123 lgamma -171.5 -> -713.14301641168481 +lgam0124 lgamma -176.5 -> -738.95247590846486 +lgam0125 lgamma -177.5 -> -744.13144651738037 +lgam0126 lgamma -178.5 -> -749.3160351186001 + +lgam0130 lgamma -1000.5 -> -5914.4377011168517 +lgam0131 lgamma -30000.5 -> -279278.6629959144 +lgam0132 lgamma -4503599627370495.5 -> -1.5782258434492883e+17 + +-- results close to 0: positive argument ... +lgam0150 lgamma 0.99999999999999989 -> 6.4083812134800075e-17 +lgam0151 lgamma 1.0000000000000002 -> -1.2816762426960008e-16 +lgam0152 lgamma 1.9999999999999998 -> -9.3876980655431170e-17 +lgam0153 lgamma 2.0000000000000004 -> 1.8775396131086244e-16 + +-- ... and negative argument +lgam0160 lgamma -2.7476826467 -> -5.2477408147689136e-11 +lgam0161 lgamma -2.457024738 -> 3.3464637541912932e-10 + + --------------------------- -- gamma: Gamma function -- --------------------------- diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index 22260231551..20524300212 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -48,6 +48,36 @@ def to_ulps(x): n = ~(n+2**63) return n +def ulps_check(expected, got, ulps=20): + """Given non-NaN floats `expected` and `got`, + check that they're equal to within the given number of ulps. + + Returns None on success and an error message on failure.""" + + ulps_error = to_ulps(got) - to_ulps(expected) + if abs(ulps_error) <= ulps: + return None + return "error = {} ulps; permitted error = {} ulps".format(ulps_error, + ulps) + +def acc_check(expected, got, rel_err=2e-15, abs_err = 5e-323): + """Determine whether non-NaN floats a and b are equal to within a + (small) rounding error. The default values for rel_err and + abs_err are chosen to be suitable for platforms where a float is + represented by an IEEE 754 double. They allow an error of between + 9 and 19 ulps.""" + + # need to special case infinities, since inf - inf gives nan + if math.isinf(expected) and got == expected: + return None + + error = got - expected + + permitted_error = max(abs_err, rel_err * abs(expected)) + if abs(error) < permitted_error: + return None + return "error = {}; permitted error = {}".format(error, + permitted_error) def parse_mtestfile(fname): """Parse a file with test values @@ -949,13 +979,23 @@ class MathTests(unittest.TestCase): except OverflowError: got = 'OverflowError' - diff_ulps = None + accuracy_failure = None if isinstance(got, float) and isinstance(expected, float): if math.isnan(expected) and math.isnan(got): continue if not math.isnan(expected) and not math.isnan(got): - diff_ulps = to_ulps(expected) - to_ulps(got) - if abs(diff_ulps) <= ALLOWED_ERROR: + # we use different closeness criteria for + # different functions. + if fn == 'gamma': + accuracy_failure = ulps_check(expected, got, 20) + elif fn == 'lgamma': + accuracy_failure = acc_check(expected, got, + rel_err = 5e-15, + abs_err = 5e-15) + else: + raise ValueError("don't know how to check accuracy " + "for this function") + if accuracy_failure is None: continue if isinstance(got, str) and isinstance(expected, str): @@ -963,8 +1003,8 @@ class MathTests(unittest.TestCase): continue fail_msg = fail_fmt.format(id, fn, arg, expected, got) - if diff_ulps is not None: - fail_msg += ' ({} ulps)'.format(diff_ulps) + if accuracy_failure is not None: + fail_msg += ' ({})'.format(accuracy_failure) failures.append(fail_msg) if failures: diff --git a/Misc/NEWS b/Misc/NEWS index 8f051b74ea3..850369e1579 100644 --- a/Misc/NEWS +++ b/Misc/NEWS @@ -443,7 +443,7 @@ Extension Modules - Issue #7078: Set struct.__doc__ from _struct.__doc__. -- Issue #3366: Add gamma function to math module. +- Issue #3366: Add gamma, lgamma functions to math module. - Issue #6877: It is now possible to link the readline extension to the libedit readline emulation on OSX 10.5 or later. diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 052fca1be31..ece68a742c8 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -321,6 +321,60 @@ m_tgamma(double x) return r; } +/* + lgamma: natural log of the absolute value of the Gamma function. + For large arguments, Lanczos' formula works extremely well here. +*/ + +static double +m_lgamma(double x) +{ + double r, absx; + + /* special cases */ + if (!Py_IS_FINITE(x)) { + if (Py_IS_NAN(x)) + return x; /* lgamma(nan) = nan */ + else + return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */ + } + + /* integer arguments */ + if (x == floor(x) && x <= 2.0) { + if (x <= 0.0) { + errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */ + return Py_HUGE_VAL; /* integers n <= 0 */ + } + else { + return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */ + } + } + + absx = fabs(x); + /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */ + if (absx < 1e-20) + return -log(absx); + + /* Lanczos' formula */ + if (x > 0.0) { + /* we could save a fraction of a ulp in accuracy by having a + second set of numerator coefficients for lanczos_sum that + absorbed the exp(-lanczos_g) term, and throwing out the + lanczos_g subtraction below; it's probably not worth it. */ + r = log(lanczos_sum(x)) - lanczos_g + + (x-0.5)*(log(x+lanczos_g-0.5)-1); + } + else { + r = log(pi) - log(fabs(sinpi(absx))) - log(absx) - + (log(lanczos_sum(absx)) - lanczos_g + + (absx-0.5)*(log(absx+lanczos_g-0.5)-1)); + } + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; +} + + /* wrapper for atan2 that deals directly with special cases before delegating to the platform libm for the remaining cases. This @@ -694,6 +748,8 @@ PyDoc_STRVAR(math_floor_doc, FUNC1A(gamma, m_tgamma, "gamma(x)\n\nGamma function at x.") +FUNC1A(lgamma, m_lgamma, + "lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.") FUNC1(log1p, log1p, 1, "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n" "The result is computed in a way which is accurate for x near zero.") @@ -1448,6 +1504,7 @@ static PyMethodDef math_methods[] = { {"isinf", math_isinf, METH_O, math_isinf_doc}, {"isnan", math_isnan, METH_O, math_isnan_doc}, {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc}, + {"lgamma", math_lgamma, METH_O, math_lgamma_doc}, {"log", math_log, METH_VARARGS, math_log_doc}, {"log1p", math_log1p, METH_O, math_log1p_doc}, {"log10", math_log10, METH_O, math_log10_doc},