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
proposes similar extension to the Python language:

    * Added a new subtype (imaginary) of the complex type.  New type
      has 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
      - Lib/test/test_fractions.py
      - Lib/test/test_socket.py
      - Lib/test/test_str.py
    * Print dot for signed zeros in the real part of repr(complex):

          >>> complex(-0.0, 1)  # was (-0+1j)
          (-0.0+1j)

    * repr(complex) naw prints real part even if it's zero.

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

1) New code allows to use complex arithmetic for implementation of
   mathematical functions without special "corner cases".  Take the inverse
   tangent as an example:

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

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

       ±a±bj = complex(±float(a), 0.0) ± complex(0.0, float(b))

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

       ±a±bj = complex(±float(a), ±float(b))

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

   For example:

       >>> -0.0+1j  # now (-0.0+1j)
       1j
       >>> float('inf')*1j  # now infj
       (nan+infj)

3) The ``eval(repr(x)) == x`` invariant was broken for the complex type:

       >>> complex(-0.0, 1.0)  # also note funny signed integer zero below
       (-0+1j)
       >>> -0+1j
       1j
       >>> complex(0.0, -cmath.inf)
       -infj
       >>> -cmath.infj
       (-0-infj)

Author: skirpichev

Doc/data/stable_abi.dat

@@ -279,7 +279,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
 
@@ -295,7 +295,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

@@ -979,11 +979,12 @@ Imaginary literals are described by the following lexical definitions:
 .. productionlist:: python-grammar
    imagnumber: (`floatnumber` | `digitpart`) ("j" | "J")
 
-An imaginary literal yields a complex number with a real part of 0.0.  Complex
-numbers are represented as a pair of floating-point numbers and have the same
-restrictions on their range.  To create a complex number with a nonzero real
-part, add a floating-point number to it, e.g., ``(3+4j)``.  Some examples of
-imaginary literals::
+An imaginary literal yields a complex number without a real part, an instance
+of :class:`imaginary`.  Complex numbers are represented as a pair of
+floating-point numbers and have the same restrictions on their range.  To
+create a complex number with a nonzero real part, add an integer or
+floating-point number to imaginary, e.g., ``3+4j``.  Some examples of imaginary
+literals::
 
    3.14j   10.j    10j     .001j   1e100j   3.14e-10j   3.14_15_93j
 

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

@@ -73,7 +73,7 @@ def _raise_malformed_node(node):
             msg += f' on line {lno}'
         raise ValueError(msg + f': {node!r}')
     def _convert_num(node):
-        if not isinstance(node, Constant) or type(node.value) not in (int, float, complex):
+        if not isinstance(node, Constant) or type(node.value) not in (int, float, imaginary, complex):
             _raise_malformed_node(node)
         return node.value
     def _convert_signed_num(node):

Lib/ast.py

@@ -102,7 +102,7 @@ def copy(x):
 
 def _copy_immutable(x):
     return x
-for t in (types.NoneType, int, float, bool, complex, str, tuple,
+for t in (types.NoneType, int, float, bool, complex, imaginary, str, tuple,
           bytes, frozenset, type, range, slice, property,
           types.BuiltinFunctionType, types.EllipsisType,
           types.NotImplementedType, types.FunctionType, types.CodeType,
@@ -173,7 +173,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

@@ -23,6 +23,9 @@ class BadComplex3:
     def __complex__(self):
         raise RuntimeError
 
+class ImaginarySubclass(imaginary):
+    pass
+
 
 class CAPIComplexTest(ComplexesAreIdenticalMixin, unittest.TestCase):
     def test_check(self):
@@ -176,7 +179,7 @@ def test_py_cr_sum(self):
         # Test _Py_cr_sum()
         _py_cr_sum = _testcapi._py_cr_sum
 
-        self.assertComplexesAreIdentical(_py_cr_sum(-0j, -0.0)[0],
+        self.assertComplexesAreIdentical(_py_cr_sum(-0.0 - 0j, -0.0)[0],
                                          complex(-0.0, -0.0))
 
     def test_py_c_diff(self):
@@ -189,7 +192,7 @@ def test_py_cr_diff(self):
         # Test _Py_cr_diff()
         _py_cr_diff = _testcapi._py_cr_diff
 
-        self.assertComplexesAreIdentical(_py_cr_diff(-0j, 0.0)[0],
+        self.assertComplexesAreIdentical(_py_cr_diff(-0.0 - 0j, 0.0)[0],
                                          complex(-0.0, -0.0))
 
     def test_py_rc_diff(self):
@@ -275,7 +278,6 @@ def test_py_c_pow(self):
         self.assertEqual(_py_c_pow(max_num, 2),
                          (complex(INF, INF), errno.ERANGE))
 
-
     def test_py_c_abs(self):
         # Test _Py_c_abs()
         _py_c_abs = _testcapi._py_c_abs
@@ -294,5 +296,41 @@ def test_py_c_abs(self):
         self.assertEqual(_py_c_abs(complex(*[DBL_MAX]*2))[1], errno.ERANGE)
 
 
+class CAPIImaginaryTest(unittest.TestCase):
+    def test_check(self):
+        # Test PyImaginary_Check()
+        check = _testlimitedcapi.imaginary_check
+
+        self.assertTrue(check(2j))
+        self.assertTrue(check(ImaginarySubclass(2)))
+        self.assertFalse(check(ComplexSubclass(1+2j)))
+        self.assertFalse(check(Complex()))
+        self.assertFalse(check(3))
+        self.assertFalse(check(3.0))
+        self.assertFalse(check(object()))
+
+        # CRASHES check(NULL)
+
+    def test_checkexact(self):
+        # PyImaginary_CheckExact()
+        checkexact = _testlimitedcapi.imaginary_checkexact
+
+        self.assertTrue(checkexact(2j))
+        self.assertFalse(checkexact(ImaginarySubclass(2)))
+        self.assertFalse(checkexact(ComplexSubclass(1+2j)))
+        self.assertFalse(checkexact(Complex()))
+        self.assertFalse(checkexact(3))
+        self.assertFalse(checkexact(3.0))
+        self.assertFalse(checkexact(object()))
+
+        # CRASHES checkexact(NULL)
+
+    def test_fromdouble(self):
+        # Test PyImaginary_FromDouble()
+        fromdouble = _testlimitedcapi.imaginary_fromdouble
+
+        self.assertEqual(fromdouble(2.0), 2.0j)
+
+
 if __name__ == "__main__":
     unittest.main()