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

[skirpichev]

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/

[hpkfft]

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)

[skirpichev]

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;
 }