diff --git a/Include/internal/pycore_pymath.h b/Include/internal/pycore_pymath.h index 532c5ceafb56395..f7c244a70cf9d62 100644 --- a/Include/internal/pycore_pymath.h +++ b/Include/internal/pycore_pymath.h @@ -8,55 +8,6 @@ extern "C" { # error "this header requires Py_BUILD_CORE define" #endif - -/* _Py_ADJUST_ERANGE1(x) - * _Py_ADJUST_ERANGE2(x, y) - * Set errno to 0 before calling a libm function, and invoke one of these - * macros after, passing the function result(s) (_Py_ADJUST_ERANGE2 is useful - * for functions returning complex results). This makes two kinds of - * adjustments to errno: (A) If it looks like the platform libm set - * errno=ERANGE due to underflow, clear errno. (B) If it looks like the - * platform libm overflowed but didn't set errno, force errno to ERANGE. In - * effect, we're trying to force a useful implementation of C89 errno - * behavior. - * Caution: - * This isn't reliable. C99 no longer requires libm to set errno under - * any exceptional condition, but does require +- HUGE_VAL return - * values on overflow. A 754 box *probably* maps HUGE_VAL to a - * double infinity, and we're cool if that's so, unless the input - * was an infinity and an infinity is the expected result. A C89 - * system sets errno to ERANGE, so we check for that too. We're - * out of luck if a C99 754 box doesn't map HUGE_VAL to +Inf, or - * if the returned result is a NaN, or if a C89 box returns HUGE_VAL - * in non-overflow cases. - */ -static inline void _Py_ADJUST_ERANGE1(double x) -{ - if (errno == 0) { - if (x == INFINITY || x == -INFINITY) { - errno = ERANGE; - } - } - else if (errno == ERANGE && x == 0.0) { - errno = 0; - } -} - -static inline void _Py_ADJUST_ERANGE2(double x, double y) -{ - if (x == INFINITY || x == -INFINITY || - y == INFINITY || y == -INFINITY) - { - if (errno == 0) { - errno = ERANGE; - } - } - else if (errno == ERANGE) { - errno = 0; - } -} - - //--- HAVE_PY_SET_53BIT_PRECISION macro ------------------------------------ // // The functions _Py_dg_strtod() and _Py_dg_dtoa() in Python/dtoa.c (which are diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index bb307191dffcc14..4ac9ba5def4544a 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -84,6 +84,10 @@ def assertClose(self, x, y, eps=1e-9): # check that relative difference < eps self.assertTrue(abs(x-y)/abs(y) < eps) + def assertSameSign(self, x, y): + if copysign(1., x) != copysign(1., y): + self.fail(f'{x!r} and {y!r} have different signs') + def check_div(self, x, y): """Compute complex z=x*y, and check that z/x==y and z/y==x.""" z = x * y @@ -446,6 +450,63 @@ def test_pow_with_small_integer_exponents(self): self.assertEqual(str(float_pow), str(int_pow)) self.assertEqual(str(complex_pow), str(int_pow)) + # Check that complex numbers with special components + # are correctly handled. + values = [complex(x, y) + for x in [5, -5, +0.0, -0.0, INF, -INF, NAN] + for y in [12, -12, +0.0, -0.0, INF, -INF, NAN]] + for c in values: + with self.subTest(value=c): + self.assertComplexesAreIdentical(c**0, complex(1, +0.0)) + self.assertComplexesAreIdentical(c**1, c) + self.assertComplexesAreIdentical(c**2, c*c) + self.assertComplexesAreIdentical(c**3, c*(c*c)) + self.assertComplexesAreIdentical(c**3, (c*c)*c) + if not c: + continue + for n in range(1, 9): + with self.subTest(exponent=-n): + self.assertComplexesAreIdentical(c**-n, 1/(c**n)) + + # Special cases for complex division. + for x in [+2, -2]: + for y in [+0.0, -0.0]: + c = complex(x, y) + with self.subTest(value=c): + self.assertComplexesAreIdentical(c**-1, complex(1/x, -y)) + c = complex(y, x) + with self.subTest(value=c): + self.assertComplexesAreIdentical(c**-1, complex(y, -1/x)) + for x in [+INF, -INF]: + for y in [+1, -1]: + c = complex(x, y) + with self.subTest(value=c): + self.assertComplexesAreIdentical(c**-1, complex(1/x, -0.0*y)) + self.assertComplexesAreIdentical(c**-2, complex(0.0, -y/x)) + c = complex(y, x) + with self.subTest(value=c): + self.assertComplexesAreIdentical(c**-1, complex(+0.0*y, -1/x)) + self.assertComplexesAreIdentical(c**-2, complex(-0.0, -y/x)) + + # Test that zeros have the same sign as small non-zero values. + eps = 1e-11 + pairs = [(complex(x, y), complex(x, copysign(0.0, y))) + for x in [+1, -1] for y in [+eps, -eps]] + pairs += [(complex(y, x), complex(copysign(0.0, y), x)) + for x in [+1, -1] for y in [+eps, -eps]] + for c1, c2 in pairs: + for n in exponents: + with self.subTest(value=c1, exponent=n): + r1 = c1**n + r2 = c2**n + self.assertClose(r1, r2) + self.assertSameSign(r1.real, r2.real) + self.assertSameSign(r1.imag, r2.imag) + self.assertNotEqual(r1.real, 0.0) + if n != 0: + self.assertNotEqual(r1.imag, 0.0) + self.assertTrue(r2.real == 0.0 or r2.imag == 0.0) + def test_boolcontext(self): for i in range(100): self.assertTrue(complex(random() + 1e-6, random() + 1e-6)) diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 3612c2699a557db..2e649f6d7849772 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -11,7 +11,6 @@ #include "pycore_freelist.h" // _Py_FREELIST_FREE(), _Py_FREELIST_POP() #include "pycore_long.h" // _PyLong_GetZero() #include "pycore_object.h" // _PyObject_Init() -#include "pycore_pymath.h" // _Py_ADJUST_ERANGE2() #define _PyComplexObject_CAST(op) ((PyComplexObject *)(op)) @@ -26,8 +25,6 @@ class complex "PyComplexObject *" "&PyComplex_Type" /* elementary operations on complex numbers */ -static Py_complex c_1 = {1., 0.}; - Py_complex _Py_c_sum(Py_complex a, Py_complex b) { @@ -166,8 +163,8 @@ _Py_rc_prod(double a, Py_complex b) #ifdef _M_ARM64 #pragma optimize("", off) #endif -Py_complex -_Py_c_quot(Py_complex a, Py_complex b) +static Py_complex +c_quot(Py_complex a, Py_complex b) { /****************************************************************** This was the original algorithm. It's grossly prone to spurious @@ -197,16 +194,11 @@ _Py_c_quot(Py_complex a, Py_complex b) if (abs_breal >= abs_bimag) { /* divide tops and bottom by b.real */ - if (abs_breal == 0.0) { - errno = EDOM; - r.real = r.imag = 0.0; - } - else { - const double ratio = b.imag / b.real; - const double denom = b.real + b.imag * ratio; - r.real = (a.real + a.imag * ratio) / denom; - r.imag = (a.imag - a.real * ratio) / denom; - } + /* If b.real == b.imag == 0.0, r.real and r.imag will be NaN. */ + const double ratio = b.imag / b.real; + const double denom = b.real + b.imag * ratio; + r.real = (a.real + a.imag * ratio) / denom; + r.imag = (a.imag - a.real * ratio) / denom; } else if (abs_bimag >= abs_breal) { /* divide tops and bottom by b.imag */ @@ -246,21 +238,30 @@ _Py_c_quot(Py_complex a, Py_complex b) } Py_complex -_Py_cr_quot(Py_complex a, double b) +_Py_c_quot(Py_complex a, Py_complex b) { - Py_complex r = a; - if (b) { - r.real /= b; - r.imag /= b; + Py_complex r; + if (b.real == 0.0 && b.imag == 0.0) { + errno = EDOM; + r.real = r.imag = 0.0; // As documented, so shall it be done. } else { - errno = EDOM; - r.real = r.imag = 0.0; + r = c_quot(a, b); } return r; } -/* an equivalent of _Py_c_quot() function, when 1st argument is real */ +Py_complex +_Py_cr_quot(Py_complex a, double b) +{ + Py_complex r = a; + // May overflow or may result in NaN. + r.real /= b; + r.imag /= b; + return r; +} + +/* an equivalent of c_quot() function, when 1st argument is real */ Py_complex _Py_rc_quot(double a, Py_complex b) { @@ -269,16 +270,10 @@ _Py_rc_quot(double a, Py_complex b) const double abs_bimag = b.imag < 0 ? -b.imag : b.imag; if (abs_breal >= abs_bimag) { - if (abs_breal == 0.0) { - errno = EDOM; - r.real = r.imag = 0.0; - } - else { - const double ratio = b.imag / b.real; - const double denom = b.real + b.imag * ratio; - r.real = a / denom; - r.imag = (-a * ratio) / denom; - } + const double ratio = b.imag / b.real; + const double denom = b.real + b.imag * ratio; + r.real = a / denom; + r.imag = (-a * ratio) / denom; } else if (abs_bimag >= abs_breal) { const double ratio = b.real / b.imag; @@ -306,8 +301,28 @@ _Py_rc_quot(double a, Py_complex b) #pragma optimize("", on) #endif -Py_complex -_Py_c_pow(Py_complex a, Py_complex b) +static Py_complex +c_powi(Py_complex x, long n) +{ + Py_complex r = x; + long absn = (n < 0) ? -n : n; + + assert(absn > 0); + if (--absn & 0x1) + r = _Py_c_prod(r, x); + while (absn >>= 1) { + x = _Py_c_prod(x, x); + if (absn & 0x1) + r = _Py_c_prod(r, x); + } + + if (n < 0) + r = _Py_rc_quot(1.0, r); + return r; +} + +static Py_complex +c_pow(Py_complex a, Py_complex b, int *e) { Py_complex r; double vabs,len,at,phase; @@ -315,11 +330,20 @@ _Py_c_pow(Py_complex a, Py_complex b) r.real = 1.; r.imag = 0.; } - else if (a.real == 0. && a.imag == 0.) { - if (b.imag != 0. || b.real < 0.) - errno = EDOM; + else if (a.real == 0. && a.imag == 0. && (b.imag != 0. || b.real < 0.)) { r.real = 0.; r.imag = 0.; + *e = EDOM; + } + else if (b.imag == 0.0 && b.real == floor(b.real) + && fabs(b.real) <= 100.0) { + // The exponent is a small integer value, so use a faster and + // more accurate algorithm. + r = c_powi(a, (long)b.real); + if ((isinf(r.real) || isinf(r.imag)) + && isfinite(a.real) && isfinite(a.imag)) { + *e = ERANGE; + } } else { vabs = hypot(a.real,a.imag); @@ -332,42 +356,30 @@ _Py_c_pow(Py_complex a, Py_complex b) } r.real = len*cos(phase); r.imag = len*sin(phase); - - _Py_ADJUST_ERANGE2(r.real, r.imag); + if ((isinf(r.real) || isinf(r.imag)) + && isfinite(a.real) && isfinite(a.imag) + && isfinite(b.real) && isfinite(b.imag)) { + *e = ERANGE; + } } return r; } -static Py_complex -c_powu(Py_complex x, long n) +Py_complex +_Py_c_pow(Py_complex a, Py_complex b) { - Py_complex r, p; - long mask = 1; - r = c_1; - p = x; - while (mask > 0 && n >= mask) { - if (n & mask) - r = _Py_c_prod(r,p); - mask <<= 1; - p = _Py_c_prod(p,p); - } + Py_complex r; + int e = 0; + int saved_errno = errno; + r = c_pow(a, b, &e); + // Sets errno = EDOM on invalid, ERANGE on overflow, else unchanged. + errno = (e) ? e : saved_errno; return r; } -static Py_complex -c_powi(Py_complex x, long n) -{ - if (n > 0) - return c_powu(x,n); - else - return _Py_c_quot(c_1, c_powu(x,-n)); - -} - -double -_Py_c_abs(Py_complex z) +static double +c_abs(Py_complex z, int *e) { - /* sets errno = ERANGE on overflow; otherwise errno = 0 */ double result; if (!isfinite(z.real) || !isfinite(z.imag)) { @@ -376,12 +388,10 @@ _Py_c_abs(Py_complex z) NaN. */ if (isinf(z.real)) { result = fabs(z.real); - errno = 0; return result; } if (isinf(z.imag)) { result = fabs(z.imag); - errno = 0; return result; } /* either the real or imaginary part is a NaN, @@ -390,12 +400,22 @@ _Py_c_abs(Py_complex z) } result = hypot(z.real, z.imag); if (!isfinite(result)) - errno = ERANGE; - else - errno = 0; + *e = ERANGE; return result; } +double +_Py_c_abs(Py_complex z) +{ + double r; + int e = 0; + int saved_errno = errno; + r = c_abs(z, &e); + // Sets errno = ERANGE on overflow, else errno is unchanged. + errno = (e) ? e : saved_errno; + return r; +} + static PyObject * complex_subtype_from_c_complex(PyTypeObject *type, Py_complex cval) { @@ -693,15 +713,17 @@ real_to_complex(PyObject **pobj, Py_complex *pc) See C11's Annex G, sections G.5.1 and G.5.2. */ -#define COMPLEX_BINOP(NAME, FUNC) \ +#define COMPLEX_BINOP(NAME, FUNC, IS_DIVISION) \ static PyObject * \ complex_##NAME(PyObject *v, PyObject *w) \ { \ Py_complex a; \ - errno = 0; \ if (PyComplex_Check(w)) { \ Py_complex b = ((PyComplexObject *)w)->cval; \ if (PyComplex_Check(v)) { \ + if (IS_DIVISION && b.real == 0.0 \ + && b.imag == 0.0) \ + goto DivByZero; \ a = ((PyComplexObject *)v)->cval; \ a = _Py_c_##FUNC(a, b); \ } \ @@ -709,6 +731,9 @@ real_to_complex(PyObject **pobj, Py_complex *pc) return v; \ } \ else { \ + if (IS_DIVISION && b.real == 0.0 \ + && b.imag == 0.0) \ + goto DivByZero; \ a = _Py_rc_##FUNC(a.real, b); \ } \ } \ @@ -721,20 +746,22 @@ real_to_complex(PyObject **pobj, Py_complex *pc) if (real_to_double(&w, &b) < 0) { \ return w; \ } \ + if (IS_DIVISION && b == 0.0) \ + goto DivByZero; \ a = _Py_cr_##FUNC(a, b); \ } \ - if (errno == EDOM) { \ - PyErr_SetString(PyExc_ZeroDivisionError, \ - "division by zero"); \ - return NULL; \ - } \ return PyComplex_FromCComplex(a); \ - } + \ + DivByZero: \ + PyErr_SetString(PyExc_ZeroDivisionError, \ + "division by zero"); \ + return NULL; \ + } -COMPLEX_BINOP(add, sum) -COMPLEX_BINOP(mul, prod) -COMPLEX_BINOP(sub, diff) -COMPLEX_BINOP(div, quot) +COMPLEX_BINOP(add, sum, 0) +COMPLEX_BINOP(mul, prod, 0) +COMPLEX_BINOP(sub, diff, 0) +COMPLEX_BINOP(div, quot, 1) static PyObject * complex_pow(PyObject *v, PyObject *w, PyObject *z) @@ -748,25 +775,17 @@ complex_pow(PyObject *v, PyObject *w, PyObject *z) PyErr_SetString(PyExc_ValueError, "complex modulo"); return NULL; } - errno = 0; - // Check whether the exponent has a small integer value, and if so use - // a faster and more accurate algorithm. - if (b.imag == 0.0 && b.real == floor(b.real) && fabs(b.real) <= 100.0) { - p = c_powi(a, (long)b.real); - _Py_ADJUST_ERANGE2(p.real, p.imag); - } - else { - p = _Py_c_pow(a, b); - } + int e = 0; + p = c_pow(a, b, &e); - if (errno == EDOM) { + if (e == EDOM) { PyErr_SetString(PyExc_ZeroDivisionError, "zero to a negative or complex power"); return NULL; } - else if (errno == ERANGE) { + else if (e == ERANGE) { PyErr_SetString(PyExc_OverflowError, - "complex exponentiation"); + "complex exponentiation result out of range"); return NULL; } return PyComplex_FromCComplex(p); @@ -796,8 +815,9 @@ static PyObject * complex_abs(PyObject *op) { PyComplexObject *v = _PyComplexObject_CAST(op); - double result = _Py_c_abs(v->cval); - if (errno == ERANGE) { + int e = 0; + double result = c_abs(v->cval, &e); + if (e == ERANGE) { PyErr_SetString(PyExc_OverflowError, "absolute value too large"); return NULL; diff --git a/Objects/floatobject.c b/Objects/floatobject.c index 17e6a729dcd83fc..e1ef3b937421d69 100644 --- a/Objects/floatobject.c +++ b/Objects/floatobject.c @@ -785,18 +785,13 @@ float_pow(PyObject *v, PyObject *w, PyObject *z) * positive and not equal to 1.0. We finally allow * the platform pow to step in and do the rest. */ - 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); + if (isinf(ix)) { + PyErr_SetString(PyExc_OverflowError, + "float exponentiation result out of range"); return NULL; } return PyFloat_FromDouble(ix); @@ -836,20 +831,12 @@ float_is_integer_impl(PyObject *self) /*[clinic end generated code: output=7112acf95a4d31ea input=311810d3f777e10d]*/ { double x = PyFloat_AsDouble(self); - PyObject *o; if (x == -1.0 && PyErr_Occurred()) return NULL; if (!isfinite(x)) Py_RETURN_FALSE; - errno = 0; - o = (floor(x) == x) ? Py_True : Py_False; - if (errno != 0) { - PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError : - PyExc_ValueError); - return NULL; - } - return Py_NewRef(o); + return PyBool_FromLong(floor(x) == x); } /*[clinic input]