pythongh-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

@@ -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
 
 

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
 }

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.

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 *

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;
 

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