Changeset 388 for python/vendor/current/Modules/mathmodule.c

Timestamp:

Mar 19, 2014, 11:11:30 AM (11 years ago)

Author:

dmik

Message:

python: Update vendor to 2.7.6.

File:

• python/vendor/current/Modules/mathmodule.c (modified) (25 diffs)

python/vendor/current/Modules/mathmodule.c

@@ -54 +54 @@
 
 #include "Python.h"
-#include "longintrepr.h" /* just for SHIFT */
+#include "_math.h"
 
 #ifdef _OSF_SOURCE
@@ -61 +61 @@
 #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)
-{
-        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;
+/*
+   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 double
+sinpi(double x)
+{
+    double y, r;
+    int n;
+    /* this function should only ever be called for finite arguments */
+    assert(Py_IS_FINITE(x));
+    y = fmod(fabs(x), 2.0);
+    n = (int)round(2.0*y);
+    assert(0 <= n && n <= 4);
+    switch (n) {
+    case 0:
+        r = sin(pi*y);
+        break;
+    case 1:
+        r = cos(pi*(y-0.5));
+        break;
+    case 2:
+        /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give
+           -0.0 instead of 0.0 when y == 1.0. */
+        r = sin(pi*(1.0-y));
+        break;
+    case 3:
+        r = -cos(pi*(y-1.5));
+        break;
+    case 4:
+        r = sin(pi*(y-2.0));
+        break;
+    default:
+        assert(0);  /* should never get here */
+        r = -1.23e200; /* silence gcc warning */
+    }
+    return copysign(1.0, x)*r;
+}
+
+/* Implementation of the real gamma function.  In extensive but non-exhaustive
+   random tests, this function proved accurate to within <= 10 ulps across the
+   entire float domain.  Note that accuracy may depend on the quality of the
+   system math functions, the pow function in particular.  Special cases
+   follow C99 annex F.  The parameters and method are tailored 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
+            return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */
+    }
+
+    /* 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 */
+    if (x > 0.0) {
+        /* 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(x)) - lanczos_g +
+            (x-0.5)*(log(x+lanczos_g-0.5)-1);
+    }
+    else {
+        r = log(pi) - log(fabs(sinpi(absx))) - log(absx) -
+            (log(lanczos_sum(absx)) - lanczos_g +
+             (absx-0.5)*(log(absx+lanczos_g-0.5)-1));
+    }
+    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;
+    }
 }
 
@@ -112 +533 @@
 m_atan2(double y, double x)
 {
-        if (Py_IS_NAN(x) || Py_IS_NAN(y))
-                return Py_NAN;
-        if (Py_IS_INFINITY(y)) {
-                if (Py_IS_INFINITY(x)) {
-                        if (copysign(1., x) == 1.)
-                                /* atan2(+-inf, +inf) == +-pi/4 */
-                                return copysign(0.25*Py_MATH_PI, y);
-                        else
-                                /* atan2(+-inf, -inf) == +-pi*3/4 */
-                                return copysign(0.75*Py_MATH_PI, y);
-                }
-                /* atan2(+-inf, x) == +-pi/2 for finite x */
-                return copysign(0.5*Py_MATH_PI, y);
-        }
-        if (Py_IS_INFINITY(x) || y == 0.) {
-                if (copysign(1., x) == 1.)
-                        /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
-                        return copysign(0., y);
-                else
-                        /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
-                        return copysign(Py_MATH_PI, y);
-        }
-        return atan2(y, x);
+    if (Py_IS_NAN(x) || Py_IS_NAN(y))
+        return Py_NAN;
+    if (Py_IS_INFINITY(y)) {
+        if (Py_IS_INFINITY(x)) {
+            if (copysign(1., x) == 1.)
+                /* atan2(+-inf, +inf) == +-pi/4 */
+                return copysign(0.25*Py_MATH_PI, y);
+            else
+                /* atan2(+-inf, -inf) == +-pi*3/4 */
+                return copysign(0.75*Py_MATH_PI, y);
+        }
+        /* atan2(+-inf, x) == +-pi/2 for finite x */
+        return copysign(0.5*Py_MATH_PI, y);
+    }
+    if (Py_IS_INFINITY(x) || y == 0.) {
+        if (copysign(1., x) == 1.)
+            /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
+            return copysign(0., y);
+        else
+            /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
+            return copysign(Py_MATH_PI, y);
+    }
+    return atan2(y, x);
 }
 
@@ -147 +568 @@
 m_log(double x)
 {
-        if (Py_IS_FINITE(x)) {
-                if (x > 0.0)
-                        return log(x);
-                errno = EDOM;
-                if (x == 0.0)
-                        return -Py_HUGE_VAL; /* log(0) = -inf */
-                else
-                        return Py_NAN; /* log(-ve) = nan */
-        }
-        else if (Py_IS_NAN(x))
-                return x; /* log(nan) = nan */
-        else if (x > 0.0)
-                return x; /* log(inf) = inf */
-        else {
-                errno = EDOM;
-                return Py_NAN; /* log(-inf) = nan */
-        }
+    if (Py_IS_FINITE(x)) {
+        if (x > 0.0)
+            return log(x);
+        errno = EDOM;
+        if (x == 0.0)
+            return -Py_HUGE_VAL; /* log(0) = -inf */
+        else
+            return Py_NAN; /* log(-ve) = nan */
+    }
+    else if (Py_IS_NAN(x))
+        return x; /* log(nan) = nan */
+    else if (x > 0.0)
+        return x; /* log(inf) = inf */
+    else {
+        errno = EDOM;
+        return Py_NAN; /* log(-inf) = nan */
+    }
 }
 
@@ -169 +590 @@
 m_log10(double x)
 {
-        if (Py_IS_FINITE(x)) {
-                if (x > 0.0)
-                        return log10(x);
-                errno = EDOM;
-                if (x == 0.0)
-                        return -Py_HUGE_VAL; /* log10(0) = -inf */
-                else
-                        return Py_NAN; /* log10(-ve) = nan */
-        }
-        else if (Py_IS_NAN(x))
-                return x; /* log10(nan) = nan */
-        else if (x > 0.0)
-                return x; /* log10(inf) = inf */
-        else {
-                errno = EDOM;
-                return Py_NAN; /* log10(-inf) = nan */
-        }
-}
-
+    if (Py_IS_FINITE(x)) {
+        if (x > 0.0)
+            return log10(x);
+        errno = EDOM;
+        if (x == 0.0)
+            return -Py_HUGE_VAL; /* log10(0) = -inf */
+        else
+            return Py_NAN; /* log10(-ve) = nan */
+    }
+    else if (Py_IS_NAN(x))
+        return x; /* log10(nan) = nan */
+    else if (x > 0.0)
+        return x; /* log10(inf) = inf */
+    else {
+        errno = EDOM;
+        return Py_NAN; /* log10(-inf) = nan */
+    }
+}
+
+
+/* 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;
+}
 
 /*
@@ -222 +683 @@
 math_1(PyObject *arg, double (*func) (double), int can_overflow)
 {
-        double x, r;
-        x = PyFloat_AsDouble(arg);
-        if (x == -1.0 && PyErr_Occurred())
-                return NULL;
-        errno = 0;
-        PyFPE_START_PROTECT("in math_1", return 0);
-        r = (*func)(x);
-        PyFPE_END_PROTECT(r);
-        if (Py_IS_NAN(r)) {
-                if (!Py_IS_NAN(x))
-                        errno = EDOM;
-                else
-                        errno = 0;
-        }
-        else if (Py_IS_INFINITY(r)) {
-                if (Py_IS_FINITE(x))
-                        errno = can_overflow ? ERANGE : EDOM;
-                else
-                        errno = 0;
-        }
-        if (errno && is_error(r))
-                return NULL;
-        else
-                return PyFloat_FromDouble(r);
+    double x, r;
+    x = PyFloat_AsDouble(arg);
+    if (x == -1.0 && PyErr_Occurred())
+        return NULL;
+    errno = 0;
+    PyFPE_START_PROTECT("in math_1", return 0);
+    r = (*func)(x);
+    PyFPE_END_PROTECT(r);
+    if (Py_IS_NAN(r)) {
+        if (!Py_IS_NAN(x))
+            errno = EDOM;
+        else
+            errno = 0;
+    }
+    else if (Py_IS_INFINITY(r)) {
+        if (Py_IS_FINITE(x))
+            errno = can_overflow ? ERANGE : EDOM;
+        else
+            errno = 0;
+    }
+    if (errno && is_error(r))
+        return NULL;
+    else
+        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);
 }
 
@@ -278 +759 @@
 math_2(PyObject *args, double (*func) (double, double), char *funcname)
 {
-        PyObject *ox, *oy;
-        double x, y, r;
-        if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
-                return NULL;
-        x = PyFloat_AsDouble(ox);
-        y = PyFloat_AsDouble(oy);
-        if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
-                return NULL;
-        errno = 0;
-        PyFPE_START_PROTECT("in math_2", return 0);
-        r = (*func)(x, y);
-        PyFPE_END_PROTECT(r);
-        if (Py_IS_NAN(r)) {
-                if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
-                        errno = EDOM;
-                else
-                        errno = 0;
-        }
-        else if (Py_IS_INFINITY(r)) {
-                if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
-                        errno = ERANGE;
-                else
-                        errno = 0;
-        }
-        if (errno && is_error(r))
-                return NULL;
-        else
-                return PyFloat_FromDouble(r);
-}
-
-#define FUNC1(funcname, func, can_overflow, docstring)                  \
-        static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
-                return math_1(args, func, can_overflow);                    \
-        }\
-        PyDoc_STRVAR(math_##funcname##_doc, docstring);
+    PyObject *ox, *oy;
+    double x, y, r;
+    if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
+        return NULL;
+    x = PyFloat_AsDouble(ox);
+    y = PyFloat_AsDouble(oy);
+    if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
+        return NULL;
+    errno = 0;
+    PyFPE_START_PROTECT("in math_2", return 0);
+    r = (*func)(x, y);
+    PyFPE_END_PROTECT(r);
+    if (Py_IS_NAN(r)) {
+        if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
+            errno = EDOM;
+        else
+            errno = 0;
+    }
+    else if (Py_IS_INFINITY(r)) {
+        if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
+            errno = ERANGE;
+        else
+            errno = 0;
+    }
+    if (errno && is_error(r))
+        return NULL;
+    else
+        return PyFloat_FromDouble(r);
+}
+
+#define FUNC1(funcname, func, can_overflow, docstring)                  \
+    static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
+        return math_1(args, func, can_overflow);                            \
+    }\
+    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); \
-        }\
-        PyDoc_STRVAR(math_##funcname##_doc, docstring);
+    static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
+        return math_2(args, func, #funcname); \
+    }\
+    PyDoc_STRVAR(math_##funcname##_doc, docstring);
 
 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,
@@ -333 +820 @@
       "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.")
 FUNC1(ceil, ceil, 0,
@@ -339 +826 @@
       "This is the smallest integral value >= x.")
 FUNC2(copysign, copysign,
-      "copysign(x,y)\n\nReturn x with the sign of y.")
+      "copysign(x, y)\n\nReturn x with the sign of y.")
 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.")
@@ -351 +846 @@
       "floor(x)\n\nReturn the floor of x as a float.\n"
       "This is the largest integral value <= x.")
-FUNC1(log1p, log1p, 1,
-      "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\
-      The result is computed in a way which is accurate for x near zero.")
+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, 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.")
 FUNC1(sin, sin, 0,
       "sin(x)\n\nReturn the sine of x (measured in radians).")
@@ -385 +884 @@
    regular doubles instead of extended long precision (80-bit) values.  This
    prevents double rounding because any addition or subtraction of two doubles
-   can be resolved exactly into double-sized hi and lo values.  As long as the 
+   can be resolved exactly into double-sized hi and lo values.  As long as the
    hi value gets forced into a double before yr and lo are computed, the extra
    bits in downstream extended precision operations (x87 for example) will be
@@ -408 +907 @@
              double  *ps,    Py_ssize_t *m_ptr)
 {
-        void *v = NULL;
-        Py_ssize_t m = *m_ptr;
-
-        m += m;  /* double */
-        if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
-                double *p = *p_ptr;
-                if (p == ps) {
-                        v = PyMem_Malloc(sizeof(double) * m);
-                        if (v != NULL)
-                                memcpy(v, ps, sizeof(double) * n);
-                }
-                else
-                        v = PyMem_Realloc(p, sizeof(double) * m);
-        }
-        if (v == NULL) {        /* size overflow or no memory */
-                PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
-                return 1;
-        }
-        *p_ptr = (double*) v;
-        *m_ptr = m;
-        return 0;
+    void *v = NULL;
+    Py_ssize_t m = *m_ptr;
+
+    m += m;  /* double */
+    if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
+        double *p = *p_ptr;
+        if (p == ps) {
+            v = PyMem_Malloc(sizeof(double) * m);
+            if (v != NULL)
+                memcpy(v, ps, sizeof(double) * n);
+        }
+        else
+            v = PyMem_Realloc(p, sizeof(double) * m);
+    }
+    if (v == NULL) {        /* size overflow or no memory */
+        PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
+        return 1;
+    }
+    *p_ptr = (double*) v;
+    *m_ptr = m;
+    return 0;
 }
 
@@ -464 +963 @@
 math_fsum(PyObject *self, PyObject *seq)
 {
-        PyObject *item, *iter, *sum = NULL;
-        Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
-        double x, y, t, ps[NUM_PARTIALS], *p = ps;
-        double xsave, special_sum = 0.0, inf_sum = 0.0;
-        volatile double hi, yr, lo;
-
-        iter = PyObject_GetIter(seq);
-        if (iter == NULL)
-                return NULL;
-
-        PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
-
-        for(;;) {           /* for x in iterable */
-                assert(0 <= n && n <= m);
-                assert((m == NUM_PARTIALS && p == ps) ||
-                       (m >  NUM_PARTIALS && p != NULL));
-
-                item = PyIter_Next(iter);
-                if (item == NULL) {
-                        if (PyErr_Occurred())
-                                goto _fsum_error;
-                        break;
-                }
-                x = PyFloat_AsDouble(item);
-                Py_DECREF(item);
-                if (PyErr_Occurred())
-                        goto _fsum_error;
-
-                xsave = x;
-                for (i = j = 0; j < n; j++) {       /* for y in partials */
-                        y = p[j];
-                        if (fabs(x) < fabs(y)) {
-                                t = x; x = y; y = t;
-                        }
-                        hi = x + y;
-                        yr = hi - x;
-                        lo = y - yr;
-                        if (lo != 0.0)
-                                p[i++] = lo;
-                        x = hi;
-                }
-
-                n = i;                              /* ps[i:] = [x] */
-                if (x != 0.0) {
-                        if (! Py_IS_FINITE(x)) {
-                                /* a nonfinite x could arise either as
-                                   a result of intermediate overflow, or
-                                   as a result of a nan or inf in the
-                                   summands */
-                                if (Py_IS_FINITE(xsave)) {
-                                        PyErr_SetString(PyExc_OverflowError,
-                                              "intermediate overflow in fsum");
-                                        goto _fsum_error;
-                                }
-                                if (Py_IS_INFINITY(xsave))
-                                        inf_sum += xsave;
-                                special_sum += xsave;
-                                /* reset partials */
-                                n = 0;
-                        }
-                        else if (n >= m && _fsum_realloc(&p, n, ps, &m))
-                                goto _fsum_error;
-                        else
-                                p[n++] = x;
-                }
-        }
-
-        if (special_sum != 0.0) {
-                if (Py_IS_NAN(inf_sum))
-                        PyErr_SetString(PyExc_ValueError,
-                                        "-inf + inf in fsum");
-                else
-                        sum = PyFloat_FromDouble(special_sum);
-                goto _fsum_error;
-        }
-
-        hi = 0.0;
-        if (n > 0) {
-                hi = p[--n];
-                /* sum_exact(ps, hi) from the top, stop when the sum becomes
-                   inexact. */
-                while (n > 0) {
-                        x = hi;
-                        y = p[--n];
-                        assert(fabs(y) < fabs(x));
-                        hi = x + y;
-                        yr = hi - x;
-                        lo = y - yr;
-                        if (lo != 0.0)
-                                break;
-                }
-                /* Make half-even rounding work across multiple partials.
-                   Needed so that sum([1e-16, 1, 1e16]) will round-up the last
-                   digit to two instead of down to zero (the 1e-16 makes the 1
-                   slightly closer to two).  With a potential 1 ULP rounding
-                   error fixed-up, math.fsum() can guarantee commutativity. */
-                if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
-                              (lo > 0.0 && p[n-1] > 0.0))) {
-                        y = lo * 2.0;
-                        x = hi + y;
-                        yr = x - hi;
-                        if (y == yr)
-                                hi = x;
-                }
-        }
-        sum = PyFloat_FromDouble(hi);
+    PyObject *item, *iter, *sum = NULL;
+    Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
+    double x, y, t, ps[NUM_PARTIALS], *p = ps;
+    double xsave, special_sum = 0.0, inf_sum = 0.0;
+    volatile double hi, yr, lo;
+
+    iter = PyObject_GetIter(seq);
+    if (iter == NULL)
+        return NULL;
+
+    PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
+
+    for(;;) {           /* for x in iterable */
+        assert(0 <= n && n <= m);
+        assert((m == NUM_PARTIALS && p == ps) ||
+               (m >  NUM_PARTIALS && p != NULL));
+
+        item = PyIter_Next(iter);
+        if (item == NULL) {
+            if (PyErr_Occurred())
+                goto _fsum_error;
+            break;
+        }
+        x = PyFloat_AsDouble(item);
+        Py_DECREF(item);
+        if (PyErr_Occurred())
+            goto _fsum_error;
+
+        xsave = x;
+        for (i = j = 0; j < n; j++) {       /* for y in partials */
+            y = p[j];
+            if (fabs(x) < fabs(y)) {
+                t = x; x = y; y = t;
+            }
+            hi = x + y;
+            yr = hi - x;
+            lo = y - yr;
+            if (lo != 0.0)
+                p[i++] = lo;
+            x = hi;
+        }
+
+        n = i;                              /* ps[i:] = [x] */
+        if (x != 0.0) {
+            if (! Py_IS_FINITE(x)) {
+                /* a nonfinite x could arise either as
+                   a result of intermediate overflow, or
+                   as a result of a nan or inf in the
+                   summands */
+                if (Py_IS_FINITE(xsave)) {
+                    PyErr_SetString(PyExc_OverflowError,
+                          "intermediate overflow in fsum");
+                    goto _fsum_error;
+                }
+                if (Py_IS_INFINITY(xsave))
+                    inf_sum += xsave;
+                special_sum += xsave;
+                /* reset partials */
+                n = 0;
+            }
+            else if (n >= m && _fsum_realloc(&p, n, ps, &m))
+                goto _fsum_error;
+            else
+                p[n++] = x;
+        }
+    }
+
+    if (special_sum != 0.0) {
+        if (Py_IS_NAN(inf_sum))
+            PyErr_SetString(PyExc_ValueError,
+                            "-inf + inf in fsum");
+        else
+            sum = PyFloat_FromDouble(special_sum);
+        goto _fsum_error;
+    }
+
+    hi = 0.0;
+    if (n > 0) {
+        hi = p[--n];
+        /* sum_exact(ps, hi) from the top, stop when the sum becomes
+           inexact. */
+        while (n > 0) {
+            x = hi;
+            y = p[--n];
+            assert(fabs(y) < fabs(x));
+            hi = x + y;
+            yr = hi - x;
+            lo = y - yr;
+            if (lo != 0.0)
+                break;
+        }
+        /* Make half-even rounding work across multiple partials.
+           Needed so that sum([1e-16, 1, 1e16]) will round-up the last
+           digit to two instead of down to zero (the 1e-16 makes the 1
+           slightly closer to two).  With a potential 1 ULP rounding
+           error fixed-up, math.fsum() can guarantee commutativity. */
+        if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
+                      (lo > 0.0 && p[n-1] > 0.0))) {
+            y = lo * 2.0;
+            x = hi + y;
+            yr = x - hi;
+            if (y == yr)
+                hi = x;
+        }
+    }
+    sum = PyFloat_FromDouble(hi);
 
 _fsum_error:
-        PyFPE_END_PROTECT(hi)
-        Py_DECREF(iter);
-        if (p != ps)
-                PyMem_Free(p);
-        return sum;
+    PyFPE_END_PROTECT(hi)
+    Py_DECREF(iter);
+    if (p != ps)
+        PyMem_Free(p);
+    return sum;
 }
 
@@ -582 +1081 @@
 
 PyDoc_STRVAR(math_fsum_doc,
-"sum(iterable)\n\n\
+"fsum(iterable)\n\n\
 Return an accurate floating point sum of values in the iterable.\n\
 Assumes IEEE-754 floating point arithmetic.");
@@ -589 +1088 @@
 math_factorial(PyObject *self, PyObject *arg)
 {
-        long i, x;
-        PyObject *result, *iobj, *newresult;
-
-        if (PyFloat_Check(arg)) {
-                double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);
-                if (dx != floor(dx)) {
-                        PyErr_SetString(PyExc_ValueError, 
-                                "factorial() only accepts integral values");
-                        return NULL;
-                }
-        }
-
-        x = PyInt_AsLong(arg);
-        if (x == -1 && PyErr_Occurred())
-                return NULL;
-        if (x < 0) {
-                PyErr_SetString(PyExc_ValueError, 
-                        "factorial() not defined for negative values");
-                return NULL;
-        }
-
-        result = (PyObject *)PyInt_FromLong(1);
-        if (result == NULL)
-                return NULL;
-        for (i=1 ; i<=x ; i++) {
-                iobj = (PyObject *)PyInt_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;
-        }
-        return result;
+    long i, x;
+    PyObject *result, *iobj, *newresult;
+
+    if (PyFloat_Check(arg)) {
+        PyObject *lx;
+        double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);
+        if (!(Py_IS_FINITE(dx) && dx == floor(dx))) {
+            PyErr_SetString(PyExc_ValueError,
+                "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 = PyInt_AsLong(arg);
+
+    if (x == -1 && PyErr_Occurred())
+        return NULL;
+    if (x < 0) {
+        PyErr_SetString(PyExc_ValueError,
+            "factorial() not defined for negative values");
+        return NULL;
+    }
+
+    result = (PyObject *)PyInt_FromLong(1);
+    if (result == NULL)
+        return NULL;
+    for (i=1 ; i<=x ; i++) {
+        iobj = (PyObject *)PyInt_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;
+    }
+    return result;
 
 error:
-        Py_DECREF(result);
-        return NULL;
+    Py_DECREF(result);
+    return NULL;
 }
 
@@ -639 +1145 @@
 math_trunc(PyObject *self, PyObject *number)
 {
-        return PyObject_CallMethod(number, "__trunc__", NULL);
+    return PyObject_CallMethod(number, "__trunc__", NULL);
 }
 
@@ -650 +1156 @@
 math_frexp(PyObject *self, PyObject *arg)
 {
-        int i;
-        double x = PyFloat_AsDouble(arg);
-        if (x == -1.0 && PyErr_Occurred())
-                return NULL;
-        /* deal with special cases directly, to sidestep platform
-           differences */
-        if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
-                i = 0;
-        }
-        else {
-                PyFPE_START_PROTECT("in math_frexp", return 0);
-                x = frexp(x, &i);
-                PyFPE_END_PROTECT(x);
-        }
-        return Py_BuildValue("(di)", x, i);
+    int i;
+    double x = PyFloat_AsDouble(arg);
+    if (x == -1.0 && PyErr_Occurred())
+        return NULL;
+    /* deal with special cases directly, to sidestep platform
+       differences */
+    if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
+        i = 0;
+    }
+    else {
+        PyFPE_START_PROTECT("in math_frexp", return 0);
+        x = frexp(x, &i);
+        PyFPE_END_PROTECT(x);
+    }
+    return Py_BuildValue("(di)", x, i);
 }
 
@@ -677 +1183 @@
 math_ldexp(PyObject *self, PyObject *args)
 {
-        double x, r;
-        PyObject *oexp;
-        long exp;
-        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;
-                        }
-                }
-        }
-        else if (PyInt_Check(oexp)) {
-                exp = PyInt_AS_LONG(oexp);
-        }
-        else {
-                PyErr_SetString(PyExc_TypeError,
-                                "Expected an int or long as second argument "
-                                "to ldexp.");
-                return NULL;
-        }
-
-        if (x == 0. || !Py_IS_FINITE(x)) {
-                /* NaNs, zeros and infinities are returned unchanged */
-                r = x;
-                errno = 0;
-        } else if (exp > INT_MAX) {
-                /* overflow */
-                r = copysign(Py_HUGE_VAL, x);
-                errno = ERANGE;
-        } else if (exp < INT_MIN) {
-                /* underflow to +-0 */
-                r = copysign(0., x);
-                errno = 0;
-        } else {
-                errno = 0;
-                PyFPE_START_PROTECT("in math_ldexp", return 0);
-                r = ldexp(x, (int)exp);
-                PyFPE_END_PROTECT(r);
-                if (Py_IS_INFINITY(r))
-                        errno = ERANGE;
-        }
-
-        if (errno && is_error(r))
-                return NULL;
-        return PyFloat_FromDouble(r);
+    double x, r;
+    PyObject *oexp;
+    long exp;
+    int overflow;
+    if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp))
+        return NULL;
+
+    if (PyLong_Check(oexp) || PyInt_Check(oexp)) {
+        /* on overflow, replace exponent with either LONG_MAX
+           or LONG_MIN, depending on the sign. */
+        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,
+                        "Expected an int or long as second argument "
+                        "to ldexp.");
+        return NULL;
+    }
+
+    if (x == 0. || !Py_IS_FINITE(x)) {
+        /* NaNs, zeros and infinities are returned unchanged */
+        r = x;
+        errno = 0;
+    } else if (exp > INT_MAX) {
+        /* overflow */
+        r = copysign(Py_HUGE_VAL, x);
+        errno = ERANGE;
+    } else if (exp < INT_MIN) {
+        /* underflow to +-0 */
+        r = copysign(0., x);
+        errno = 0;
+    } else {
+        errno = 0;
+        PyFPE_START_PROTECT("in math_ldexp", return 0);
+        r = ldexp(x, (int)exp);
+        PyFPE_END_PROTECT(r);
+        if (Py_IS_INFINITY(r))
+            errno = ERANGE;
+    }
+
+    if (errno && is_error(r))
+        return NULL;
+    return PyFloat_FromDouble(r);
 }
 
 PyDoc_STRVAR(math_ldexp_doc,
-"ldexp(x, i) -> x * (2**i)");
+"ldexp(x, i)\n\n\
+Return x * (2**i).");
 
 static PyObject *
 math_modf(PyObject *self, PyObject *arg)
 {
-        double y, x = PyFloat_AsDouble(arg);
-        if (x == -1.0 && PyErr_Occurred())
-                return NULL;
-        /* some platforms don't do the right thing for NaNs and
-           infinities, so we take care of special cases directly. */
-        if (!Py_IS_FINITE(x)) {
-                if (Py_IS_INFINITY(x))
-                        return Py_BuildValue("(dd)", copysign(0., x), x);
-                else if (Py_IS_NAN(x))
-                        return Py_BuildValue("(dd)", x, x);
-        }          
-
-        errno = 0;
-        PyFPE_START_PROTECT("in math_modf", return 0);
-        x = modf(x, &y);
-        PyFPE_END_PROTECT(x);
-        return Py_BuildValue("(dd)", x, y);
+    double y, x = PyFloat_AsDouble(arg);
+    if (x == -1.0 && PyErr_Occurred())
+        return NULL;
+    /* some platforms don't do the right thing for NaNs and
+       infinities, so we take care of special cases directly. */
+    if (!Py_IS_FINITE(x)) {
+        if (Py_IS_INFINITY(x))
+            return Py_BuildValue("(dd)", copysign(0., x), x);
+        else if (Py_IS_NAN(x))
+            return Py_BuildValue("(dd)", x, x);
+    }
+
+    errno = 0;
+    PyFPE_START_PROTECT("in math_modf", return 0);
+    x = modf(x, &y);
+    PyFPE_END_PROTECT(x);
+    return Py_BuildValue("(dd)", x, y);
 }
 
@@ -772 +1266 @@
 /* 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) {
-                        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);
-        }
-
-        /* Else let libm handle it by itself. */
-        return math_1(arg, func, 0);
+    /* If it is long, do it ourselves. */
+    if (PyLong_Check(arg)) {
+        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;
+        }
+
+        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. */
+    return math_1(arg, func, 0);
 }
 
@@ -806 +1314 @@
 math_log(PyObject *self, PyObject *args)
 {
-        PyObject *arg;
-        PyObject *base = NULL;
-        PyObject *num, *den;
-        PyObject *ans;
-
-        if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
-                return NULL;
-
-        num = loghelper(arg, m_log, "log");
-        if (num == NULL || base == NULL)
-                return num;
-
-        den = loghelper(base, m_log, "log");
-        if (den == NULL) {
-                Py_DECREF(num);
-                return NULL;
-        }
-
-        ans = PyNumber_Divide(num, den);
-        Py_DECREF(num);
-        Py_DECREF(den);
-        return ans;
+    PyObject *arg;
+    PyObject *base = NULL;
+    PyObject *num, *den;
+    PyObject *ans;
+
+    if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
+        return NULL;
+
+    num = loghelper(arg, m_log, "log");
+    if (num == NULL || base == NULL)
+        return num;
+
+    den = loghelper(base, m_log, "log");
+    if (den == NULL) {
+        Py_DECREF(num);
+        return NULL;
+    }
+
+    ans = PyNumber_Divide(num, den);
+    Py_DECREF(num);
+    Py_DECREF(den);
+    return ans;
 }
 
 PyDoc_STRVAR(math_log_doc,
-"log(x[, base]) -> the logarithm of x to the given base.\n\
+"log(x[, base])\n\n\
+Return the logarithm of x to the given base.\n\
 If the base not specified, returns the natural logarithm (base e) of x.");
 
@@ -837 +1346 @@
 math_log10(PyObject *self, PyObject *arg)
 {
-        return loghelper(arg, m_log10, "log10");
+    return loghelper(arg, m_log10, "log10");
 }
 
 PyDoc_STRVAR(math_log10_doc,
-"log10(x) -> the base 10 logarithm of x.");
+"log10(x)\n\nReturn the base 10 logarithm of x.");
 
 static PyObject *
 math_fmod(PyObject *self, PyObject *args)
 {
-        PyObject *ox, *oy;
-        double r, x, y;
-        if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
-                return NULL;
-        x = PyFloat_AsDouble(ox);
-        y = PyFloat_AsDouble(oy);
-        if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
-                return NULL;
-        /* fmod(x, +/-Inf) returns x for finite x. */
-        if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
-                return PyFloat_FromDouble(x);
-        errno = 0;
-        PyFPE_START_PROTECT("in math_fmod", return 0);
-        r = fmod(x, y);
-        PyFPE_END_PROTECT(r);
-        if (Py_IS_NAN(r)) {
-                if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
-                        errno = EDOM;
-                else
-                        errno = 0;
-        }
-        if (errno && is_error(r))
-                return NULL;
-        else
-                return PyFloat_FromDouble(r);
+    PyObject *ox, *oy;
+    double r, x, y;
+    if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
+        return NULL;
+    x = PyFloat_AsDouble(ox);
+    y = PyFloat_AsDouble(oy);
+    if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
+        return NULL;
+    /* fmod(x, +/-Inf) returns x for finite x. */
+    if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
+        return PyFloat_FromDouble(x);
+    errno = 0;
+    PyFPE_START_PROTECT("in math_fmod", return 0);
+    r = fmod(x, y);
+    PyFPE_END_PROTECT(r);
+    if (Py_IS_NAN(r)) {
+        if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
+            errno = EDOM;
+        else
+            errno = 0;
+    }
+    if (errno && is_error(r))
+        return NULL;
+    else
+        return PyFloat_FromDouble(r);
 }
 
 PyDoc_STRVAR(math_fmod_doc,
-"fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
+"fmod(x, y)\n\nReturn fmod(x, y), according to platform C."
 "  x % y may differ.");
 
@@ -880 +1389 @@
 math_hypot(PyObject *self, PyObject *args)
 {
-        PyObject *ox, *oy;
-        double r, x, y;
-        if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
-                return NULL;
-        x = PyFloat_AsDouble(ox);
-        y = PyFloat_AsDouble(oy);
-        if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
-                return NULL;
-        /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
-        if (Py_IS_INFINITY(x))
-                return PyFloat_FromDouble(fabs(x));
-        if (Py_IS_INFINITY(y))
-                return PyFloat_FromDouble(fabs(y));
-        errno = 0;
-        PyFPE_START_PROTECT("in math_hypot", return 0);
-        r = hypot(x, y);
-        PyFPE_END_PROTECT(r);
-        if (Py_IS_NAN(r)) {
-                if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
-                        errno = EDOM;
-                else
-                        errno = 0;
-        }
-        else if (Py_IS_INFINITY(r)) {
-                if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
-                        errno = ERANGE;
-                else
-                        errno = 0;
-        }
-        if (errno && is_error(r))
-                return NULL;
-        else
-                return PyFloat_FromDouble(r);
+    PyObject *ox, *oy;
+    double r, x, y;
+    if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
+        return NULL;
+    x = PyFloat_AsDouble(ox);
+    y = PyFloat_AsDouble(oy);
+    if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
+        return NULL;
+    /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
+    if (Py_IS_INFINITY(x))
+        return PyFloat_FromDouble(fabs(x));
+    if (Py_IS_INFINITY(y))
+        return PyFloat_FromDouble(fabs(y));
+    errno = 0;
+    PyFPE_START_PROTECT("in math_hypot", return 0);
+    r = hypot(x, y);
+    PyFPE_END_PROTECT(r);
+    if (Py_IS_NAN(r)) {
+        if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
+            errno = EDOM;
+        else
+            errno = 0;
+    }
+    else if (Py_IS_INFINITY(r)) {
+        if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
+            errno = ERANGE;
+        else
+            errno = 0;
+    }
+    if (errno && is_error(r))
+        return NULL;
+    else
+        return PyFloat_FromDouble(r);
 }
 
 PyDoc_STRVAR(math_hypot_doc,
-"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
+"hypot(x, y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
 
 /* pow can't use math_2, but needs its own wrapper: the problem is
@@ -927 +1436 @@
 math_pow(PyObject *self, PyObject *args)
 {
-        PyObject *ox, *oy;
-        double r, x, y;
-        int odd_y;
-
-        if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
-                return NULL;
-        x = PyFloat_AsDouble(ox);
-        y = PyFloat_AsDouble(oy);
-        if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
-                return NULL;
-
-        /* deal directly with IEEE specials, to cope with problems on various
-           platforms whose semantics don't exactly match C99 */
-        r = 0.; /* silence compiler warning */
-        if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
-                errno = 0;
-                if (Py_IS_NAN(x))
-                        r = y == 0. ? 1. : x; /* NaN**0 = 1 */
-                else if (Py_IS_NAN(y))
-                        r = x == 1. ? 1. : y; /* 1**NaN = 1 */
-                else if (Py_IS_INFINITY(x)) {
-                        odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
-                        if (y > 0.)
-                                r = odd_y ? x : fabs(x);
-                        else if (y == 0.)
-                                r = 1.;
-                        else /* y < 0. */
-                                r = odd_y ? copysign(0., x) : 0.;
-                }
-                else if (Py_IS_INFINITY(y)) {
-                        if (fabs(x) == 1.0)
-                                r = 1.;
-                        else if (y > 0. && fabs(x) > 1.0)
-                                r = y;
-                        else if (y < 0. && fabs(x) < 1.0) {
-                                r = -y; /* result is +inf */
-                                if (x == 0.) /* 0**-inf: divide-by-zero */
-                                        errno = EDOM;
-                        }
-                        else
-                                r = 0.;
-                }
-        }
-        else {
-                /* let libm handle finite**finite */
-                errno = 0;
-                PyFPE_START_PROTECT("in math_pow", return 0);
-                r = pow(x, y);
-                PyFPE_END_PROTECT(r);
-                /* a NaN result should arise only from (-ve)**(finite
-                   non-integer); in this case we want to raise ValueError. */
-                if (!Py_IS_FINITE(r)) {
-                        if (Py_IS_NAN(r)) {
-                                errno = EDOM;
-                        }
-                        /* 
-                           an infinite result here arises either from:
-                           (A) (+/-0.)**negative (-> divide-by-zero)
-                           (B) overflow of x**y with x and y finite
-                        */
-                        else if (Py_IS_INFINITY(r)) {
-                                if (x == 0.)
-                                        errno = EDOM;
-                                else
-                                        errno = ERANGE;
-                        }
-                }
-        }
-
-        if (errno && is_error(r))
-                return NULL;
-        else
-                return PyFloat_FromDouble(r);
+    PyObject *ox, *oy;
+    double r, x, y;
+    int odd_y;
+
+    if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
+        return NULL;
+    x = PyFloat_AsDouble(ox);
+    y = PyFloat_AsDouble(oy);
+    if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
+        return NULL;
+
+    /* deal directly with IEEE specials, to cope with problems on various
+       platforms whose semantics don't exactly match C99 */
+    r = 0.; /* silence compiler warning */
+    if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
+        errno = 0;
+        if (Py_IS_NAN(x))
+            r = y == 0. ? 1. : x; /* NaN**0 = 1 */
+        else if (Py_IS_NAN(y))
+            r = x == 1. ? 1. : y; /* 1**NaN = 1 */
+        else if (Py_IS_INFINITY(x)) {
+            odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
+            if (y > 0.)
+                r = odd_y ? x : fabs(x);
+            else if (y == 0.)
+                r = 1.;
+            else /* y < 0. */
+                r = odd_y ? copysign(0., x) : 0.;
+        }
+        else if (Py_IS_INFINITY(y)) {
+            if (fabs(x) == 1.0)
+                r = 1.;
+            else if (y > 0. && fabs(x) > 1.0)
+                r = y;
+            else if (y < 0. && fabs(x) < 1.0) {
+                r = -y; /* result is +inf */
+                if (x == 0.) /* 0**-inf: divide-by-zero */
+                    errno = EDOM;
+            }
+            else
+                r = 0.;
+        }
+    }
+    else {
+        /* let libm handle finite**finite */
+        errno = 0;
+        PyFPE_START_PROTECT("in math_pow", return 0);
+        r = pow(x, y);
+        PyFPE_END_PROTECT(r);
+        /* a NaN result should arise only from (-ve)**(finite
+           non-integer); in this case we want to raise ValueError. */
+        if (!Py_IS_FINITE(r)) {
+            if (Py_IS_NAN(r)) {
+                errno = EDOM;
+            }
+            /*
+               an infinite result here arises either from:
+               (A) (+/-0.)**negative (-> divide-by-zero)
+               (B) overflow of x**y with x and y finite
+            */
+            else if (Py_IS_INFINITY(r)) {
+                if (x == 0.)
+                    errno = EDOM;
+                else
+                    errno = ERANGE;
+            }
+        }
+    }
+
+    if (errno && is_error(r))
+        return NULL;
+    else
+        return PyFloat_FromDouble(r);
 }
 
 PyDoc_STRVAR(math_pow_doc,
-"pow(x,y)\n\nReturn x**y (x to the power of y).");
+"pow(x, y)\n\nReturn x**y (x to the power of y).");
 
 static const double degToRad = Py_MATH_PI / 180.0;
@@ -1011 +1520 @@
 math_degrees(PyObject *self, PyObject *arg)
 {
-        double x = PyFloat_AsDouble(arg);
-        if (x == -1.0 && PyErr_Occurred())
-                return NULL;
-        return PyFloat_FromDouble(x * radToDeg);
+    double x = PyFloat_AsDouble(arg);
+    if (x == -1.0 && PyErr_Occurred())
+        return NULL;
+    return PyFloat_FromDouble(x * radToDeg);
 }
 
 PyDoc_STRVAR(math_degrees_doc,
-"degrees(x) -> converts angle x from radians to degrees");
+"degrees(x)\n\n\
+Convert angle x from radians to degrees.");
 
 static PyObject *
 math_radians(PyObject *self, PyObject *arg)
 {
-        double x = PyFloat_AsDouble(arg);
-        if (x == -1.0 && PyErr_Occurred())
-                return NULL;
-        return PyFloat_FromDouble(x * degToRad);
+    double x = PyFloat_AsDouble(arg);
+    if (x == -1.0 && PyErr_Occurred())
+        return NULL;
+    return PyFloat_FromDouble(x * degToRad);
 }
 
 PyDoc_STRVAR(math_radians_doc,
-"radians(x) -> converts angle x from degrees to radians");
+"radians(x)\n\n\
+Convert angle x from degrees to radians.");
 
 static PyObject *
 math_isnan(PyObject *self, PyObject *arg)
 {
-        double x = PyFloat_AsDouble(arg);
-        if (x == -1.0 && PyErr_Occurred())
-                return NULL;
-        return PyBool_FromLong((long)Py_IS_NAN(x));
+    double x = PyFloat_AsDouble(arg);
+    if (x == -1.0 && PyErr_Occurred())
+        return NULL;
+    return PyBool_FromLong((long)Py_IS_NAN(x));
 }
 
 PyDoc_STRVAR(math_isnan_doc,
-"isnan(x) -> bool\n\
-Checks if float x is not a number (NaN)");
+"isnan(x) -> bool\n\n\
+Check if float x is not a number (NaN).");
 
 static PyObject *
 math_isinf(PyObject *self, PyObject *arg)
 {
-        double x = PyFloat_AsDouble(arg);
-        if (x == -1.0 && PyErr_Occurred())
-                return NULL;
-        return PyBool_FromLong((long)Py_IS_INFINITY(x));
+    double x = PyFloat_AsDouble(arg);
+    if (x == -1.0 && PyErr_Occurred())
+        return NULL;
+    return PyBool_FromLong((long)Py_IS_INFINITY(x));
 }
 
 PyDoc_STRVAR(math_isinf_doc,
-"isinf(x) -> bool\n\
-Checks if float x is infinite (positive or negative)");
+"isinf(x) -> bool\n\n\
+Check if float x is infinite (positive or negative).");
 
 static PyMethodDef math_methods[] = {
-        {"acos",        math_acos,      METH_O,         math_acos_doc},
-        {"acosh",       math_acosh,     METH_O,         math_acosh_doc},
-        {"asin",        math_asin,      METH_O,         math_asin_doc},
-        {"asinh",       math_asinh,     METH_O,         math_asinh_doc},
-        {"atan",        math_atan,      METH_O,         math_atan_doc},
-        {"atan2",       math_atan2,     METH_VARARGS,   math_atan2_doc},
-        {"atanh",       math_atanh,     METH_O,         math_atanh_doc},
-        {"ceil",        math_ceil,      METH_O,         math_ceil_doc},
-        {"copysign",    math_copysign,  METH_VARARGS,   math_copysign_doc},
-        {"cos",         math_cos,       METH_O,         math_cos_doc},
-        {"cosh",        math_cosh,      METH_O,         math_cosh_doc},
-        {"degrees",     math_degrees,   METH_O,         math_degrees_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},
-        {"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},
-        {"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},
-        {"log",         math_log,       METH_VARARGS,   math_log_doc},
-        {"log1p",       math_log1p,     METH_O,         math_log1p_doc},
-        {"log10",       math_log10,     METH_O,         math_log10_doc},
-        {"modf",        math_modf,      METH_O,         math_modf_doc},
-        {"pow",         math_pow,       METH_VARARGS,   math_pow_doc},
-        {"radians",     math_radians,   METH_O,         math_radians_doc},
-        {"sin",         math_sin,       METH_O,         math_sin_doc},
-        {"sinh",        math_sinh,      METH_O,         math_sinh_doc},
-        {"sqrt",        math_sqrt,      METH_O,         math_sqrt_doc},
-        {"tan",         math_tan,       METH_O,         math_tan_doc},
-        {"tanh",        math_tanh,      METH_O,         math_tanh_doc},
-        {"trunc",       math_trunc,     METH_O,         math_trunc_doc},
-        {NULL,          NULL}           /* sentinel */
+    {"acos",            math_acos,      METH_O,         math_acos_doc},
+    {"acosh",           math_acosh,     METH_O,         math_acosh_doc},
+    {"asin",            math_asin,      METH_O,         math_asin_doc},
+    {"asinh",           math_asinh,     METH_O,         math_asinh_doc},
+    {"atan",            math_atan,      METH_O,         math_atan_doc},
+    {"atan2",           math_atan2,     METH_VARARGS,   math_atan2_doc},
+    {"atanh",           math_atanh,     METH_O,         math_atanh_doc},
+    {"ceil",            math_ceil,      METH_O,         math_ceil_doc},
+    {"copysign",        math_copysign,  METH_VARARGS,   math_copysign_doc},
+    {"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},
+    {"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},
+    {"modf",            math_modf,      METH_O,         math_modf_doc},
+    {"pow",             math_pow,       METH_VARARGS,   math_pow_doc},
+    {"radians",         math_radians,   METH_O,         math_radians_doc},
+    {"sin",             math_sin,       METH_O,         math_sin_doc},
+    {"sinh",            math_sinh,      METH_O,         math_sinh_doc},
+    {"sqrt",            math_sqrt,      METH_O,         math_sqrt_doc},
+    {"tan",             math_tan,       METH_O,         math_tan_doc},
+    {"tanh",            math_tanh,      METH_O,         math_tanh_doc},
+    {"trunc",           math_trunc,     METH_O,         math_trunc_doc},
+    {NULL,              NULL}           /* sentinel */
 };
 
@@ -1105 +1621 @@
 initmath(void)
 {
-        PyObject *m;
-
-        m = Py_InitModule3("math", math_methods, module_doc);
-        if (m == NULL)
-                goto finally;
-
-        PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
-        PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
+    PyObject *m;
+
+    m = Py_InitModule3("math", math_methods, module_doc);
+    if (m == NULL)
+        goto finally;
+
+    PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
+    PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
 
     finally:
-        return;
-}
+    return;
+}