Index: PCbuild/pythoncore.vcproj
===================================================================
--- PCbuild/pythoncore.vcproj (revision 74666)
+++ PCbuild/pythoncore.vcproj (working copy)
@@ -1091,6 +1091,14 @@
>
+
+
+
+
Index: setup.py
===================================================================
--- setup.py (revision 74666)
+++ 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_cextra.c'],
libraries=math_libs) )
# fast string operations implemented in C
exts.append( Extension('strop', ['stropmodule.c']) )
Index: PC/VS7.1/pythoncore.vcproj
===================================================================
--- PC/VS7.1/pythoncore.vcproj (revision 74666)
+++ 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 74666)
+++ PC/VC6/pythoncore.dsp (working copy)
@@ -519,6 +519,10 @@
# End Source File
# Begin Source File
+SOURCE=..\..\Modules\math_cextra.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 74666)
+++ Doc/library/math.rst (working copy)
@@ -308,6 +308,30 @@
Return the hyperbolic tangent of *x*.
+Special functions
+-----------------
+
+.. function:: gamma(x)
+
+ Return the Gamma function at *x*.
+
+
+.. function:: lgamma(x)
+
+ The lgamma function computes the natural log of the absolute value
+ of the Gamma function.
+
+
+.. function:: erf(x)
+
+ Return the error function at *x*.
+
+
+.. function:: erfc(x)
+
+ Return the complementary error function at *x*.
+
+
Constants
---------
Index: Lib/test/test_math.py
===================================================================
--- Lib/test/test_math.py (revision 74666)
+++ 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')
@@ -19,16 +20,54 @@
HAVE_DOUBLE_ROUNDING = (x + y == 1e16 + 4)
# locate file with test values
-if __name__ == '__main__':
- file = sys.argv[0]
-else:
- file = __file__
-test_dir = os.path.dirname(file) or os.curdir
-test_file = os.path.join(test_dir, 'cmath_testcases.txt')
+test_dir = os.path.dirname(__file__) or os.curdir
+math_testcases = os.path.join(test_dir, 'math_testcases.txt')
+cmath_testcases = os.path.join(test_dir, 'cmath_testcases.txt')
-def parse_testfile(fname):
+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_ctestfile(fname):
+ """Parse a file with test values
+
Empty lines or lines starting with -- are ignored
yields id, fn, arg_real, arg_imag, exp_real, exp_imag
"""
@@ -860,10 +899,10 @@
else:
self.fail("sqrt(-1) didn't raise ValueError")
- def test_testfile(self):
- if not float.__getformat__("double").startswith("IEEE"):
- return
- for id, fn, ar, ai, er, ei, flags in parse_testfile(test_file):
+ @unittest.skipUnless(float.__getformat__("double").startswith("IEEE"),
+ "test requires IEEE 754 doubles")
+ def test_ctestfile(self):
+ for id, fn, ar, ai, er, ei, flags in parse_ctestfile(cmath_testcases):
# Skip if either the input or result is complex, or if
# flags is nonempty
if ai != 0. or ei != 0. or flags:
@@ -884,6 +923,50 @@
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 = 100 # 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,308 @@
+-- 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.
+
+-------------------------
+-- erf: error function --
+-------------------------
+
+erf0000 erf 0.0 -> 0.0
+erf0001 erf -0.0 -> -0.0
+erf0002 erf inf -> 1.0
+erf0003 erf -inf -> -1.0
+erf0004 erf nan -> nan
+
+-- tiny values
+erf0010 erf 1e-308 -> 1.1283791670955125e-308
+erf0011 erf 5e-324 -> 4.9406564584124654e-324
+erf0012 erf 1e-10 -> 1.1283791670955126e-10
+
+-- small integers
+erf0020 erf 1 -> 0.84270079294971489
+erf0021 erf 2 -> 0.99532226501895271
+erf0022 erf 3 -> 0.99997790950300136
+erf0023 erf 4 -> 0.99999998458274209
+erf0024 erf 5 -> 0.99999999999846256
+erf0025 erf 6 -> 1.0
+
+erf0030 erf -1 -> -0.84270079294971489
+erf0031 erf -2 -> -0.99532226501895271
+erf0032 erf -3 -> -0.99997790950300136
+erf0033 erf -4 -> -0.99999998458274209
+erf0034 erf -5 -> -0.99999999999846256
+erf0035 erf -6 -> -1.0
+
+-- huge values should all go to +/-1, depending on sign
+erf0040 erf -40 -> -1.0
+erf0041 erf 1e16 -> 1.0
+erf0042 erf -1e150 -> -1.0
+erf0043 erf 1.7e308 -> 1.0
+
+
+----------------------------------------
+-- erfc: complementary error function --
+----------------------------------------
+
+erfc0000 erfc 0.0 -> 1.0
+erfc0001 erfc -0.0 -> 1.0
+erfc0002 erfc inf -> 0.0
+erfc0003 erfc -inf -> 2.0
+erfc0004 erfc nan -> nan
+
+-- tiny values
+erfc0010 erfc 1e-308 -> 1.0
+erfc0011 erfc 5e-324 -> 1.0
+erfc0012 erfc 1e-10 -> 0.99999999988716204
+
+-- small integers
+erfc0020 erfc 1 -> 0.15729920705028513
+erfc0021 erfc 2 -> 0.0046777349810472662
+erfc0022 erfc 3 -> 2.2090496998585441e-05
+erfc0023 erfc 4 -> 1.541725790028002e-08
+erfc0024 erfc 5 -> 1.5374597944280349e-12
+erfc0025 erfc 6 -> 2.1519736712498913e-17
+
+erfc0030 erfc -1 -> 1.8427007929497148
+erfc0031 erfc -2 -> 1.9953222650189528
+erfc0032 erfc -3 -> 1.9999779095030015
+erfc0033 erfc -4 -> 1.9999999845827421
+erfc0034 erfc -5 -> 1.9999999999984626
+erfc0035 erfc -6 -> 2.0
+
+-- as x -> infinity, erfc(x) behaves like exp(-x*x)/x/sqrt(pi)
+erfc0040 erfc 20 -> 5.3958656116079012e-176
+erfc0041 erfc 25 -> 8.3001725711965228e-274
+erfc0042 erfc 27 -> 5.2370464393526292e-319
+erfc0043 erfc 28 -> 0.0
+
+-- huge values
+erfc0050 erfc -40 -> 2.0
+erfc0051 erfc 1e16 -> 0.0
+erfc0052 erfc -1e150 -> 2.0
+erfc0053 erfc 1.7e308 -> 0.0
+
+
+
+---------------------------------------------------------
+-- 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
+
+-- 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 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
+
+---------------------------
+-- 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
+
Index: Lib/test/test_cmath.py
===================================================================
--- Lib/test/test_cmath.py (revision 74666)
+++ Lib/test/test_cmath.py (working copy)
@@ -1,5 +1,5 @@
from test.test_support import run_unittest
-from test.test_math import parse_testfile, test_file
+from test.test_math import parse_ctestfile, cmath_testcases
import unittest
import os, sys
import cmath, math
@@ -88,7 +88,7 @@
test_functions.append(lambda x : cmath.log(14.-27j, x))
def setUp(self):
- self.test_values = open(test_file)
+ self.test_values = open(cmath_testcases)
def tearDown(self):
self.test_values.close()
@@ -310,7 +310,7 @@
two floats."""
return complex(*polar(z))
- for id, fn, ar, ai, er, ei, flags in parse_testfile(test_file):
+ for id, fn, ar, ai, er, ei, flags in parse_ctestfile(cmath_testcases):
arg = complex(ar, ai)
expected = complex(er, ei)
if fn == 'rect':
Index: Modules/math_cextra.c
===================================================================
--- Modules/math_cextra.c (revision 0)
+++ Modules/math_cextra.c (revision 0)
@@ -0,0 +1,656 @@
+/* This file defines substitutes for the C99 gamma, lgamma, erf and
+ erfc functions, for those platforms whose math library doesn't
+ already have them.
+
+ Exceptional cases: ideally, the functions defined in the file
+ should behave as follows (on IEEE 754 machines)
+
+ - 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 set errno
+
+ */
+
+#include "Python.h"
+
+#if defined(HAVE_LGAMMA) && defined(HAVE_TGAMMA) && \
+ defined(HAVE_ERF) && defined(HAVE_ERFC) && !defined(TEST_GAMMA)
+#define HAVE_SPECIAL_FUNCTIONS
+#endif
+
+#ifdef HAVE_SPECIAL_FUNCTIONS
+
+double py_erf(double x) {
+ return erf(x);
+}
+
+double py_erfc(double x) {
+ return erfc(x);
+}
+
+double py_lgamma(double x) {
+ /* deal with special cases, then delegate to the libm version */
+ double r;
+ if (!Py_IS_FINITE(x)) {
+ if (Py_IS_NAN(x))
+ return x; /* lgamma(nan) = nan */
+ else
+ return Py_HUGE_VAL; /* lgamma(+-inf) = inf */
+ }
+
+ if (x <= 0.0 && x == floor(x)) {
+ errno = EDOM;
+ return Py_HUGE_VAL; /* lgamma(x) = inf, divide-by-zero for integer x <= 0 */
+ }
+ r = lgamma(x);
+ /* make sure that we set errno on overflow */
+ if (Py_IS_INFINITY(r))
+ errno = ERANGE;
+ return r;
+}
+
+double py_tgamma(double x) {
+ /* deal with special cases, then delegate to the libm version */
+ double r;
+ if (!Py_IS_FINITE(x)) {
+ if (Py_IS_NAN(x) || x > 0.0) /* tgamma(nan) = nan, tgamma(inf) = inf */
+ return x;
+ else {
+ errno = EDOM;
+ return Py_NAN; /* tgamma(-inf) = nan, invalid */
+ }
+ }
+
+ if (x <= 0.0) {
+ if (x == 0.0) {
+ errno = EDOM;
+ return 1.0/x; /* tgamma(+-0) = +-inf, divide-by-zero */
+ }
+ else if (x == floor(x)) {
+ errno = EDOM;
+ return Py_NAN; /* tgamma(x) = nan, invalid for integer x < 0 */
+ }
+ }
+ r = tgamma(x);
+ /* make sure to set errno on overflow */
+ if (Py_IS_INFINITY(r))
+ errno = ERANGE;
+ return r;
+}
+
+#else /* HAVE_SPECIAL_FUNCTIONS */
+/* erf(x)
+ * erfc(x)
+ * original code from: http://sourceware.org/cgi-bin/cvsweb.cgi/~checkout~/src/newlib/libm/math/s_erf.c?rev=1.1.1.1&cvsroot=src
+ * x
+ * 2 |\
+ * erf(x) = --------- | exp(-t*t)dt
+ * sqrt(pi) \|
+ * 0
+ *
+ * erfc(x) = 1-erf(x)
+ * Note that
+ * erf(-x) = -erf(x)
+ * erfc(-x) = 2 - erfc(x)
+ *
+ * Method:
+ * 1. For |x| in [0, 0.84375]
+ * erf(x) = x + x*R(x^2)
+ * erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
+ * = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
+ * where R = P/Q where P is an odd poly of degree 8 and
+ * Q is an odd poly of degree 10.
+ * -57.90
+ * | R - (erf(x)-x)/x | <= 2
+ *
+ *
+ * Remark. The formula is derived by noting
+ * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
+ * and that
+ * 2/sqrt(pi) = 1.128379167095512573896158903121545171688
+ * is close to one. The interval is chosen because the fix
+ * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
+ * near 0.6174), and by some experiment, 0.84375 is chosen to
+ * guarantee the error is less than one ulp for erf.
+ *
+ * 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
+ * c = 0.84506291151 rounded to single (24 bits)
+ * erf(x) = sign(x) * (c + P1(s)/Q1(s))
+ * erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
+ * 1+(c+P1(s)/Q1(s)) if x < 0
+ * |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
+ * Remark: here we use the taylor series expansion at x=1.
+ * erf(1+s) = erf(1) + s*Poly(s)
+ * = 0.845.. + P1(s)/Q1(s)
+ * That is, we use rational approximation to approximate
+ * erf(1+s) - (c = (single)0.84506291151)
+ * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
+ * where
+ * P1(s) = degree 6 poly in s
+ * Q1(s) = degree 6 poly in s
+ *
+ * 3. For x in [1.25,1/0.35(~2.857143)],
+ * erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
+ * erf(x) = 1 - erfc(x)
+ * where
+ * R1(z) = degree 7 poly in z, (z=1/x^2)
+ * S1(z) = degree 8 poly in z
+ *
+ * 4. For x in [1/0.35,28]
+ * erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
+ * = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6 x >= 28
+ * erf(x) = sign(x) *(1 - tiny) (raise inexact)
+ * erfc(x) = tiny*tiny (raise underflow) if x > 0
+ * = 2 - tiny if x<0
+ *
+ * 7. Special case:
+ * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
+ * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
+ * erfc/erf(NaN) is NaN
+ */
+
+static const double
+tiny = 1e-300,
+half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
+one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
+ /* c = (float)0.84506291151 */
+erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
+/*
+ * Coefficients for approximation to erf on [0,0.84375]
+ */
+efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
+efx8= 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
+pp0 = 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
+pp1 = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
+pp2 = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
+pp3 = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
+pp4 = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
+qq1 = 3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
+qq2 = 6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
+qq3 = 5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
+qq4 = 1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
+qq5 = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
+/*
+ * Coefficients for approximation to erf in [0.84375,1.25]
+ */
+pa0 = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
+pa1 = 4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
+pa2 = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
+pa3 = 3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
+pa4 = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
+pa5 = 3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
+pa6 = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
+qa1 = 1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
+qa2 = 5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
+qa3 = 7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
+qa4 = 1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
+qa5 = 1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
+qa6 = 1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
+/*
+ * Coefficients for approximation to erfc in [1.25,1/0.35]
+ */
+ra0 = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
+ra1 = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
+ra2 = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
+ra3 = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
+ra4 = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
+ra5 = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
+ra6 = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
+ra7 = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
+sa1 = 1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
+sa2 = 1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
+sa3 = 4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
+sa4 = 6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
+sa5 = 4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
+sa6 = 1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
+sa7 = 6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
+sa8 = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
+/*
+ * Coefficients for approximation to erfc in [1/.35,28]
+ */
+rb0 = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
+rb1 = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
+rb2 = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
+rb3 = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
+rb4 = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
+rb5 = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
+rb6 = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
+sb1 = 3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
+sb2 = 3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
+sb3 = 1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
+sb4 = 3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
+sb5 = 2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
+sb6 = 4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
+sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
+
+double
+py_erf(double x)
+{
+ int hx,ix,i;
+ double R,S,P,Q,s,y,z,r,f;
+ f = frexp(x,&hx);
+ ix = hx&0x7fffffff;
+ if(ix>=1023) { /* erf(nan)=nan */
+ i = ((unsigned int)hx>>31)<<1;
+ return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */
+ }
+
+ if(fabs(x) < 0.84375) { /* |x|<0.84375 */
+ if(fabs(hx) > 28) { /* |x|<2**-28 */
+ if (fabs(hx) > 57)
+ return 0.125*(8.0*x+efx8*x); /*avoid underflow */
+ return x + efx*x;
+ }
+ z = x*x;
+ r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
+ s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
+ y = r/s;
+ return x + x*y;
+ }
+ if(fabs(x) < 1.25) { /* 0.84375 <= |x| < 1.25 */
+ s = fabs(x)-one;
+ P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
+ Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
+ if(f>0) return erx + P/Q; else return -erx - P/Q;
+ }
+ if (fabs(x) >= 6) { /* inf>|x|>=6 */
+ if(f>0) return one-tiny; else return tiny-one;
+ }
+ x = fabs(x);
+ s = one/(x*x);
+ if(fabs(x)<1/0.35) { /* |x| < 1/0.35 */
+ R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
+ ra5+s*(ra6+s*ra7))))));
+ S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
+ sa5+s*(sa6+s*(sa7+s*sa8)))))));
+ } else { /* |x| >= 1/0.35 */
+ R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
+ rb5+s*rb6)))));
+ S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
+ sb5+s*(sb6+s*sb7))))));
+ }
+ z = ldexp(x,0);
+ r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
+ if(f>0) return one-r/x; else return r/x-one;
+}
+
+double
+py_erfc(double x)
+{
+ int hx,ix;
+ double R,S,P,Q,s,y,z,r,f;
+ f = frexp(x,&hx);
+ ix = hx&0x7fffffff;
+ if(ix>=0x7ff00000) { /* erfc(nan)=nan */
+ /* erfc(+-inf)=0,2 */
+ return (double)(((unsigned int)hx>>31)<<1)+one/x;
+ }
+
+ if(fabs(x) < 0.84375) { /* |x|<0.84375 */
+ if(fabs(hx) > 56) /* |x|<2**-56 */
+ return one-x;
+ z = x*x;
+ r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
+ s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
+ y = r/s;
+ if(x < 1/4.) { /* x<1/4 */
+ return one-(x+x*y);
+ } else {
+ r = x*y;
+ r += (x-half);
+ return half - r ;
+ }
+ }
+ if(fabs(x) < 1.25) { /* 0.84375 <= |x| < 1.25 */
+ s = fabs(x)-one;
+ P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
+ Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
+ if( f>=0) {
+ z = one-erx; return z - P/Q;
+ } else {
+ z = erx+P/Q; return one+z;
+ }
+ }
+ if (fabs(x)<28) { /* |x|<28 */
+ x = fabs(x);
+ s = one/(x*x);
+ if(fabs(x)<1/.35) { /* |x| < 1/.35 ~ 2.857143*/
+ R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
+ ra5+s*(ra6+s*ra7))))));
+ S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
+ sa5+s*(sa6+s*(sa7+s*sa8)))))));
+ } else { /* |x| >= 1/.35 ~ 2.857143 */
+ if(x<-6) return two-tiny;/* x < -6 */
+ R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
+ rb5+s*rb6)))));
+ S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
+ sb5+s*(sb6+s*sb7))))));
+ }
+ z = ldexp(x,0);
+ r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
+ if(f>0) return r/x; else return two-r/x;
+ } else {
+ if(f>0) return tiny*tiny; else return two-tiny;
+ }
+}
+
+/* lgamma(x)
+ * gamma(x)
+ * original code from: http://www.sourceware.org/cgi-bin/cvsweb.cgi/~checkout~/src/newlib/libm/mathfp/er_lgamma.c?rev=1.6&content-type=text/plain&cvsroot=src
+ *
+ * Method:
+ * 1. Argument Reduction for 0 < x <= 8
+ * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
+ * reduce x to a number in [1.5,2.5] by
+ * lgamma(1+s) = log(s) + lgamma(s)
+ * for example,
+ * lgamma(7.3) = log(6.3) + lgamma(6.3)
+ * = log(6.3*5.3) + lgamma(5.3)
+ * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
+ * 2. Polynomial approximation of lgamma around its
+ * minimun ymin=1.461632144968362245 to maintain monotonicity.
+ * On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
+ * Let z = x-ymin;
+ * lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
+ * where
+ * poly(z) is a 14 degree polynomial.
+ * 2. Rational approximation in the primary interval [2,3]
+ * We use the following approximation:
+ * s = x-2.0;
+ * lgamma(x) = 0.5*s + s*P(s)/Q(s)
+ * with accuracy
+ * |P/Q - (lgamma(x)-0.5s)| < 2**-61.71
+ * Our algorithms are based on the following observation
+ *
+ * zeta(2)-1 2 zeta(3)-1 3
+ * lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ...
+ * 2 3
+ *
+ * where Euler = 0.5771... is the Euler constant, which is very
+ * close to 0.5.
+ *
+ * 3. For x>=8, we have
+ * lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
+ * (better formula:
+ * lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
+ * Let z = 1/x, then we approximation
+ * f(z) = lgamma(x) - (x-0.5)(log(x)-1)
+ * by
+ * 3 5 11
+ * w = w0 + w1*z + w2*z + w3*z + ... + w6*z
+ * where
+ * |w - f(z)| < 2**-58.74
+ *
+ * 4. For negative x, since (G is gamma function)
+ * -x*G(-x)*G(x) = pi/sin(pi*x),
+ * we have
+ * G(x) = pi/(sin(pi*x)*(-x)*G(-x))
+ * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
+ * Hence, for x<0, signgam = sign(sin(pi*x)) and
+ * lgamma(x) = log(|Gamma(x)|)
+ * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
+ * Note: one should avoid compute pi*(-x) directly in the
+ * computation of sin(pi*(-x)).
+ *
+ * 5. Special Cases
+ * lgamma(2+s) ~ s*(1-Euler) for tiny s
+ * lgamma(1)=lgamma(2)=0
+ * lgamma(x) ~ -log(x) for tiny x
+ * lgamma(0) = lgamma(inf) = inf
+ * lgamma(-integer) = +-inf
+ *
+ */
+
+static const double
+two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
+/* half= 5.00000000000000000000e-01, */ /* 0x3FE00000, 0x00000000 */
+/* one = 1.00000000000000000000e+00, */ /* 0x3FF00000, 0x00000000 */
+pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
+a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
+a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
+a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */
+a3 = 2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */
+a4 = 7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */
+a5 = 2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */
+a6 = 1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */
+a7 = 5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */
+a8 = 2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */
+a9 = 1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */
+a10 = 2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */
+a11 = 4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */
+tc = 1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */
+tf = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */
+/* tt = -(tail of tf) */
+tt = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */
+t0 = 4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */
+t1 = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */
+t2 = 6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */
+t3 = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */
+t4 = 1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */
+t5 = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */
+t6 = 6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */
+t7 = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */
+t8 = 2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */
+t9 = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */
+t10 = 8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */
+t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */
+t12 = 3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */
+t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */
+t14 = 3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */
+u0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
+u1 = 6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */
+u2 = 1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */
+u3 = 9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */
+u4 = 2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */
+u5 = 1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */
+v1 = 2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */
+v2 = 2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */
+v3 = 7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */
+v4 = 1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */
+v5 = 3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */
+s0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
+s1 = 2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */
+s2 = 3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */
+s3 = 1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */
+s4 = 2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */
+s5 = 1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */
+s6 = 3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */
+r1 = 1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */
+r2 = 7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */
+r3 = 1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */
+r4 = 1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */
+r5 = 7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */
+r6 = 7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */
+w0 = 4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */
+w1 = 8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */
+w2 = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */
+w3 = 7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */
+w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
+w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
+w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
+
+static double zero = 0.0;
+
+static double sin_pi(double x)
+{
+ double y,z,f;
+ int n,ix;
+ f = frexp(x, &ix);
+
+ if(fabs(x)<0.25) return -sin(pi*x);
+ y = -x; /* x is assume negative */
+
+ /*
+ * argument reduction, make sure inexact flag not raised if input
+ * is an integer
+ */
+ z = floor(y);
+ if(z!=y) { /* inexact anyway */
+ y *= 0.5;
+ y = 2.0*(y - floor(y)); /* y = |x| mod 2.0 */
+ n = (int) (y*4.0);
+ } else {
+ if(ix>=52) {
+ y = zero; n = 0; /* y must be even */
+ } else {
+ if(ix<52) z = y+two52; /* exact */
+ f = frexp(z, &n);
+ n &= 1;
+ y = n;
+ n<<= 2;
+ }
+ }
+ switch (n) {
+ case 0: y = sin(pi*y); break;
+ case 1:
+ case 2: y = cos(pi*(0.5-y)); break;
+ case 3:
+ case 4: y = sin(pi*(one-y)); break;
+ case 5:
+ case 6: y = -cos(pi*(y-1.5)); break;
+ default: y = sin(pi*(y-2.0)); break;
+ }
+ z = cos(pi*(z+1.0));
+ return -y*z;
+}
+
+double
+py_lgamma(double x)
+{
+ double t,y,z,nadj,p,p1,p2,p3,q,r,w,f;
+ int i,hx,ix,signgamp;
+
+ nadj = 0;
+
+ f = frexp(x, &ix);
+
+ /* purge off +-inf, NaN, +-0, and negative arguments */
+ signgamp = 1;
+ hx = ix&0x7fffffff;
+ if(hx>=0x7ff00000) return x*x;
+ if((f==0.0)&&(ix==0)) return one/zero;
+ if(ix<-71) { /* |x|<2**-70, return -log(|x|) */
+ if(x<0) {
+ signgamp = -1;
+ return -log(-x);
+ } else return -log(x);
+ }
+ if(f<0) {
+ if(ix>=52) /* |x|>=2**52, must be -integer */
+ return one/zero;
+ t = sin_pi(x);
+ if(t==zero) return one/zero; /* -integer */
+ nadj = log(pi/fabs(t*x));
+ if(t=0.7316) {y = one-x; i= 0;}
+ else if(x>=0.23164) {y= x-(tc-one); i=1;}
+ else {y = x; i=2;}
+ } else {
+ r = zero;
+ if(x>=1.7316) {y=2.0-x;i=0;} /* [1.7316,2] */
+ else if(x>=1.23164) {y=x-tc;i=1;} /* [1.23,1.73] */
+ else {y=x-one;i=2;}
+ }
+ switch(i) {
+ case 0:
+ z = y*y;
+ p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
+ p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
+ p = y*p1+p2;
+ r += (p-0.5*y); break;
+ case 1:
+ z = y*y;
+ w = z*y;
+ p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */
+ p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
+ p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
+ p = z*p1-(tt-w*(p2+y*p3));
+ r += (tf + p); break;
+ case 2:
+ p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
+ p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
+ r += (-0.5*y + p1/p2);
+ }
+ /* x < 8.0 */
+ } else if(ix<4) {
+ i = (int)x;
+ t = zero;
+ y = x-(double)i;
+ p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
+ q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
+ r = half*y+p/q;
+ z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
+ switch(i) {
+ case 7: z *= (y+6.0); /* FALLTHRU */
+ case 6: z *= (y+5.0); /* FALLTHRU */
+ case 5: z *= (y+4.0); /* FALLTHRU */
+ case 4: z *= (y+3.0); /* FALLTHRU */
+ case 3: z *= (y+2.0); /* FALLTHRU */
+ r += log(z); break;
+ }
+ /* 8.0 <= x < 2**58 */
+ } else if(ix < 58) {
+ t = log(x);
+ z = one/x;
+ y = z*z;
+ w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
+ r = (x-half)*(t-one)+w;
+ /* 2**58 <= x <= inf */
+ } else
+ r = x*(log(x)-one);
+ if(f<0) r = nadj - r;
+ return signgamp*r;
+}
+
+double
+py_tgamma(double x)
+{
+ double s;
+ s = 1.0;
+ if (x<0) {
+ s = cos(pi*floor(x));
+ }
+ return s*exp(py_lgamma(x));
+}
+
+#endif /* HAVE_SPECIAL_FUNCTIONS */
Index: Modules/math_cextra.h
===================================================================
--- Modules/math_cextra.h (revision 0)
+++ Modules/math_cextra.h (revision 0)
@@ -0,0 +1,4 @@
+double py_erf(double);
+double py_erfc(double);
+double py_lgamma(double);
+double py_tgamma(double);
Index: Modules/mathmodule.c
===================================================================
--- Modules/mathmodule.c (revision 74666)
+++ Modules/mathmodule.c (working copy)
@@ -54,6 +54,7 @@
#include "Python.h"
#include "longintrepr.h" /* just for SHIFT */
+#include "math_cextra.h"
#ifdef _OSF_SOURCE
/* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
@@ -247,6 +248,25 @@
return PyFloat_FromDouble(r);
}
+/* variant of math_1, to be used when the function being wrapped is known
+ to set errno properly. */
+
+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 +333,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); \
@@ -343,6 +369,10 @@
"cos(x)\n\nReturn the cosine of x (measured in radians).")
FUNC1(cosh, cosh, 1,
"cosh(x)\n\nReturn the hyperbolic cosine of x.")
+FUNC1(erf, py_erf, 0,
+ "erf(x)\n\nError function at x.")
+FUNC1(erfc, py_erfc, 0,
+ "erfc(x)\n\nComplementary error function at x.")
FUNC1(exp, exp, 1,
"exp(x)\n\nReturn e raised to the power of x.")
FUNC1(fabs, fabs, 0,
@@ -350,6 +380,10 @@
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.")
+FUNC1A(lgamma, py_lgamma,
+ "lgamma(x)\n\nNatural log of the absolute value of the gamma function.")
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.")
@@ -1070,6 +1104,8 @@
{"cos", math_cos, METH_O, math_cos_doc},
{"cosh", math_cosh, METH_O, math_cosh_doc},
{"degrees", math_degrees, METH_O, math_degrees_doc},
+ {"erf", math_erf, METH_O, math_erf_doc},
+ {"erfc", math_erfc, METH_O, math_erfc_doc},
{"exp", math_exp, METH_O, math_exp_doc},
{"fabs", math_fabs, METH_O, math_fabs_doc},
{"factorial", math_factorial, METH_O, math_factorial_doc},
@@ -1077,10 +1113,12 @@
{"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},
{"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},
Index: Modules/Setup.dist
===================================================================
--- Modules/Setup.dist (revision 74666)
+++ 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_cextra.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