Issue #3366: Add lgamma function to math module.

This commit is contained in:
Mark Dickinson 2009-12-11 17:29:33 +00:00
parent 5cc4e2a040
commit 9be87bc992
5 changed files with 216 additions and 6 deletions

View File

@ -318,6 +318,14 @@ Special functions
.. versionadded:: 2.7 .. versionadded:: 2.7
.. function:: lgamma(x)
Return the natural logarithm of the absolute value of the Gamma
function at *x*.
.. versionadded:: 2.7
Constants Constants
--------- ---------

View File

@ -47,6 +47,111 @@
-- MPFR homepage at http://www.mpfr.org for more information about the -- MPFR homepage at http://www.mpfr.org for more information about the
-- MPFR project. -- 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 -- -- gamma: Gamma function --
--------------------------- ---------------------------

View File

@ -48,6 +48,36 @@ def to_ulps(x):
n = ~(n+2**63) n = ~(n+2**63)
return n 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): def parse_mtestfile(fname):
"""Parse a file with test values """Parse a file with test values
@ -952,13 +982,23 @@ class MathTests(unittest.TestCase):
except OverflowError: except OverflowError:
got = 'OverflowError' got = 'OverflowError'
diff_ulps = None accuracy_failure = None
if isinstance(got, float) and isinstance(expected, float): if isinstance(got, float) and isinstance(expected, float):
if math.isnan(expected) and math.isnan(got): if math.isnan(expected) and math.isnan(got):
continue continue
if not math.isnan(expected) and not math.isnan(got): if not math.isnan(expected) and not math.isnan(got):
diff_ulps = to_ulps(expected) - to_ulps(got) # we use different closeness criteria for
if abs(diff_ulps) <= ALLOWED_ERROR: # 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 continue
if isinstance(got, str) and isinstance(expected, str): if isinstance(got, str) and isinstance(expected, str):
@ -966,8 +1006,8 @@ class MathTests(unittest.TestCase):
continue continue
fail_msg = fail_fmt.format(id, fn, arg, expected, got) fail_msg = fail_fmt.format(id, fn, arg, expected, got)
if diff_ulps is not None: if accuracy_failure is not None:
fail_msg += ' ({} ulps)'.format(diff_ulps) fail_msg += ' ({})'.format(accuracy_failure)
failures.append(fail_msg) failures.append(fail_msg)
if failures: if failures:

View File

@ -1654,7 +1654,7 @@ Extension Modules
- Issue #7078: Set struct.__doc__ from _struct.__doc__. - 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 #6823: Allow time.strftime() to accept a tuple with a isdst field - Issue #6823: Allow time.strftime() to accept a tuple with a isdst field
outside of the range of [-1, 1] by normalizing the value to within that outside of the range of [-1, 1] by normalizing the value to within that

View File

@ -321,6 +321,60 @@ m_tgamma(double x)
return r; 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 wrapper for atan2 that deals directly with special cases before
delegating to the platform libm for the remaining cases. This delegating to the platform libm for the remaining cases. This
@ -639,6 +693,8 @@ FUNC1(floor, floor, 0,
"This is the largest integral value <= x.") "This is the largest integral value <= x.")
FUNC1A(gamma, m_tgamma, FUNC1A(gamma, m_tgamma,
"gamma(x)\n\nGamma function at x.") "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, FUNC1(log1p, log1p, 1,
"log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n" "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.") "The result is computed in a way which is accurate for x near zero.")
@ -1375,6 +1431,7 @@ static PyMethodDef math_methods[] = {
{"isinf", math_isinf, METH_O, math_isinf_doc}, {"isinf", math_isinf, METH_O, math_isinf_doc},
{"isnan", math_isnan, METH_O, math_isnan_doc}, {"isnan", math_isnan, METH_O, math_isnan_doc},
{"ldexp", math_ldexp, METH_VARARGS, math_ldexp_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}, {"log", math_log, METH_VARARGS, math_log_doc},
{"log1p", math_log1p, METH_O, math_log1p_doc}, {"log1p", math_log1p, METH_O, math_log1p_doc},
{"log10", math_log10, METH_O, math_log10_doc}, {"log10", math_log10, METH_O, math_log10_doc},