Modules/mathmodule.c

cpython 3.14 @ ab2d84fe1023/Modules/mathmodule.c ↗

The math module. A thin Go-side port is pending; the module currently falls back to the pure-Python equivalent for any function not yet bridged.

The file is organized around three implementation strategies:

• libm delegation. Most trig, exponential, and logarithm functions call the C standard library directly (sin, cos, exp, log, etc.) through a uniform math_1 / math_2 dispatch shim that converts Python floats to double, checks errno, and converts NaN/inf results to Python exceptions or special values.
• Pure-integer algorithms. factorial, comb, and perm operate on Python integers without floating-point rounding. factorial uses a divide-and-conquer prime-swing or partial-product method; comb uses the multiplicative formula with intermediate reduction.
• Compensated summation. fsum maintains an adaptive list of partial sums (Shewchuk's algorithm) to return an exactly-rounded result even for inputs that would lose precision in a naive accumulation.

math.floor, math.ceil, and math.trunc are special: they try the nb_floor, nb_ceil, and __trunc__ special methods first (to support Decimal and Fraction), and only call libm on plain float objects.

Map

Lines	Symbol	Role	gopy
1-120	includes, math_1, math_2, is_error	Uniform libm-call dispatch and error-checking shims.	module/math/module.go:math1
120-400	math_sin, math_cos, math_tan, math_asin, math_acos, math_atan, math_atan2, math_exp, math_expm1, math_log, math_log1p, math_log2, math_log10, math_pow, math_sqrt	libm-delegated trig, exponential, and log functions.	module/math/module.go:Sin … Sqrt
400-900	math_fsum	Exact floating-point summation via partial sums list (Shewchuk).	module/math/module.go:Fsum
900-1500	math_factorial, factorial_partial_product, factorial_odd_part, math_comb, math_perm	Integer combinatorics. No floating-point is used.	module/math/module.go:Factorial
1500-1800	math_floor, math_ceil, math_trunc	Delegate to nb_floor / nb_ceil / __trunc__ first, then libm.	module/math/module.go:Floor
1800-2200	math_isfinite, math_isinf, math_isnan, math_isclose, math_copysign, math_fabs, math_fmod, math_modf, math_frexp, math_ldexp	Numeric predicates and low-level float manipulation.	module/math/module.go:IsFinite
2200-2700	math_hypot, vector_norm	Multi-argument hypot with compensated accumulation to avoid intermediate overflow.	module/math/module.go:Hypot
2700-2950	math_prod, math_gcd, math_lcm	Integer-aware product, GCD (Euclidean), and LCM.	module/math/module.go:Prod
2950-3200	math_methods table, constant definitions, PyInit_math	Module definition. Registers pi, e, tau, inf, nan.	module/math/module.go:Module

Reading

fsum partial sums algorithm (lines 400 to 900)

cpython 3.14 @ ab2d84fe1023/Modules/mathmodule.c#L400-900 ↗

math.fsum returns the exact sum of a finite iterable of floats, rounded once to the nearest representable double. The core is Shewchuk's expansion-summation algorithm, which maintains a growing list of partial sums such that their exact mathematical sum equals the true sum of all inputs processed so far.

static PyObject *
math_fsum_impl(PyObject *module, PyObject *seq)
{
    PyObject *item, *iter;
    Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
    double x, y, t, *p = ps;
    double xsave, special_sum = 0.0, inf_sum = 0.0;
    ...
    iter = PyObject_GetIter(seq);
    while ((item = PyIter_Next(iter))) {
        x = PyFloat_AsDouble(item);
        ...
        /* Merge x into the partials list */
        for (i = j = 0; j < n; j++) {
            y = p[j];
            if (fabs(x) < fabs(y)) { double tmp = x; x = y; y = tmp; }
            /* Two-sum: hi + lo = x + y exactly */
            t = x + y;
            if (fabs(t) >= fabs(x)) {
                y = y - (t - x);
            } else {
                y = x - (t - y);
                x = t;
            }
            if (y != 0.0)
                p[i++] = y;
            x = t;
        }
        n = i;
        if (x != 0.0) {
            if (n >= m) { /* grow partials array */ ... }
            p[n++] = x;
        }
    }
    /* Round-to-nearest from the partials list */
    ...
}

The two-sum trick (t = x + y plus the error term) is the key: it splits each addition into a "hi" part stored back as x and a "lo" error term stored as a new partial. At the end, the partials list is collapsed right-to-left using the same careful rounding to produce a correctly-rounded final result.

The result is identical to computing the sum in arbitrary precision and then rounding once: IEEE 754 guarantees this if the partial list is maintained correctly. The list grows dynamically but is pre-allocated to 32 slots (NUM_PARTIALS) to avoid allocation for short iterables.

factorial divide-and-conquer (lines 900 to 1200)

cpython 3.14 @ ab2d84fe1023/Modules/mathmodule.c#L900-1200 ↗

math.factorial(n) raises ValueError for negative inputs and OverflowError for inputs that do not fit in a C long. For small n (up to ~20) a simple loop suffices. For large n, CPython uses a split-radix recursive algorithm:

static PyObject *
factorial_partial_product(unsigned long start, unsigned long stop,
                          unsigned long step)
{
    unsigned long midpoint, num_operands;
    PyObject *left = NULL, *right = NULL, *result = NULL;
    ...
    num_operands = (stop - start) / step;
    if (num_operands <= 8) {
        /* Base case: multiply directly */
        result = PyLong_FromUnsignedLong(start);
        for (u = start + step; u < stop; u += step) {
            factor = PyLong_FromUnsignedLong(u);
            tmp = PyNumber_Multiply(result, factor);
            Py_DECREF(result); Py_DECREF(factor);
            result = tmp;
        }
        return result;
    }
    midpoint = start + step * (num_operands / 2);
    left = factorial_partial_product(start, midpoint, step);
    right = factorial_partial_product(midpoint, stop, step);
    result = PyNumber_Multiply(left, right);
    ...
    return result;
}

Multiplying two numbers of similar size is faster than repeated multiplication of one growing accumulator against small factors, because Python's big-integer multiply uses Karatsuba (O(n^1.585)) and the balanced split keeps both operands near the same bit-length throughout.

math_comb and math_perm use the multiplicative formula n * (n-1) * ... * (n-k+1) / k! with GCD reduction at each step to keep intermediate values small and avoid a final big-integer division.

math.hypot multi-argument compensated norm (lines 2200 to 2700)

cpython 3.14 @ ab2d84fe1023/Modules/mathmodule.c#L2200-2700 ↗

Since Python 3.8, math.hypot accepts any number of arguments and computes the Euclidean norm sqrt(sum(x**2 for x in args)) without intermediate overflow or underflow.

static double
vector_norm(Py_ssize_t n, double *vec, double max, int found_nan)
{
    double x, scale, oldcmax, csum = 1.0, frac = 0.0;
    ...
    if (Py_IS_INFINITY(max)) return max;
    if (found_nan) return Py_NAN;
    if (max == 0.0 || n <= 1) return max;
    /* Scale all inputs by 1/max so that max^2 == 1, preventing overflow */
    scale = 1.0 / max;
    for (i = 0; i < n; i++) {
        x = vec[i];
        x *= scale;      /* x = x_i / max */
        x *= x;          /* x = (x_i/max)^2 */
        csum += x;       /* compensated accumulation */
        ...
    }
    return max * sqrt(csum);
}

The "scale by max" trick ensures (x_i / max)^2 is always in [0, 1], so no intermediate value exceeds 1. The final result is max * sqrt(csum), which recovers the true magnitude. A separate compensated accumulation (similar to Kahan summation) is used for csum so that the squared terms do not lose low bits.

math.floor / math.ceil / math.trunc dispatch (lines 1500 to 1800)

cpython 3.14 @ ab2d84fe1023/Modules/mathmodule.c#L1500-1800 ↗

static PyObject *
math_floor(PyObject *module, PyObject *number)
{
    PyObject *method;
    ...
    method = _PyObject_LookupSpecial(number, &_Py_ID(__floor__));
    if (method != NULL)
        return PyObject_CallNoArgs(method);
    /* Fall back to C floor() for plain float */
    ...
    x = PyFloat_AsDouble(number);
    ...
    return PyLong_FromDouble(floor(x));
}

For Decimal, Fraction, and any other type that defines __floor__, __ceil__, or __trunc__, the method is called directly. This is the mechanism that lets math.floor(Decimal('3.7')) return Decimal('3') rather than int(3). For bare float, the libm function is called and the result is converted to int (CPython's math.floor always returns int for float inputs; it is the slot that may return the original type).

gopy mirror

module/math/module.go (pending). The libm-delegation functions are straightforward wrappers around Go's math package. Fsum requires a careful port of the Shewchuk partial-sum loop; Go's float64 is IEEE 754 double, so the algorithm translates directly. Factorial will use math/big.Int for inputs above the native integer threshold.

CPython 3.14 changes

math.hypot gained the multi-argument form in 3.8. math.prod was added in 3.8. math.comb and math.perm were added in 3.8. math.isclose was added in 3.5. math.gcd accepted multiple arguments starting in 3.9, and math.lcm was added in 3.9. The fsum implementation has been stable since 2.6. In 3.12, math.factorial began raising TypeError (not ValueError) for non-integer float inputs to align with the rest of the numeric tower.