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

[skirpichev]

"Generally, mixed-mode arithmetic combining real and complex variables should be performed directly, not by first coercing the real to complex, lest the sign of zero be rendered uninformative; the same goes for combinations of pure imaginary quantities with complex variables." (c) Kahan, W: Branch cuts for complex elementary functions.

This patch implements mixed-mode arithmetic rules, combining real and complex variables as specified by C standards since C99. Most C compilers implementing C99+ Annex G have only these special rules (without support for imaginary type, which is going to be deprecated in C2y).

New rules allow to use complex arithmetic for implementation of mathematical functions without "corner cases" for special numbers (like signed zero or infinities). Well, at least it works now for more cases;)

Examples:

>>> 2-0j  # was (2+0j)
(2-0j)
>>> z = complex(-0.0, 2)
>>> cmath.asinh(z)
(-1.3169578969248166+1.5707963267948966j)
>>> cmath.log(z + cmath.sqrt(1 + z*z))  # real component had wrong sign before
(-1.3169578969248164+1.5707963267948966j)
>>> (cmath.inf+1j)*2  # imaginary component was nan before
(inf+2j)

That doesn't work (requires support for imaginary type):

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

---

Notes for reviewers

• Maybe it worth add missing (as noted in the commit message) case for the true division (i.e. x/(u + vj)==(x*u + (-x*v)j)/(u**2 + v**2)).
• Maybe we can also fix repr(complex(-0.0, 1)), which currently prints funny negative integer zero; see this commit. Probably this doesn't make sense alone.
• JFR, full-fledged implementation of mixed-mode complex arithmetic (with new type): Imaginary type and IEC 60559-compatible complex arithmetic skirpichev/cpython#1
• For some advocacy of mixed-mode rules (and imaginary type in particular), see Rationale for C99 (Annex G), WG14 document N3215 and Augmenting a Programming Language with Complex Arithmetic
• On simple benchmarks (like ./python -m timeit -s 'c=1+1j;d=1.2' 'c*d') I got a measurable performance boost for mixed arithmetic (~10-12%, except for _Py_dc_quot() case) and a performance degradation for complex arithmetic (~4-5%).

• Issue: Arithmetics with complex infinities is inconsistent with C/C++ #69639

---

📚 Documentation preview 📚: https://cpython-previews--124829.org.readthedocs.build/

[serhiy-storchaka]

LGTM in general. Very LGTM! I only left a few suggestions for style.

[skirpichev]

@serhiy-storchaka, what do you think on fixing repr(complex(-0.0, 1)) or adding special case for the true division (this is what you did in #124243 or I did in skirpichev#1)?

I also worry that this is a half-way solution: this neither fixes all eval(repr) round-trip issues or allows to use any analytical identity from textbook, e.g.:

>>> -0.0+1j
1j
>>> z = complex(-0.0, 2)
>>> 1j*(cmath.log(1 - 1j*z) - cmath.log(1 + 1j*z))/2
(1.5707963267948966+0.5493061443340549j)
>>> cmath.atan(z)
(-1.5707963267948966+0.5493061443340549j)

Proper fix seems to be only something like skirpichev#1. Does it look complicated for you? (Note, that arithmetic methods on the complex type are mostly unchanged wrt the current pr.)

[skirpichev]

I also was thinking about using additional C-API helper functions like _Py_c_sum_real() (as GSL does), instead of inlining that stuff. Does it make sense for you?

[serhiy-storchaka]

I think that the idea of special handling of mixed complex-real arithmetic is much easier to "sell" than the idea of the imaginary number class. And it solves a half of problems. It will help to convince in necessarily to support pure imaginaries. Let's eat the elephant piece by piece.

adding special case for the true division

I was surprised you did not include it. What does C99+ say about it? It looks natural and is useful in some cases, so I would implement it even if the C standard omits it.

I also was thinking about using additional C-API helper functions like _Py_c_sum_real() (as GSL does), instead of inlining that stuff. Does it make sense for you?

It makes sense to me. The status of functions like _Py_c_sum() is unclear -- they are private, but documented. Some of them are used outside of complexobject.c, and some of them are non-trivial, -- this is likely the reason of exposing them. I think that this all would also be true for the new functions.

Of course, the changes in arithmetic and the new C API (even if it is private) should be documented in many places.

[skirpichev]

And it solves a half of problems.

And that is a problem.

I was surprised you did not include it. What does C99+ say about it?

It seems, there is no such special version in the C standard. Not sure why. And such case miss in implementations, e.g. clang:
https://github.com/llvm/llvm-project/blob/496187b3b81ea76a8b67d796609d7f09992cf96d/libcxx/include/complex#L835-L841

On another hand, GSL or MPC libraries implement such case.

I would implement it even if the C standard omits it.

Ok, I'll add this with naming scheme like in your pr.

Of course, the changes in arithmetic and the new C API (even if it is private) should be documented in many places.

New arithmetic rules documented in stdtypes.rst, C API stuff will go to Doc/c-api/complex.rst. Did I miss something else?

[serhiy-storchaka]

New arithmetic rules documented in stdtypes.rst, C API stuff will go to Doc/c-api/complex.rst. Did I miss something else?

In Doc/reference/expressions.rst: "If either argument is a complex number, the other is converted to complex". There may be other leftovers.

You missed a What's New entry and versionchanged directives.

[picnixz]

One last missing modification and LGTM! Thanks for your patience Sergey!

[skirpichev]

@serhiy-storchaka, are you ok with current wording in docs?

[serhiy-storchaka]

LGTM. 👍

[skirpichev]

Thanks to all reviewers, especially Serhiy!

What's next, @serhiy-storchaka? This pr, as per description, address #69639 only partially. It's still easy to smash sign of zero or make some component to be nan, e.g.:

>>> -0.0+1j
1j
>>> float('inf')*1j
(nan+infj)
>>> z = complex(-0.0, 2)
>>> import cmath
>>> cmath.atan(z)
(-1.5707963267948966+0.5493061443340549j)
>>> 1j*(cmath.log(1 - 1j*z) - cmath.log(1 + 1j*z))/2
(1.5707963267948966+0.5493061443340549j)

and we still have funny integer zeros in repr:

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

Probably, it doesn't make sense to change repr alone. And I think that the previous discussion showed - this set of problems can be solved only with new imaginary type.

Should I prepare a PEP draft, would you like to sponsor (or be PEP co-author) such idea? Or it's better to introduce this first on the d.p.o/ideas?

PS:
skirpichev#1 rebased and adjusted. (BTW it can be made smaller: some tests in the main just use imaginary literals like 1j e.g. to trigger errors, where 1+1j will work just as fine.)

[vstinner]

This change adds 6 private functions and documents them:

• _Py_cr_diff()
• _Py_cr_prod()
• _Py_cr_quot()
• _Py_cr_sum()
• _Py_rc_diff()
• _Py_rc_quot()

I would prefer to move these functions to the internal C API, or to make them public. I dislike adding private functions (functions with a name prefixed by _Py).

[skirpichev]

@vstinner, these functions follows existing convention for comlexobject.c.

move these functions to the internal C API

I don't think it's a good idea. These functions are extensively used in the cmath module and outside of the CPython it happens too.

make them public

This does make sense for me. We could just drop _ prefixes or do you suggest some more severe renaming? Should I open issue in the C-API repo or this seems to be a minor issue for you?

[vstinner]

Previous discussions about _Py_c_xxx() complex functions:

• _Py_c_sum, _Py_c_diff, etc - were removed from Python.h #112019
• C API: Remove private C API functions (move them to the internal C API) #106320 (comment)

[vstinner]

We could just drop _ prefixes or do you suggest some more severe renaming?

I would prefer to rename _Py_c_sum() to PyComplex_Add() and pass PyComplex struct by reference (pointer) instead of passing it by value. So it would be a new C API.

Note: _Py_c_abs() is not documented.

Should I open issue in the C-API repo or this seems to be a minor issue for you?

If you agree with the idea of adding a public C API, you can open a new issue in this project (cpython) first.

[skirpichev]

FYI: #128813