]> granicus.if.org Git - python/commitdiff
Merged revisions 75117 via svnmerge from
authorMark Dickinson <dickinsm@gmail.com>
Mon, 28 Sep 2009 19:21:11 +0000 (19:21 +0000)
committerMark Dickinson <dickinsm@gmail.com>
Mon, 28 Sep 2009 19:21:11 +0000 (19:21 +0000)
svn+ssh://pythondev@svn.python.org/python/trunk

........
  r75117 | mark.dickinson | 2009-09-28 19:54:55 +0100 (Mon, 28 Sep 2009) | 3 lines

  Issue #3366:  Add gamma function to math module.
  (lgamma, erf and erfc to follow).
........

Doc/library/math.rst
Lib/test/math_testcases.txt [new file with mode: 0644]
Lib/test/test_math.py
Misc/NEWS
Modules/mathmodule.c

index 566603d4c3f27aad2bec05d4d24c05255c7848b1..270c25f16392ccc02d9acd7b8e70996dce60a857 100644 (file)
@@ -278,6 +278,16 @@ Hyperbolic functions
    Return the hyperbolic tangent of *x*.
 
 
+Special functions
+-----------------
+
+.. function:: gamma(x)
+
+   Return the Gamma function at *x*.
+
+   .. versionadded:: 2.7
+
+
 Constants
 ---------
 
diff --git a/Lib/test/math_testcases.txt b/Lib/test/math_testcases.txt
new file mode 100644 (file)
index 0000000..764c0e0
--- /dev/null
@@ -0,0 +1,146 @@
+-- Testcases for functions in math.
+--
+-- Each line takes the form:
+--
+-- <testid> <function> <input_value> -> <output_value> <flags>
+--
+-- where:
+--
+--   <testid> is a short name identifying the test,
+--
+--   <function> is the function to be tested (exp, cos, asinh, ...),
+--
+--   <input_value> is a string representing a floating-point value
+--
+--   <output_value> is the expected (ideal) output value, again
+--     represented as a string.
+--
+--   <flags> is a list of the floating-point flags required by C99
+--
+-- The possible flags are:
+--
+--   divide-by-zero : raised when a finite input gives a
+--     mathematically infinite result.
+--
+--   overflow : raised when a finite input gives a finite result that
+--     is too large to fit in the usual range of an IEEE 754 double.
+--
+--   invalid : raised for invalid inputs (e.g., sqrt(-1))
+--
+--   ignore-sign : indicates that the sign of the result is
+--     unspecified; e.g., if the result is given as inf,
+--     then both -inf and inf should be accepted as correct.
+--
+-- Flags may appear in any order.
+--
+-- Lines beginning with '--' (like this one) start a comment, and are
+-- ignored.  Blank lines, or lines containing only whitespace, are also
+-- ignored.
+
+-- Many of the values below were computed with the help of
+-- version 2.4 of the MPFR library for multiple-precision
+-- floating-point computations with correct rounding.  All output
+-- values in this file are (modulo yet-to-be-discovered bugs)
+-- correctly rounded, provided that each input and output decimal
+-- floating-point value below is interpreted as a representation of
+-- the corresponding nearest IEEE 754 double-precision value.  See the
+-- MPFR homepage at http://www.mpfr.org for more information about the
+-- MPFR project.
+
+---------------------------
+-- gamma: Gamma function --
+---------------------------
+
+-- special values
+gam0000 gamma 0.0 -> inf        divide-by-zero
+gam0001 gamma -0.0 -> -inf      divide-by-zero
+gam0002 gamma inf -> inf
+gam0003 gamma -inf -> nan       invalid
+gam0004 gamma nan -> nan
+
+-- negative integers inputs are invalid
+gam0010 gamma -1 -> nan         invalid
+gam0011 gamma -2 -> nan         invalid
+gam0012 gamma -1e16 -> nan      invalid
+gam0013 gamma -1e300 -> nan     invalid
+
+-- small positive integers give factorials
+gam0020 gamma 1 -> 1
+gam0021 gamma 2 -> 1
+gam0022 gamma 3 -> 2
+gam0023 gamma 4 -> 6
+gam0024 gamma 5 -> 24
+gam0025 gamma 6 -> 120
+
+-- half integers
+gam0030 gamma 0.5 -> 1.7724538509055161
+gam0031 gamma 1.5 -> 0.88622692545275805
+gam0032 gamma 2.5 -> 1.3293403881791370
+gam0033 gamma 3.5 -> 3.3233509704478426
+gam0034 gamma -0.5 -> -3.5449077018110322
+gam0035 gamma -1.5 -> 2.3632718012073548
+gam0036 gamma -2.5 -> -0.94530872048294190
+gam0037 gamma -3.5 -> 0.27008820585226911
+
+-- values near 0
+gam0040 gamma 0.1 -> 9.5135076986687306
+gam0041 gamma 0.01 -> 99.432585119150602
+gam0042 gamma 1e-8 -> 99999999.422784343
+gam0043 gamma 1e-16 -> 10000000000000000
+gam0044 gamma 1e-30 -> 9.9999999999999988e+29
+gam0045 gamma 1e-160 -> 1.0000000000000000e+160
+gam0046 gamma 1e-308 -> 1.0000000000000000e+308
+gam0047 gamma 5.6e-309 -> 1.7857142857142848e+308
+gam0048 gamma 5.5e-309 -> inf   overflow
+gam0049 gamma 1e-309 -> inf     overflow
+gam0050 gamma 1e-323 -> inf     overflow
+gam0051 gamma 5e-324 -> inf     overflow
+gam0060 gamma -0.1 -> -10.686287021193193
+gam0061 gamma -0.01 -> -100.58719796441078
+gam0062 gamma -1e-8 -> -100000000.57721567
+gam0063 gamma -1e-16 -> -10000000000000000
+gam0064 gamma -1e-30 -> -9.9999999999999988e+29
+gam0065 gamma -1e-160 -> -1.0000000000000000e+160
+gam0066 gamma -1e-308 -> -1.0000000000000000e+308
+gam0067 gamma -5.6e-309 -> -1.7857142857142848e+308
+gam0068 gamma -5.5e-309 -> -inf overflow
+gam0069 gamma -1e-309 -> -inf   overflow
+gam0070 gamma -1e-323 -> -inf   overflow
+gam0071 gamma -5e-324 -> -inf   overflow
+
+-- values near negative integers
+gam0080 gamma -0.99999999999999989 -> -9007199254740992.0
+gam0081 gamma -1.0000000000000002 -> 4503599627370495.5
+gam0082 gamma -1.9999999999999998 -> 2251799813685248.5
+gam0083 gamma -2.0000000000000004 -> -1125899906842623.5
+gam0084 gamma -100.00000000000001 -> -7.5400833348831090e-145
+gam0085 gamma -99.999999999999986 -> 7.5400833348840962e-145
+
+-- large inputs
+gam0100 gamma 170 -> 4.2690680090047051e+304
+gam0101 gamma 171 -> 7.2574156153079990e+306
+gam0102 gamma 171.624 -> 1.7942117599248104e+308
+gam0103 gamma 171.625 -> inf    overflow
+gam0104 gamma 172 -> inf        overflow
+gam0105 gamma 2000 -> inf       overflow
+gam0106 gamma 1.7e308 -> inf    overflow
+
+-- inputs for which gamma(x) is tiny
+gam0120 gamma -100.5 -> -3.3536908198076787e-159
+gam0121 gamma -160.5 -> -5.2555464470078293e-286
+gam0122 gamma -170.5 -> -3.3127395215386074e-308
+gam0123 gamma -171.5 -> 1.9316265431711902e-310
+gam0124 gamma -176.5 -> -1.1956388629358166e-321
+gam0125 gamma -177.5 -> 4.9406564584124654e-324
+gam0126 gamma -178.5 -> -0.0
+gam0127 gamma -179.5 -> 0.0
+gam0128 gamma -201.0001 -> 0.0
+gam0129 gamma -202.9999 -> -0.0
+gam0130 gamma -1000.5 -> -0.0
+gam0131 gamma -1000000000.3 -> -0.0
+gam0132 gamma -4503599627370495.5 -> 0.0
+
+-- inputs that cause problems for the standard reflection formula,
+-- thanks to loss of accuracy in 1-x
+gam0140 gamma -63.349078729022985 -> 4.1777971677761880e-88
+gam0141 gamma -127.45117632943295 -> 1.1831110896236810e-214
index 6f1578285962a82910e98d1bc1403a03d61cbb10..f29bddd8dec5eda484fda35a8b55113f1009375d 100644 (file)
@@ -7,6 +7,7 @@ import math
 import os
 import sys
 import random
+import struct
 
 eps = 1E-05
 NAN = float('nan')
@@ -29,8 +30,50 @@ if __name__ == '__main__':
 else:
     file = __file__
 test_dir = os.path.dirname(file) or os.curdir
+math_testcases = os.path.join(test_dir, 'math_testcases.txt')
 test_file = os.path.join(test_dir, 'cmath_testcases.txt')
 
+def to_ulps(x):
+    """Convert a non-NaN float x to an integer, in such a way that
+    adjacent floats are converted to adjacent integers.  Then
+    abs(ulps(x) - ulps(y)) gives the difference in ulps between two
+    floats.
+
+    The results from this function will only make sense on platforms
+    where C doubles are represented in IEEE 754 binary64 format.
+
+    """
+    n = struct.unpack('q', struct.pack('<d', x))[0]
+    if n < 0:
+        n = ~(n+2**63)
+    return n
+
+
+def parse_mtestfile(fname):
+    """Parse a file with test values
+
+    -- starts a comment
+    blank lines, or lines containing only a comment, are ignored
+    other lines are expected to have the form
+      id fn arg -> expected [flag]*
+
+    """
+    with open(fname) as fp:
+        for line in fp:
+            # strip comments, and skip blank lines
+            if '--' in line:
+                line = line[:line.index('--')]
+            if not line.strip():
+                continue
+
+            lhs, rhs = line.split('->')
+            id, fn, arg = lhs.split()
+            rhs_pieces = rhs.split()
+            exp = rhs_pieces[0]
+            flags = rhs_pieces[1:]
+
+            yield (id, fn, float(arg), float(exp), flags)
+
 def parse_testfile(fname):
     """Parse a file with test values
 
@@ -884,6 +927,51 @@ class MathTests(unittest.TestCase):
                 self.fail(message)
             self.ftest("%s:%s(%r)" % (id, fn, ar), result, er)
 
+    @unittest.skipUnless(float.__getformat__("double").startswith("IEEE"),
+                         "test requires IEEE 754 doubles")
+    def test_mtestfile(self):
+        ALLOWED_ERROR = 20  # permitted error, in ulps
+        fail_fmt = "{}:{}({!r}): expected {!r}, got {!r}"
+
+        failures = []
+        for id, fn, arg, expected, flags in parse_mtestfile(math_testcases):
+            func = getattr(math, fn)
+
+            if 'invalid' in flags or 'divide-by-zero' in flags:
+                expected = 'ValueError'
+            elif 'overflow' in flags:
+                expected = 'OverflowError'
+
+            try:
+                got = func(arg)
+            except ValueError:
+                got = 'ValueError'
+            except OverflowError:
+                got = 'OverflowError'
+
+            diff_ulps = 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 diff_ulps <= ALLOWED_ERROR:
+                        continue
+
+            if isinstance(got, str) and isinstance(expected, str):
+                if got == expected:
+                    continue
+
+            fail_msg = fail_fmt.format(id, fn, arg, expected, got)
+            if diff_ulps is not None:
+                fail_msg += ' ({} ulps)'.format(diff_ulps)
+            failures.append(fail_msg)
+
+        if failures:
+            self.fail('Failures in test_mtestfile:\n  ' +
+                      '\n  '.join(failures))
+
+
 def test_main():
     from doctest import DocFileSuite
     suite = unittest.TestSuite()
index e48d8f2e75dcaf45952c49cfc4f60f2598f49a2a..93b8b6e146fa7d088be0ee4e0478ee211143129c 100644 (file)
--- a/Misc/NEWS
+++ b/Misc/NEWS
@@ -201,6 +201,9 @@ Library
 
 Extension Modules
 -----------------
+
+- Issue #3366: Add gamma function 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.
 
index 38d214a5897a6f78fcaf1f3e9db96543163f1c56..249c22765e802f13076ca3c68900411ccf5de91a 100644 (file)
@@ -60,44 +60,265 @@ raised for division by zero and mod by zero.
 extern double copysign(double, double);
 #endif
 
-/* Call is_error when errno != 0, and where x is the result libm
- * returned.  is_error will usually set up an exception and return
- * true (1), but may return false (0) without setting up an exception.
- */
-static int
-is_error(double x)
+/*
+   sin(pi*x), giving accurate results for all finite x (especially x
+   integral or close to an integer).  This is here for use in the
+   reflection formula for the gamma function.  It conforms to IEEE
+   754-2008 for finite arguments, but not for infinities or nans.
+*/
+
+static const double pi = 3.141592653589793238462643383279502884197;
+
+static double
+sinpi(double x)
 {
-       int result = 1; /* presumption of guilt */
-       assert(errno);  /* non-zero errno is a precondition for calling */
-       if (errno == EDOM)
-               PyErr_SetString(PyExc_ValueError, "math domain error");
+       double y, r;
+       int n;
+       /* this function should only ever be called for finite arguments */
+       assert(Py_IS_FINITE(x));
+       y = fmod(fabs(x), 2.0);
+       n = (int)round(2.0*y);
+       assert(0 <= n && n <= 4);
+       switch (n) {
+       case 0:
+               r = sin(pi*y);
+               break;
+       case 1:
+               r = cos(pi*(y-0.5));
+               break;
+       case 2:
+               /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give
+                  -0.0 instead of 0.0 when y == 1.0. */
+               r = sin(pi*(1.0-y));
+               break;
+       case 3:
+               r = -cos(pi*(y-1.5));
+               break;
+       case 4:
+               r = sin(pi*(y-2.0));
+               break;
+       default:
+               assert(0);  /* should never get here */
+               r = -1.23e200; /* silence gcc warning */
+       }
+       return copysign(1.0, x)*r;
+}
 
-       else if (errno == ERANGE) {
-               /* ANSI C generally requires libm functions to set ERANGE
-                * on overflow, but also generally *allows* them to set
-                * ERANGE on underflow too.  There's no consistency about
-                * the latter across platforms.
-                * Alas, C99 never requires that errno be set.
-                * Here we suppress the underflow errors (libm functions
-                * should return a zero on underflow, and +- HUGE_VAL on
-                * overflow, so testing the result for zero suffices to
-                * distinguish the cases).
-                *
-                * On some platforms (Ubuntu/ia64) it seems that errno can be
-                * set to ERANGE for subnormal results that do *not* underflow
-                * to zero.  So to be safe, we'll ignore ERANGE whenever the
-                * function result is less than one in absolute value.
-                */
-               if (fabs(x) < 1.0)
-                       result = 0;
-               else
-                       PyErr_SetString(PyExc_OverflowError,
-                                       "math range error");
+/* Implementation of the real gamma function.  In extensive but non-exhaustive
+   random tests, this function proved accurate to within <= 10 ulps across the
+   entire float domain.  Note that accuracy may depend on the quality of the
+   system math functions, the pow function in particular.  Special cases
+   follow C99 annex F.  The parameters and method are tailored to platforms
+   whose double format is the IEEE 754 binary64 format.
+
+   Method: for x > 0.0 we use the Lanczos approximation with parameters N=13
+   and g=6.024680040776729583740234375; these parameters are amongst those
+   used by the Boost library.  Following Boost (again), we re-express the
+   Lanczos sum as a rational function, and compute it that way.  The
+   coefficients below were computed independently using MPFR, and have been
+   double-checked against the coefficients in the Boost source code.
+
+   For x < 0.0 we use the reflection formula.
+
+   There's one minor tweak that deserves explanation: Lanczos' formula for
+   Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5).  For many x
+   values, x+g-0.5 can be represented exactly.  However, in cases where it
+   can't be represented exactly the small error in x+g-0.5 can be magnified
+   significantly by the pow and exp calls, especially for large x.  A cheap
+   correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error
+   involved in the computation of x+g-0.5 (that is, e = computed value of
+   x+g-0.5 - exact value of x+g-0.5).  Here's the proof:
+
+   Correction factor
+   -----------------
+   Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754
+   double, and e is tiny.  Then:
+
+     pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e)
+     = pow(y, x-0.5)/exp(y) * C,
+
+   where the correction_factor C is given by
+
+     C = pow(1-e/y, x-0.5) * exp(e)
+
+   Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so:
+
+     C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y
+
+   But y-(x-0.5) = g+e, and g+e ~ g.  So we get C ~ 1 + e*g/y, and
+
+     pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y),
+
+   Note that for accuracy, when computing r*C it's better to do
+
+     r + e*g/y*r;
+
+   than
+
+     r * (1 + e*g/y);
+
+   since the addition in the latter throws away most of the bits of
+   information in e*g/y.
+*/
+
+#define LANCZOS_N 13
+static const double lanczos_g = 6.024680040776729583740234375;
+static const double lanczos_g_minus_half = 5.524680040776729583740234375;
+static const double lanczos_num_coeffs[LANCZOS_N] = {
+       23531376880.410759688572007674451636754734846804940,
+       42919803642.649098768957899047001988850926355848959,
+       35711959237.355668049440185451547166705960488635843,
+       17921034426.037209699919755754458931112671403265390,
+       6039542586.3520280050642916443072979210699388420708,
+       1439720407.3117216736632230727949123939715485786772,
+       248874557.86205415651146038641322942321632125127801,
+       31426415.585400194380614231628318205362874684987640,
+       2876370.6289353724412254090516208496135991145378768,
+       186056.26539522349504029498971604569928220784236328,
+       8071.6720023658162106380029022722506138218516325024,
+       210.82427775157934587250973392071336271166969580291,
+       2.5066282746310002701649081771338373386264310793408
+};
+
+/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */
+static const double lanczos_den_coeffs[LANCZOS_N] = {
+       0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
+       13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};
+
+/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */
+#define NGAMMA_INTEGRAL 23
+static const double gamma_integral[NGAMMA_INTEGRAL] = {
+       1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
+       3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
+       1307674368000.0, 20922789888000.0, 355687428096000.0,
+       6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
+       51090942171709440000.0, 1124000727777607680000.0,
+};
+
+/* Lanczos' sum L_g(x), for positive x */
+
+static double
+lanczos_sum(double x)
+{
+       double num = 0.0, den = 0.0;
+       int i;
+       assert(x > 0.0);
+       /* evaluate the rational function lanczos_sum(x).  For large
+          x, the obvious algorithm risks overflow, so we instead
+          rescale the denominator and numerator of the rational
+          function by x**(1-LANCZOS_N) and treat this as a
+          rational function in 1/x.  This also reduces the error for
+          larger x values.  The choice of cutoff point (5.0 below) is
+          somewhat arbitrary; in tests, smaller cutoff values than
+          this resulted in lower accuracy. */
+       if (x < 5.0) {
+               for (i = LANCZOS_N; --i >= 0; ) {
+                       num = num * x + lanczos_num_coeffs[i];
+                       den = den * x + lanczos_den_coeffs[i];
+               }
        }
-       else
-                /* Unexpected math error */
-               PyErr_SetFromErrno(PyExc_ValueError);
-       return result;
+       else {
+               for (i = 0; i < LANCZOS_N; i++) {
+                       num = num / x + lanczos_num_coeffs[i];
+                       den = den / x + lanczos_den_coeffs[i];
+               }
+       }
+       return num/den;
+}
+
+static double
+m_tgamma(double x)
+{
+       double absx, r, y, z, sqrtpow;
+
+       /* special cases */
+       if (!Py_IS_FINITE(x)) {
+               if (Py_IS_NAN(x) || x > 0.0)
+                       return x;  /* tgamma(nan) = nan, tgamma(inf) = inf */
+               else {
+                       errno = EDOM;
+                       return Py_NAN;  /* tgamma(-inf) = nan, invalid */
+               }
+       }
+       if (x == 0.0) {
+               errno = EDOM;
+               return 1.0/x; /* tgamma(+-0.0) = +-inf, divide-by-zero */
+       }
+
+       /* integer arguments */
+       if (x == floor(x)) {
+               if (x < 0.0) {
+                       errno = EDOM;  /* tgamma(n) = nan, invalid for */
+                       return Py_NAN; /* negative integers n */
+               }
+               if (x <= NGAMMA_INTEGRAL)
+                       return gamma_integral[(int)x - 1];
+       }
+       absx = fabs(x);
+
+       /* tiny arguments:  tgamma(x) ~ 1/x for x near 0 */
+       if (absx < 1e-20) {
+               r = 1.0/x;
+               if (Py_IS_INFINITY(r))
+                       errno = ERANGE;
+               return r;
+       }
+
+       /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for
+          x > 200, and underflows to +-0.0 for x < -200, not a negative
+          integer. */
+       if (absx > 200.0) {
+               if (x < 0.0) {
+                       return 0.0/sinpi(x);
+               }
+               else {
+                       errno = ERANGE;
+                       return Py_HUGE_VAL;
+               }
+       }
+
+       y = absx + lanczos_g_minus_half;
+       /* compute error in sum */
+       if (absx > lanczos_g_minus_half) {
+               /* note: the correction can be foiled by an optimizing
+                  compiler that (incorrectly) thinks that an expression like
+                  a + b - a - b can be optimized to 0.0.  This shouldn't
+                  happen in a standards-conforming compiler. */
+               double q = y - absx;
+               z = q - lanczos_g_minus_half;
+       }
+       else {
+               double q = y - lanczos_g_minus_half;
+               z = q - absx;
+       }
+       z = z * lanczos_g / y;
+       if (x < 0.0) {
+               r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx);
+               r -= z * r;
+               if (absx < 140.0) {
+                       r /= pow(y, absx - 0.5);
+               }
+               else {
+                       sqrtpow = pow(y, absx / 2.0 - 0.25);
+                       r /= sqrtpow;
+                       r /= sqrtpow;
+               }
+       }
+       else {
+               r = lanczos_sum(absx) / exp(y);
+               r += z * r;
+               if (absx < 140.0) {
+                       r *= pow(y, absx - 0.5);
+               }
+               else {
+                       sqrtpow = pow(y, absx / 2.0 - 0.25);
+                       r *= sqrtpow;
+                       r *= sqrtpow;
+               }
+       }
+       if (Py_IS_INFINITY(r))
+               errno = ERANGE;
+       return r;
 }
 
 /*
@@ -188,6 +409,46 @@ m_log10(double x)
 }
 
 
+/* Call is_error when errno != 0, and where x is the result libm
+ * returned.  is_error will usually set up an exception and return
+ * true (1), but may return false (0) without setting up an exception.
+ */
+static int
+is_error(double x)
+{
+       int result = 1; /* presumption of guilt */
+       assert(errno);  /* non-zero errno is a precondition for calling */
+       if (errno == EDOM)
+               PyErr_SetString(PyExc_ValueError, "math domain error");
+
+       else if (errno == ERANGE) {
+               /* ANSI C generally requires libm functions to set ERANGE
+                * on overflow, but also generally *allows* them to set
+                * ERANGE on underflow too.  There's no consistency about
+                * the latter across platforms.
+                * Alas, C99 never requires that errno be set.
+                * Here we suppress the underflow errors (libm functions
+                * should return a zero on underflow, and +- HUGE_VAL on
+                * overflow, so testing the result for zero suffices to
+                * distinguish the cases).
+                *
+                * On some platforms (Ubuntu/ia64) it seems that errno can be
+                * set to ERANGE for subnormal results that do *not* underflow
+                * to zero.  So to be safe, we'll ignore ERANGE whenever the
+                * function result is less than one in absolute value.
+                */
+               if (fabs(x) < 1.0)
+                       result = 0;
+               else
+                       PyErr_SetString(PyExc_OverflowError,
+                                       "math range error");
+       }
+       else
+                /* Unexpected math error */
+               PyErr_SetFromErrno(PyExc_ValueError);
+       return result;
+}
+
 /*
    math_1 is used to wrap a libm function f that takes a double
    arguments and returns a double.
@@ -252,6 +513,26 @@ math_1_to_whatever(PyObject *arg, double (*func) (double),
        return (*from_double_func)(r);
 }
 
+/* variant of math_1, to be used when the function being wrapped is known to
+   set errno properly (that is, errno = EDOM for invalid or divide-by-zero,
+   errno = ERANGE for overflow). */
+
+static PyObject *
+math_1a(PyObject *arg, double (*func) (double))
+{
+       double x, r;
+       x = PyFloat_AsDouble(arg);
+       if (x == -1.0 && PyErr_Occurred())
+               return NULL;
+       errno = 0;
+       PyFPE_START_PROTECT("in math_1a", return 0);
+       r = (*func)(x);
+       PyFPE_END_PROTECT(r);
+       if (errno && is_error(r))
+               return NULL;
+       return PyFloat_FromDouble(r);
+}
+
 /*
    math_2 is used to wrap a libm function f that takes two double
    arguments and returns a double.
@@ -330,6 +611,12 @@ math_2(PyObject *args, double (*func) (double, double), char *funcname)
        }\
         PyDoc_STRVAR(math_##funcname##_doc, docstring);
 
+#define FUNC1A(funcname, func, docstring)                              \
+       static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
+               return math_1a(args, func);                             \
+       }\
+        PyDoc_STRVAR(math_##funcname##_doc, docstring);
+
 #define FUNC2(funcname, func, docstring) \
        static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
                return math_2(args, func, #funcname); \
@@ -405,6 +692,8 @@ PyDoc_STRVAR(math_floor_doc,
             "floor(x)\n\nReturn the floor of x as an int.\n"
             "This is the largest integral value <= x.");
 
+FUNC1A(gamma, m_tgamma,
+      "gamma(x)\n\nGamma 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.")
@@ -1150,6 +1439,7 @@ static PyMethodDef math_methods[] = {
        {"fmod",        math_fmod,      METH_VARARGS,   math_fmod_doc},
        {"frexp",       math_frexp,     METH_O,         math_frexp_doc},
        {"fsum",        math_fsum,      METH_O,         math_fsum_doc},
+       {"gamma",       math_gamma,     METH_O,         math_gamma_doc},
        {"hypot",       math_hypot,     METH_VARARGS,   math_hypot_doc},
        {"isinf",       math_isinf,     METH_O,         math_isinf_doc},
        {"isnan",       math_isnan,     METH_O,         math_isnan_doc},