Index: PCbuild/pythoncore.vcproj
===================================================================
--- PCbuild/pythoncore.vcproj (revision 74993)
+++ PCbuild/pythoncore.vcproj (working copy)
@@ -1091,6 +1091,14 @@
>
+
+
+
+
Index: setup.py
===================================================================
--- setup.py (revision 74993)
+++ setup.py (working copy)
@@ -414,7 +414,7 @@
libraries=math_libs) )
# math library functions, e.g. sin()
- exts.append( Extension('math', ['mathmodule.c'],
+ exts.append( Extension('math', ['mathmodule.c', 'math_extra.c'],
libraries=math_libs) )
# fast string operations implemented in C
exts.append( Extension('strop', ['stropmodule.c']) )
Index: Misc/NEWS
===================================================================
--- Misc/NEWS (revision 74993)
+++ Misc/NEWS (working copy)
@@ -1335,6 +1335,8 @@
Extension Modules
-----------------
+- Issue #3366: Add Gamma function to math module (math.gamma).
+
- Issue #6877: Make it possible to link the readline extension to libedit
on OSX.
Index: PC/VS7.1/pythoncore.vcproj
===================================================================
--- PC/VS7.1/pythoncore.vcproj (revision 74993)
+++ PC/VS7.1/pythoncore.vcproj (working copy)
@@ -659,6 +659,12 @@
RelativePath="..\..\Python\marshal.c">
+
+
+
+
+
+
+
+
Index: PC/VC6/pythoncore.dsp
===================================================================
--- PC/VC6/pythoncore.dsp (revision 74993)
+++ PC/VC6/pythoncore.dsp (working copy)
@@ -519,6 +519,10 @@
# End Source File
# Begin Source File
+SOURCE=..\..\Modules\math_extra.c
+# End Source File
+# Begin Source File
+
SOURCE=..\..\Modules\mathmodule.c
# End Source File
# Begin Source File
Index: Doc/library/math.rst
===================================================================
--- Doc/library/math.rst (revision 74993)
+++ Doc/library/math.rst (working copy)
@@ -308,6 +308,16 @@
Return the hyperbolic tangent of *x*.
+Special functions
+-----------------
+
+.. function:: gamma(x)
+
+ Return the Gamma function at *x*.
+
+ .. versionadded:: 2.7
+
+
Constants
---------
Index: Lib/test/test_math.py
===================================================================
--- Lib/test/test_math.py (revision 74993)
+++ Lib/test/test_math.py (working copy)
@@ -7,6 +7,7 @@
import os
import sys
import random
+import struct
eps = 1E-05
NAN = float('nan')
@@ -29,8 +30,50 @@
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(' 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
@@ -887,6 +930,51 @@
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: Lib/test/math_testcases.txt
===================================================================
--- Lib/test/math_testcases.txt (revision 0)
+++ Lib/test/math_testcases.txt (revision 0)
@@ -0,0 +1,146 @@
+-- Testcases for functions in math.
+--
+-- Each line takes the form:
+--
+-- ->
+--
+-- where:
+--
+-- is a short name identifying the test,
+--
+-- is the function to be tested (exp, cos, asinh, ...),
+--
+-- is a string representing a floating-point value
+--
+-- is the expected (ideal) output value, again
+-- represented as a string.
+--
+-- 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: Modules/math_extra.h
===================================================================
--- Modules/math_extra.h (revision 0)
+++ Modules/math_extra.h (revision 0)
@@ -0,0 +1,9 @@
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+double py_tgamma(double);
+
+#ifdef __cplusplus
+}
+#endif
Index: Modules/mathmodule.c
===================================================================
--- Modules/mathmodule.c (revision 74993)
+++ Modules/mathmodule.c (working copy)
@@ -54,6 +54,7 @@
#include "Python.h"
#include "longintrepr.h" /* just for SHIFT */
+#include "math_extra.h"
#ifdef _OSF_SOURCE
/* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
@@ -247,6 +248,26 @@
return PyFloat_FromDouble(r);
}
+/* variant of math_1, to be used when the function being wrapped is known
+ to set errno properly. In particular, math_1a should be used to wrap
+ functions defined in the Modules/math_extra.c file. */
+
+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.
@@ -313,6 +334,12 @@
}\
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); \
@@ -350,6 +377,8 @@
FUNC1(floor, floor, 0,
"floor(x)\n\nReturn the floor of x as a float.\n"
"This is the largest integral value <= x.")
+FUNC1A(gamma, py_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.")
@@ -1077,6 +1106,7 @@
{"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},
Index: Modules/Setup.dist
===================================================================
--- Modules/Setup.dist (revision 74993)
+++ Modules/Setup.dist (working copy)
@@ -169,7 +169,7 @@
#array arraymodule.c # array objects
#cmath cmathmodule.c # -lm # complex math library functions
-#math mathmodule.c # -lm # math library functions, e.g. sin()
+#math mathmodule.c math_extra.c # -lm # math library functions, e.g. sin()
#_struct _struct.c # binary structure packing/unpacking
#time timemodule.c # -lm # time operations and variables
#operator operator.c # operator.add() and similar goodies
Index: Modules/math_extra.c
===================================================================
--- Modules/math_extra.c (revision 0)
+++ Modules/math_extra.c (revision 0)
@@ -0,0 +1,269 @@
+/* This file contains substitutes for various C99 math functions, for
+ those platforms whose math library doesn't define these functions. */
+
+/* Note on exceptional cases: ideally, the functions defined in this file
+ should behave as follows on IEEE 754 hardware:
+
+ - where C99 Annex G recommends signalling divide-by-zero,
+ return +-infinity and set errno = EDOM
+ - where C99 Annex G recommends signalling invalid,
+ return NaN and set errno = EDOM
+ - where C99 Annex G recommends signalling overflow,
+ return +-infinity and set errno = ERANGE
+ - in all other cases, don't touch errno
+
+ In general, Python needs both the return value *and*
+ the errno value to determine what to do.
+
+ */
+
+#include "Python.h"
+
+static const double pi = 3.141592653589793238462643383279502884197;
+
+/* sin(pi*x), giving accurate results for all x (especially x
+ integral or close to an integer). */
+static double
+sinpi(double x)
+{
+ 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;
+}
+
+/* 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 for platforms
+ whose float type 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 are taken directly from the Boost source code. For x <
+ 0.0 we use the usual 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_NTERMS 13
+static const double lanczos_g = 6.024680040776729583740234375;
+static const double g_minus_half = 5.524680040776729583740234375;
+
+static const double lanczos_num_coeffs[LANCZOS_NTERMS] = {
+ 23531376880.41075968857200767445163675473,
+ 42919803642.64909876895789904700198885093,
+ 35711959237.35566804944018545154716670596,
+ 17921034426.03720969991975575445893111267,
+ 6039542586.352028005064291644307297921070,
+ 1439720407.311721673663223072794912393972,
+ 248874557.8620541565114603864132294232163,
+ 31426415.58540019438061423162831820536287,
+ 2876370.628935372441225409051620849613599,
+ 186056.2653952234950402949897160456992822,
+ 8071.672002365816210638002902272250613822,
+ 210.8242777515793458725097339207133627117,
+ 2.506628274631000270164908177133837338626
+};
+
+/* denominator is x*(x+1)*...*(x+11) */
+static const double lanczos_den_coeffs[LANCZOS_NTERMS] = {
+ 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 Lg(double x) {
+ double num = 0.0, den = 0.0;
+ int i;
+ assert(x > 0.0);
+ /* evaluate the rational function Lg(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_NTERMS) 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_NTERMS; --i >= 0; ) {
+ num = num * x + lanczos_num_coeffs[i];
+ den = den * x + lanczos_den_coeffs[i];
+ }
+ }
+ else {
+ for (i = 0; i < LANCZOS_NTERMS; i++) {
+ num = num / x + lanczos_num_coeffs[i];
+ den = den / x + lanczos_den_coeffs[i];
+ }
+ }
+ return num/den;
+}
+
+double py_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;
+ return Py_NAN; /* tgamma(n) = nan, invalid for neg. integer 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 + g_minus_half;
+ /* compute error in sum */
+ if (absx > 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 standard-conforming
+ compiler. */
+ double q = y - absx;
+ z = q - g_minus_half;
+ }
+ else {
+ double q = y - g_minus_half;
+ z = q - absx;
+ }
+ z = z * lanczos_g / y;
+ if (x < 0.0) {
+ r = -pi / sinpi(absx) / absx * exp(y) / Lg(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 = Lg(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;
+}