Skip to content

Don't rely on errno values, coming from platform's math functions #156145

Open
Open
Don't rely on errno values, coming from platform's math functions#156145
Labels
extension-modulesC modules in the Modules dirinterpreter-core(Objects, Python, Grammar, and Parser dirs)type-featureA feature request or enhancement

Description

@skirpichev
skirpichev
opened on Aug 21, 2026

Feature or enhancement

Proposal:

In the cmath module most functions don't check errno values. See cosh() example below, which set nonzero errno in some cases, but don't read errno (e.g. to detect overflow in libm's cosh/sinh()). Here we essentially assume, that platform implements Annex F of the C standard:

cpython/Modules/cmathmodule.c

Lines 441 to 490 in 04242c0

static Py_complex
cmath_cosh_impl(PyObject *module, Py_complex z)
/*[clinic end generated code: output=2e969047da601bdb input=d6b66339e9cc332b]*/
{
Py_complex r;
double x_minus_one;
/* special treatment for cosh(+/-inf + iy) if y is not a NaN */
if (!isfinite(z.real) || !isfinite(z.imag)) {
if (isinf(z.real) && isfinite(z.imag) &&
(z.imag != 0.)) {
if (z.real > 0) {
r.real = copysign(INF, cos(z.imag));
r.imag = copysign(INF, sin(z.imag));
}
else {
r.real = copysign(INF, cos(z.imag));
r.imag = -copysign(INF, sin(z.imag));
}
}
else {
r = cosh_special_values[special_type(z.real)]
[special_type(z.imag)];
}
/* need to set errno = EDOM if y is +/- infinity and x is not
a NaN */
if (isinf(z.imag) && !isnan(z.real))
errno = EDOM;
else
errno = 0;
return r;
}
if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
/* deal correctly with cases where cosh(z.real) overflows but
cosh(z) does not. */
x_minus_one = z.real - copysign(1., z.real);
r.real = cos(z.imag) * cosh(x_minus_one) * Py_MATH_E;
r.imag = sin(z.imag) * sinh(x_minus_one) * Py_MATH_E;
} else {
r.real = cos(z.imag) * cosh(z.real);
r.imag = sin(z.imag) * sinh(z.real);
}
/* detect overflow, and set errno accordingly */
if (isinf(r.real) || isinf(r.imag))
errno = ERANGE;
else
errno = 0;
return r;
}

What if we extend this approach to the rest of CPython's floating-point arithmetics? Now this seems natural, as we require that platforms C double has IEEE 754 binary64 format. On practice, that means much more, i.e. that platform supports Annex F (modulo bugs). And we can utilize that, assuming that for invalid input - NaN's produced, for overflows - infinities.

For example, instead of current code of float_pow() that handle finite input

cpython/Objects/floatobject.c

Lines 788 to 801 in 04242c0

errno = 0;
ix = pow(iv, iw);
_Py_ADJUST_ERANGE1(ix);
if (negate_result)
ix = -ix;
if (errno != 0) {
/* We don't expect any errno value other than ERANGE, but
* the range of libm bugs appears unbounded.
*/
PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
PyExc_ValueError);
return NULL;
}

we could use something much more simple:

    ix = pow(iv, iw);
    if (negate_result)
        ix = -ix;
    if (isinf(ix)) {
        PyErr_SetString(PyExc_OverflowError,
                        "float exponentiation result out of range");
        return NULL;
    }

(_Py_ADJUST_ERANGE1 and _Py_ADJUST_ERANGE2 helpers will be removed.)

See recent d.p.o thread for illustration of subtle issues, that can be introduced trying to fix errno, coming from platforms functions. We also have occurring bugreports, when platform libm set errno wrongly, e.g. #153144.

Disclaimer: this issue filled, based on my humble understanding of the proposal by @mdickinson. It seems, there are volunteers to work on this (CC @hpkfft), hence issue was opened.

Has this already been discussed elsewhere?

I have already discussed this feature proposal on Discourse

Links to previous discussion of this feature:

https://discuss.python.org/t/108539/

Linked PRs

Activity

  1. skirpichev
    added
    type-featureA feature request or enhancement
    interpreter-core(Objects, Python, Grammar, and Parser dirs)
    extension-modulesC modules in the Modules dir
    on Aug 21, 2026

  2. hpkfft commented on Aug 21, 2026

    @hpkfft
    hpkfft
    on Aug 21, 2026
    Contributor

    Thank you, Sergey. Note that errno is thread-local storage which is defined in a different shared library (i.e., not libpython3.16.so), so there is a performance cost to reading/writing it. I've started exploring avoiding its use in compleobject.c. (This is more than just not relying on the value coming from the platform's math library--it's avoiding errno as a communication channel between Python functions in the same translation unit.) Early results for complex abs(z) and result = z + 1.0 on an Intel Xeon using GCC 14.2.0 are as follows:

    Benchmark main patch
    abs(z) 29.5 ns 27.6 ns: 1.07x faster
    inc(z) 16.4 ns 15.0 ns: 1.09x faster
    pyperf benchmark code 
    #!/usr/bin/env ./python
import time
import pyperf

z = complex(3.140625, 1.0)

def bench_inc(loops, z):
    range_it = range(loops)
    t0 = time.perf_counter()
    for _ in range_it:
        result = z +  1.0
        result = z +  2.0
        result = z +  3.0
        result = z +  4.0
        result = z +  5.0
        result = z +  6.0
        result = z +  7.0
        result = z +  8.0
        result = z +  9.0
        result = z + 10.0
        result = z + 11.0
        result = z + 12.0
        result = z + 13.0
        result = z + 14.0
        result = z + 15.0
        result = z + 16.0
        result = z + 17.0
        result = z + 18.0
        result = z + 19.0
        result = z + 20.0
    return time.perf_counter() - t0

runner = pyperf.Runner(values=10, processes=24)
runner.bench_func("abs(z)", abs, z)
runner.bench_time_func("inc(z)", bench_inc, z, inner_loops=20)
  3. hpkfft
    added a commit that references this issue on Aug 21, 2026

    Code as benchmarked for comment in Issue python#156145

    eb6f38e

  4. skirpichev commented on Aug 21, 2026

    @skirpichev
    skirpichev
    on Aug 21, 2026
    MemberAuthor

    Note that errno is thread-local storage which is defined in a different shared library (i.e., not libpython3.16.so), so there is a performance cost to reading/writing it.

    This will be irrelevant in some cases. For example, patch below replaces errno using in the cmath module by a global variable (err). That's incorrect. Valid version would use module state, etc. But difference already barely noticeable:

    Mean +- std dev: [ref] 556 ns +- 2 ns -> [patch] 539 ns +- 5 ns: 1.03x faster
    Details 
    import pyperf
from cmath import asinh

z = complex(3.140625, 1.0)

runner = pyperf.Runner()
runner.bench_func("asinh(z)", asinh, z)
    diff --git a/Modules/cmathmodule.c b/Modules/cmathmodule.c
index f6e1475b00e..53545586c5a 100644
--- a/Modules/cmathmodule.c
+++ b/Modules/cmathmodule.c
@@ -16,6 +16,8 @@
 /* For _Py_log1p with workarounds for buggy handling of zeros. */
 #include "_math.h"
 
+static int err = 0;
+
 #include "clinic/cmathmodule.c.h"
 /*[clinic input]
 module cmath
@@ -25,7 +27,7 @@ module cmath
 /*[python input]
 class Py_complex_protected_converter(Py_complex_converter):
     def modify(self):
-        return 'errno = 0;'
+        return 'err = 0;'
 
 
 class Py_complex_protected_return_converter(CReturnConverter):
@@ -34,11 +36,11 @@ class Py_complex_protected_return_converter(CReturnConverter):
     def render(self, function, data):
         self.declare(data)
         data.return_conversion.append("""
-if (errno == EDOM) {
+if (err == EDOM) {
     PyErr_SetString(PyExc_ValueError, "math domain error");
     goto exit;
 }
-else if (errno == ERANGE) {
+else if (err == ERANGE) {
     PyErr_SetString(PyExc_OverflowError, "math range error");
     goto exit;
 }
@@ -47,7 +49,7 @@ else {
 }
 """.strip())
 [python start generated code]*/
-/*[python end generated code: output=da39a3ee5e6b4b0d input=8b27adb674c08321]*/
+/*[python end generated code: output=da39a3ee5e6b4b0d input=b4c607d4f07dc83c]*/
 
 #if (FLT_RADIX != 2 && FLT_RADIX != 16)
 #error "Modules/cmathmodule.c expects FLT_RADIX to be 2 or 16"
@@ -141,7 +143,7 @@ special_type(double d)
 
 #define SPECIAL_VALUE(z, table)                       \
     if (!isfinite((z).real) || !isfinite((z).imag)) { \
-        errno = 0;                                    \
+        err = 0;                                      \
         return table[special_type((z).real)]          \
                     [special_type((z).imag)];         \
     }
@@ -156,10 +158,10 @@ special_type(double d)
 
 /* First, the C functions that do the real work.  Each of the c_*
    functions computes and returns the C99 Annex G recommended result
-   and also sets errno as follows: errno = 0 if no floating-point
-   exception is associated with the result; errno = EDOM if C99 Annex
+   and also sets err as follows: err = 0 if no floating-point
+   exception is associated with the result; err = EDOM if C99 Annex
    G recommends raising divide-by-zero or invalid for this result; and
-   errno = ERANGE where the overflow floating-point signal should be
+   err = ERANGE where the overflow floating-point signal should be
    raised.
 */
 
@@ -204,7 +206,7 @@ cmath_acos_impl(PyObject *module, Py_complex z)
         r.real = 2.*atan2(s1.real, s2.real);
         r.imag = asinh(s2.real*s1.imag - s2.imag*s1.real);
     }
-    errno = 0;
+    err = 0;
     return r;
 }
 
@@ -247,7 +249,7 @@ cmath_acosh_impl(PyObject *module, Py_complex z)
         r.real = asinh(s1.real*s2.real + s1.imag*s2.imag);
         r.imag = 2.*atan2(s1.imag, s2.real);
     }
-    errno = 0;
+    err = 0;
     return r;
 }
 
@@ -315,7 +317,7 @@ cmath_asinh_impl(PyObject *module, Py_complex z)
         r.real = asinh(s1.real*s2.imag-s2.real*s1.imag);
         r.imag = atan2(z.imag, s1.real*s2.real-s1.imag*s2.imag);
     }
-    errno = 0;
+    err = 0;
     return r;
 }
 
@@ -381,22 +383,22 @@ cmath_atanh_impl(PyObject *module, Py_complex z)
         h = hypot(z.real/2., z.imag/2.);  /* safe from overflow */
         r.real = z.real/4./h/h;
         r.imag = copysign(Py_MATH_PI/2., z.imag);
-        errno = 0;
+        err = 0;
     } else if (z.real == 1. && ay < CM_SQRT_DBL_MIN) {
         /* C99 standard says:  atanh(1+/-0.) should be inf +/- 0i */
         if (ay == 0.) {
             r.real = INF;
             r.imag = z.imag;
-            errno = EDOM;
+            err = EDOM;
         } else {
             r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.)));
             r.imag = copysign(atan2(2., -ay)/2, z.imag);
-            errno = 0;
+            err = 0;
         }
     } else {
         r.real = m_log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.;
         r.imag = -atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.;
-        errno = 0;
+        err = 0;
     }
     return r;
 }
@@ -462,12 +464,12 @@ cmath_cosh_impl(PyObject *module, Py_complex z)
             r = cosh_special_values[special_type(z.real)]
                                    [special_type(z.imag)];
         }
-        /* need to set errno = EDOM if y is +/- infinity and x is not
+        /* need to set err = EDOM if y is +/- infinity and x is not
            a NaN */
         if (isinf(z.imag) && !isnan(z.real))
-            errno = EDOM;
+            err = EDOM;
         else
-            errno = 0;
+            err = 0;
         return r;
     }
 
@@ -483,9 +485,9 @@ cmath_cosh_impl(PyObject *module, Py_complex z)
     }
     /* detect overflow, and set errno accordingly */
     if (isinf(r.real) || isinf(r.imag))
-        errno = ERANGE;
+        err = ERANGE;
     else
-        errno = 0;
+        err = 0;
     return r;
 }
 
@@ -531,14 +533,14 @@ cmath_exp_impl(PyObject *module, Py_complex z)
             r = exp_special_values[special_type(z.real)]
                                   [special_type(z.imag)];
         }
-        /* need to set errno = EDOM if y is +/- infinity and x is not
+        /* need to set err = EDOM if y is +/- infinity and x is not
            a NaN and not -infinity */
         if (isinf(z.imag) &&
             (isfinite(z.real) ||
              (isinf(z.real) && z.real > 0)))
-            errno = EDOM;
+            err = EDOM;
         else
-            errno = 0;
+            err = 0;
         return r;
     }
 
@@ -551,11 +553,11 @@ cmath_exp_impl(PyObject *module, Py_complex z)
         r.real = l*cos(z.imag);
         r.imag = l*sin(z.imag);
     }
-    /* detect overflow, and set errno accordingly */
+    /* detect overflow, and set err accordingly */
     if (isinf(r.real) || isinf(r.imag))
-        errno = ERANGE;
+        err = ERANGE;
     else
-        errno = 0;
+        err = 0;
     return r;
 }
 
@@ -595,7 +597,7 @@ c_log(Py_complex z)
        (4) z = 0.  The simplest thing to do here is to call the
        floating-point log with an argument of 0, and let its behaviour
        (returning -infinity, signaling a floating-point exception, setting
-       errno, or whatever) determine that of c_log.  So the usual formula
+       err, or whatever) determine that of c_log.  So the usual formula
        is fine here.
 
      */
@@ -620,7 +622,7 @@ c_log(Py_complex z)
             /* log(+/-0. +/- 0i) */
             r.real = -INF;
             r.imag = atan2(z.imag, z.real);
-            errno = EDOM;
+            err = EDOM;
             return r;
         }
     } else {
@@ -634,7 +636,7 @@ c_log(Py_complex z)
         }
     }
     r.imag = atan2(z.imag, z.real);
-    errno = 0;
+    err = 0;
     return r;
 }
@@ -650,13 +652,10 @@ cmath_log10_impl(PyObject *module, Py_complex z)
 /*[clinic end generated code: output=2922779a7c38cbe1 input=cff5644f73c1519c]*/
 {
     Py_complex r;
-    int errno_save;
 
     r = c_log(z);
-    errno_save = errno; /* just in case the divisions affect errno */
     r.real = r.real / M_LN10;
     r.imag = r.imag / M_LN10;
-    errno = errno_save;
     return r;
 }
 
@@ -724,12 +723,12 @@ cmath_sinh_impl(PyObject *module, Py_complex z)
             r = sinh_special_values[special_type(z.real)]
                                    [special_type(z.imag)];
         }
-        /* need to set errno = EDOM if y is +/- infinity and x is not
+        /* need to set err = EDOM if y is +/- infinity and x is not
            a NaN */
         if (isinf(z.imag) && !isnan(z.real))
-            errno = EDOM;
+            err = EDOM;
         else
-            errno = 0;
+            err = 0;
         return r;
     }
 
@@ -743,9 +742,9 @@ cmath_sinh_impl(PyObject *module, Py_complex z)
     }
     /* detect overflow, and set errno accordingly */
     if (isinf(r.real) || isinf(r.imag))
-        errno = ERANGE;
+        err = ERANGE;
     else
-        errno = 0;
+        err = 0;
     return r;
 }
@@ -830,7 +829,7 @@ cmath_sqrt_impl(PyObject *module, Py_complex z)
         r.real = d;
         r.imag = copysign(s, z.imag);
     }
-    errno = 0;
+    err = 0;
     return r;
 }
 
@@ -912,12 +911,12 @@ cmath_tanh_impl(PyObject *module, Py_complex z)
             r = tanh_special_values[special_type(z.real)]
                                    [special_type(z.imag)];
         }
-        /* need to set errno = EDOM if z.imag is +/-infinity and
+        /* need to set err = EDOM if z.imag is +/-infinity and
            z.real is finite */
         if (isinf(z.imag) && isfinite(z.real))
-            errno = EDOM;
+            err = EDOM;
         else
-            errno = 0;
+            err = 0;
         return r;
     }
 
@@ -934,7 +933,7 @@ cmath_tanh_impl(PyObject *module, Py_complex z)
         r.real = tx*(1.+ty*ty)/denom;
         r.imag = ((ty/denom)*cx)*cx;
     }
-    errno = 0;
+    err = 0;
     return r;
 }
 
@@ -958,7 +957,7 @@ cmath_log_impl(PyObject *module, Py_complex x, PyObject *y_obj)
 {
     Py_complex y;
 
-    errno = 0;
+    err = 0;
     x = c_log(x);
     if (y_obj != NULL) {
         y = PyComplex_AsCComplex(y_obj);
@@ -968,8 +967,10 @@ cmath_log_impl(PyObject *module, Py_complex x, PyObject *y_obj)
         y = c_log(y);
         x = _Py_c_quot(x, y);
     }
-    if (errno != 0)
+    if (err != 0) {
+        errno = err;
         return math_error();
+    }
     return PyComplex_FromCComplex(x);
 }
 
@@ -1074,7 +1075,7 @@ cmath_rect_impl(PyObject *module, double r, double phi)
 /*[clinic end generated code: output=74ff3d17585f3388 input=50e60c5d28c834e6]*/
 {
     Py_complex z;
-    errno = 0;
+    err = 0;
 
     /* deal with special values */
     if (!isfinite(r) || !isfinite(phi)) {
@@ -1099,21 +1100,21 @@ cmath_rect_impl(PyObject *module, double r, double phi)
         /* need to set errno = EDOM if r is a nonzero number and phi
            is infinite */
         if (r != 0. && !isnan(r) && isinf(phi))
-            errno = EDOM;
+            err = EDOM;
         else
-            errno = 0;
+            err = 0;
     }
     else if (phi == 0.0) {
         /* Workaround for buggy results with phi=-0.0 on OS X 10.8.  See
            bugs.python.org/issue18513. */
         z.real = r;
         z.imag = r * phi;
-        errno = 0;
+        err = 0;
     }
     else {
         z.real = r * cos(phi);
         z.imag = r * sin(phi);
-        errno = 0;
+        err = 0;
     }
     return z;
 }

    Edit:
    Unfortunately, using floating-point environment seems to be a poor alternative to testing errno. For abs() benchmark below I got:

    Mean +- std dev: [ref] 259 ns +- 3 ns -> [patch] 388 ns +- 2 ns: 1.50x slower
    Details 
    diff --git a/Objects/complexobject.c b/Objects/complexobject.c
index 3612c2699a5..01078a121d7 100644
--- a/Objects/complexobject.c
+++ b/Objects/complexobject.c
@@ -5,6 +5,7 @@
 /* Submitted by Jim Hugunin */
 
 #include "Python.h"
+#include "fenv.h"
 #include "pycore_call.h"          // _PyObject_CallNoArgs()
 #include "pycore_complexobject.h" // _PyComplex_FormatAdvancedWriter()
 #include "pycore_floatobject.h"   // _Py_convert_int_to_double()
@@ -796,8 +797,11 @@ static PyObject *
 complex_abs(PyObject *op)
 {
     PyComplexObject *v = _PyComplexObject_CAST(op);
-    double result = _Py_c_abs(v->cval);
-    if (errno == ERANGE) {
+    double result;
+
+    feclearexcept(FE_ALL_EXCEPT);
+    result = hypot(v->cval.real, v->cval.imag);
+    if (fetestexcept(FE_OVERFLOW)) {
         PyErr_SetString(PyExc_OverflowError,
                         "absolute value too large");
         return NULL;
    # bench.py
import pyperf

z = complex(3.140625, 1.0)

runner = pyperf.Runner()
runner.bench_func("abs(z)", abs, z)

    Note, that in some cases we will need to examine either errno or floating-point exceptions. For example, lgamma() of nonpositive integer.

  5. skirpichev commented on Aug 29, 2026

    @skirpichev
    skirpichev
    on Aug 29, 2026
    MemberAuthor

    @hpkfft, let me know if you are still interested in working on this.

    No need to do everything in one shot. You can start with changing floatobject.c, for example.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    extension-modulesC modules in the Modules dirinterpreter-core(Objects, Python, Grammar, and Parser dirs)type-featureA feature request or enhancement

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions