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

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

1) First, new code allows to use complex arithmetic for implementation of
   mathematical functions without special "corner cases".  Take the inverse
   hyperbolic sine as an example:

       >>> z = complex(-0.0, 2)
       >>> cmath.asinh(z)
       (-1.3169578969248166+1.5707963267948966j)
       >>> # naive textbook formula doesn't work:
       >>> cmath.log(z + cmath.sqrt(1 + z*z))
       (1.3169578969248166+1.5707963267948966j)
       >>> # "fixed" version does:
       >>> cmath.log(z + cmath.sqrt(complex(1 + (z*z).real, (z*z).imag)))
       (-1.3169578969248164+1.5707963267948966j)

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.

   Things are worse, because the CPython docs sometimes asserts that the
   rectangular form is used and that some simple invariants holds.  For
   example, sphinx docs for the complex class says: "complex(real=0,
   imag=0) ...  Return a complex number with the value real + imag*1j ...".
   But:

       >>> complex(0.0, cmath.inf)
       infj
       >>> 0.0 + cmath.inf*1j
       (nan+infj)

3) The ``eval(repr(x)) == x`` invariant was broken for the complex type.  Below
   are simple examples with signed zero:

       >>> complex(-0.0, 1.0)  # also note funny signed integer zero below
       (-0+1j)
       >>> -0+1j
       1j
       >> -(0.0-1j)  # "correct" input for above with Python numeric literals
       (-0+1j)
       >>> -(0-1j)  # also "correct"; integer 0 matters!
       (-0+1j)
       >>> complex(1.0, -0.0)
       (1-0j)
       >>> 1-0j
       (1+0j)
       >>> -(-1 + 0j)
       (1-0j)

   Similar is true for complex numbers with other special components:

       >>> complex(0.0, -cmath.inf)
       -infj
       >>> -cmath.infj
       (-0-infj)

Author: skirpichev

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"

Include/internal/pycore_floatobject.h

@@ -55,6 +55,8 @@ extern PyObject* _Py_string_to_number_with_underscores(
 
 extern double _Py_parse_inf_or_nan(const char *p, char **endptr);
 
+extern int convert_to_double(PyObject **v, double *dbl);
+
 
 #ifdef __cplusplus
 }

Include/marshal.h

@@ -13,7 +13,7 @@ PyAPI_FUNC(PyObject *) PyMarshal_ReadObjectFromString(const char *,
                                                       Py_ssize_t);
 PyAPI_FUNC(PyObject *) PyMarshal_WriteObjectToString(PyObject *, int);
 
-#define Py_MARSHAL_VERSION 4
+#define Py_MARSHAL_VERSION 5
 
 PyAPI_FUNC(long) PyMarshal_ReadLongFromFile(FILE *);
 PyAPI_FUNC(int) PyMarshal_ReadShortFromFile(FILE *);

Lib/ast.py

@@ -74,7 +74,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/copy.py

@@ -101,7 +101,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,
@@ -178,6 +178,7 @@ def _deepcopy_atomic(x, memo):
 d[float] = _deepcopy_atomic
 d[bool] = _deepcopy_atomic
 d[complex] = _deepcopy_atomic
+d[imaginary] = _deepcopy_atomic
 d[bytes] = _deepcopy_atomic
 d[str] = _deepcopy_atomic
 d[types.CodeType] = _deepcopy_atomic

Lib/test/test_complex.py

@@ -21,6 +21,23 @@
     (1, 0+0j),
 )
 
+class float2:
+    def __init__(self, value):
+        self.value = value
+    def __float__(self):
+        return self.value
+
+class MyIndex:
+    def __init__(self, value):
+        self.value = value
+    def __index__(self):
+        return self.value
+
+class MyInt:
+    def __int__(self):
+        return 42
+
+
 class ComplexTest(unittest.TestCase):
 
     def assertAlmostEqual(self, a, b):
@@ -73,6 +90,10 @@ def assertFloatsAreIdentical(self, x, y):
                 msg += ': zeros have different signs'
         self.fail(msg.format(x, y))
 
+    def assertComplexesAreIdentical(self, x, y):
+        self.assertFloatsAreIdentical(x.real, y.real)
+        self.assertFloatsAreIdentical(x.imag, y.imag)
+
     def assertClose(self, x, y, eps=1e-9):
         """Return true iff complexes x and y "are close"."""
         self.assertCloseAbs(x.real, y.real, eps)
@@ -118,6 +139,26 @@ def test_truediv(self):
             self.assertTrue(isnan(z.real))
             self.assertTrue(isnan(z.imag))
 
+        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):
@@ -188,24 +229,83 @@ def check(n, deltas, is_equal, imag = 0.0):
     def test_add(self):
         self.assertEqual(1j + int(+1), complex(+1, 1))
         self.assertEqual(1j + int(-1), complex(-1, 1))
+        self.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))
+        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)
         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):
@@ -338,6 +438,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):
         class NS:
@@ -481,31 +582,15 @@ def __complex__(self):
 
         self.assertRaises(EvilExc, complex, evilcomplex())
 
-        class float2:
-            def __init__(self, value):
-                self.value = value
-            def __float__(self):
-                return self.value
-
         self.assertAlmostEqual(complex(float2(42.)), 42)
         self.assertAlmostEqual(complex(real=float2(17.), imag=float2(23.)), 17+23j)
         self.assertRaises(TypeError, complex, float2(None))
 
-        class MyIndex:
-            def __init__(self, value):
-                self.value = value
-            def __index__(self):
-                return self.value
-
         self.assertAlmostEqual(complex(MyIndex(42)), 42.0+0.0j)
         self.assertAlmostEqual(complex(123, MyIndex(42)), 123.0+42.0j)
         self.assertRaises(OverflowError, complex, MyIndex(2**2000))
         self.assertRaises(OverflowError, complex, 123, MyIndex(2**2000))
 
-        class MyInt:
-            def __int__(self):
-                return 42
-
         self.assertRaises(TypeError, complex, MyInt())
         self.assertRaises(TypeError, complex, 123, MyInt())
 
@@ -532,6 +617,23 @@ def __complex__(self):
             self.assertEqual(complex(complex1(1j)), 2j)
         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(float2(42.)), 42j)
+        self.assertRaises(TypeError, imaginary, float2(None))
+
+        self.assertEqual(imaginary(MyIndex(42)), 42j)
+        self.assertRaises(OverflowError, imaginary, MyIndex(2**2000))
+
+        self.assertRaises(TypeError, imaginary, MyInt())
+        self.assertRaises(TypeError, imaginary, 123, MyInt())
+
+        self.assertRaises(TypeError, imaginary, complex())
+        self.assertRaises(TypeError, imaginary, object())
+
     def test___complex__(self):
         z = 3 + 4j
         self.assertEqual(z.__complex__(), z)
@@ -647,19 +749,29 @@ def test(v, expected, test_fn=self.assertEqual):
     def test_pos(self):
         class ComplexSubclass(complex):
             pass
+        class ImaginarySubclass(imaginary):
+            pass
 
         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))
 
@@ -675,15 +787,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

Lib/test/test_fractions.py

@@ -1610,7 +1610,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,

Lib/test/test_socket.py

@@ -931,7 +931,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:
@@ -943,7 +943,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:

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, '%c requires int or char', operator.mod, '%c', pi)
 
         class RaisingNumber:

Modules/cmathmodule.c

@@ -1231,15 +1231,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;
     }