mirror of https://github.com/python/cpython.git
Issue #3366: Add gamma function to math module.
(lgamma, erf and erfc to follow).
This commit is contained in:
parent
ddfb6cdc2b
commit
b93fff0a57
|
@ -308,6 +308,16 @@ Hyperbolic functions
|
||||||
Return the hyperbolic tangent of *x*.
|
Return the hyperbolic tangent of *x*.
|
||||||
|
|
||||||
|
|
||||||
|
Special functions
|
||||||
|
-----------------
|
||||||
|
|
||||||
|
.. function:: gamma(x)
|
||||||
|
|
||||||
|
Return the Gamma function at *x*.
|
||||||
|
|
||||||
|
.. versionadded:: 2.7
|
||||||
|
|
||||||
|
|
||||||
Constants
|
Constants
|
||||||
---------
|
---------
|
||||||
|
|
||||||
|
|
|
@ -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
|
|
@ -7,6 +7,7 @@
|
||||||
import os
|
import os
|
||||||
import sys
|
import sys
|
||||||
import random
|
import random
|
||||||
|
import struct
|
||||||
|
|
||||||
eps = 1E-05
|
eps = 1E-05
|
||||||
NAN = float('nan')
|
NAN = float('nan')
|
||||||
|
@ -29,8 +30,50 @@
|
||||||
else:
|
else:
|
||||||
file = __file__
|
file = __file__
|
||||||
test_dir = os.path.dirname(file) or os.curdir
|
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')
|
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):
|
def parse_testfile(fname):
|
||||||
"""Parse a file with test values
|
"""Parse a file with test values
|
||||||
|
|
||||||
|
@ -887,6 +930,51 @@ def test_testfile(self):
|
||||||
self.fail(message)
|
self.fail(message)
|
||||||
self.ftest("%s:%s(%r)" % (id, fn, ar), result, er)
|
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():
|
def test_main():
|
||||||
from doctest import DocFileSuite
|
from doctest import DocFileSuite
|
||||||
suite = unittest.TestSuite()
|
suite = unittest.TestSuite()
|
||||||
|
|
|
@ -1341,6 +1341,8 @@ C-API
|
||||||
Extension Modules
|
Extension Modules
|
||||||
-----------------
|
-----------------
|
||||||
|
|
||||||
|
- Issue #3366: Add gamma function 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
|
||||||
range.
|
range.
|
||||||
|
|
|
@ -60,44 +60,265 @@ raised for division by zero and mod by zero.
|
||||||
extern double copysign(double, double);
|
extern double copysign(double, double);
|
||||||
#endif
|
#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
|
sin(pi*x), giving accurate results for all finite x (especially x
|
||||||
* true (1), but may return false (0) without setting up an exception.
|
integral or close to an integer). This is here for use in the
|
||||||
*/
|
reflection formula for the gamma function. It conforms to IEEE
|
||||||
static int
|
754-2008 for finite arguments, but not for infinities or nans.
|
||||||
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) {
|
static const double pi = 3.141592653589793238462643383279502884197;
|
||||||
/* ANSI C generally requires libm functions to set ERANGE
|
|
||||||
* on overflow, but also generally *allows* them to set
|
static double
|
||||||
* ERANGE on underflow too. There's no consistency about
|
sinpi(double x)
|
||||||
* the latter across platforms.
|
{
|
||||||
* Alas, C99 never requires that errno be set.
|
double y, r;
|
||||||
* Here we suppress the underflow errors (libm functions
|
int n;
|
||||||
* should return a zero on underflow, and +- HUGE_VAL on
|
/* this function should only ever be called for finite arguments */
|
||||||
* overflow, so testing the result for zero suffices to
|
assert(Py_IS_FINITE(x));
|
||||||
* distinguish the cases).
|
y = fmod(fabs(x), 2.0);
|
||||||
*
|
n = (int)round(2.0*y);
|
||||||
* On some platforms (Ubuntu/ia64) it seems that errno can be
|
assert(0 <= n && n <= 4);
|
||||||
* set to ERANGE for subnormal results that do *not* underflow
|
switch (n) {
|
||||||
* to zero. So to be safe, we'll ignore ERANGE whenever the
|
case 0:
|
||||||
* function result is less than one in absolute value.
|
r = sin(pi*y);
|
||||||
*/
|
break;
|
||||||
if (fabs(x) < 1.0)
|
case 1:
|
||||||
result = 0;
|
r = cos(pi*(y-0.5));
|
||||||
else
|
break;
|
||||||
PyErr_SetString(PyExc_OverflowError,
|
case 2:
|
||||||
"math range error");
|
/* 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 */
|
||||||
}
|
}
|
||||||
else
|
return copysign(1.0, x)*r;
|
||||||
/* Unexpected math error */
|
}
|
||||||
PyErr_SetFromErrno(PyExc_ValueError);
|
|
||||||
return result;
|
/* 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 {
|
||||||
|
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
|
math_1 is used to wrap a libm function f that takes a double
|
||||||
arguments and returns a double.
|
arguments and returns a double.
|
||||||
|
@ -247,6 +508,26 @@ math_1(PyObject *arg, double (*func) (double), int can_overflow)
|
||||||
return PyFloat_FromDouble(r);
|
return PyFloat_FromDouble(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
|
math_2 is used to wrap a libm function f that takes two double
|
||||||
arguments and returns a double.
|
arguments and returns a double.
|
||||||
|
@ -313,6 +594,12 @@ math_2(PyObject *args, double (*func) (double, double), char *funcname)
|
||||||
}\
|
}\
|
||||||
PyDoc_STRVAR(math_##funcname##_doc, docstring);
|
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) \
|
#define FUNC2(funcname, func, docstring) \
|
||||||
static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
|
static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
|
||||||
return math_2(args, func, #funcname); \
|
return math_2(args, func, #funcname); \
|
||||||
|
@ -350,6 +637,8 @@ FUNC1(fabs, fabs, 0,
|
||||||
FUNC1(floor, floor, 0,
|
FUNC1(floor, floor, 0,
|
||||||
"floor(x)\n\nReturn the floor of x as a float.\n"
|
"floor(x)\n\nReturn the floor of x as a float.\n"
|
||||||
"This is the largest integral value <= x.")
|
"This is the largest integral value <= x.")
|
||||||
|
FUNC1A(gamma, m_tgamma,
|
||||||
|
"gamma(x)\n\nGamma 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.")
|
||||||
|
@ -1077,6 +1366,7 @@ static PyMethodDef math_methods[] = {
|
||||||
{"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
|
{"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
|
||||||
{"frexp", math_frexp, METH_O, math_frexp_doc},
|
{"frexp", math_frexp, METH_O, math_frexp_doc},
|
||||||
{"fsum", math_fsum, METH_O, math_fsum_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},
|
{"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
|
||||||
{"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},
|
||||||
|
|
Loading…
Reference in New Issue