From 850e1151eaaa9903e802e8addb7a37fbb7fd4f08 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Tue, 23 Apr 2024 18:00:55 +0300 Subject: [PATCH] gh-69639: add mixed-mode rules for complex arithmetic (float case) "Generally, mixed-mode arithmetic combining real and complex variables should be performed directly, not by first coercing the real to complex, lest the sign of zero be rendered uninformative; the same goes for combinations of pure imaginary quantities with complex variables." (c) Kahan, W: Branch cuts for complex elementary functions. This patch implements mixed-mode arithmetic rules, combining real and complex variables as specified by C standards since C99. Most C compilers implementing C99+ Annex G have only these special rules (without support for imaginary type, which is going to be deprecated in C2y). --- Doc/library/cmath.rst | 12 +- Include/internal/pycore_floatobject.h | 2 + ...4-08-03-14-02-27.gh-issue-69639.mW3iKq.rst | 2 + Objects/complexobject.c | 164 +++++++++++++----- Objects/floatobject.c | 6 +- Python/bltinmodule.c | 2 - 6 files changed, 138 insertions(+), 50 deletions(-) create mode 100644 Misc/NEWS.d/next/Core and Builtins/2024-08-03-14-02-27.gh-issue-69639.mW3iKq.rst diff --git a/Doc/library/cmath.rst b/Doc/library/cmath.rst index 381a8332f4b187..98519a11c95eb6 100644 --- a/Doc/library/cmath.rst +++ b/Doc/library/cmath.rst @@ -24,17 +24,17 @@ the function is then applied to the result of the conversion. imaginary axis we look at the sign of the real part. For example, the :func:`cmath.sqrt` function has a branch cut along the - negative real axis. An argument of ``complex(-2.0, -0.0)`` is treated as + negative real axis. An argument of ``-2-0j`` is treated as though it lies *below* the branch cut, and so gives a result on the negative imaginary axis:: - >>> cmath.sqrt(complex(-2.0, -0.0)) + >>> cmath.sqrt(-2-0j) -1.4142135623730951j - But an argument of ``complex(-2.0, 0.0)`` is treated as though it lies above + But an argument of ``-2+0j`` is treated as though it lies above the branch cut:: - >>> cmath.sqrt(complex(-2.0, 0.0)) + >>> cmath.sqrt(-2+0j) 1.4142135623730951j @@ -63,9 +63,9 @@ rectangular coordinates to polar coordinates and back. along the negative real axis. The sign of the result is the same as the sign of ``x.imag``, even when ``x.imag`` is zero:: - >>> phase(complex(-1.0, 0.0)) + >>> phase(-1+0j) 3.141592653589793 - >>> phase(complex(-1.0, -0.0)) + >>> phase(-1-0j) -3.141592653589793 diff --git a/Include/internal/pycore_floatobject.h b/Include/internal/pycore_floatobject.h index be1c6cc97720d2..fd52fe006392b2 100644 --- a/Include/internal/pycore_floatobject.h +++ b/Include/internal/pycore_floatobject.h @@ -54,6 +54,8 @@ extern PyObject* _Py_string_to_number_with_underscores( extern double _Py_parse_inf_or_nan(const char *p, char **endptr); +extern int _Py_convert_to_double(PyObject **v, double *dbl); + #ifdef __cplusplus } diff --git a/Misc/NEWS.d/next/Core and Builtins/2024-08-03-14-02-27.gh-issue-69639.mW3iKq.rst b/Misc/NEWS.d/next/Core and Builtins/2024-08-03-14-02-27.gh-issue-69639.mW3iKq.rst new file mode 100644 index 00000000000000..0581bcd8c13681 --- /dev/null +++ b/Misc/NEWS.d/next/Core and Builtins/2024-08-03-14-02-27.gh-issue-69639.mW3iKq.rst @@ -0,0 +1,2 @@ +Implement mixed-mode arithmetic rules combining real and complex variables +as specified by C standards since C99. Patch by Sergey B Kirpichev. diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 787235c63a6be1..7c8e64a96ebab4 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -8,6 +8,7 @@ #include "Python.h" #include "pycore_call.h" // _PyObject_CallNoArgs() #include "pycore_complexobject.h" // _PyComplex_FormatAdvancedWriter() +#include "pycore_floatobject.h" // _Py_convert_to_double() #include "pycore_long.h" // _PyLong_GetZero() #include "pycore_object.h" // _PyObject_Init() #include "pycore_pymath.h" // _Py_ADJUST_ERANGE2() @@ -481,75 +482,160 @@ complex_hash(PyComplexObject *v) return (obj) static int -to_complex(PyObject **pobj, Py_complex *pc) +to_float(PyObject **pobj, double *dbl) { PyObject *obj = *pobj; - pc->real = pc->imag = 0.0; - if (PyLong_Check(obj)) { - pc->real = PyLong_AsDouble(obj); - if (pc->real == -1.0 && PyErr_Occurred()) { - *pobj = NULL; - return -1; - } - return 0; - } if (PyFloat_Check(obj)) { - pc->real = PyFloat_AsDouble(obj); - return 0; + *dbl = PyFloat_AS_DOUBLE(obj); + } + else if (_Py_convert_to_double(pobj, dbl) < 0) { + return -1; } - *pobj = Py_NewRef(Py_NotImplemented); - return -1; + return 0; } +static int +to_complex(PyObject **pobj, Py_complex *pc) +{ + pc->imag = 0.0; + return to_float(pobj, &(pc->real)); +} + +/* Complex arithmetic rules implement special mixed-mode case: combining + pure-real (float's or int's) value and complex value performed directly, not + by first coercing the real value to complex. + + Lets consider the addition as an example, assuming that int's are implicitly + converted to float's. We have following rules (up to variants with changed + order of operands): + + complex(x, y) + complex(u, v) = complex(x + u, y + v) (1) + float(x) + complex(u, v) = complex(x + u, v) (2) + + Similar rules are implemented for subtraction, multiplication and division. + See C11's Annex G, sections G.5.1 and G.5.2. + */ static PyObject * complex_add(PyObject *v, PyObject *w) { - Py_complex result; - Py_complex a, b; - TO_COMPLEX(v, a); - TO_COMPLEX(w, b); - result = _Py_c_sum(a, b); - return PyComplex_FromCComplex(result); + if (PyComplex_Check(w)) { + PyObject *tmp = v; + v = w; + w = tmp; + } + + Py_complex a = ((PyComplexObject *)(v))->cval; + double b; + + if (PyComplex_Check(w)) { + Py_complex b = ((PyComplexObject *)(w))->cval; + a = _Py_c_sum(a, b); + } + else if (to_float(&w, &b) < 0) { + return w; + } + else { + a.real += b; + } + + return PyComplex_FromCComplex(a); } static PyObject * complex_sub(PyObject *v, PyObject *w) { - Py_complex result; - Py_complex a, b; - TO_COMPLEX(v, a); - TO_COMPLEX(w, b); - result = _Py_c_diff(a, b); - return PyComplex_FromCComplex(result); + Py_complex a; + + if (PyComplex_Check(w)) { + Py_complex b = ((PyComplexObject *)(w))->cval; + + if (PyComplex_Check(v)) { + a = ((PyComplexObject *)(v))->cval; + errno = 0; + a = _Py_c_diff(a, b); + } + else if (to_float(&v, &a.real) < 0) { + return v; + } + else { + a = (Py_complex) {a.real, -b.imag}; + a.real -= b.real; + } + } + else { + a = ((PyComplexObject *)(v))->cval; + double b; + + if (to_float(&w, &b) < 0) { + return w; + } + a.real -= b; + } + + return PyComplex_FromCComplex(a); } static PyObject * complex_mul(PyObject *v, PyObject *w) { - Py_complex result; - Py_complex a, b; - TO_COMPLEX(v, a); - TO_COMPLEX(w, b); - result = _Py_c_prod(a, b); - return PyComplex_FromCComplex(result); + if (PyComplex_Check(w)) { + PyObject *tmp = v; + v = w; + w = tmp; + } + + Py_complex a = ((PyComplexObject *)(v))->cval; + double b; + + if (PyComplex_Check(w)) { + Py_complex b = ((PyComplexObject *)(w))->cval; + a = _Py_c_prod(a, b); + } + else if (to_float(&w, &b) < 0) { + return w; + } + else { + a.real *= b; + a.imag *= b; + } + + return PyComplex_FromCComplex(a); } static PyObject * complex_div(PyObject *v, PyObject *w) { - Py_complex quot; - Py_complex a, b; - TO_COMPLEX(v, a); - TO_COMPLEX(w, b); - errno = 0; - quot = _Py_c_quot(a, b); + Py_complex a; + + if (PyComplex_Check(w)) { + Py_complex b = ((PyComplexObject *)(w))->cval; + TO_COMPLEX(v, a); + errno = 0; + a = _Py_c_quot(a, b); + } + else { + double b; + + if (to_float(&w, &b) < 0) { + return w; + } + if (b) { + a = ((PyComplexObject *)(v))->cval; + a.real /= b; + a.imag /= b; + } + else { + errno = EDOM; + } + } + if (errno == EDOM) { PyErr_SetString(PyExc_ZeroDivisionError, "division by zero"); return NULL; } - return PyComplex_FromCComplex(quot); + return PyComplex_FromCComplex(a); } static PyObject * diff --git a/Objects/floatobject.c b/Objects/floatobject.c index 68fd3e54632950..5a559298e9abf7 100644 --- a/Objects/floatobject.c +++ b/Objects/floatobject.c @@ -314,13 +314,13 @@ PyFloat_AsDouble(PyObject *op) #define CONVERT_TO_DOUBLE(obj, dbl) \ if (PyFloat_Check(obj)) \ dbl = PyFloat_AS_DOUBLE(obj); \ - else if (convert_to_double(&(obj), &(dbl)) < 0) \ + else if (_Py_convert_to_double(&(obj), &(dbl)) < 0) \ return obj; /* Methods */ -static int -convert_to_double(PyObject **v, double *dbl) +int +_Py_convert_to_double(PyObject **v, double *dbl) { PyObject *obj = *v; diff --git a/Python/bltinmodule.c b/Python/bltinmodule.c index f87f942cc76258..b2a18855d04d7f 100644 --- a/Python/bltinmodule.c +++ b/Python/bltinmodule.c @@ -2742,7 +2742,6 @@ builtin_sum_impl(PyObject *module, PyObject *iterable, PyObject *start) double value = PyLong_AsDouble(item); if (value != -1.0 || !PyErr_Occurred()) { re_sum = cs_add(re_sum, value); - im_sum.hi += 0.0; Py_DECREF(item); continue; } @@ -2755,7 +2754,6 @@ builtin_sum_impl(PyObject *module, PyObject *iterable, PyObject *start) if (PyFloat_Check(item)) { double value = PyFloat_AS_DOUBLE(item); re_sum = cs_add(re_sum, value); - im_sum.hi += 0.0; _Py_DECREF_SPECIALIZED(item, _PyFloat_ExactDealloc); continue; }