Merged revisions 76755 via svnmerge from

svn+ssh://pythondev@svn.python.org/python/trunk

........
  r76755 | mark.dickinson | 2009-12-11 17:29:33 +0000 (Fri, 11 Dec 2009) | 2 lines

  Issue #3366:  Add lgamma function to math module.
........

5 files changed, 216 insertions, 6 deletions

diff --git a/Doc/library/math.rst b/Doc/library/math.rst
index fdd3a1f..e903c5f 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 764c0e0..2b1fc69 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 2226023..2052430 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 8f051b7..850369e 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 052fca1..ece68a7 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -322,6 +322,60 @@ m_tgamma(double x)
 }
 
 /*
+   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
    is necessary to get consistent behaviour across platforms.
@@ -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},