pythongh-69639: add mixed-mode rules for complex arithmetic (C-like)

"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 (in particular,
there is no special version for the true division with real lhs
operand).  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).

Author: skirpichev

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
 
 

Doc/library/stdtypes.rst

@@ -243,6 +243,9 @@ numeric literal yields an imaginary number (a complex number with a zero real
 part) which you can add to an integer or float to get a complex number with real
 and imaginary parts.
 
+The constructors :func:`int`, :func:`float`, and
+:func:`complex` can be used to produce numbers of a specific type.
+
 .. index::
    single: arithmetic
    pair: built-in function; int
@@ -262,12 +265,15 @@ and imaginary parts.
 
 Python fully supports mixed arithmetic: when a binary arithmetic operator has
 operands of different numeric types, the operand with the "narrower" type is
-widened to that of the other, where integer is narrower than floating point,
-which is narrower than complex. A comparison between numbers of different types
-behaves as though the exact values of those numbers were being compared. [2]_
+widened to that of the other, where integer is narrower than floating point.
+Arithmetic with complex and real operands is defined by the usual mathematical
+formula, for example::
 
-The constructors :func:`int`, :func:`float`, and
-:func:`complex` can be used to produce numbers of a specific type.
+    x + complex(u, v) = complex(x + u, v)
+    x * complex(u, v) = complex(x * u, x * v)
+
+A comparison between numbers of different types behaves as though the exact
+values of those numbers were being compared. [2]_
 
 All numeric types (except complex) support the following operations (for priorities of
 the operations, see :ref:`operator-summary`):

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
 }

Lib/test/test_complex.py

@@ -127,6 +127,9 @@ def test_truediv(self):
             self.assertTrue(isnan(z.real))
             self.assertTrue(isnan(z.imag))
 
+        self.assertComplexesAreIdentical(complex(INF, NAN) / 2,
+                                         complex(INF, NAN))
+
         self.assertComplexesAreIdentical(complex(INF, 1)/(0.0+1j),
                                          complex(NAN, -INF))
 
@@ -224,20 +227,32 @@ def check(n, deltas, is_equal, imag = 0.0):
     def test_add(self):
         self.assertEqual(1j + int(+1), complex(+1, 1))
         self.assertEqual(1j + int(-1), complex(-1, 1))
+        self.assertComplexesAreIdentical(complex(-0.0, -0.0) + (-0.0),
+                                         complex(-0.0, -0.0))
+        self.assertComplexesAreIdentical((-0.0) + complex(-0.0, -0.0),
+                                         complex(-0.0, -0.0))
         self.assertRaises(OverflowError, operator.add, 1j, 10**1000)
         self.assertRaises(TypeError, operator.add, 1j, None)
         self.assertRaises(TypeError, operator.add, None, 1j)
 
     def test_sub(self):
         self.assertEqual(1j - int(+1), complex(-1, 1))
         self.assertEqual(1j - int(-1), complex(1, 1))
+        self.assertComplexesAreIdentical(complex(-0.0, -0.0) - 0.0,
+                                         complex(-0.0, -0.0))
+        self.assertComplexesAreIdentical(-0.0 - complex(0.0, 0.0),
+                                         complex(-0.0, -0.0))
         self.assertRaises(OverflowError, operator.sub, 1j, 10**1000)
         self.assertRaises(TypeError, operator.sub, 1j, None)
         self.assertRaises(TypeError, operator.sub, None, 1j)
 
     def test_mul(self):
         self.assertEqual(1j * int(20), complex(0, 20))
         self.assertEqual(1j * int(-1), complex(0, -1))
+        self.assertComplexesAreIdentical(complex(INF, NAN) * 2,
+                                         complex(INF, NAN))
+        self.assertComplexesAreIdentical(2 * complex(INF, NAN),
+                                         complex(INF, NAN))
         self.assertRaises(OverflowError, operator.mul, 1j, 10**1000)
         self.assertRaises(TypeError, operator.mul, 1j, None)
         self.assertRaises(TypeError, operator.mul, None, 1j)

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,163 @@ 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)
+       float(x) + complex(u, v) = complex(x + u, v)
+
+   Similar rules are implemented for subtraction and multiplication.  See C11's
+   Annex G, sections G.5.1 and G.5.2.  The true division is special:
+
+       complex(x, y) / float(u) = complex(x/u, y/u)
+       float(x) / complex(u, v) = complex(x, 0) / complex(u, v)
+ */
 
 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

@@ -309,13 +309,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;
             }