Imaginary type and IEC 60559-compatible complex arithmetic

"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.

That's why C standards since C99 introduce imaginary types.  This patch
implements similar extension to the Python language.

Lets consider (actually interrelated) problems, which will be solved on
this way.

1) Now complex arithmetic could be used for implementation of
   mathematical functions without special "corner cases", with textbooks
   formulae.  Take the inverse tangent as an example:

       >>> z = complex(-0.0, 2)
       >>> cmath.atan(z)
       (-1.5707963267948966+0.5493061443340549j)
       >>> atan = lambda z: 1j*(cmath.log(1 - 1j*z) - cmath.log(1 + 1j*z))/2
       >>> atan(z)  # real part had wrong sign before
       (-1.5707963267948966+0.5493061443340549j)

2) Previously, we have unsigned imaginary literals with the following
   semantics:

       a±bj = complex(a ± 0.0, ±b)
       complex(a, ±b)*cj = complex(a*0.0 ∓ b*c, a*c ± b*0.0)

   While this behaviour was well documented, most users would expect
   instead here:

       a±bj = complex(a, ±b)
       complex(a, ±b)*cj = complex(∓b*c, a*c)

   i.e. that it follows to the rectangular notation for complex numbers.

   For example:

       >>> -0.0+1j  # was 1j
       (-0.0+1j)
       >>> float('inf')*1j  # was (nan+infj)
       infj
       >>> -0.0+1j  # was 1j
       (-0.0+1j)
       >>> complex(-0.0, 1)  # was (-0+1j), not funny signed integer zero
       (-0.0+1j)

3) The ``eval(repr(x)) == x`` invariant now holds for the complex type.

What's changed:

    * Added a new subtype (imaginary) of the complex type with
      few overloaded methods (conjugate() and __getnewargs__()).
    * Complex and imaginary types implement IEC 60559-compatible complex
      arithmetic (as specified by C11 Annex G).
    * Imaginary literals now produce instances of imaginary type.
    * cmath.infj/nanj were changed to be of imaginary type.
    * Modules ast, code, copy, marshal got support for imaginary type.
    * Few tests adapted to use complex, instead of imaginary literals
    * Print dot for signed zeros in the real part of repr(complex)
    * repr(complex) now prints real part even if it's zero.

Author: skirpichev

Doc/c-api/complex.rst

@@ -47,6 +47,14 @@ pointers.  This is consistent throughout the API.
       This function is :term:`soft deprecated`.
 
 
+.. c:function:: Py_complex _Py_ci_sum(Py_complex left, double right)
+
+   Return the sum of a complex number and an imaginary number, using the C
+   :c:type:`Py_complex` representation.
+
+   .. versionadded:: next
+
+
 .. c:function:: Py_complex _Py_c_diff(Py_complex left, Py_complex right)
 
    Return the difference between two complex numbers, using the C
@@ -56,6 +64,22 @@ pointers.  This is consistent throughout the API.
       This function is :term:`soft deprecated`.
 
 
+.. c:function:: Py_complex _Py_ci_diff(Py_complex left, double right)
+
+   Return the difference between a complex number and an imaginary number,
+   using the C :c:type:`Py_complex` representation.
+
+   .. versionadded:: next
+
+
+.. c:function:: Py_complex _Py_ic_diff(double left, Py_complex right)
+
+   Return the difference between an imaginary number and a complex number,
+   using the C :c:type:`Py_complex` representation.
+
+   .. versionadded:: next
+
+
 .. c:function:: Py_complex _Py_c_neg(Py_complex num)
 
    Return the negation of the complex number *num*, using the C
@@ -74,6 +98,14 @@ pointers.  This is consistent throughout the API.
       This function is :term:`soft deprecated`.
 
 
+.. c:function:: Py_complex _Py_ci_prod(Py_complex left, double right)
+
+   Return the product of a complex number and an imaginary number, using the C
+   :c:type:`Py_complex` representation.
+
+   .. versionadded:: next
+
+
 .. c:function:: Py_complex _Py_c_quot(Py_complex dividend, Py_complex divisor)
 
    Return the quotient of two complex numbers, using the C :c:type:`Py_complex`
@@ -86,6 +118,28 @@ pointers.  This is consistent throughout the API.
       This function is :term:`soft deprecated`.
 
 
+.. c:function:: Py_complex _Py_ci_quot(Py_complex dividend, double divisor)
+
+   Return the quotient of a complex number and an imaginary number, using the C
+   :c:type:`Py_complex` representation.
+
+   If *divisor* is zero, this method returns zero and sets
+   :c:data:`errno` to :c:macro:`!EDOM`.
+
+   .. versionadded:: next
+
+
+.. c:function:: Py_complex _Py_ic_quot(double dividend, Py_complex divisor)
+
+   Return the quotient of an imaginary number and a complex number, using the C
+   :c:type:`Py_complex` representation.
+
+   If *divisor* is zero, this method returns zero and sets
+   :c:data:`errno` to :c:macro:`!EDOM`.
+
+   .. versionadded:: next
+
+
 .. c:function:: Py_complex _Py_c_pow(Py_complex num, Py_complex exp)
 
    Return the exponentiation of *num* by *exp*, using the C :c:type:`Py_complex`

Doc/data/stable_abi.dat

@@ -330,7 +330,7 @@ Constants
 .. data:: infj
 
    Complex number with zero real part and positive infinity imaginary
-   part. Equivalent to ``complex(0.0, float('inf'))``.
+   part. Equivalent to ``float('inf')*1j``.
 
    .. versionadded:: 3.6
 
@@ -346,7 +346,7 @@ Constants
 .. data:: nanj
 
    Complex number with zero real part and NaN imaginary part. Equivalent to
-   ``complex(0.0, float('nan'))``.
+   ``float('nan')*1j``.
 
    .. versionadded:: 3.6
 

Doc/library/cmath.rst

@@ -33,10 +33,11 @@ are always available.  They are listed here in alphabetical order.
 | |  :func:`complex`      | |                     | |  **P**              | |  **V**                |
 | |                       | |  **I**              | |  :func:`pow`        | |  :func:`vars`         |
 | |  **D**                | |  :func:`id`         | |  :func:`print`      | |                       |
-| |  :func:`delattr`      | |  :func:`input`      | |  :func:`property`   | |  **Z**                |
-| |  |func-dict|_         | |  :func:`int`        | |                     | |  :func:`zip`          |
-| |  :func:`dir`          | |  :func:`isinstance` | |                     | |                       |
-| |  :func:`divmod`       | |  :func:`issubclass` | |                     | |  **_**                |
+| |  :func:`delattr`      | |  :func:`imaginary`  | |  :func:`property`   | |                       |
+| |  |func-dict|_         | |  :func:`input`      | |                     | |  **Z**                |
+| |  :func:`dir`          | |  :func:`int`        | |                     | |  :func:`zip`          |
+| |  :func:`divmod`       | |  :func:`isinstance` | |                     | |                       |
+| |                       | |  :func:`issubclass` | |                     | |  **_**                |
 | |                       | |  :func:`iter`       | |                     | |  :func:`__import__`   |
 +-------------------------+-----------------------+-----------------------+-------------------------+
 
@@ -388,7 +389,7 @@ are always available.  They are listed here in alphabetical order.
       >>> complex('+1.23')
       (1.23+0j)
       >>> complex('-4.5j')
-      -4.5j
+      (0.0-4.5j)
       >>> complex('-1.23+4.5j')
       (-1.23+4.5j)
       >>> complex('\t( -1.23+4.5J )\n')
@@ -398,7 +399,7 @@ are always available.  They are listed here in alphabetical order.
       >>> complex(1.23)
       (1.23+0j)
       >>> complex(imag=-4.5)
-      -4.5j
+      (0.0-4.5j)
       >>> complex(-1.23, 4.5)
       (-1.23+4.5j)
 
@@ -442,7 +443,7 @@ are always available.  They are listed here in alphabetical order.
 
    See also :meth:`complex.from_number` which only accepts a single numeric argument.
 
-   If all arguments are omitted, returns ``0j``.
+   If all arguments are omitted, returns ``0.0+0j``.
 
    The complex type is described in :ref:`typesnumeric`.
 
@@ -970,6 +971,14 @@ are always available.  They are listed here in alphabetical order.
    .. audit-event:: builtins.id id id
 
 
+.. class:: imaginary(x=0.0)
+
+   Return an imaginary number with the value ``float(x)*1j``.  If argument is
+   omitted, returns ``0j``.
+
+   .. versionadded:: next
+
+
 .. function:: input()
               input(prompt)
 

Doc/library/functions.rst

@@ -1297,8 +1297,8 @@ Imaginary literals
 
 Python has :ref:`complex number <typesnumeric>` objects, but no complex
 literals.
-Instead, *imaginary literals* denote complex numbers with a zero
-real part.
+Instead, *imaginary literals* denote complex numbers without a real part, an instance
+of :class:`imaginary`.
 
 For example, in math, the complex number 3+4.2\ *i* is written
 as the real number 3 added to the imaginary number 4.2\ *i*.

Doc/reference/lexical_analysis.rst

@@ -9,15 +9,21 @@ extern "C" {
 /* Complex object interface */
 
 PyAPI_DATA(PyTypeObject) PyComplex_Type;
+PyAPI_DATA(PyTypeObject) PyImaginary_Type;
 
 #define PyComplex_Check(op) PyObject_TypeCheck((op), &PyComplex_Type)
 #define PyComplex_CheckExact(op) Py_IS_TYPE((op), &PyComplex_Type)
 
+#define PyImaginary_Check(op) PyObject_TypeCheck((op), &PyImaginary_Type)
+#define PyImaginary_CheckExact(op) Py_IS_TYPE((op), &PyImaginary_Type)
+
 PyAPI_FUNC(PyObject *) PyComplex_FromDoubles(double real, double imag);
 
 PyAPI_FUNC(double) PyComplex_RealAsDouble(PyObject *op);
 PyAPI_FUNC(double) PyComplex_ImagAsDouble(PyObject *op);
 
+PyAPI_FUNC(PyObject *) PyImaginary_FromDouble(double imag);
+
 #ifndef Py_LIMITED_API
 #  define Py_CPYTHON_COMPLEXOBJECT_H
 #  include "cpython/complexobject.h"

Include/complexobject.h

@@ -10,10 +10,22 @@ typedef struct {
 /* Operations on complex numbers (soft deprecated
    since Python 3.15). */
 PyAPI_FUNC(Py_complex) _Py_c_sum(Py_complex, Py_complex);
+PyAPI_FUNC(Py_complex) _Py_cr_sum(Py_complex, double);
+PyAPI_FUNC(Py_complex) _Py_ci_sum(Py_complex, double);
 PyAPI_FUNC(Py_complex) _Py_c_diff(Py_complex, Py_complex);
+PyAPI_FUNC(Py_complex) _Py_cr_diff(Py_complex, double);
+PyAPI_FUNC(Py_complex) _Py_rc_diff(double, Py_complex);
+PyAPI_FUNC(Py_complex) _Py_ci_diff(Py_complex, double);
+PyAPI_FUNC(Py_complex) _Py_ic_diff(double, Py_complex);
 PyAPI_FUNC(Py_complex) _Py_c_neg(Py_complex);
 PyAPI_FUNC(Py_complex) _Py_c_prod(Py_complex, Py_complex);
+PyAPI_FUNC(Py_complex) _Py_cr_prod(Py_complex, double);
+PyAPI_FUNC(Py_complex) _Py_ci_prod(Py_complex, double);
 PyAPI_FUNC(Py_complex) _Py_c_quot(Py_complex, Py_complex);
+PyAPI_FUNC(Py_complex) _Py_cr_quot(Py_complex, double);
+PyAPI_FUNC(Py_complex) _Py_rc_quot(double, Py_complex);
+PyAPI_FUNC(Py_complex) _Py_ci_quot(Py_complex, double);
+PyAPI_FUNC(Py_complex) _Py_ic_quot(double, Py_complex);
 PyAPI_FUNC(Py_complex) _Py_c_pow(Py_complex, Py_complex);
 PyAPI_FUNC(double) _Py_c_abs(Py_complex);
 

Include/cpython/complexobject.h

@@ -89,7 +89,7 @@ def _convert_literal(node):
         isinstance(node, UnaryOp)
         and isinstance(node.op, (UAdd, USub))
         and isinstance(node.operand, Constant)
-        and type(operand := node.operand.value) in (int, float, complex)
+        and type(operand := node.operand.value) in (int, float, imaginary)
     ):
         if isinstance(node.op, UAdd):
             return + operand
@@ -101,7 +101,7 @@ def _convert_literal(node):
         and isinstance(node.left, (Constant, UnaryOp))
         and isinstance(node.right, Constant)
         and type(left := _convert_literal(node.left)) in (int, float)
-        and type(right := _convert_literal(node.right)) is complex
+        and type(right := _convert_literal(node.right)) is imaginary
     ):
         if isinstance(node.op, Add):
             return left + right

Lib/ast.py

@@ -101,7 +101,7 @@ def copy(x):
 
 
 _copy_atomic_types = {types.NoneType, int, float, bool, complex, str, tuple,
-          bytes, frozenset, type, range, slice, property,
+          imaginary, bytes, frozenset, type, range, slice, property,
           types.BuiltinFunctionType, types.EllipsisType,
           types.NotImplementedType, types.FunctionType, types.CodeType,
           weakref.ref, super}
@@ -164,7 +164,8 @@ def deepcopy(x, memo=None, _nil=[]):
 
 _atomic_types =  {types.NoneType, types.EllipsisType, types.NotImplementedType,
           int, float, bool, complex, bytes, str, types.CodeType, type, range,
-          types.BuiltinFunctionType, types.FunctionType, weakref.ref, property}
+          types.BuiltinFunctionType, types.FunctionType, weakref.ref,
+          property, imaginary}
 
 _deepcopy_dispatch = d = {}
 

Lib/copy.py

@@ -785,12 +785,12 @@ def __int__(self):
             ('%X format: an integer is required, not float', b'%X', 2.11),
             ('%o format: an integer is required, not float', b'%o', 1.79),
             ('%x format: an integer is required, not PseudoFloat', b'%x', pi),
-            ('%x format: an integer is required, not complex', b'%x', 3j),
-            ('%X format: an integer is required, not complex', b'%X', 2j),
-            ('%o format: an integer is required, not complex', b'%o', 1j),
-            ('%u format: a real number is required, not complex', b'%u', 3j),
-            ('%i format: a real number is required, not complex', b'%i', 2j),
-            ('%d format: a real number is required, not complex', b'%d', 2j),
+            ('%x format: an integer is required, not complex', b'%x', 1+3j),
+            ('%X format: an integer is required, not complex', b'%X', 1+2j),
+            ('%o format: an integer is required, not complex', b'%o', 1+1j),
+            ('%u format: a real number is required, not complex', b'%u', 1+3j),
+            ('%i format: a real number is required, not complex', b'%i', 1+2j),
+            ('%d format: a real number is required, not complex', b'%d', 1+2j),
             (
                 r'%c requires an integer in range\(256\)'
                 r' or a single byte, not .*\.PseudoFloat',