diff options
Diffstat (limited to 'Modules/mathmodule.c')
| -rw-r--r-- | Modules/mathmodule.c | 986 |
1 files changed, 862 insertions, 124 deletions
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 05bad5de6e..29c32a30d8 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -53,54 +53,471 @@ raised for division by zero and mod by zero. */ #include "Python.h" -#include "longintrepr.h" /* just for SHIFT */ +#include "_math.h" #ifdef _OSF_SOURCE /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */ extern double copysign(double, double); #endif -/* Call is_error when errno != 0, and where x is the result libm - * returned. is_error will usually set up an exception and return - * true (1), but may return false (0) without setting up an exception. - */ -static int -is_error(double x) +/* + sin(pi*x), giving accurate results for all finite x (especially x + integral or close to an integer). This is here for use in the + reflection formula for the gamma function. It conforms to IEEE + 754-2008 for finite arguments, but not for infinities or nans. +*/ + +static const double pi = 3.141592653589793238462643383279502884197; +static const double sqrtpi = 1.772453850905516027298167483341145182798; +static const double logpi = 1.144729885849400174143427351353058711647; + +static double +sinpi(double x) { - int result = 1; /* presumption of guilt */ - assert(errno); /* non-zero errno is a precondition for calling */ - if (errno == EDOM) - PyErr_SetString(PyExc_ValueError, "math domain error"); + double y, r; + int n; + /* this function should only ever be called for finite arguments */ + assert(Py_IS_FINITE(x)); + y = fmod(fabs(x), 2.0); + n = (int)round(2.0*y); + assert(0 <= n && n <= 4); + switch (n) { + case 0: + r = sin(pi*y); + break; + case 1: + r = cos(pi*(y-0.5)); + break; + case 2: + /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give + -0.0 instead of 0.0 when y == 1.0. */ + r = sin(pi*(1.0-y)); + break; + case 3: + r = -cos(pi*(y-1.5)); + break; + case 4: + r = sin(pi*(y-2.0)); + break; + default: + assert(0); /* should never get here */ + r = -1.23e200; /* silence gcc warning */ + } + return copysign(1.0, x)*r; +} - else if (errno == ERANGE) { - /* ANSI C generally requires libm functions to set ERANGE - * on overflow, but also generally *allows* them to set - * ERANGE on underflow too. There's no consistency about - * the latter across platforms. - * Alas, C99 never requires that errno be set. - * Here we suppress the underflow errors (libm functions - * should return a zero on underflow, and +- HUGE_VAL on - * overflow, so testing the result for zero suffices to - * distinguish the cases). - * - * On some platforms (Ubuntu/ia64) it seems that errno can be - * set to ERANGE for subnormal results that do *not* underflow - * to zero. So to be safe, we'll ignore ERANGE whenever the - * function result is less than one in absolute value. - */ - if (fabs(x) < 1.0) - result = 0; +/* 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; +} + +/* + lgamma: natural log of the absolute value of the Gamma function. + For large arguments, Lanczos' formula works extremely well here. +*/ + +static double +m_lgamma(double x) +{ + double r, absx; + + /* special cases */ + if (!Py_IS_FINITE(x)) { + if (Py_IS_NAN(x)) + return x; /* lgamma(nan) = nan */ else - PyErr_SetString(PyExc_OverflowError, - "math range error"); + return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */ } - else - /* Unexpected math error */ - PyErr_SetFromErrno(PyExc_ValueError); + + /* integer arguments */ + if (x == floor(x) && x <= 2.0) { + if (x <= 0.0) { + errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */ + return Py_HUGE_VAL; /* integers n <= 0 */ + } + else { + return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */ + } + } + + absx = fabs(x); + /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */ + if (absx < 1e-20) + return -log(absx); + + /* Lanczos' formula. We could save a fraction of a ulp in accuracy by + having a second set of numerator coefficients for lanczos_sum that + absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g + subtraction below; it's probably not worth it. */ + r = log(lanczos_sum(absx)) - lanczos_g; + r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1); + if (x < 0.0) + /* Use reflection formula to get value for negative x. */ + r = logpi - log(fabs(sinpi(absx))) - log(absx) - r; + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; +} + +/* + Implementations of the error function erf(x) and the complementary error + function erfc(x). + + Method: following 'Numerical Recipes' by Flannery, Press et. al. (2nd ed., + Cambridge University Press), we use a series approximation for erf for + small x, and a continued fraction approximation for erfc(x) for larger x; + combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x), + this gives us erf(x) and erfc(x) for all x. + + The series expansion used is: + + erf(x) = x*exp(-x*x)/sqrt(pi) * [ + 2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...] + + The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k). + This series converges well for smallish x, but slowly for larger x. + + The continued fraction expansion used is: + + erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - ) + 3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...] + + after the first term, the general term has the form: + + k*(k-0.5)/(2*k+0.5 + x**2 - ...). + + This expansion converges fast for larger x, but convergence becomes + infinitely slow as x approaches 0.0. The (somewhat naive) continued + fraction evaluation algorithm used below also risks overflow for large x; + but for large x, erfc(x) == 0.0 to within machine precision. (For + example, erfc(30.0) is approximately 2.56e-393). + + Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and + continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) < + ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the + numbers of terms to use for the relevant expansions. */ + +#define ERF_SERIES_CUTOFF 1.5 +#define ERF_SERIES_TERMS 25 +#define ERFC_CONTFRAC_CUTOFF 30.0 +#define ERFC_CONTFRAC_TERMS 50 + +/* + Error function, via power series. + + Given a finite float x, return an approximation to erf(x). + Converges reasonably fast for small x. +*/ + +static double +m_erf_series(double x) +{ + double x2, acc, fk, result; + int i, saved_errno; + + x2 = x * x; + acc = 0.0; + fk = (double)ERF_SERIES_TERMS + 0.5; + for (i = 0; i < ERF_SERIES_TERMS; i++) { + acc = 2.0 + x2 * acc / fk; + fk -= 1.0; + } + /* Make sure the exp call doesn't affect errno; + see m_erfc_contfrac for more. */ + saved_errno = errno; + result = acc * x * exp(-x2) / sqrtpi; + errno = saved_errno; return result; } /* + Complementary error function, via continued fraction expansion. + + Given a positive float x, return an approximation to erfc(x). Converges + reasonably fast for x large (say, x > 2.0), and should be safe from + overflow if x and nterms are not too large. On an IEEE 754 machine, with x + <= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller + than the smallest representable nonzero float. */ + +static double +m_erfc_contfrac(double x) +{ + double x2, a, da, p, p_last, q, q_last, b, result; + int i, saved_errno; + + if (x >= ERFC_CONTFRAC_CUTOFF) + return 0.0; + + x2 = x*x; + a = 0.0; + da = 0.5; + p = 1.0; p_last = 0.0; + q = da + x2; q_last = 1.0; + for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) { + double temp; + a += da; + da += 2.0; + b = da + x2; + temp = p; p = b*p - a*p_last; p_last = temp; + temp = q; q = b*q - a*q_last; q_last = temp; + } + /* Issue #8986: On some platforms, exp sets errno on underflow to zero; + save the current errno value so that we can restore it later. */ + saved_errno = errno; + result = p / q * x * exp(-x2) / sqrtpi; + errno = saved_errno; + return result; +} + +/* Error function erf(x), for general x */ + +static double +m_erf(double x) +{ + double absx, cf; + + if (Py_IS_NAN(x)) + return x; + absx = fabs(x); + if (absx < ERF_SERIES_CUTOFF) + return m_erf_series(x); + else { + cf = m_erfc_contfrac(absx); + return x > 0.0 ? 1.0 - cf : cf - 1.0; + } +} + +/* Complementary error function erfc(x), for general x. */ + +static double +m_erfc(double x) +{ + double absx, cf; + + if (Py_IS_NAN(x)) + return x; + absx = fabs(x); + if (absx < ERF_SERIES_CUTOFF) + return 1.0 - m_erf_series(x); + else { + cf = m_erfc_contfrac(absx); + return x > 0.0 ? cf : 2.0 - cf; + } +} + +/* wrapper for atan2 that deals directly with special cases before delegating to the platform libm for the remaining cases. This is necessary to get consistent behaviour across platforms. @@ -188,6 +605,46 @@ m_log10(double x) } +/* Call is_error when errno != 0, and where x is the result libm + * returned. is_error will usually set up an exception and return + * true (1), but may return false (0) without setting up an exception. + */ +static int +is_error(double x) +{ + int result = 1; /* presumption of guilt */ + assert(errno); /* non-zero errno is a precondition for calling */ + if (errno == EDOM) + PyErr_SetString(PyExc_ValueError, "math domain error"); + + else if (errno == ERANGE) { + /* ANSI C generally requires libm functions to set ERANGE + * on overflow, but also generally *allows* them to set + * ERANGE on underflow too. There's no consistency about + * the latter across platforms. + * Alas, C99 never requires that errno be set. + * Here we suppress the underflow errors (libm functions + * should return a zero on underflow, and +- HUGE_VAL on + * overflow, so testing the result for zero suffices to + * distinguish the cases). + * + * On some platforms (Ubuntu/ia64) it seems that errno can be + * set to ERANGE for subnormal results that do *not* underflow + * to zero. So to be safe, we'll ignore ERANGE whenever the + * function result is less than one in absolute value. + */ + if (fabs(x) < 1.0) + result = 0; + else + PyErr_SetString(PyExc_OverflowError, + "math range error"); + } + else + /* Unexpected math error */ + PyErr_SetFromErrno(PyExc_ValueError); + return result; +} + /* math_1 is used to wrap a libm function f that takes a double arguments and returns a double. @@ -252,6 +709,26 @@ math_1_to_whatever(PyObject *arg, double (*func) (double), return (*from_double_func)(r); } +/* variant of math_1, to be used when the function being wrapped is known to + set errno properly (that is, errno = EDOM for invalid or divide-by-zero, + errno = ERANGE for overflow). */ + +static PyObject * +math_1a(PyObject *arg, double (*func) (double)) +{ + double x, r; + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + errno = 0; + PyFPE_START_PROTECT("in math_1a", return 0); + r = (*func)(x); + PyFPE_END_PROTECT(r); + if (errno && is_error(r)) + return NULL; + return PyFloat_FromDouble(r); +} + /* math_2 is used to wrap a libm function f that takes two double arguments and returns a double. @@ -330,6 +807,12 @@ math_2(PyObject *args, double (*func) (double, double), char *funcname) }\ PyDoc_STRVAR(math_##funcname##_doc, docstring); +#define FUNC1A(funcname, func, docstring) \ + static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ + return math_1a(args, func); \ + }\ + PyDoc_STRVAR(math_##funcname##_doc, docstring); + #define FUNC2(funcname, func, docstring) \ static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ return math_2(args, func, #funcname); \ @@ -338,35 +821,33 @@ math_2(PyObject *args, double (*func) (double, double), char *funcname) FUNC1(acos, acos, 0, "acos(x)\n\nReturn the arc cosine (measured in radians) of x.") -FUNC1(acosh, acosh, 0, +FUNC1(acosh, m_acosh, 0, "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.") FUNC1(asin, asin, 0, "asin(x)\n\nReturn the arc sine (measured in radians) of x.") -FUNC1(asinh, asinh, 0, +FUNC1(asinh, m_asinh, 0, "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.") FUNC1(atan, atan, 0, "atan(x)\n\nReturn the arc tangent (measured in radians) of x.") FUNC2(atan2, m_atan2, "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n" "Unlike atan(y/x), the signs of both x and y are considered.") -FUNC1(atanh, atanh, 0, +FUNC1(atanh, m_atanh, 0, "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.") static PyObject * math_ceil(PyObject *self, PyObject *number) { static PyObject *ceil_str = NULL; - PyObject *method; + PyObject *method, *result; - if (ceil_str == NULL) { - ceil_str = PyUnicode_InternFromString("__ceil__"); - if (ceil_str == NULL) + method = _PyObject_LookupSpecial(number, "__ceil__", &ceil_str); + if (method == NULL) { + if (PyErr_Occurred()) return NULL; - } - - method = _PyType_Lookup(Py_TYPE(number), ceil_str); - if (method == NULL) return math_1_to_int(number, ceil, 0); - else - return PyObject_CallFunction(method, "O", number); + } + result = PyObject_CallFunctionObjArgs(method, NULL); + Py_DECREF(method); + return result; } PyDoc_STRVAR(math_ceil_doc, @@ -379,33 +860,43 @@ FUNC1(cos, cos, 0, "cos(x)\n\nReturn the cosine of x (measured in radians).") FUNC1(cosh, cosh, 1, "cosh(x)\n\nReturn the hyperbolic cosine of x.") +FUNC1A(erf, m_erf, + "erf(x)\n\nError function at x.") +FUNC1A(erfc, m_erfc, + "erfc(x)\n\nComplementary error function at x.") FUNC1(exp, exp, 1, "exp(x)\n\nReturn e raised to the power of x.") +FUNC1(expm1, m_expm1, 1, + "expm1(x)\n\nReturn exp(x)-1.\n" + "This function avoids the loss of precision involved in the direct " + "evaluation of exp(x)-1 for small x.") FUNC1(fabs, fabs, 0, "fabs(x)\n\nReturn the absolute value of the float x.") static PyObject * math_floor(PyObject *self, PyObject *number) { static PyObject *floor_str = NULL; - PyObject *method; + PyObject *method, *result; - if (floor_str == NULL) { - floor_str = PyUnicode_InternFromString("__floor__"); - if (floor_str == NULL) + method = _PyObject_LookupSpecial(number, "__floor__", &floor_str); + if (method == NULL) { + if (PyErr_Occurred()) return NULL; - } - - method = _PyType_Lookup(Py_TYPE(number), floor_str); - if (method == NULL) return math_1_to_int(number, floor, 0); - else - return PyObject_CallFunction(method, "O", number); + } + result = PyObject_CallFunctionObjArgs(method, NULL); + Py_DECREF(method); + return result; } PyDoc_STRVAR(math_floor_doc, "floor(x)\n\nReturn the floor of x as an int.\n" "This is the largest integral value <= x."); -FUNC1(log1p, log1p, 1, +FUNC1A(gamma, m_tgamma, + "gamma(x)\n\nGamma function at x.") +FUNC1A(lgamma, m_lgamma, + "lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.") +FUNC1(log1p, m_log1p, 0, "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.") FUNC1(sin, sin, 0, @@ -640,49 +1131,276 @@ PyDoc_STRVAR(math_fsum_doc, Return an accurate floating point sum of values in the iterable.\n\ Assumes IEEE-754 floating point arithmetic."); +/* Return the smallest integer k such that n < 2**k, or 0 if n == 0. + * Equivalent to floor(lg(x))+1. Also equivalent to: bitwidth_of_type - + * count_leading_zero_bits(x) + */ + +/* XXX: This routine does more or less the same thing as + * bits_in_digit() in Objects/longobject.c. Someday it would be nice to + * consolidate them. On BSD, there's a library function called fls() + * that we could use, and GCC provides __builtin_clz(). + */ + +static unsigned long +bit_length(unsigned long n) +{ + unsigned long len = 0; + while (n != 0) { + ++len; + n >>= 1; + } + return len; +} + +static unsigned long +count_set_bits(unsigned long n) +{ + unsigned long count = 0; + while (n != 0) { + ++count; + n &= n - 1; /* clear least significant bit */ + } + return count; +} + +/* Divide-and-conquer factorial algorithm + * + * Based on the formula and psuedo-code provided at: + * http://www.luschny.de/math/factorial/binarysplitfact.html + * + * Faster algorithms exist, but they're more complicated and depend on + * a fast prime factorization algorithm. + * + * Notes on the algorithm + * ---------------------- + * + * factorial(n) is written in the form 2**k * m, with m odd. k and m are + * computed separately, and then combined using a left shift. + * + * The function factorial_odd_part computes the odd part m (i.e., the greatest + * odd divisor) of factorial(n), using the formula: + * + * factorial_odd_part(n) = + * + * product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j + * + * Example: factorial_odd_part(20) = + * + * (1) * + * (1) * + * (1 * 3 * 5) * + * (1 * 3 * 5 * 7 * 9) + * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19) + * + * Here i goes from large to small: the first term corresponds to i=4 (any + * larger i gives an empty product), and the last term corresponds to i=0. + * Each term can be computed from the last by multiplying by the extra odd + * numbers required: e.g., to get from the penultimate term to the last one, + * we multiply by (11 * 13 * 15 * 17 * 19). + * + * To see a hint of why this formula works, here are the same numbers as above + * but with the even parts (i.e., the appropriate powers of 2) included. For + * each subterm in the product for i, we multiply that subterm by 2**i: + * + * factorial(20) = + * + * (16) * + * (8) * + * (4 * 12 * 20) * + * (2 * 6 * 10 * 14 * 18) * + * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19) + * + * The factorial_partial_product function computes the product of all odd j in + * range(start, stop) for given start and stop. It's used to compute the + * partial products like (11 * 13 * 15 * 17 * 19) in the example above. It + * operates recursively, repeatedly splitting the range into two roughly equal + * pieces until the subranges are small enough to be computed using only C + * integer arithmetic. + * + * The two-valuation k (i.e., the exponent of the largest power of 2 dividing + * the factorial) is computed independently in the main math_factorial + * function. By standard results, its value is: + * + * two_valuation = n//2 + n//4 + n//8 + .... + * + * It can be shown (e.g., by complete induction on n) that two_valuation is + * equal to n - count_set_bits(n), where count_set_bits(n) gives the number of + * '1'-bits in the binary expansion of n. + */ + +/* factorial_partial_product: Compute product(range(start, stop, 2)) using + * divide and conquer. Assumes start and stop are odd and stop > start. + * max_bits must be >= bit_length(stop - 2). */ + +static PyObject * +factorial_partial_product(unsigned long start, unsigned long stop, + unsigned long max_bits) +{ + unsigned long midpoint, num_operands; + PyObject *left = NULL, *right = NULL, *result = NULL; + + /* If the return value will fit an unsigned long, then we can + * multiply in a tight, fast loop where each multiply is O(1). + * Compute an upper bound on the number of bits required to store + * the answer. + * + * Storing some integer z requires floor(lg(z))+1 bits, which is + * conveniently the value returned by bit_length(z). The + * product x*y will require at most + * bit_length(x) + bit_length(y) bits to store, based + * on the idea that lg product = lg x + lg y. + * + * We know that stop - 2 is the largest number to be multiplied. From + * there, we have: bit_length(answer) <= num_operands * + * bit_length(stop - 2) + */ + + num_operands = (stop - start) / 2; + /* The "num_operands <= 8 * SIZEOF_LONG" check guards against the + * unlikely case of an overflow in num_operands * max_bits. */ + if (num_operands <= 8 * SIZEOF_LONG && + num_operands * max_bits <= 8 * SIZEOF_LONG) { + unsigned long j, total; + for (total = start, j = start + 2; j < stop; j += 2) + total *= j; + return PyLong_FromUnsignedLong(total); + } + + /* find midpoint of range(start, stop), rounded up to next odd number. */ + midpoint = (start + num_operands) | 1; + left = factorial_partial_product(start, midpoint, + bit_length(midpoint - 2)); + if (left == NULL) + goto error; + right = factorial_partial_product(midpoint, stop, max_bits); + if (right == NULL) + goto error; + result = PyNumber_Multiply(left, right); + + error: + Py_XDECREF(left); + Py_XDECREF(right); + return result; +} + +/* factorial_odd_part: compute the odd part of factorial(n). */ + +static PyObject * +factorial_odd_part(unsigned long n) +{ + long i; + unsigned long v, lower, upper; + PyObject *partial, *tmp, *inner, *outer; + + inner = PyLong_FromLong(1); + if (inner == NULL) + return NULL; + outer = inner; + Py_INCREF(outer); + + upper = 3; + for (i = bit_length(n) - 2; i >= 0; i--) { + v = n >> i; + if (v <= 2) + continue; + lower = upper; + /* (v + 1) | 1 = least odd integer strictly larger than n / 2**i */ + upper = (v + 1) | 1; + /* Here inner is the product of all odd integers j in the range (0, + n/2**(i+1)]. The factorial_partial_product call below gives the + product of all odd integers j in the range (n/2**(i+1), n/2**i]. */ + partial = factorial_partial_product(lower, upper, bit_length(upper-2)); + /* inner *= partial */ + if (partial == NULL) + goto error; + tmp = PyNumber_Multiply(inner, partial); + Py_DECREF(partial); + if (tmp == NULL) + goto error; + Py_DECREF(inner); + inner = tmp; + /* Now inner is the product of all odd integers j in the range (0, + n/2**i], giving the inner product in the formula above. */ + + /* outer *= inner; */ + tmp = PyNumber_Multiply(outer, inner); + if (tmp == NULL) + goto error; + Py_DECREF(outer); + outer = tmp; + } + + goto done; + + error: + Py_DECREF(outer); + done: + Py_DECREF(inner); + return outer; +} + +/* Lookup table for small factorial values */ + +static const unsigned long SmallFactorials[] = { + 1, 1, 2, 6, 24, 120, 720, 5040, 40320, + 362880, 3628800, 39916800, 479001600, +#if SIZEOF_LONG >= 8 + 6227020800, 87178291200, 1307674368000, + 20922789888000, 355687428096000, 6402373705728000, + 121645100408832000, 2432902008176640000 +#endif +}; + static PyObject * math_factorial(PyObject *self, PyObject *arg) { - long i, x; - PyObject *result, *iobj, *newresult; + long x; + PyObject *result, *odd_part, *two_valuation; if (PyFloat_Check(arg)) { + PyObject *lx; double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg); - if (dx != floor(dx)) { + if (!(Py_IS_FINITE(dx) && dx == floor(dx))) { PyErr_SetString(PyExc_ValueError, - "factorial() only accepts integral values"); + "factorial() only accepts integral values"); return NULL; } + lx = PyLong_FromDouble(dx); + if (lx == NULL) + return NULL; + x = PyLong_AsLong(lx); + Py_DECREF(lx); } + else + x = PyLong_AsLong(arg); - x = PyLong_AsLong(arg); if (x == -1 && PyErr_Occurred()) return NULL; if (x < 0) { PyErr_SetString(PyExc_ValueError, - "factorial() not defined for negative values"); + "factorial() not defined for negative values"); return NULL; } - result = (PyObject *)PyLong_FromLong(1); - if (result == NULL) + /* use lookup table if x is small */ + if (x < (long)(sizeof(SmallFactorials)/sizeof(SmallFactorials[0]))) + return PyLong_FromUnsignedLong(SmallFactorials[x]); + + /* else express in the form odd_part * 2**two_valuation, and compute as + odd_part << two_valuation. */ + odd_part = factorial_odd_part(x); + if (odd_part == NULL) + return NULL; + two_valuation = PyLong_FromLong(x - count_set_bits(x)); + if (two_valuation == NULL) { + Py_DECREF(odd_part); return NULL; - for (i=1 ; i<=x ; i++) { - iobj = (PyObject *)PyLong_FromLong(i); - if (iobj == NULL) - goto error; - newresult = PyNumber_Multiply(result, iobj); - Py_DECREF(iobj); - if (newresult == NULL) - goto error; - Py_DECREF(result); - result = newresult; } + result = PyNumber_Lshift(odd_part, two_valuation); + Py_DECREF(two_valuation); + Py_DECREF(odd_part); return result; - -error: - Py_DECREF(result); - return NULL; } PyDoc_STRVAR(math_factorial_doc, @@ -694,27 +1412,24 @@ static PyObject * math_trunc(PyObject *self, PyObject *number) { static PyObject *trunc_str = NULL; - PyObject *trunc; + PyObject *trunc, *result; if (Py_TYPE(number)->tp_dict == NULL) { if (PyType_Ready(Py_TYPE(number)) < 0) return NULL; } - if (trunc_str == NULL) { - trunc_str = PyUnicode_InternFromString("__trunc__"); - if (trunc_str == NULL) - return NULL; - } - - trunc = _PyType_Lookup(Py_TYPE(number), trunc_str); + trunc = _PyObject_LookupSpecial(number, "__trunc__", &trunc_str); if (trunc == NULL) { - PyErr_Format(PyExc_TypeError, - "type %.100s doesn't define __trunc__ method", - Py_TYPE(number)->tp_name); + if (!PyErr_Occurred()) + PyErr_Format(PyExc_TypeError, + "type %.100s doesn't define __trunc__ method", + Py_TYPE(number)->tp_name); return NULL; } - return PyObject_CallFunctionObjArgs(trunc, number, NULL); + result = PyObject_CallFunctionObjArgs(trunc, NULL); + Py_DECREF(trunc); + return result; } PyDoc_STRVAR(math_trunc_doc, @@ -755,28 +1470,18 @@ math_ldexp(PyObject *self, PyObject *args) double x, r; PyObject *oexp; long exp; + int overflow; if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp)) return NULL; if (PyLong_Check(oexp)) { /* on overflow, replace exponent with either LONG_MAX or LONG_MIN, depending on the sign. */ - exp = PyLong_AsLong(oexp); - if (exp == -1 && PyErr_Occurred()) { - if (PyErr_ExceptionMatches(PyExc_OverflowError)) { - if (Py_SIZE(oexp) < 0) { - exp = LONG_MIN; - } - else { - exp = LONG_MAX; - } - PyErr_Clear(); - } - else { - /* propagate any unexpected exception */ - return NULL; - } - } + exp = PyLong_AsLongAndOverflow(oexp, &overflow); + if (exp == -1 && PyErr_Occurred()) + return NULL; + if (overflow) + exp = overflow < 0 ? LONG_MIN : LONG_MAX; } else { PyErr_SetString(PyExc_TypeError, @@ -845,31 +1550,45 @@ PyDoc_STRVAR(math_modf_doc, /* A decent logarithm is easy to compute even for huge longs, but libm can't do that by itself -- loghelper can. func is log or log10, and name is - "log" or "log10". Note that overflow isn't possible: a long can contain - no more than INT_MAX * SHIFT bits, so has value certainly less than - 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is + "log" or "log10". Note that overflow of the result isn't possible: a long + can contain no more than INT_MAX * SHIFT bits, so has value certainly less + than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is small enough to fit in an IEEE single. log and log10 are even smaller. -*/ + However, intermediate overflow is possible for a long if the number of bits + in that long is larger than PY_SSIZE_T_MAX. */ static PyObject* loghelper(PyObject* arg, double (*func)(double), char *funcname) { /* If it is long, do it ourselves. */ if (PyLong_Check(arg)) { - double x; - int e; - x = _PyLong_AsScaledDouble(arg, &e); - if (x <= 0.0) { + double x, result; + Py_ssize_t e; + + /* Negative or zero inputs give a ValueError. */ + if (Py_SIZE(arg) <= 0) { PyErr_SetString(PyExc_ValueError, "math domain error"); return NULL; } - /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~= - log(x) + log(2) * e * PyLong_SHIFT. - CAUTION: e*PyLong_SHIFT may overflow using int arithmetic, - so force use of double. */ - x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0); - return PyFloat_FromDouble(x); + + x = PyLong_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + if (!PyErr_ExceptionMatches(PyExc_OverflowError)) + return NULL; + /* Here the conversion to double overflowed, but it's possible + to compute the log anyway. Clear the exception and continue. */ + PyErr_Clear(); + x = _PyLong_Frexp((PyLongObject *)arg, &e); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */ + result = func(x) + func(2.0) * e; + } + else + /* Successfully converted x to a double. */ + result = func(x); + return PyFloat_FromDouble(result); } /* Else let libm handle it by itself. */ @@ -1109,6 +1828,19 @@ PyDoc_STRVAR(math_radians_doc, Convert angle x from degrees to radians."); static PyObject * +math_isfinite(PyObject *self, PyObject *arg) +{ + double x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + return PyBool_FromLong((long)Py_IS_FINITE(x)); +} + +PyDoc_STRVAR(math_isfinite_doc, +"isfinite(x) -> bool\n\n\ +Return True if x is neither an infinity nor a NaN, and False otherwise."); + +static PyObject * math_isnan(PyObject *self, PyObject *arg) { double x = PyFloat_AsDouble(arg); @@ -1119,7 +1851,7 @@ math_isnan(PyObject *self, PyObject *arg) PyDoc_STRVAR(math_isnan_doc, "isnan(x) -> bool\n\n\ -Check if float x is not a number (NaN)."); +Return True if x is a NaN (not a number), and False otherwise."); static PyObject * math_isinf(PyObject *self, PyObject *arg) @@ -1132,7 +1864,7 @@ math_isinf(PyObject *self, PyObject *arg) PyDoc_STRVAR(math_isinf_doc, "isinf(x) -> bool\n\n\ -Check if float x is infinite (positive or negative)."); +Return True if x is a positive or negative infinity, and False otherwise."); static PyMethodDef math_methods[] = { {"acos", math_acos, METH_O, math_acos_doc}, @@ -1147,17 +1879,23 @@ static PyMethodDef math_methods[] = { {"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}, + {"expm1", math_expm1, METH_O, math_expm1_doc}, {"fabs", math_fabs, METH_O, math_fabs_doc}, {"factorial", math_factorial, METH_O, math_factorial_doc}, {"floor", math_floor, METH_O, math_floor_doc}, {"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}, + {"isfinite", math_isfinite, METH_O, math_isfinite_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}, |