From d062169eaeb93ee1ca87fd5e373365be47225859 Mon Sep 17 00:00:00 2001
From: Sergey B Kirpichev <skirpichev@gmail.com>
Date: Tue, 23 Apr 2024 18:00:55 +0300
Subject: [PATCH 1/2] Imaginary type and IEC 60559-compatible complex
 arithmetic
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

"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)
---
 Doc/data/stable_abi.dat                     |   2 +
 Doc/library/cmath.rst                       |   4 +-
 Doc/library/functions.rst                   |  23 +-
 Doc/reference/lexical_analysis.rst          |  11 +-
 Include/complexobject.h                     |   6 +
 Lib/ast.py                                  |   2 +-
 Lib/copy.py                                 |   5 +-
 Lib/test/test_capi/test_complex.py          |  44 ++-
 Lib/test/test_complex.py                    | 184 +++++++++---
 Lib/test/test_format.py                     |  13 +-
 Lib/test/test_fractions.py                  |   2 +-
 Lib/test/test_marshal.py                    |  11 +-
 Lib/test/test_socket.py                     |   4 +-
 Lib/test/test_stable_abi_ctypes.py          |   2 +
 Lib/test/test_str.py                        |  12 +-
 Misc/stable_abi.toml                        |   5 +
 Modules/_testlimitedcapi/complex.c          |  23 ++
 Modules/cmathmodule.c                       |   6 +-
 Objects/clinic/complexobject.c.h            |  99 ++++++-
 Objects/codeobject.c                        |  14 +
 Objects/complexobject.c                     | 293 +++++++++++++++++++-
 Objects/object.c                            |   1 +
 PC/python3dll.c                             |   2 +
 Parser/action_helpers.c                     |   5 +-
 Parser/pegen.c                              |   9 +-
 Python/ast.c                                |   5 +-
 Python/ast_unparse.c                        |   2 +-
 Python/bltinmodule.c                        |   6 +
 Python/formatter_unicode.c                  |  15 +-
 Python/marshal.c                            |  32 +++
 Tools/c-analyzer/cpython/globals-to-fix.tsv |   1 +
 31 files changed, 743 insertions(+), 100 deletions(-)

diff --git a/Doc/data/stable_abi.dat b/Doc/data/stable_abi.dat
index 6f9d27297e8f65..a4335311d8bd7d 100644
--- a/Doc/data/stable_abi.dat
+++ b/Doc/data/stable_abi.dat
@@ -304,6 +304,8 @@ func,PyGILState_Release,3.2,,
 type,PyGILState_STATE,3.2,,
 type,PyGetSetDef,3.2,,full-abi
 data,PyGetSetDescr_Type,3.2,,
+func,PyImaginary_FromDouble,3.14,,
+data,PyImaginary_Type,3.14,,
 func,PyImport_AddModule,3.2,,
 func,PyImport_AddModuleObject,3.7,,
 func,PyImport_AddModuleRef,3.13,,
diff --git a/Doc/library/cmath.rst b/Doc/library/cmath.rst
index e7c027dd4d0c22..ef613863e34496 100644
--- a/Doc/library/cmath.rst
+++ b/Doc/library/cmath.rst
@@ -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
 
diff --git a/Doc/library/functions.rst b/Doc/library/functions.rst
index a7549b9bce76e2..03e9c86b51795c 100644
--- a/Doc/library/functions.rst
+++ b/Doc/library/functions.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)
 
diff --git a/Doc/reference/lexical_analysis.rst b/Doc/reference/lexical_analysis.rst
index f7167032ad7df9..6279a6c8a85e8b 100644
--- a/Doc/reference/lexical_analysis.rst
+++ b/Doc/reference/lexical_analysis.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
 
diff --git a/Include/complexobject.h b/Include/complexobject.h
index ebe49a832f7414..2ff40ecfe52128 100644
--- a/Include/complexobject.h
+++ b/Include/complexobject.h
@@ -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"
diff --git a/Lib/ast.py b/Lib/ast.py
index 154d2c8c1f9ebb..0991927a58a819 100644
--- a/Lib/ast.py
+++ b/Lib/ast.py
@@ -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):
diff --git a/Lib/copy.py b/Lib/copy.py
index f27e109973cfb7..333491ad901c1c 100644
--- a/Lib/copy.py
+++ b/Lib/copy.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 = {}
 
diff --git a/Lib/test/test_capi/test_complex.py b/Lib/test/test_capi/test_complex.py
index 97e0eb3f043080..1f30fd24307d99 100644
--- a/Lib/test/test_capi/test_complex.py
+++ b/Lib/test/test_capi/test_complex.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()
diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py
index 179556f57e884f..25c2b2983d21ab 100644
--- a/Lib/test/test_complex.py
+++ b/Lib/test/test_complex.py
@@ -37,6 +37,9 @@ def __float__(self):
 class ComplexSubclass(complex):
     pass
 
+class ImaginarySubclass(imaginary):
+    pass
+
 class OtherComplexSubclass(complex):
     pass
 
@@ -193,6 +196,26 @@ def test_truediv(self):
         self.assertComplexesAreIdentical(float(INF)/complex(NAN, INF),
                                          complex(NAN, NAN))
 
+        self.assertEqual((1+1j)/float(2), 0.5+0.5j)
+        self.assertRaises(TypeError, operator.truediv, None, 1+1j)
+        self.assertRaises(TypeError, operator.truediv, 1+1j, None)
+
+        self.assertEqual(1j/imaginary(2), 0.5)
+        self.assertEqual((1+1j)/imaginary(2), 0.5-0.5j)
+        self.assertEqual(1j/complex(1, 1), 0.5+0.5j)
+        self.assertEqual(1j/float(2), 0.5j)
+        self.assertEqual(float(1)/(1+2j), 0.2-0.4j)
+        self.assertEqual(float(1)/(-1+2j), -0.2-0.4j)
+        self.assertEqual(float(1)/(1-2j), 0.2+0.4j)
+
+        z = float(1)/(NAN+2j)
+        self.assertTrue(isnan(z.real))
+        self.assertTrue(isnan(z.imag))
+
+        self.assertRaises(ZeroDivisionError, operator.truediv, 1j, 0j)
+        self.assertRaises(ZeroDivisionError, operator.truediv, 1j, complex(0, 0))
+        self.assertRaises(ZeroDivisionError, operator.truediv, 1j, 0.0)
+
     def test_truediv_zero_division(self):
         for a, b in ZERO_DIVISION:
             with self.assertRaises(ZeroDivisionError):
@@ -267,10 +290,28 @@ def test_add(self):
                                          complex(-0.0, -0.0))
         self.assertComplexesAreIdentical((-0.0) + complex(-0.0, -0.0),
                                          complex(-0.0, -0.0))
+        self.assertEqual(1j + imaginary(1), imaginary(2))
         self.assertRaises(OverflowError, operator.add, 1j, 10**1000)
         self.assertRaises(TypeError, operator.add, 1j, None)
+        self.assertRaises(TypeError, operator.add, 1+1j, None)
         self.assertRaises(TypeError, operator.add, None, 1j)
 
+        self.assertComplexesAreIdentical(float(0.0) + 0j, complex(0, 0))
+        self.assertComplexesAreIdentical(0j + float(0.0), complex(0, 0))
+        self.assertComplexesAreIdentical(float(-0.0) + 0j, complex(-0.0, 0))
+        self.assertComplexesAreIdentical(0j + float(-0.0), complex(-0.0, 0))
+        self.assertComplexesAreIdentical((-0.0+0j) + float(0.0), complex(0, 0))
+        self.assertComplexesAreIdentical(float(0.0) + (-0.0+0j), complex(0, 0))
+        self.assertComplexesAreIdentical((1+0j) + complex(1-0j), complex(2, 0))
+        self.assertComplexesAreIdentical(0j + complex(-0.0-0j), complex(-0.0, 0))
+        self.assertComplexesAreIdentical(0j + complex(-0j), complex(0, 0))
+        self.assertComplexesAreIdentical((1+0j) + complex(-0.0-0j), complex(1, 0))
+        self.assertComplexesAreIdentical(complex(-0.0+0j) + (-0j), complex(-0.0, 0))
+        self.assertComplexesAreIdentical((1+0j) + float(1.0), complex(2, 0))
+        self.assertComplexesAreIdentical(float(1.0) + (1+0j), complex(2, 0))
+        self.assertComplexesAreIdentical((1-0j) + float(1.0), complex(2, -0.0))
+        self.assertComplexesAreIdentical(float(1.0) + (1-0j), complex(2, -0.0))
+
     def test_sub(self):
         self.assertEqual(1j - int(+1), complex(-1, 1))
         self.assertEqual(1j - int(-1), complex(1, 1))
@@ -282,13 +323,34 @@ def test_sub(self):
                                          complex(-1, 1))
         self.assertComplexesAreIdentical(complex(2, 1) - complex(1, 2),
                                          complex(1, -1))
+        self.assertEqual(1j - imaginary(2), imaginary(-1))
         self.assertRaises(OverflowError, operator.sub, 1j, 10**1000)
         self.assertRaises(TypeError, operator.sub, 1j, None)
+        self.assertRaises(TypeError, operator.sub, 1+1j, None)
+        self.assertRaises(TypeError, operator.sub, None, 1+1j)
         self.assertRaises(TypeError, operator.sub, None, 1j)
+        self.assertRaises(TypeError, operator.sub, 1j, None)
+
+        self.assertComplexesAreIdentical(float(0.0) - 0j, complex(0, -0.0))
+        self.assertComplexesAreIdentical(0j - float(0.0), complex(-0.0, 0))
+        self.assertComplexesAreIdentical(float(-0.0) - 0j, complex(-0.0, -0.0))
+        self.assertComplexesAreIdentical(0j - float(-0.0), complex(0, 0))
+        self.assertComplexesAreIdentical((-0.0+0j) - float(0.0), complex(-0.0, 0))
+        self.assertComplexesAreIdentical(float(0.0) - (-0.0+0j), complex(0, -0.0))
+        self.assertComplexesAreIdentical((1+0j) - complex(1-0j), complex(0, 0))
+        self.assertComplexesAreIdentical(0j - complex(-0.0-0j), complex(0, 0))
+        self.assertComplexesAreIdentical(0j - complex(-0j), complex(-0.0, 0))
+        self.assertComplexesAreIdentical((1+0j) - complex(-0.0-0j), complex(1, 0))
+        self.assertComplexesAreIdentical(complex(-0.0+0j) - (-0j), complex(-0.0, 0))
+        self.assertComplexesAreIdentical((1+0j) - float(1.0), complex(0, 0))
+        self.assertComplexesAreIdentical(float(1.0) - (1+0j), complex(0, -0.0))
+        self.assertComplexesAreIdentical((1-0j) - float(1.0), complex(0, -0.0))
+        self.assertComplexesAreIdentical(float(1.0) - (1-0j), complex(0, 0))
 
     def test_mul(self):
         self.assertEqual(1j * int(20), complex(0, 20))
         self.assertEqual(1j * int(-1), complex(0, -1))
+        self.assertEqual(2j * imaginary(3), -6.0)
         for c, r in [(2, complex(INF, 2)), (INF, complex(INF, INF)),
                      (0, complex(NAN, 0)), (-0.0, complex(NAN, -0.0)),
                      (NAN, complex(NAN, NAN))]:
@@ -297,8 +359,27 @@ def test_mul(self):
                 self.assertComplexesAreIdentical(c * complex(INF, 1), r)
         self.assertRaises(OverflowError, operator.mul, 1j, 10**1000)
         self.assertRaises(TypeError, operator.mul, 1j, None)
+        self.assertRaises(TypeError, operator.mul, 1+1j, None)
         self.assertRaises(TypeError, operator.mul, None, 1j)
 
+        self.assertComplexesAreIdentical(float(0.0) * 0j, complex(0, 0))
+        self.assertComplexesAreIdentical(0j * float(0.0), complex(0.0, 0))
+        self.assertComplexesAreIdentical(float(-0.0) * 0j, complex(0.0, -0.0))
+        self.assertComplexesAreIdentical(0j * float(-0.0), complex(0, -0.0))
+        self.assertComplexesAreIdentical((-0.0+0j) * float(0.0), complex(-0.0, 0))
+        self.assertComplexesAreIdentical(float(0.0) * (-0.0+0j), complex(-0.0, 0))
+        self.assertComplexesAreIdentical((-0.0+0j) * float(-0.0), complex(0, -0.0))
+        self.assertComplexesAreIdentical(float(-0.0) * (-0.0+0j), complex(0, -0.0))
+        self.assertComplexesAreIdentical((-0.0-0j) * float(-0.0), complex(0, 0))
+        self.assertComplexesAreIdentical(float(-0.0) * (-0.0-0j), complex(0, 0))
+        self.assertComplexesAreIdentical((1+0j) * complex(1-0j), complex(1, 0))
+        self.assertComplexesAreIdentical(0j * complex(-0.0-0j), complex(0.0, -0.0))
+        self.assertComplexesAreIdentical(0j * complex(-0j), complex(0, 0))
+        self.assertComplexesAreIdentical((1+0j) * complex(-0.0-0j), complex(0, -0.0))
+        self.assertComplexesAreIdentical((-0.0+0j) * complex(-0j), complex(0, 0))
+        self.assertComplexesAreIdentical((1+0j) * float(1.0), complex(1, 0))
+        self.assertComplexesAreIdentical(float(1.0) * (1+0j), complex(1, 0))
+
     def test_mod(self):
         # % is no longer supported on complex numbers
         with self.assertRaises(TypeError):
@@ -436,6 +517,7 @@ def test_boolcontext(self):
 
     def test_conjugate(self):
         self.assertClose(complex(5.3, 9.8).conjugate(), 5.3-9.8j)
+        self.assertEqual(1j.conjugate(), -1j)
 
     def test_constructor(self):
         def check(z, x, y):
@@ -469,10 +551,10 @@ def check(z, x, y):
                 "argument 'real' must be a real number, not .*WithComplex"):
             check(complex(WithComplex(4.25+0j), 0), 4.25, 0.0)
         with self.assertWarnsRegex(DeprecationWarning,
-                "argument 'real' must be a real number, not complex"):
+                "argument 'real' must be a real number, not imaginary"):
             check(complex(4.25j, 0), 0.0, 4.25)
         with self.assertWarnsRegex(DeprecationWarning,
-                "argument 'real' must be a real number, not complex"):
+                "argument 'real' must be a real number, not imaginary"):
             check(complex(0j, 4.25), 0.0, 4.25)
         with self.assertWarnsRegex(DeprecationWarning,
                 "argument 'imag' must be a real number, not complex"):
@@ -484,19 +566,19 @@ def check(z, x, y):
                 "argument 'imag' must be a real number, not .*WithComplex"):
             complex(0, WithComplex(4.25+0j))
         with self.assertWarnsRegex(DeprecationWarning,
-                "argument 'imag' must be a real number, not complex"):
+                "argument 'imag' must be a real number, not imaginary"):
             check(complex(0.0, 4.25j), -4.25, 0.0)
         with self.assertWarnsRegex(DeprecationWarning,
                 "argument 'real' must be a real number, not complex"):
             check(complex(4.25+0j, 0j), 4.25, 0.0)
         with self.assertWarnsRegex(DeprecationWarning,
-                "argument 'real' must be a real number, not complex"):
+                "argument 'real' must be a real number, not imaginary"):
             check(complex(4.25j, 0j), 0.0, 4.25)
         with self.assertWarnsRegex(DeprecationWarning,
-                "argument 'real' must be a real number, not complex"):
+                "argument 'real' must be a real number, not imaginary"):
             check(complex(0j, 4.25+0j), 0.0, 4.25)
         with self.assertWarnsRegex(DeprecationWarning,
-                "argument 'real' must be a real number, not complex"):
+                "argument 'real' must be a real number, not imaginary"):
             check(complex(0j, 4.25j), -4.25, 0.0)
 
         check(complex(real=4.25), 4.25, 0.0)
@@ -580,18 +662,14 @@ def __complex__(self):
         self.assertRaises(TypeError, complex, WithIndex(None), 1.5)
         self.assertRaises(TypeError, complex, 1.5, WithIndex(None))
 
-        class MyInt:
-            def __int__(self):
-                return 42
-
-        self.assertRaises(TypeError, complex, MyInt())
-        self.assertRaises(TypeError, complex, MyInt(), 1.5)
-        self.assertRaises(TypeError, complex, 1.5, MyInt())
+        self.assertRaises(TypeError, complex, MyInt(42))
+        self.assertRaises(TypeError, complex, MyInt(42), 1.5)
+        self.assertRaises(TypeError, complex, 1.5, MyInt(42))
 
         class complex0(complex):
             """Test usage of __complex__() when inheriting from 'complex'"""
             def __complex__(self):
-                return 42j
+                return 1+42j
 
         class complex1(complex):
             """Test usage of __complex__() with a __new__() method"""
@@ -606,11 +684,30 @@ class complex2(complex):
             def __complex__(self):
                 return None
 
-        check(complex(complex0(1j)), 0.0, 42.0)
+        check(complex(complex0(1j)), 1.0, 42.0)
         with self.assertWarns(DeprecationWarning):
             check(complex(complex1(1j)), 0.0, 2.0)
         self.assertRaises(TypeError, complex, complex2(1j))
 
+    def test_imaginary_constructor(self):
+        self.assertEqual(imaginary(), 0j)
+        self.assertEqual(imaginary(-2), -2j)
+        self.assertEqual(imaginary(1.25), 1.25j)
+        self.assertEqual(imaginary("1.25"), 1.25j)
+
+        self.assertEqual(imaginary(WithFloat(42.)), 42j)
+        self.assertRaises(TypeError, imaginary, WithFloat(None))
+
+        self.assertEqual(imaginary(WithIndex(42)), 42j)
+        self.assertRaises(OverflowError, imaginary, WithIndex(2**2000))
+
+        self.assertRaises(TypeError, imaginary, MyInt(42))
+        self.assertRaises(TypeError, imaginary, 123, MyInt(42))
+        self.assertRaises(ValueError, imaginary, "1.25j")
+
+        self.assertRaises(TypeError, imaginary, complex())
+        self.assertRaises(TypeError, imaginary, object())
+
     def test___complex__(self):
         z = 3 + 4j
         self.assertEqual(z.__complex__(), z)
@@ -784,9 +881,13 @@ def test(v, expected, test_fn=self.assertEqual):
         test(complex(NAN, NAN), "(nan+nanj)")
         test(complex(-NAN, -NAN), "(nan+nanj)")
 
-        test(complex(0, INF), "infj")
-        test(complex(0, -INF), "-infj")
-        test(complex(0, NAN), "nanj")
+        test(complex(0, INF), "(0.0+infj)")
+        test(complex(0, -INF), "(0.0-infj)")
+        test(complex(0, NAN), "(0.0+nanj)")
+
+        test(imaginary(INF), "infj")
+        test(imaginary(-INF), "-infj")
+        test(imaginary(NAN), "nanj")
 
         self.assertEqual(1-6j,complex(repr(1-6j)))
         self.assertEqual(1+6j,complex(repr(1+6j)))
@@ -799,29 +900,43 @@ def test(v, expected, test_fn=self.assertEqual):
             test_fn(repr(v), expected)
             test_fn(str(v), expected)
 
-        test(complex(0., 1.),   "1j")
-        test(complex(-0., 1.),  "(-0+1j)")
-        test(complex(0., -1.),  "-1j")
-        test(complex(-0., -1.), "(-0-1j)")
+        test(complex(0., 1.),   "(0.0+1j)")
+        test(complex(-0., 1.),  "(-0.0+1j)")
+        test(complex(0., -1.),  "(0.0-1j)")
+        test(complex(-0., -1.), "(-0.0-1j)")
 
-        test(complex(0., 0.),   "0j")
-        test(complex(0., -0.),  "-0j")
-        test(complex(-0., 0.),  "(-0+0j)")
-        test(complex(-0., -0.), "(-0-0j)")
+        test(imaginary(+1.), "1j")
+        test(imaginary(-1.), "-1j")
+
+        test(complex(0., 0.),   "(0.0+0j)")
+        test(complex(0., -0.),  "(0.0-0j)")
+        test(complex(-0., 0.),  "(-0.0+0j)")
+        test(complex(-0., -0.), "(-0.0-0j)")
+
+        test(imaginary(+0.0), "0j")
+        test(imaginary(-0.0), "-0j")
 
     def test_pos(self):
         self.assertEqual(+(1+6j), 1+6j)
         self.assertEqual(+ComplexSubclass(1, 6), 1+6j)
         self.assertIs(type(+ComplexSubclass(1, 6)), complex)
 
+        self.assertEqual(+1j, 1j)
+        self.assertEqual(+ImaginarySubclass(1), 1j)
+        self.assertIs(type(+ImaginarySubclass(1)), imaginary)
+
     def test_neg(self):
         self.assertEqual(-(1+6j), -1-6j)
+        self.assertComplexesAreIdentical(-0j, complex(0, -0.0))
+        self.assertComplexesAreIdentical(-complex(-0.0+0j), complex(0, -0.0))
 
     def test_getnewargs(self):
         self.assertEqual((1+2j).__getnewargs__(), (1.0, 2.0))
         self.assertEqual((1-2j).__getnewargs__(), (1.0, -2.0))
-        self.assertEqual((2j).__getnewargs__(), (0.0, 2.0))
-        self.assertEqual((-0j).__getnewargs__(), (0.0, -0.0))
+        self.assertEqual((0.0+2j).__getnewargs__(), (0.0, 2.0))
+        self.assertEqual((2j).__getnewargs__(), (2.0,))
+        self.assertEqual((0.0-0j).__getnewargs__(), (0.0, -0.0))
+        self.assertEqual((-0j).__getnewargs__(), (-0.0,))
         self.assertEqual(complex(0, INF).__getnewargs__(), (0.0, INF))
         self.assertEqual(complex(INF, 0).__getnewargs__(), (INF, 0.0))
 
@@ -837,15 +952,11 @@ def test_negated_imaginary_literal(self):
         z0 = -0j
         z1 = -7j
         z2 = -1e1000j
-        # Note: In versions of Python < 3.2, a negated imaginary literal
-        # accidentally ended up with real part 0.0 instead of -0.0, thanks to a
-        # modification during CST -> AST translation (see issue #9011).  That's
-        # fixed in Python 3.2.
-        self.assertFloatsAreIdentical(z0.real, -0.0)
+        self.assertFloatsAreIdentical(z0.real, +0.0)
         self.assertFloatsAreIdentical(z0.imag, -0.0)
-        self.assertFloatsAreIdentical(z1.real, -0.0)
+        self.assertFloatsAreIdentical(z1.real, +0.0)
         self.assertFloatsAreIdentical(z1.imag, -7.0)
-        self.assertFloatsAreIdentical(z2.real, -0.0)
+        self.assertFloatsAreIdentical(z2.real, +0.0)
         self.assertFloatsAreIdentical(z2.imag, -INF)
 
     @support.requires_IEEE_754
@@ -911,7 +1022,8 @@ def test_format(self):
         self.assertEqual(format(z, '3'), str(z))
 
         self.assertEqual(format(1+3j, 'g'), '1+3j')
-        self.assertEqual(format(3j, 'g'), '0+3j')
+        self.assertEqual(format(0+3j, 'g'), '0.0+3j')
+        self.assertEqual(format(3j, 'g'), '3j')
         self.assertEqual(format(1.5+3.5j, 'g'), '1.5+3.5j')
 
         self.assertEqual(format(1.5+3.5j, '+g'), '+1.5+3.5j')
diff --git a/Lib/test/test_format.py b/Lib/test/test_format.py
index 9dde63e40d06db..7d65dcd3e5cbe4 100644
--- a/Lib/test/test_format.py
+++ b/Lib/test/test_format.py
@@ -595,10 +595,15 @@ def test_negative_zero(self):
         self.assertEqual(f"{-1.:+z.0f}", "-1")
         self.assertEqual(f"{-1.:-z.0f}", "-1")
 
-        self.assertEqual(f"{0.j:z.1f}", "0.0+0.0j")
-        self.assertEqual(f"{-0.j:z.1f}", "0.0+0.0j")
-        self.assertEqual(f"{.01j:z.1f}", "0.0+0.0j")
-        self.assertEqual(f"{-.01j:z.1f}", "0.0+0.0j")
+        self.assertEqual(f"{0.0+0.j:z.1f}", "0.0+0.0j")
+        self.assertEqual(f"{0.0-0.j:z.1f}", "0.0+0.0j")
+        self.assertEqual(f"{0.0+.01j:z.1f}", "0.0+0.0j")
+        self.assertEqual(f"{0.0-.01j:z.1f}", "0.0+0.0j")
+
+        self.assertEqual(f"{0.j:z.1f}", "0.0j")
+        self.assertEqual(f"{-0.j:z.1f}", "0.0j")
+        self.assertEqual(f"{.01j:z.1f}", "0.0j")
+        self.assertEqual(f"{-.01j:z.1f}", "0.0j")
 
         self.assertEqual(f"{-0.:z>6.1f}", "zz-0.0")  # test fill, esp. 'z' fill
         self.assertEqual(f"{-0.:z>z6.1f}", "zzz0.0")
diff --git a/Lib/test/test_fractions.py b/Lib/test/test_fractions.py
index 98dccbec9566ac..dfbc289070adc7 100644
--- a/Lib/test/test_fractions.py
+++ b/Lib/test/test_fractions.py
@@ -1679,7 +1679,7 @@ def test_complex_handling(self):
         # See issue gh-102840 for more details.
 
         a = F(1, 2)
-        b = 1j
+        b = 0.0+1j
         message = "unsupported operand type(s) for %s: '%s' and '%s'"
         # test forward
         self.assertRaisesMessage(TypeError,
diff --git a/Lib/test/test_marshal.py b/Lib/test/test_marshal.py
index 93b8684c725d24..eec950cc5a2efb 100644
--- a/Lib/test/test_marshal.py
+++ b/Lib/test/test_marshal.py
@@ -318,7 +318,8 @@ def test_exact_type_match(self):
         #   >>> class Int(int): pass
         #   >>> type(loads(dumps(Int())))
         #   <type 'int'>
-        for typ in (int, float, complex, tuple, list, dict, set, frozenset):
+        for typ in (int, float, complex, tuple, list, dict, set, frozenset,
+                    imaginary):
             # Note: str subclasses are not tested because they get handled
             # by marshal's routines for objects supporting the buffer API.
             subtyp = type('subtyp', (typ,), {})
@@ -372,7 +373,7 @@ def readinto(self, buf):
                 if n is not None and n > 4:
                     n += 10**6
                 return n
-        for value in (1.0, 1j, b'0123456789', '0123456789'):
+        for value in (1.0, 1j, 1+1j, b'0123456789', '0123456789'):
             self.assertRaises(ValueError, marshal.load,
                               BadReader(marshal.dumps(value)))
 
@@ -628,7 +629,7 @@ def test_write_long_to_file(self):
             self.assertEqual(data, b'\x78\x56\x34\x12')
 
     def test_write_object_to_file(self):
-        obj = ('\u20ac', b'abc', 123, 45.6, 7+8j, 'long line '*1000)
+        obj = ('\u20ac', b'abc', 123, 45.6, 7+8j, 122j, 'long line '*1000)
         for v in range(marshal.version + 1):
             _testcapi.pymarshal_write_object_to_file(obj, os_helper.TESTFN, v)
             with open(os_helper.TESTFN, 'rb') as f:
@@ -665,7 +666,7 @@ def test_read_long_from_file(self):
         os_helper.unlink(os_helper.TESTFN)
 
     def test_read_last_object_from_file(self):
-        obj = ('\u20ac', b'abc', 123, 45.6, 7+8j)
+        obj = ('\u20ac', b'abc', 123, 45.6, 7+8j, 666j)
         for v in range(marshal.version + 1):
             data = marshal.dumps(obj, v)
             with open(os_helper.TESTFN, 'wb') as f:
@@ -681,7 +682,7 @@ def test_read_last_object_from_file(self):
             os_helper.unlink(os_helper.TESTFN)
 
     def test_read_object_from_file(self):
-        obj = ('\u20ac', b'abc', 123, 45.6, 7+8j)
+        obj = ('\u20ac', b'abc', 123, 45.6, 7+8j, 1j)
         for v in range(marshal.version + 1):
             data = marshal.dumps(obj, v)
             with open(os_helper.TESTFN, 'wb') as f:
diff --git a/Lib/test/test_socket.py b/Lib/test/test_socket.py
index 7b3914f30e5f52..a8b2dc906f748d 100644
--- a/Lib/test/test_socket.py
+++ b/Lib/test/test_socket.py
@@ -940,7 +940,7 @@ def testSendtoErrors(self):
         self.assertEqual(str(cm.exception),
                          "a bytes-like object is required, not 'str'")
         with self.assertRaises(TypeError) as cm:
-            s.sendto(5j, sockname)
+            s.sendto(1+5j, sockname)
         self.assertEqual(str(cm.exception),
                          "a bytes-like object is required, not 'complex'")
         with self.assertRaises(TypeError) as cm:
@@ -952,7 +952,7 @@ def testSendtoErrors(self):
         self.assertEqual(str(cm.exception),
                          "a bytes-like object is required, not 'str'")
         with self.assertRaises(TypeError) as cm:
-            s.sendto(5j, 0, sockname)
+            s.sendto(1+5j, 0, sockname)
         self.assertEqual(str(cm.exception),
                          "a bytes-like object is required, not 'complex'")
         with self.assertRaises(TypeError) as cm:
diff --git a/Lib/test/test_stable_abi_ctypes.py b/Lib/test/test_stable_abi_ctypes.py
index fa08dc6a25b0ea..0708e9c4ef71eb 100644
--- a/Lib/test/test_stable_abi_ctypes.py
+++ b/Lib/test/test_stable_abi_ctypes.py
@@ -341,6 +341,8 @@ def test_windows_feature_macros(self):
     "PyGILState_GetThisThreadState",
     "PyGILState_Release",
     "PyGetSetDescr_Type",
+    "PyImaginary_FromDouble",
+    "PyImaginary_Type",
     "PyImport_AddModule",
     "PyImport_AddModuleObject",
     "PyImport_AddModuleRef",
diff --git a/Lib/test/test_str.py b/Lib/test/test_str.py
index 4de6c1cba152bd..042da6a0f5266d 100644
--- a/Lib/test/test_str.py
+++ b/Lib/test/test_str.py
@@ -1572,12 +1572,12 @@ def __int__(self):
         self.assertRaisesRegex(TypeError, '%X format: an integer is required, not float', operator.mod, '%X', 2.11)
         self.assertRaisesRegex(TypeError, '%o format: an integer is required, not float', operator.mod, '%o', 1.79)
         self.assertRaisesRegex(TypeError, '%x format: an integer is required, not PseudoFloat', operator.mod, '%x', pi)
-        self.assertRaisesRegex(TypeError, '%x format: an integer is required, not complex', operator.mod, '%x', 3j)
-        self.assertRaisesRegex(TypeError, '%X format: an integer is required, not complex', operator.mod, '%X', 2j)
-        self.assertRaisesRegex(TypeError, '%o format: an integer is required, not complex', operator.mod, '%o', 1j)
-        self.assertRaisesRegex(TypeError, '%u format: a real number is required, not complex', operator.mod, '%u', 3j)
-        self.assertRaisesRegex(TypeError, '%i format: a real number is required, not complex', operator.mod, '%i', 2j)
-        self.assertRaisesRegex(TypeError, '%d format: a real number is required, not complex', operator.mod, '%d', 1j)
+        self.assertRaisesRegex(TypeError, '%x format: an integer is required, not complex', operator.mod, '%x', 1+3j)
+        self.assertRaisesRegex(TypeError, '%X format: an integer is required, not complex', operator.mod, '%X', 1+2j)
+        self.assertRaisesRegex(TypeError, '%o format: an integer is required, not complex', operator.mod, '%o', 1+1j)
+        self.assertRaisesRegex(TypeError, '%u format: a real number is required, not complex', operator.mod, '%u', 1+3j)
+        self.assertRaisesRegex(TypeError, '%i format: a real number is required, not complex', operator.mod, '%i', 1+2j)
+        self.assertRaisesRegex(TypeError, '%d format: a real number is required, not complex', operator.mod, '%d', 1+1j)
         self.assertRaisesRegex(TypeError, r'%c requires an int or a unicode character, not .*\.PseudoFloat', operator.mod, '%c', pi)
 
         class RaisingNumber:
diff --git a/Misc/stable_abi.toml b/Misc/stable_abi.toml
index f9e51f0683c965..491cc8dc15cc98 100644
--- a/Misc/stable_abi.toml
+++ b/Misc/stable_abi.toml
@@ -2506,6 +2506,11 @@
 [function.PyEval_GetFrameLocals]
     added = '3.13'
 
+[function.PyImaginary_FromDouble]
+    added = '3.14'
+[data.PyImaginary_Type]
+    added = '3.14'
+
 [function.Py_TYPE]
     added = '3.14'
 [function.Py_REFCNT]
diff --git a/Modules/_testlimitedcapi/complex.c b/Modules/_testlimitedcapi/complex.c
index e4c244e5c88d06..a71502ad9959b2 100644
--- a/Modules/_testlimitedcapi/complex.c
+++ b/Modules/_testlimitedcapi/complex.c
@@ -16,6 +16,20 @@ complex_checkexact(PyObject *Py_UNUSED(module), PyObject *obj)
     return PyLong_FromLong(PyComplex_CheckExact(obj));
 }
 
+static PyObject *
+imaginary_check(PyObject *Py_UNUSED(module), PyObject *obj)
+{
+    NULLABLE(obj);
+    return PyLong_FromLong(PyImaginary_Check(obj));
+}
+
+static PyObject *
+imaginary_checkexact(PyObject *Py_UNUSED(module), PyObject *obj)
+{
+    NULLABLE(obj);
+    return PyLong_FromLong(PyImaginary_CheckExact(obj));
+}
+
 static PyObject *
 complex_fromdoubles(PyObject *Py_UNUSED(module), PyObject *args)
 {
@@ -28,6 +42,12 @@ complex_fromdoubles(PyObject *Py_UNUSED(module), PyObject *args)
     return PyComplex_FromDoubles(real, imag);
 }
 
+static PyObject *
+imaginary_fromdouble(PyObject *Py_UNUSED(module), PyObject *obj)
+{
+    return PyImaginary_FromDouble(PyFloat_AsDouble(obj));
+}
+
 static PyObject *
 complex_realasdouble(PyObject *Py_UNUSED(module), PyObject *obj)
 {
@@ -65,6 +85,9 @@ static PyMethodDef test_methods[] = {
     {"complex_fromdoubles", complex_fromdoubles, METH_VARARGS},
     {"complex_realasdouble", complex_realasdouble, METH_O},
     {"complex_imagasdouble", complex_imagasdouble, METH_O},
+    {"imaginary_check", imaginary_check, METH_O},
+    {"imaginary_checkexact", imaginary_checkexact, METH_O},
+    {"imaginary_fromdouble", imaginary_fromdouble, METH_O},
     {NULL},
 };
 
diff --git a/Modules/cmathmodule.c b/Modules/cmathmodule.c
index 81cbf0d554de3c..cbd4b8ef451f55 100644
--- a/Modules/cmathmodule.c
+++ b/Modules/cmathmodule.c
@@ -1188,15 +1188,13 @@ cmath_exec(PyObject *mod)
         return -1;
     }
 
-    Py_complex infj = {0.0, Py_INFINITY};
-    if (PyModule_Add(mod, "infj", PyComplex_FromCComplex(infj)) < 0) {
+    if (PyModule_Add(mod, "infj", PyImaginary_FromDouble(Py_INFINITY)) < 0) {
         return -1;
     }
     if (PyModule_Add(mod, "nan", PyFloat_FromDouble(fabs(Py_NAN))) < 0) {
         return -1;
     }
-    Py_complex nanj = {0.0, fabs(Py_NAN)};
-    if (PyModule_Add(mod, "nanj", PyComplex_FromCComplex(nanj)) < 0) {
+    if (PyModule_Add(mod, "nanj", PyImaginary_FromDouble(Py_NAN)) < 0) {
         return -1;
     }
 
diff --git a/Objects/clinic/complexobject.c.h b/Objects/clinic/complexobject.c.h
index 3c3d1071b6eec3..259a5bec24f1e7 100644
--- a/Objects/clinic/complexobject.c.h
+++ b/Objects/clinic/complexobject.c.h
@@ -170,4 +170,101 @@ PyDoc_STRVAR(complex_from_number__doc__,
 
 #define COMPLEX_FROM_NUMBER_METHODDEF    \
     {"from_number", (PyCFunction)complex_from_number, METH_O|METH_CLASS, complex_from_number__doc__},
-/*[clinic end generated code: output=8c49a41c5a7f0aee input=a9049054013a1b77]*/
+
+PyDoc_STRVAR(imaginary_new__doc__,
+"imaginary(x=0)\n"
+"--\n"
+"\n"
+"Create an imaginary number from a real number or string.\n"
+"\n"
+"This is equivalent of float(x)*1j.");
+
+static PyObject *
+imaginary_new_impl(PyTypeObject *type, PyObject *x);
+
+static PyObject *
+imaginary_new(PyTypeObject *type, PyObject *args, PyObject *kwargs)
+{
+    PyObject *return_value = NULL;
+    #if defined(Py_BUILD_CORE) && !defined(Py_BUILD_CORE_MODULE)
+
+    #define NUM_KEYWORDS 1
+    static struct {
+        PyGC_Head _this_is_not_used;
+        PyObject_VAR_HEAD
+        PyObject *ob_item[NUM_KEYWORDS];
+    } _kwtuple = {
+        .ob_base = PyVarObject_HEAD_INIT(&PyTuple_Type, NUM_KEYWORDS)
+        .ob_item = { _Py_LATIN1_CHR('x'), },
+    };
+    #undef NUM_KEYWORDS
+    #define KWTUPLE (&_kwtuple.ob_base.ob_base)
+
+    #else  // !Py_BUILD_CORE
+    #  define KWTUPLE NULL
+    #endif  // !Py_BUILD_CORE
+
+    static const char * const _keywords[] = {"x", NULL};
+    static _PyArg_Parser _parser = {
+        .keywords = _keywords,
+        .fname = "imaginary",
+        .kwtuple = KWTUPLE,
+    };
+    #undef KWTUPLE
+    PyObject *argsbuf[1];
+    PyObject * const *fastargs;
+    Py_ssize_t nargs = PyTuple_GET_SIZE(args);
+    Py_ssize_t noptargs = nargs + (kwargs ? PyDict_GET_SIZE(kwargs) : 0) - 0;
+    PyObject *x = NULL;
+
+    fastargs = _PyArg_UnpackKeywords(_PyTuple_CAST(args)->ob_item, nargs, kwargs, NULL, &_parser,
+            /*minpos*/ 0, /*maxpos*/ 1, /*minkw*/ 0, /*varpos*/ 0, argsbuf);
+    if (!fastargs) {
+        goto exit;
+    }
+    if (!noptargs) {
+        goto skip_optional_pos;
+    }
+    x = fastargs[0];
+skip_optional_pos:
+    return_value = imaginary_new_impl(type, x);
+
+exit:
+    return return_value;
+}
+
+PyDoc_STRVAR(imaginary_conjugate__doc__,
+"conjugate($self, /)\n"
+"--\n"
+"\n"
+"Return the complex conjugate of its argument. (-4j).conjugate() == 4j.");
+
+#define IMAGINARY_CONJUGATE_METHODDEF    \
+    {"conjugate", (PyCFunction)imaginary_conjugate, METH_NOARGS, imaginary_conjugate__doc__},
+
+static PyObject *
+imaginary_conjugate_impl(PyComplexObject *self);
+
+static PyObject *
+imaginary_conjugate(PyComplexObject *self, PyObject *Py_UNUSED(ignored))
+{
+    return imaginary_conjugate_impl(self);
+}
+
+PyDoc_STRVAR(imaginary___getnewargs____doc__,
+"__getnewargs__($self, /)\n"
+"--\n"
+"\n");
+
+#define IMAGINARY___GETNEWARGS___METHODDEF    \
+    {"__getnewargs__", (PyCFunction)imaginary___getnewargs__, METH_NOARGS, imaginary___getnewargs____doc__},
+
+static PyObject *
+imaginary___getnewargs___impl(PyComplexObject *self);
+
+static PyObject *
+imaginary___getnewargs__(PyComplexObject *self, PyObject *Py_UNUSED(ignored))
+{
+    return imaginary___getnewargs___impl(self);
+}
+/*[clinic end generated code: output=90feb509dccfb970 input=a9049054013a1b77]*/
diff --git a/Objects/codeobject.c b/Objects/codeobject.c
index 148350cc4b9195..adfcdcd82e2055 100644
--- a/Objects/codeobject.c
+++ b/Objects/codeobject.c
@@ -2434,6 +2434,15 @@ _PyCode_ConstantKey(PyObject *op)
         else
             key = PyTuple_Pack(2, Py_TYPE(op), op);
     }
+    else if (PyImaginary_CheckExact(op)) {
+        double d = ((PyComplexObject *)(op))->cval.imag;
+        if (d == 0.0 && copysign(1.0, d) < 0.0) {
+            key = PyTuple_Pack(3, Py_TYPE(op), op, Py_None);
+        }
+        else {
+            key = PyTuple_Pack(2, Py_TYPE(op), op);
+        }
+    }
     else if (PyComplex_CheckExact(op)) {
         Py_complex z;
         int real_negzero, imag_negzero;
@@ -2659,6 +2668,11 @@ compare_constants(const void *key1, const void *key2) {
         Py_complex c2 = ((PyComplexObject *)op2)->cval;
         return memcmp(&c1, &c2, sizeof(Py_complex)) == 0;
     }
+    else if (PyImaginary_CheckExact(op1)) {
+        double i1 = ((PyComplexObject *)op1)->cval.imag;
+        double i2 = ((PyComplexObject *)op2)->cval.imag;
+        return memcmp(&i1, &i2, sizeof(double)) == 0;
+    }
     _Py_FatalErrorFormat("unexpected type in compare_constants: %s",
                          Py_TYPE(op1)->tp_name);
     return 0;
diff --git a/Objects/complexobject.c b/Objects/complexobject.c
index 8fbca3cb02d80a..8ef0f6e65c7bfc 100644
--- a/Objects/complexobject.c
+++ b/Objects/complexobject.c
@@ -17,8 +17,9 @@
 
 /*[clinic input]
 class complex "PyComplexObject *" "&PyComplex_Type"
+class imaginary "PyComplexObject *" "&PyComplex_Type"
 [clinic start generated code]*/
-/*[clinic end generated code: output=da39a3ee5e6b4b0d input=819e057d2d10f5ec]*/
+/*[clinic end generated code: output=da39a3ee5e6b4b0d input=041bba3f29a299d0]*/
 
 #include "clinic/complexobject.c.h"
 
@@ -518,7 +519,7 @@ complex_repr(PyComplexObject *v)
     const char *lead = "";
     const char *tail = "";
 
-    if (v->cval.real == 0. && copysign(1.0, v->cval.real)==1.0) {
+    if (PyImaginary_Check(v)) {
         /* Real part is +0: just output the imaginary part and do not
            include parens. */
         re = "";
@@ -532,7 +533,8 @@ complex_repr(PyComplexObject *v)
         /* Format imaginary part with sign, real part without. Include
            parens in the result. */
         pre = PyOS_double_to_string(v->cval.real, format_code,
-                                    precision, 0, NULL);
+                                    precision,
+                                    v->cval.real? 0 : Py_DTSF_ADD_DOT_0, NULL);
         if (!pre) {
             PyErr_NoMemory();
             goto done;
@@ -607,19 +609,29 @@ real_to_complex(PyObject **pobj, Py_complex *pc)
 }
 
 /* Complex arithmetic rules implement special mixed-mode case where combining
-   a pure-real (float or int) value and a complex value is performed directly
-   without first coercing the real value to a complex value.
+   a pure-real (float or int) value and a complex value or a pure-imaginary
+   value (instances of the imaginary type, like 1j) and a complex value are
+   performed directly, without first coercing operands to a complex value.
 
    Let us consider the addition as an example, assuming that ints are implicitly
    converted to floats.  We have the following rules (up to variants with changed
    order of operands):
 
-       complex(a, b) + complex(c, d) = complex(a + c, b + d)
-       float(a) + complex(b, c) = complex(a + b, c)
+       complex(a, b) + complex(c, d) = complex(a + c, b + d)   (1)
+       float(a) + complex(b, c) = complex(a + b, c)            (2)
+       imaginary(a) + complex(b, c) = complex(b, a + c)        (3)
+       imaginary(a) + imaginary(b) = imaginary(a + b)          (4)
+       float(a) + imaginary(b) = complex(a, b)                 (5)
 
    Similar rules are implemented for subtraction, multiplication and division.
    See C11's Annex G, sections G.5.1 and G.5.2.
- */
+
+   Note, that the imaginary type implemented as a subtype of the complex
+   type.  So, complex_op() functions below must implement only one type of
+   mixed-mode operations, i.e.  when one argument is a float (2).  The rest
+   (respectively: complex op imaginary (3), imaginary op imaginary (4) and
+   imaginary op float (5)) should be implemented by an appropriate
+   imaginary_op() function. */
 
 #define COMPLEX_BINOP(NAME, FUNC)                           \
     static PyObject *                                       \
@@ -1364,3 +1376,268 @@ PyTypeObject PyComplex_Type = {
     PyObject_Free,                              /* tp_free */
     .tp_version_tag = _Py_TYPE_VERSION_COMPLEX,
 };
+
+/*[clinic input]
+@classmethod
+imaginary.__new__ as imaginary_new
+    x: object(c_default="NULL") = 0
+
+Create an imaginary number from a real number or string.
+
+This is equivalent of float(x)*1j.
+[clinic start generated code]*/
+
+static PyObject *
+imaginary_new_impl(PyTypeObject *type, PyObject *x)
+/*[clinic end generated code: output=07242ea06eb219f6 input=0ec29c2535687795]*/
+{
+    if (x == NULL) {
+        x = _PyLong_GetZero();
+    }
+
+    PyNumberMethods *nbx = Py_TYPE(x)->tp_as_number;
+
+    if (nbx == NULL || (nbx->nb_float == NULL && nbx->nb_index == NULL
+                        && !PyFloat_Check(x) && !PyUnicode_Check(x)))
+    {
+        PyErr_Format(PyExc_TypeError,
+                     "imaginary() first argument must be a real number, "
+                     "not '%.200s'", Py_TYPE(x)->tp_name);
+            return NULL;
+    }
+
+    PyObject* tmp = PyNumber_Float(x);
+
+    if (tmp == NULL) {
+        return NULL;
+    }
+
+    PyObject *ret = type->tp_alloc(type, 0);
+
+    if (ret != NULL) {
+        ((PyComplexObject *)ret)->cval = (Py_complex) {0.0, PyFloat_AS_DOUBLE(tmp)};
+    }
+    Py_DECREF(tmp);
+    return ret;
+}
+
+PyObject *
+PyImaginary_FromDouble(double imag)
+{
+    /* Inline PyObject_New */
+    PyComplexObject *op = PyObject_Malloc(sizeof(PyComplexObject));
+
+    if (op == NULL) {
+        return PyErr_NoMemory();
+    }
+    _PyObject_Init((PyObject*)op, &PyImaginary_Type);
+    op->cval = (Py_complex){0.0, imag};
+    return (PyObject *) op;
+}
+
+static PyObject *
+imaginary_neg(PyComplexObject *v)
+{
+    return PyImaginary_FromDouble(-v->cval.imag);
+}
+
+static PyObject *
+imaginary_pos(PyComplexObject *v)
+{
+    if (PyImaginary_CheckExact(v)) {
+        return Py_NewRef(v);
+    }
+    return PyImaginary_FromDouble(v->cval.imag);
+}
+
+/* Imaginary type arithmetic.  Binary operations below
+   must support also complex and float (or int) operands. */
+
+static PyObject *
+imaginary_add(PyObject *v, PyObject *w)
+{
+    if (!PyImaginary_Check(v)) {
+        PyObject *tmp = v;
+
+        v = w;
+        w = tmp;
+    }
+
+    double a = ((PyComplexObject *)(v))->cval.imag;
+    double b;
+
+    if (PyComplex_Check(w)) {
+        a += ((PyComplexObject *)(w))->cval.imag;
+        b = ((PyComplexObject *)(w))->cval.real;
+    }
+    else if (real_to_double(&w, &b) < 0) {
+        return w;
+    }
+    return PyComplex_FromDoubles(b, a);
+}
+
+static PyObject *
+imaginary_sub(PyObject *v, PyObject *w)
+{
+    if (PyImaginary_Check(w)) {
+        double b = -((PyComplexObject *)(w))->cval.imag;
+        double a;
+
+        if (PyImaginary_Check(v)) {
+            b += ((PyComplexObject *)(v))->cval.imag;
+            return PyImaginary_FromDouble(b);
+        }
+        if (PyComplex_Check(v)) {
+            a = ((PyComplexObject *)(v))->cval.real;
+            b += ((PyComplexObject *)(v))->cval.imag;
+        }
+        else if (real_to_double(&v, &a) < 0) {
+            return v;
+        }
+        return PyComplex_FromDoubles(a, b);
+    }
+
+    double a = ((PyComplexObject *)(v))->cval.imag;
+    double b;
+
+    if (PyComplex_Check(w)) {
+        a -= ((PyComplexObject *)(w))->cval.imag;
+        b = ((PyComplexObject *)(w))->cval.real;
+    }
+    else if (real_to_double(&w, &b) < 0) {
+        return w;
+    }
+    return PyComplex_FromDoubles(-b, a);
+}
+
+static PyObject *
+imaginary_mul(PyObject *v, PyObject *w)
+{
+    if (!PyImaginary_Check(v)) {
+        PyObject *tmp = v;
+
+        v = w;
+        w = tmp;
+    }
+
+    double a = ((PyComplexObject *)(v))->cval.imag;
+    double b;
+
+    if (PyComplex_Check(w)) {
+        b = -a * ((PyComplexObject *)(w))->cval.imag;
+        if (PyImaginary_Check(w)) {
+            return PyFloat_FromDouble(b);
+        }
+        a *= ((PyComplexObject *)(w))->cval.real;
+        return PyComplex_FromDoubles(b, a);
+    }
+    else if (real_to_double(&w, &b) < 0) {
+        return w;
+    }
+    return PyImaginary_FromDouble(a*b);
+}
+
+static PyObject *
+imaginary_div(PyObject *v, PyObject *w)
+{
+    if (PyImaginary_Check(w)) {
+        double b = ((PyComplexObject *)(w))->cval.imag;
+        double a;
+
+        if (b) {
+            if (PyImaginary_Check(v)) {
+                double a = ((PyComplexObject *)(v))->cval.imag;
+
+                return PyFloat_FromDouble(a/b);
+            }
+            if (PyComplex_Check(v)) {
+                Py_complex a = ((PyComplexObject *)(v))->cval;
+
+                return PyComplex_FromDoubles(a.imag/b, -a.real/b);
+            }
+            if (real_to_double(&v, &a) < 0) {
+                return v;
+            }
+            return PyImaginary_FromDouble(-a/b);
+        }
+    }
+    else {
+        double a = ((PyComplexObject *)(v))->cval.imag;
+        double b;
+
+        if (PyComplex_Check(w)) {
+            Py_complex b = ((PyComplexObject *)(w))->cval;
+
+            errno = 0;
+            b = _Py_rc_quot(a, b);
+            b = (Py_complex){-b.imag, b.real};
+            if (!errno) {
+                return PyComplex_FromDoubles(b.real, b.imag);
+            }
+        }
+        else if (real_to_double(&w, &b) < 0) {
+            return w;
+        }
+        else if (b) {
+            return PyImaginary_FromDouble(a/b);
+        }
+    }
+    PyErr_SetString(PyExc_ZeroDivisionError, "complex division by zero");
+    return NULL;
+}
+
+/*[clinic input]
+imaginary.conjugate
+
+Return the complex conjugate of its argument. (-4j).conjugate() == 4j.
+[clinic start generated code]*/
+
+static PyObject *
+imaginary_conjugate_impl(PyComplexObject *self)
+/*[clinic end generated code: output=247bb3742efd769d input=8dcd1c93a873c492]*/
+{
+    Py_complex c = self->cval;
+
+    c.imag = -c.imag;
+    return PyImaginary_FromDouble(c.imag);
+}
+
+/*[clinic input]
+imaginary.__getnewargs__
+
+[clinic start generated code]*/
+
+static PyObject *
+imaginary___getnewargs___impl(PyComplexObject *self)
+/*[clinic end generated code: output=a587156eda821f7d input=da4b7e53915987d5]*/
+{
+    Py_complex c = self->cval;
+
+    return Py_BuildValue("(d)", c.imag);
+}
+
+static PyMethodDef imaginary_methods[] = {
+    IMAGINARY_CONJUGATE_METHODDEF
+    IMAGINARY___GETNEWARGS___METHODDEF
+    {NULL}  /* sentinel */
+};
+
+static PyNumberMethods imaginary_as_number = {
+    .nb_add = (binaryfunc)imaginary_add,
+    .nb_subtract = (binaryfunc)imaginary_sub,
+    .nb_multiply = (binaryfunc)imaginary_mul,
+    .nb_true_divide = (binaryfunc)imaginary_div,
+    .nb_negative = (unaryfunc)imaginary_neg,
+    .nb_positive = (unaryfunc)imaginary_pos,
+};
+
+PyTypeObject PyImaginary_Type = {
+    PyVarObject_HEAD_INIT(NULL, 0)
+    .tp_name = "imaginary",
+    .tp_as_number = &imaginary_as_number,
+    .tp_flags = Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,
+    .tp_doc = imaginary_new__doc__,
+    .tp_base = &PyComplex_Type,
+    .tp_new = imaginary_new,
+    .tp_methods = imaginary_methods,
+};
diff --git a/Objects/object.c b/Objects/object.c
index 8868fa29066404..95b5ade6d03b3c 100644
--- a/Objects/object.c
+++ b/Objects/object.c
@@ -2381,6 +2381,7 @@ static PyTypeObject* static_types[] = {
     // subclasses: _PyTypes_FiniTypes() deallocates them before their base
     // class
     &PyBool_Type,         // base=&PyLong_Type
+    &PyImaginary_Type,    // base=&PyComplex_Type
     &PyCMethod_Type,      // base=&PyCFunction_Type
     &PyODictItems_Type,   // base=&PyDictItems_Type
     &PyODictKeys_Type,    // base=&PyDictKeys_Type
diff --git a/PC/python3dll.c b/PC/python3dll.c
index 8657ddb9fa5155..bebe4bbe14fbc8 100755
--- a/PC/python3dll.c
+++ b/PC/python3dll.c
@@ -297,6 +297,7 @@ EXPORT_FUNC(PyGC_IsEnabled)
 EXPORT_FUNC(PyGILState_Ensure)
 EXPORT_FUNC(PyGILState_GetThisThreadState)
 EXPORT_FUNC(PyGILState_Release)
+EXPORT_FUNC(PyImaginary_FromDouble)
 EXPORT_FUNC(PyImport_AddModule)
 EXPORT_FUNC(PyImport_AddModuleObject)
 EXPORT_FUNC(PyImport_AddModuleRef)
@@ -900,6 +901,7 @@ EXPORT_DATA(PyFilter_Type)
 EXPORT_DATA(PyFloat_Type)
 EXPORT_DATA(PyFrozenSet_Type)
 EXPORT_DATA(PyGetSetDescr_Type)
+EXPORT_DATA(PyImaginary_Type)
 EXPORT_DATA(PyList_Type)
 EXPORT_DATA(PyListIter_Type)
 EXPORT_DATA(PyListRevIter_Type)
diff --git a/Parser/action_helpers.c b/Parser/action_helpers.c
index 5ac1dd7813689c..2063e623addf44 100644
--- a/Parser/action_helpers.c
+++ b/Parser/action_helpers.c
@@ -837,7 +837,7 @@ _PyPegen_seq_delete_starred_exprs(Parser *p, asdl_seq *kwargs)
 expr_ty
 _PyPegen_ensure_imaginary(Parser *p, expr_ty exp)
 {
-    if (exp->kind != Constant_kind || !PyComplex_CheckExact(exp->v.Constant.value)) {
+    if (exp->kind != Constant_kind || !PyImaginary_CheckExact(exp->v.Constant.value)) {
         RAISE_SYNTAX_ERROR_KNOWN_LOCATION(exp, "imaginary number required in complex literal");
         return NULL;
     }
@@ -847,7 +847,8 @@ _PyPegen_ensure_imaginary(Parser *p, expr_ty exp)
 expr_ty
 _PyPegen_ensure_real(Parser *p, expr_ty exp)
 {
-    if (exp->kind != Constant_kind || PyComplex_CheckExact(exp->v.Constant.value)) {
+    if (exp->kind != Constant_kind || (PyComplex_CheckExact(exp->v.Constant.value)
+                                       || PyImaginary_CheckExact(exp->v.Constant.value))) {
         RAISE_SYNTAX_ERROR_KNOWN_LOCATION(exp, "real number required in complex literal");
         return NULL;
     }
diff --git a/Parser/pegen.c b/Parser/pegen.c
index bb98e7b184a4dc..e04c2afa09e073 100644
--- a/Parser/pegen.c
+++ b/Parser/pegen.c
@@ -614,7 +614,7 @@ parsenumber_raw(const char *s)
     const char *end;
     long x;
     double dx;
-    Py_complex compl;
+    double imag;
     int imflag;
 
     assert(s != NULL);
@@ -638,12 +638,11 @@ parsenumber_raw(const char *s)
     }
     /* XXX Huge floats may silently fail */
     if (imflag) {
-        compl.real = 0.;
-        compl.imag = PyOS_string_to_double(s, (char **)&end, NULL);
-        if (compl.imag == -1.0 && PyErr_Occurred()) {
+        imag = PyOS_string_to_double(s, (char **)&end, NULL);
+        if (imag == -1.0 && PyErr_Occurred()) {
             return NULL;
         }
-        return PyComplex_FromCComplex(compl);
+        return PyImaginary_FromDouble(imag);
     }
     dx = PyOS_string_to_double(s, NULL, NULL);
     if (dx == -1.0 && PyErr_Occurred()) {
diff --git a/Python/ast.c b/Python/ast.c
index bf1ff5f3ec18ba..e95474141c870c 100644
--- a/Python/ast.c
+++ b/Python/ast.c
@@ -174,6 +174,7 @@ validate_constant(struct validator *state, PyObject *value)
     if (PyLong_CheckExact(value)
             || PyFloat_CheckExact(value)
             || PyComplex_CheckExact(value)
+            || PyImaginary_CheckExact(value)
             || PyBool_Check(value)
             || PyUnicode_CheckExact(value)
             || PyBytes_CheckExact(value))
@@ -419,7 +420,7 @@ ensure_literal_number(expr_ty exp, bool allow_real, bool allow_imaginary)
     PyObject *value = exp->v.Constant.value;
     return (allow_real && PyFloat_CheckExact(value)) ||
            (allow_real && PyLong_CheckExact(value)) ||
-           (allow_imaginary && PyComplex_CheckExact(value));
+           (allow_imaginary && PyImaginary_CheckExact(value));
 }
 
 static int
@@ -499,7 +500,7 @@ validate_pattern_match_value(struct validator *state, expr_ty exp)
             PyObject *literal = exp->v.Constant.value;
             if (PyLong_CheckExact(literal) || PyFloat_CheckExact(literal) ||
                 PyBytes_CheckExact(literal) || PyComplex_CheckExact(literal) ||
-                PyUnicode_CheckExact(literal)) {
+                PyImaginary_CheckExact(literal) || PyUnicode_CheckExact(literal)) {
                 return 1;
             }
             PyErr_SetString(PyExc_ValueError,
diff --git a/Python/ast_unparse.c b/Python/ast_unparse.c
index 8017cfc7fcf268..a97c24c15f8908 100644
--- a/Python/ast_unparse.c
+++ b/Python/ast_unparse.c
@@ -73,7 +73,7 @@ append_repr(_PyUnicodeWriter *writer, PyObject *obj)
     }
 
     if ((PyFloat_CheckExact(obj) && isinf(PyFloat_AS_DOUBLE(obj))) ||
-        PyComplex_CheckExact(obj))
+        PyComplex_CheckExact(obj) || PyImaginary_CheckExact(obj))
     {
         _Py_DECLARE_STR(str_replace_inf, "1e309");  // evaluates to inf
         PyObject *new_repr = PyUnicode_Replace(
diff --git a/Python/bltinmodule.c b/Python/bltinmodule.c
index 17df9208f224f4..e8e42cc2f84692 100644
--- a/Python/bltinmodule.c
+++ b/Python/bltinmodule.c
@@ -2818,6 +2818,11 @@ builtin_sum_impl(PyObject *module, PyObject *iterable, PyObject *start)
                 return PyComplex_FromDoubles(cs_to_double(re_sum),
                                              cs_to_double(im_sum));
             }
+            if (PyImaginary_CheckExact(item)) {
+                double value = PyComplex_ImagAsDouble(item);
+                im_sum = cs_add(im_sum, value);
+                continue;
+            }
             if (PyComplex_CheckExact(item)) {
                 z = PyComplex_AsCComplex(item);
                 re_sum = cs_add(re_sum, z.real);
@@ -3345,6 +3350,7 @@ _PyBuiltin_Init(PyInterpreterState *interp)
     SETBUILTIN("float",                 &PyFloat_Type);
     SETBUILTIN("frozenset",             &PyFrozenSet_Type);
     SETBUILTIN("property",              &PyProperty_Type);
+    SETBUILTIN("imaginary",             &PyImaginary_Type);
     SETBUILTIN("int",                   &PyLong_Type);
     SETBUILTIN("list",                  &PyList_Type);
     SETBUILTIN("map",                   &PyMap_Type);
diff --git a/Python/formatter_unicode.c b/Python/formatter_unicode.c
index 16f711184990ac..2bd130a3c66bfe 100644
--- a/Python/formatter_unicode.c
+++ b/Python/formatter_unicode.c
@@ -1271,7 +1271,7 @@ format_complex_internal(PyObject *value,
         /* Omitted type specifier. Should be like str(self). */
         type = 'r';
         default_precision = 0;
-        if (re == 0.0 && copysign(1.0, re) == 1.0)
+        if (PyImaginary_Check(value))
             skip_re = 1;
         else
             add_parens = 1;
@@ -1290,8 +1290,17 @@ format_complex_internal(PyObject *value,
     /* Cast "type", because if we're in unicode we need to pass an
        8-bit char. This is safe, because we've restricted what "type"
        can be. */
-    re_buf = PyOS_double_to_string(re, (char)type, precision, flags,
-                                   &re_float_type);
+    if (re == 0.0) {
+        if (PyImaginary_Check(value)) {
+            skip_re = 1;
+        }
+        re_buf = PyOS_double_to_string(re, (char)type, precision,
+                                       flags | Py_DTSF_ADD_DOT_0,
+                                       &re_float_type);
+    }
+    else
+        re_buf = PyOS_double_to_string(re, (char)type, precision, flags,
+                                       &re_float_type);
     if (re_buf == NULL)
         goto done;
     im_buf = PyOS_double_to_string(im, (char)type, precision, flags,
diff --git a/Python/marshal.c b/Python/marshal.c
index 72afa4ff89432c..965dc808a701b4 100644
--- a/Python/marshal.c
+++ b/Python/marshal.c
@@ -59,6 +59,8 @@ module marshal
 #define TYPE_ELLIPSIS           '.'
 #define TYPE_BINARY_FLOAT       'g'  // Version 0 uses TYPE_FLOAT instead.
 #define TYPE_BINARY_COMPLEX     'y'  // Version 0 uses TYPE_COMPLEX instead.
+#define TYPE_IMAGINARY          'j'
+#define TYPE_BINARY_IMAGINARY   'J'
 #define TYPE_LONG               'l'  // See also TYPE_INT.
 #define TYPE_STRING             's'  // Bytes. (Name comes from Python 2.)
 #define TYPE_TUPLE              '('  // See also TYPE_SMALL_TUPLE.
@@ -451,6 +453,16 @@ w_complex_object(PyObject *v, char flag, WFILE *p)
             w_float_str(PyComplex_ImagAsDouble(v), p);
         }
     }
+    else if (PyImaginary_CheckExact(v)) {
+        if (p->version > 1) {
+            W_TYPE(TYPE_BINARY_IMAGINARY, p);
+            w_float_bin(PyComplex_ImagAsDouble(v), p);
+        }
+        else {
+            W_TYPE(TYPE_IMAGINARY, p);
+            w_float_str(PyComplex_ImagAsDouble(v), p);
+        }
+    }
     else if (PyBytes_CheckExact(v)) {
         W_TYPE(TYPE_STRING, p);
         w_pstring(PyBytes_AS_STRING(v), PyBytes_GET_SIZE(v), p);
@@ -1153,6 +1165,26 @@ r_object(RFILE *p)
             break;
         }
 
+    case TYPE_IMAGINARY:
+        {
+            double i = r_float_str(p);
+            if (i == -1.0 && PyErr_Occurred())
+                break;
+            retval = PyImaginary_FromDouble(i);
+            R_REF(retval);
+            break;
+        }
+
+    case TYPE_BINARY_IMAGINARY:
+        {
+            double i = r_float_bin(p);
+            if (i == -1.0 && PyErr_Occurred())
+                break;
+            retval = PyImaginary_FromDouble(i);
+            R_REF(retval);
+            break;
+        }
+
     case TYPE_STRING:
         {
             const char *ptr;
diff --git a/Tools/c-analyzer/cpython/globals-to-fix.tsv b/Tools/c-analyzer/cpython/globals-to-fix.tsv
index badd7b79102310..6dfe4cc110a643 100644
--- a/Tools/c-analyzer/cpython/globals-to-fix.tsv
+++ b/Tools/c-analyzer/cpython/globals-to-fix.tsv
@@ -20,6 +20,7 @@ Objects/classobject.c	-	PyInstanceMethod_Type	-
 Objects/classobject.c	-	PyMethod_Type	-
 Objects/codeobject.c	-	PyCode_Type	-
 Objects/complexobject.c	-	PyComplex_Type	-
+Objects/complexobject.c	-	PyImaginary_Type	-
 Objects/descrobject.c	-	PyClassMethodDescr_Type	-
 Objects/descrobject.c	-	PyDictProxy_Type	-
 Objects/descrobject.c	-	PyGetSetDescr_Type	-

From c0dcb3a639147324cb2279b7887f63f4fe078dc2 Mon Sep 17 00:00:00 2001
From: Sergey B Kirpichev <skirpichev@gmail.com>
Date: Sat, 30 Nov 2024 06:12:15 +0300
Subject: [PATCH 2/2] + _Py_ci/ic* functions

---
 Doc/c-api/complex.rst              |  54 ++++++++
 Include/cpython/complexobject.h    |   6 +
 Lib/test/test_capi/test_complex.py |  42 ++++++
 Modules/_testcapi/complex.c        |  80 +++++++-----
 Objects/complexobject.c            | 203 ++++++++++++++++++-----------
 5 files changed, 272 insertions(+), 113 deletions(-)

diff --git a/Doc/c-api/complex.rst b/Doc/c-api/complex.rst
index d1f5d8eda676ef..3d4b3bbe99c069 100644
--- a/Doc/c-api/complex.rst
+++ b/Doc/c-api/complex.rst
@@ -52,6 +52,14 @@ pointers.  This is consistent throughout the API.
    .. versionadded:: 3.14
 
 
+.. 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
@@ -74,6 +82,22 @@ pointers.  This is consistent throughout the API.
    .. versionadded:: 3.14
 
 
+.. 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
@@ -94,6 +118,14 @@ pointers.  This is consistent throughout the API.
    .. versionadded:: 3.14
 
 
+.. 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`
@@ -125,6 +157,28 @@ pointers.  This is consistent throughout the API.
    .. versionadded:: 3.14
 
 
+.. 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`
diff --git a/Include/cpython/complexobject.h b/Include/cpython/complexobject.h
index 28576afad0b6b5..492bd8dfa7c165 100644
--- a/Include/cpython/complexobject.h
+++ b/Include/cpython/complexobject.h
@@ -10,15 +10,21 @@ typedef struct {
 // Operations on complex numbers.
 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);
 
diff --git a/Lib/test/test_capi/test_complex.py b/Lib/test/test_capi/test_complex.py
index 1f30fd24307d99..3fa59f109f4ee5 100644
--- a/Lib/test/test_capi/test_complex.py
+++ b/Lib/test/test_capi/test_complex.py
@@ -182,6 +182,13 @@ def test_py_cr_sum(self):
         self.assertComplexesAreIdentical(_py_cr_sum(-0.0 - 0j, -0.0)[0],
                                          complex(-0.0, -0.0))
 
+    def test_py_ci_sum(self):
+        # Test _Py_cr_sum()
+        _py_ci_sum = _testcapi._py_ci_sum
+
+        self.assertComplexesAreIdentical(_py_ci_sum(-0.0 - 0j, -0.0)[0],
+                                         complex(-0.0, -0.0))
+
     def test_py_c_diff(self):
         # Test _Py_c_diff()
         _py_c_diff = _testcapi._py_c_diff
@@ -195,6 +202,13 @@ def test_py_cr_diff(self):
         self.assertComplexesAreIdentical(_py_cr_diff(-0.0 - 0j, 0.0)[0],
                                          complex(-0.0, -0.0))
 
+    def test_py_ci_diff(self):
+        # Test _Py_ci_diff()
+        _py_ci_diff = _testcapi._py_ci_diff
+
+        self.assertComplexesAreIdentical(_py_ci_diff(-0.0 - 0j, 0.0)[0],
+                                         complex(-0.0, -0.0))
+
     def test_py_rc_diff(self):
         # Test _Py_rc_diff()
         _py_rc_diff = _testcapi._py_rc_diff
@@ -202,6 +216,13 @@ def test_py_rc_diff(self):
         self.assertComplexesAreIdentical(_py_rc_diff(-0.0, 0j)[0],
                                          complex(-0.0, -0.0))
 
+    def test_py_ic_diff(self):
+        # Test _Py_ic_diff()
+        _py_ic_diff = _testcapi._py_ic_diff
+
+        self.assertComplexesAreIdentical(_py_ic_diff(-0.0, 0j)[0],
+                                         complex(-0.0, -0.0))
+
     def test_py_c_neg(self):
         # Test _Py_c_neg()
         _py_c_neg = _testcapi._py_c_neg
@@ -221,6 +242,13 @@ def test_py_cr_prod(self):
         self.assertComplexesAreIdentical(_py_cr_prod(complex('inf+1j'), INF)[0],
                                                      complex('inf+infj'))
 
+    def test_py_ci_prod(self):
+        # Test _Py_ci_prod()
+        _py_ci_prod = _testcapi._py_ci_prod
+
+        self.assertComplexesAreIdentical(_py_ci_prod(complex('inf+1j'), INF)[0],
+                                                     complex('-inf+infj'))
+
     def test_py_c_quot(self):
         # Test _Py_c_quot()
         _py_c_quot = _testcapi._py_c_quot
@@ -250,6 +278,13 @@ def test_py_cr_quot(self):
         self.assertComplexesAreIdentical(_py_cr_quot(complex('inf+1j'), 2**1000)[0],
                                          INF + 2**-1000*1j)
 
+    def test_py_ci_quot(self):
+        # Test _Py_ci_quot()
+        _py_ci_quot = _testcapi._py_ci_quot
+
+        self.assertComplexesAreIdentical(_py_ci_quot(complex('1+infj'), 2**1000)[0],
+                                         INF - 2**-1000*1j)
+
     def test_py_rc_quot(self):
         # Test _Py_rc_quot()
         _py_rc_quot = _testcapi._py_rc_quot
@@ -257,6 +292,13 @@ def test_py_rc_quot(self):
         self.assertComplexesAreIdentical(_py_rc_quot(1.0, complex('nan-infj'))[0],
                                          0j)
 
+    def test_py_ic_quot(self):
+        # Test _Py_ic_quot()
+        _py_ic_quot = _testcapi._py_ic_quot
+
+        self.assertComplexesAreIdentical(_py_ic_quot(1.0, complex('inf-nanj'))[0],
+                                         -0.0)
+
     def test_py_c_pow(self):
         # Test _Py_c_pow()
         _py_c_pow = _testcapi._py_c_pow
diff --git a/Modules/_testcapi/complex.c b/Modules/_testcapi/complex.c
index b726cd3236f179..0fa1e065ea106a 100644
--- a/Modules/_testcapi/complex.c
+++ b/Modules/_testcapi/complex.c
@@ -57,48 +57,54 @@ _py_c_neg(PyObject *Py_UNUSED(module), PyObject *num)
         return Py_BuildValue("Di", &res, errno);                 \
     };
 
-#define _PY_CR_FUNC2(suffix)                                     \
-    static PyObject *                                            \
-    _py_cr_##suffix(PyObject *Py_UNUSED(module), PyObject *args) \
-    {                                                            \
-        Py_complex a, res;                                       \
-        double b;                                                \
-                                                                 \
-        if (!PyArg_ParseTuple(args, "Dd", &a, &b)) {             \
-            return NULL;                                         \
-        }                                                        \
-                                                                 \
-        errno = 0;                                               \
-        res = _Py_cr_##suffix(a, b);                             \
-        return Py_BuildValue("Di", &res, errno);                 \
+#define _PY_CX_FUNC2(suffix, prefix)                                     \
+    static PyObject *                                                    \
+    _py_##prefix##_##suffix(PyObject *Py_UNUSED(module), PyObject *args) \
+    {                                                                    \
+        Py_complex a, res;                                               \
+        double b;                                                        \
+                                                                         \
+        if (!PyArg_ParseTuple(args, "Dd", &a, &b)) {                     \
+            return NULL;                                                 \
+        }                                                                \
+                                                                         \
+        errno = 0;                                                       \
+        res = _Py_##prefix##_##suffix(a, b);                             \
+        return Py_BuildValue("Di", &res, errno);                         \
     };
 
-#define _PY_RC_FUNC2(suffix)                                     \
-    static PyObject *                                            \
-    _py_rc_##suffix(PyObject *Py_UNUSED(module), PyObject *args) \
-    {                                                            \
-        Py_complex b, res;                                       \
-        double a;                                                \
-                                                                 \
-        if (!PyArg_ParseTuple(args, "dD", &a, &b)) {             \
-            return NULL;                                         \
-        }                                                        \
-                                                                 \
-        errno = 0;                                               \
-        res = _Py_rc_##suffix(a, b);                             \
-        return Py_BuildValue("Di", &res, errno);                 \
+#define _PY_XC_FUNC2(suffix, prefix)                                     \
+    static PyObject *                                                    \
+    _py_##prefix##_##suffix(PyObject *Py_UNUSED(module), PyObject *args) \
+    {                                                                    \
+        Py_complex b, res;                                               \
+        double a;                                                        \
+                                                                         \
+        if (!PyArg_ParseTuple(args, "dD", &a, &b)) {                     \
+            return NULL;                                                 \
+        }                                                                \
+                                                                         \
+        errno = 0;                                                       \
+        res = _Py_##prefix##_##suffix(a, b);                             \
+        return Py_BuildValue("Di", &res, errno);                         \
     };
 
 _PY_C_FUNC2(sum)
-_PY_CR_FUNC2(sum)
+_PY_CX_FUNC2(sum, cr)
+_PY_CX_FUNC2(sum, ci)
 _PY_C_FUNC2(diff)
-_PY_CR_FUNC2(diff)
-_PY_RC_FUNC2(diff)
+_PY_CX_FUNC2(diff, cr)
+_PY_CX_FUNC2(diff, ci)
+_PY_XC_FUNC2(diff, rc)
+_PY_XC_FUNC2(diff, ic)
 _PY_C_FUNC2(prod)
-_PY_CR_FUNC2(prod)
+_PY_CX_FUNC2(prod, cr)
+_PY_CX_FUNC2(prod, ci)
 _PY_C_FUNC2(quot)
-_PY_CR_FUNC2(quot)
-_PY_RC_FUNC2(quot)
+_PY_CX_FUNC2(quot, cr)
+_PY_CX_FUNC2(quot, ci)
+_PY_XC_FUNC2(quot, rc)
+_PY_XC_FUNC2(quot, ic)
 _PY_C_FUNC2(pow)
 
 static PyObject*
@@ -125,15 +131,21 @@ static PyMethodDef test_methods[] = {
     {"complex_asccomplex", complex_asccomplex, METH_O},
     {"_py_c_sum", _py_c_sum, METH_VARARGS},
     {"_py_cr_sum", _py_cr_sum, METH_VARARGS},
+    {"_py_ci_sum", _py_ci_sum, METH_VARARGS},
     {"_py_c_diff", _py_c_diff, METH_VARARGS},
     {"_py_cr_diff", _py_cr_diff, METH_VARARGS},
+    {"_py_ci_diff", _py_ci_diff, METH_VARARGS},
     {"_py_rc_diff", _py_rc_diff, METH_VARARGS},
+    {"_py_ic_diff", _py_ic_diff, METH_VARARGS},
     {"_py_c_neg", _py_c_neg, METH_O},
     {"_py_c_prod", _py_c_prod, METH_VARARGS},
     {"_py_cr_prod", _py_cr_prod, METH_VARARGS},
+    {"_py_ci_prod", _py_ci_prod, METH_VARARGS},
     {"_py_c_quot", _py_c_quot, METH_VARARGS},
     {"_py_cr_quot", _py_cr_quot, METH_VARARGS},
+    {"_py_ci_quot", _py_ci_quot, METH_VARARGS},
     {"_py_rc_quot", _py_rc_quot, METH_VARARGS},
+    {"_py_ic_quot", _py_ic_quot, METH_VARARGS},
     {"_py_c_pow", _py_c_pow, METH_VARARGS},
     {"_py_c_abs", _py_c_abs, METH_O},
     {NULL},
diff --git a/Objects/complexobject.c b/Objects/complexobject.c
index 8ef0f6e65c7bfc..9a202d1b37aaed 100644
--- a/Objects/complexobject.c
+++ b/Objects/complexobject.c
@@ -50,6 +50,20 @@ _Py_rc_sum(double a, Py_complex b)
     return _Py_cr_sum(b, a);
 }
 
+Py_complex
+_Py_ci_sum(Py_complex a, double b)
+{
+    Py_complex r = a;
+    r.imag += b;
+    return r;
+}
+
+static inline Py_complex
+_Py_ic_sum(double a, Py_complex b)
+{
+    return _Py_ci_sum(b, a);
+}
+
 Py_complex
 _Py_c_diff(Py_complex a, Py_complex b)
 {
@@ -76,6 +90,23 @@ _Py_rc_diff(double a, Py_complex b)
     return r;
 }
 
+Py_complex
+_Py_ci_diff(Py_complex a, double b)
+{
+    Py_complex r = a;
+    r.imag -= b;
+    return r;
+}
+
+Py_complex
+_Py_ic_diff(double a, Py_complex b)
+{
+    Py_complex r;
+    r.real = -b.real;
+    r.imag = a - b.imag;
+    return r;
+}
+
 Py_complex
 _Py_c_neg(Py_complex a)
 {
@@ -109,6 +140,21 @@ _Py_rc_prod(double a, Py_complex b)
     return _Py_cr_prod(b, a);
 }
 
+Py_complex
+_Py_ci_prod(Py_complex a, double b)
+{
+    Py_complex r;
+    r.real = -a.imag*b;
+    r.imag = a.real*b;
+    return r;
+}
+
+static inline Py_complex
+_Py_ic_prod(double a, Py_complex b)
+{
+    return _Py_ci_prod(b, a);
+}
+
 /* Avoid bad optimization on Windows ARM64 until the compiler is fixed */
 #ifdef _M_ARM64
 #pragma optimize("", off)
@@ -249,6 +295,28 @@ _Py_rc_quot(double a, Py_complex b)
 
     return r;
 }
+
+Py_complex
+_Py_ci_quot(Py_complex a, double b)
+{
+    Py_complex r = a;
+    if (b) {
+        r.real = a.imag / b;
+        r.imag = -a.real / b;
+    }
+    else {
+        errno = EDOM;
+        r.real = r.imag = 0.0;
+    }
+    return r;
+}
+
+Py_complex
+_Py_ic_quot(double a, Py_complex b)
+{
+    Py_complex r = _Py_rc_quot(a, b);
+    return (Py_complex){-r.imag, r.real};
+}
 #ifdef _M_ARM64
 #pragma optimize("", on)
 #endif
@@ -1435,6 +1503,10 @@ PyImaginary_FromDouble(double imag)
     return (PyObject *) op;
 }
 
+#define CVAL(op) ((PyComplexObject *)(op))->cval
+#define CREAL(op) ((PyComplexObject *)(op))->cval.real
+#define CIMAG(op) ((PyComplexObject *)(op))->cval.imag
+
 static PyObject *
 imaginary_neg(PyComplexObject *v)
 {
@@ -1453,62 +1525,40 @@ imaginary_pos(PyComplexObject *v)
 /* Imaginary type arithmetic.  Binary operations below
    must support also complex and float (or int) operands. */
 
-static PyObject *
-imaginary_add(PyObject *v, PyObject *w)
-{
-    if (!PyImaginary_Check(v)) {
-        PyObject *tmp = v;
-
-        v = w;
-        w = tmp;
-    }
-
-    double a = ((PyComplexObject *)(v))->cval.imag;
-    double b;
-
-    if (PyComplex_Check(w)) {
-        a += ((PyComplexObject *)(w))->cval.imag;
-        b = ((PyComplexObject *)(w))->cval.real;
-    }
-    else if (real_to_double(&w, &b) < 0) {
-        return w;
-    }
-    return PyComplex_FromDoubles(b, a);
-}
-
-static PyObject *
-imaginary_sub(PyObject *v, PyObject *w)
-{
-    if (PyImaginary_Check(w)) {
-        double b = -((PyComplexObject *)(w))->cval.imag;
-        double a;
-
-        if (PyImaginary_Check(v)) {
-            b += ((PyComplexObject *)(v))->cval.imag;
-            return PyImaginary_FromDouble(b);
-        }
-        if (PyComplex_Check(v)) {
-            a = ((PyComplexObject *)(v))->cval.real;
-            b += ((PyComplexObject *)(v))->cval.imag;
-        }
-        else if (real_to_double(&v, &a) < 0) {
-            return v;
-        }
-        return PyComplex_FromDoubles(a, b);
-    }
-
-    double a = ((PyComplexObject *)(v))->cval.imag;
-    double b;
-
-    if (PyComplex_Check(w)) {
-        a -= ((PyComplexObject *)(w))->cval.imag;
-        b = ((PyComplexObject *)(w))->cval.real;
-    }
-    else if (real_to_double(&w, &b) < 0) {
-        return w;
-    }
-    return PyComplex_FromDoubles(-b, a);
-}
+#define IMAGINARY_ADDITIVE_BINOP(NAME, FUNC, OP)                       \
+    static PyObject *                                                  \
+    imaginary_##NAME(PyObject *v, PyObject *w)                         \
+    {                                                                  \
+        Py_complex a;                                                  \
+        if (PyImaginary_Check(w)) {                                    \
+            if (PyImaginary_Check(v)) {                                \
+                return PyImaginary_FromDouble(CIMAG(v) OP CIMAG(w));   \
+            }                                                          \
+            if (PyComplex_Check(v)) {                                  \
+                a = _Py_ci_##FUNC(CVAL(v), CIMAG(w));                  \
+            }                                                          \
+            else if (real_to_double(&v, &a.real) < 0) {                \
+                return v;                                              \
+            }                                                          \
+            else {                                                     \
+                a.imag = OP CIMAG(w);                                  \
+            }                                                          \
+        }                                                              \
+        else if (PyComplex_Check(w)) {                                 \
+            a = _Py_ic_##FUNC(CIMAG(v), CVAL(w));                      \
+        }                                                              \
+        else if (real_to_double(&w, &a.real) < 0) {                    \
+            return w;                                                  \
+        }                                                              \
+        else {                                                         \
+            a.real = OP a.real;                                        \
+            a.imag = CIMAG(v);                                         \
+        }                                                              \
+        return PyComplex_FromCComplex(a);                              \
+    }
+
+IMAGINARY_ADDITIVE_BINOP(add, sum, +)
+IMAGINARY_ADDITIVE_BINOP(sub, diff, -)
 
 static PyObject *
 imaginary_mul(PyObject *v, PyObject *w)
@@ -1520,41 +1570,42 @@ imaginary_mul(PyObject *v, PyObject *w)
         w = tmp;
     }
 
-    double a = ((PyComplexObject *)(v))->cval.imag;
     double b;
 
     if (PyComplex_Check(w)) {
-        b = -a * ((PyComplexObject *)(w))->cval.imag;
         if (PyImaginary_Check(w)) {
-            return PyFloat_FromDouble(b);
+            return PyFloat_FromDouble(-CIMAG(v)*CIMAG(w));
         }
-        a *= ((PyComplexObject *)(w))->cval.real;
-        return PyComplex_FromDoubles(b, a);
+
+        Py_complex a = _Py_ic_prod(CIMAG(v), CVAL(w));
+
+        return PyComplex_FromCComplex(a);
     }
     else if (real_to_double(&w, &b) < 0) {
         return w;
     }
-    return PyImaginary_FromDouble(a*b);
+    return PyImaginary_FromDouble(CIMAG(v)*b);
 }
 
 static PyObject *
 imaginary_div(PyObject *v, PyObject *w)
 {
-    if (PyImaginary_Check(w)) {
-        double b = ((PyComplexObject *)(w))->cval.imag;
-        double a;
+    double b;
 
+    if (PyImaginary_Check(w)) {
+        b = CIMAG(w);
         if (b) {
             if (PyImaginary_Check(v)) {
-                double a = ((PyComplexObject *)(v))->cval.imag;
-
-                return PyFloat_FromDouble(a/b);
+                return PyFloat_FromDouble(CIMAG(v)/b);
             }
             if (PyComplex_Check(v)) {
-                Py_complex a = ((PyComplexObject *)(v))->cval;
+                Py_complex a = _Py_ci_quot(CVAL(v), b);
 
-                return PyComplex_FromDoubles(a.imag/b, -a.real/b);
+                return PyComplex_FromCComplex(a);
             }
+
+            double a;
+
             if (real_to_double(&v, &a) < 0) {
                 return v;
             }
@@ -1562,24 +1613,18 @@ imaginary_div(PyObject *v, PyObject *w)
         }
     }
     else {
-        double a = ((PyComplexObject *)(v))->cval.imag;
-        double b;
-
         if (PyComplex_Check(w)) {
-            Py_complex b = ((PyComplexObject *)(w))->cval;
+            Py_complex a = _Py_ic_quot(CIMAG(v), CVAL(w));
 
-            errno = 0;
-            b = _Py_rc_quot(a, b);
-            b = (Py_complex){-b.imag, b.real};
             if (!errno) {
-                return PyComplex_FromDoubles(b.real, b.imag);
+                return PyComplex_FromCComplex(a);
             }
         }
         else if (real_to_double(&w, &b) < 0) {
             return w;
         }
         else if (b) {
-            return PyImaginary_FromDouble(a/b);
+            return PyImaginary_FromDouble(CIMAG(v)/b);
         }
     }
     PyErr_SetString(PyExc_ZeroDivisionError, "complex division by zero");