1

This is a follow up of the old discussion. It’s outcome was, in short, that problems with complex arithmetic in Python can be solved only with a new imaginary type.

I’ll remind those problems with the latest implementation (CPython 3.13) of complex arithmetic on few examples:

>>> 2-0j  # (1)
(2+0j)
>>> z = complex(-0.0, 2)
>>> cmath.asinh(z)
(-1.3169578969248166+1.5707963267948966j)
>>> cmath.log(z + cmath.sqrt(1 + z*z))  # (2)
(1.3169578969248166+1.5707963267948966j)
>>> (cmath.inf+1j)*2  # (3)
(inf+nanj)
>>> -0.0+1j  # (4)
1j
>>> cmath.atan(z)
(-1.5707963267948966+0.5493061443340549j)
>>> 1j*(cmath.log(1 - 1j*z) - cmath.log(1 + 1j*z))/2  # (5)
(1.5707963267948966+0.5493061443340549j)
>>> repr(complex(-0.0, 1))  # (6)
'(-0+1j)'

Examples (1) and (4) show, that complex numbers in the CPython doesn’t follow to the rectangular notation (e.g. Complex number), i.e. a+bj != complex(a, b) in general. Examples (2), (3) and (5) show that we can’t use complex arithmetic say for implementation of mathematical functions without a special treatment for complex numbers with special components (zero, signed zero or infinities). Even if the function is defined by simple formula like 1j/z or 2*z as in example (3). (You can take look on the mpmath sources to see that it’s a recurring problem for real-world applications.) Finally, example (6) has funny negative integer zero, that also breaks eval(repr) round-trip for complex numbers.

There is a known solution for all above issues from numerical experts: “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.

Recently my pr added support for 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. So, in the CPython 3.14 examples (1)-(3) work correctly.

Lets discuss if we can fix the rest. That will add support for a new builtin type, say imaginary, and imaginary literals will create instances of this type (no real part!). The patch itself is not too big, here is a draft to play with: Imaginary type and IEC 60559-compatible complex arithmetic by skirpichev · Pull Request #1 · skirpichev/cpython · GitHub.

For some further 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. Both the GNU GSL and the GNU MPC have special routines to do arithmetic with operands of different types, including pure-imaginary.

Edit: alternative approach is using special constructor forms. Say, complex(real=x, pure=True) (a copy of real line) and complex(imag=y, pure=True) (for pure-imaginary numbers). Imaginary literals will be translated to the second form. Also, doing complex arithmetic, floats will be automatically coerced to the first form. The rest (i.e. arithmetic rules) is same as in the proposed version. But I suspect that this implementation will be more complicated internally.

3 Likes
2

I’m aware that I’m probably the person who’s brought this idea up the most, historically. From a purity point of view it seems as though it has a lot of advantages, but I’m afraid I remain skeptical that this would work well in practice in current Python.

My impression is that the imaginary type initiative introduced in C99 (via the _Imaginary keyword) is essentially a failure, in that very few compilers ended up adopting and supporting it. It might be good to understand why adoption was so poor, in case the reasons cross over to Python land.

My main concern would be the knock-on effect on 3rd party packages (NumPy, SymPy, many more).

3

I suspect we all know the obvious answer, but I think it’s rather hard to verify that it’s the right answer. For example, not all major compilers (e.g. msvc) even support Annex G. Does it mean there is something wrong with the standard?

Another answer might be, that Annex G still has small user base. For example, the GNU GSL supports (incorrectly) native complex numbers since 2021, while it’s API has all primitives for mixed-mode arithmetic.

Also, it seems most (all?) current compilers, supporting Annex G, implement mixed-mode arithmetic for real operands. Should we interpret this, that only that part of the approach, advocated by Kahan, is important? For me, it seems that removal message from the WG 14 just ignores all mathematical arguments for this feature (e.g. WG’s documents, quoted in my post).

BTW, in principle, I think we can avoid new type at cost of adding some special forms for the complex constructor (mostly invisible for users). Say, complex(real=x, pure=True) (a copy of real line) and complex(imag=y, pure=True) (for pure-imaginary numbers). Imaginary literals will be translated to the second form. Also, doing complex arithmetic, floats will be automatically coerced to the first form. The rest (i.e. arithmetic rules) is same as in the proposed version. But I suspect that this implementation will be more complicated internally.

Maybe positive? What’s bad in replacing incorrect or meaningless answers like (cmath.inf+1j)*2 currently?

For the mpmath it eventually allows to drop some boilerplate code and write instead:

@defun_wrapped
def acot(ctx, z):
    return ctx.atan(ctx.one / z)

@defun_wrapped
def asec(ctx, z):
    return ctx.acos(ctx.one / z)

@defun_wrapped
def acsc(ctx, z):
    return ctx.asin(ctx.one / z)

(And silently fix all functions, that now miss required checks for special complex numbers :-))

For the gmpy2, mixed-mode rules could be implemented as for CPython. As it was noted before, the MPC library has all primitives to support this.

4

Do we? I’m not sure that I do, and it’s not clear to me which “obvious” answer you’re referring to. If I had to guess, my guess would be that imaginary types are not supported in most C compilers because there just isn’t much demand for them. That’ll be at least partly because we’re talking about corner-case details of what’s already a fairly niche use-case (even in the domain of scientific computing, which is where I spend most of my time, complex numbers don’t show up that often). But all of this applies to Python-land, too. Anyway, I’m just guessing; some research would be valuable. (E.g., are there related gcc or clang feature requests, and what was the response to those requests?)

Nice; I hadn’t seen that document. Thanks for the link.

I think that will have its own set of problems. For the PEP, it might be helpful to pick one of the two options (new imaginary type versus altered complex type) and focus the discussion on that option, to avoid an overly-wide-ranging discussion.

Perhaps, in the long term (say 5-10 years from now and beyond). Then again, perhaps not. But I’m more concerned about the short term. How would you address backwards compatibility concerns? Changing the type of a built-in expression isn’t a minor change, and it would be astonishing if it didn’t break at least some third-party code. And at that point you have to ask whether the long-term benefits of the change really outweigh the cost of the code churn required across the Python ecosystem.

Suppose, for the sake of argument, that this change went into Python 3.14. Would NumPy and SymPy (for example) be forced to change their codebases if they wanted to remain compatible with Python 3.14, or would there be a backwards compatibility path that allowed them to update at their leisure?

As a data point, have you tried running the SymPy and NumPy test suites against your PR branch?

2 Likes
5

That’s exactly what I meant as “obvious answer”. However, mixed-mode arithmetic rules are supported, despite they have much less sense, having implemented in a half-way manner.

Perhaps, this depends on application. Quantum mechanics, for example?

Probably those issues are relevant (gcc issue mention practical problems with implementation, but I’m not certain about it’s severity):

Well, new type implemented as a subtype of complex. I would expect that most breakages will be related to type(z) is complex-like checks (not sure how common they are wrt isinstance checks, see below for tests). I doubt this can be workarounded in a backward-compatible way. Alternative approach, using new special constructor forms, shouldn’t suffer from this issue.

Another source of troubles - some changes in the repr:

  • repr(complex(-0.0, x)) → (-0.0+xj)
  • repr(complex(+0.0, x)) → (0.0+xj) (to distinguish from repr(imaginary(x)) → xj)

These changes in formatting could be introduced gradually.

I don’t expect too much troubles for SymPy (of course, this will be tested), rather they could come from the mpmath dependency (which uses complex type for the fp context, and imaginary literals of course). So I run it’s test suite and got 2 failures: both related to changes in the repr output.

1 Like
6

When I had to extensively fight with branch cuts for a project years ago, all issues came to be via (the lack of) signed zeros. Is that true more generally? If so, is special-casing zeros in the arithmetic routines not a viable option?

7

I am a great supporter of the imaginary type. I have only one concern – that it is a quite large change in the Python core (ast, marshal, C API, …) for a solution of rather marginal problems. But this is the only way, we can only delay the decision.

In hindsight, there was no need to add support of complex numbers in the Python core. It could be supported in extension, like decimals and fractions.

4 Likes
8

I don’t know what you were fighting extensively with but in my experience signed zeros cause more problems than they solve. If you want to handle branch cuts in general then you need more than two types of zero:

  • -0.0
  • +0.0
  • 0.0 (actual zero)

The suggestion to add an imaginary type here amounts to wanting to be able to have a complex number whose real part is the actual zero that somehow doesn’t exist. I don’t propose to change the use of signed zeros in Python since it is based on standards implemented in hardware and exists in other languages etc but in general I think we would all be better off without the signed zeros and with finding better ways to handle underflow.

I do think that it is good if the stdlib and libraries like numpy are consistent though e.g.:

In [11]: np.float64(-0.0) + 1j
Out[11]: np.complex128(1j)

In [12]: float(-0.0) + 1j
Out[12]: 1j

In [13]: np.complex128(-0.0, 1)
Out[13]: np.complex128(-0+1j)

In [14]: complex(-0.0, 1)
Out[14]: (-0+1j)

I’m not sure if this behaviour in numpy is just designed to match the behaviour of Python core but it would be good if whatever is considered to be the right behaviour here can be consistent.

Other libraries may not have signed zeros e.g. gmpy2 does but mpmath, sympy and python-flint (arb) do not because they handle overflow/underflow differently e.g.:

In [20]: 1 / flint.arb(10**10).exp()
Out[20]: [9.27858442032487e-4342944820 +/- 6.08e-4342944835]

In [21]: 1 / flint.arb(10**100).exp()
Out[21]: nan

A nan is often more useful than a result based on signed zeros. Note that arb can compute this if given more precision but otherwise it won’t just underflow to any signed or unsigned zero and keeps a clear distinction between a small possibly nonzero number vs an exact zero:

In [5]: z = flint.arb(-1e100)

In [6]: z.exp()
Out[6]: [+/- 1.69e-102435199438739363750012109250103232700]

In [7]: z - z
Out[7]: 0

It isn’t possible to make things consistent across types that do or don’t have signed zeros but if the signed zeros are just ignored then it seems that the imaginary type doesn’t affect anything.

1 Like
9

Any special components are affected. Consider just one example from my post with a linear map z\to 2 z. In CPython 3.13 you will get nans, unless add a special case for infinite components.

I would say that changes outside of the complexobject.c seems small and more or less trivial.

Probably there is a variant with new special forms for the complex constructor (say, complex(imag=y, pure=True) to serve role of imaginary numbers). Maybe it’s better from compatibility point of view, but it looks for me more complicated internally.

That’s good. On another hand, that’s bad, because it’s not consistent with C. You can start with my (1) example, it’s equal to 2.0 + CMPLX(-0.0, -0.0).

It does, e.g.:

>>> complex('inf+1j')*1j  # should be `(-1+infj)`
(nan+infj)
10

The special case can be for zero components rather than infinite components:

>>> (cmath.inf+1j)*2  # (3)
(inf+nanj)

The obvious thing to do here is recognise that the imaginary part of the second operand is zero and handle that exactly:

>>> complex(z.real * 2, z.imag * 2)
(inf+2j)

That is precisely what arb does:

In [7]: flint.acb('inf', 1) * 2
Out[7]: [+/- inf] + 2.00000000000000j

In [8]: flint.acb('inf', 1) * (2 + 1e-100j)
Out[8]: [+/- inf] + [+/- inf]j

This is another example in which the obvious thing to do is to observe that there is a zero component so that the full complex multiplication formula can be bypassed:

>>> complex(-(1.0*z.imag), 1.0*z.real)
(-1+infj)

You cite what you consider the correct answers here that seem to be based on handling the zero components but for some reason don’t want to implement this by handling the zero components. Is that because it wouldn’t handle signed zeros in the expected way?

I don’t think many sensible calculations are done with things like inf+1j so I also don’t see a big problem with having nans in cases like this. As I said before though I would like things to be consistent in so far as possible. Consistency with C is useful but consistency within the Python ecosystem is more important.

1 Like
11

I would say - both.

Exactly! “The obvious thing” is to forget about complex arithmetic in such corner cases and work by hand, doing component-wise computations.

I’m not sure it does special cases for mixed-mode arithmetic with arb (like the mpmath does for mpf). At least, this is not a part of the public API.

Is this a positive infinity?

Edit: Oh, I see - that means the whole extended real line in the real component. Somewhat better that just nan, but apparently this was obtained without using mixed-mode arithmetic. Such is possible if your components are interval-like objects, like arb, but not with usual floating-point arithmetic.

Sure) “Not use complex arithmetic.” This is an issue I’m trying to address.

Infinities also for you “cause more problems than solve”? :slight_smile: Well, maybe…

But I feel this goes slightly off-topic. We all agreed that the IEC 60559 is something, that worth to follow for floats. There is a way to follow to this standard also in implementation of the complex arithmetic. Is this worth or not?

Unfortunately, in other cases you might have non-nans:

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

That’s much harder to catch than nans, isn’t? Do you have a plan?

The Python is often serve for prototyping algorithms, that will be implemented later in C. I don’t see any problem to achieve a consistency wrt NumPy as well.

BTW, I would appreciate opinion of NumPy’s people. So far it’s yet possible even revert https://github.com/python/cpython/pull/124829, if they feel it’s too much and/or against consistency.

12

There are functions like acb_mul_arb but generally the expectation is that you would work with the complex type and the functions will keep track of whether or not the real or imaginary parts are zero. There is no need for an imaginary type because it is just a complex number with a zero real part. The code in acb_mul is just:

def acb_mul(a, b):
    if is_zero(a.real):
        ...
    elif is_zero(a.imag):
        ...
    elif is_zero(b.real):
        ...
1 Like
13

Not only an impression. _Imaginary type is removed in C2y draft (by N3274).

Up to now we have no information about an up-to-date conforming implementation that would have imaginary types. Nonetheless, for most of it such a support would not make a hypothetical implementation behave badly for C2y. It would simply transition through imaginary types on rare occasions, but the user-visible results of complex type would be the same.

When writing C code, it’s not really needed.

  1. Imaginary part operations (with other imaginary part) can be done via normal floating point type. Since _Imaginary had poor suport, compilers can assume code would prefer this approach.
  2. Conversion to _Complex can be done via bit casting. From C spec, 6.2.5 Types:

    Each complex type has the same representation and alignment requirements as an array type containing exactly two elements of the corresponding real type; the first element is equal to the real part, and the second element to the imaginary part, of the complex number.

14

Same is true for real types, of course. Yet in this case compilers usually follow Annex G and implement mixed-mode operations for complex and real types. While it’s also possible to do component-wise, of course.

Why compilers (rather, it’s authors) “assume” that this part of the standard could be violated? That’s something Mark asked to explain rather than show us as a matter of fact.

15

I meant documented API (see reference above). Perhaps, it’s a bug;)

Unfortunately, it’s not an option for python complex type. We can’t implement mixed-mode arithmetic this way.

It can be done for type like mpmath’s mpf. BTW, here is just another example where it could be useful: polylog at negative floating point infinity · Issue #473 · mpmath/mpmath · GitHub Existing mpmath code is correct for infinite z too, it’s just assumes complex arithmetic “as in textbook”, where e.g. z * i = -\Im z + i \Re z. But with current complex arithmetic in the CPython - proper fix will include a special case…

16

It is documented. Your link was to the docs for acf (arbitrary complex float) rather than acb (arbitrary complex ball). The difference between these is like the difference between arf and arb in the real case.

17

Ok, I will elaborate.

1. Operations on imaginary part.

//example with imaginary
#include <complex.h>
#ifndef _Imaginary_I
#  error "No imaginary support"
#endif
double imaginary a, b;
double imaginary c = a + b;

//example without imaginary
double a, b;
double c = a + b;

Note, example without imaginary is not less readable, but it’s easier on the compiler. There is no point in requiring imaginary types, for calculations on them.

2. Converting to complex

#include <complex.h>
#ifndef CMPLX //possible before C11
#  ifdef _Imaginary_I
#    define CMPLX(x, y) ((double complex)( \
     (double)(x) + (double)(y) * _Imaginary_I ))
#  else
#    define CMPLX(x, y) (*(double complex*)(double[2]){x,y})
#  endif
#endif

Note, if defined like that, it can’t be used were constant expression is required (cause of 2nd variant). This includes static initializer, but why would anyone need negative zero in static?

3. Standard violation

No, there is no violation. Annex G is fully optional. Compiler not implementing Annex G can be fully standard complaint, if it doesn’t pre-define __STDC_IEC_559_COMPLEX__ (since C99, until C23), or __STDC_IEC_60559_COMPLEX__ (since C23).

Now, according to cppreference, before C11 some compilers did pre-define __STDC_IEC_559_COMPLEX__ without implementing imaginary types. So, “POSIX recommends checking if the macro _Imaginary_I is defined […]”. I don’t know why those compilers did this, but it won’t tell us why more compilers didn’t implement imaginary types within 24 years.

Btw, also from cppreference:

To date, only Oracle C compiler is known to have implemented imaginary types.

Why would compilers not support this optional feature?

First, I need to say, I have an assumption that it’s not a choice out of negligence. I think it’s deliberate choice, for simplicity and optimization. My reasoning is that, they are implementing other less documented extensions. For example, (because of current PEP 505 discussion) the oldest coalescing operator I know about: the binary ternary operator. (For pointers a and b,
expression a?:b will return a only if it’s not null.)

It needs to be understood, that C creation process is fundamentally different from (for example) Python’s. First lets look at standard creation:

  1. ISO WG14 drafts update to standard.
  2. Some compiler maintainers will see those and implement them, or comment why those changes are bad idea.
  3. ISO updates the standard.

One might think the language will only grow, but no it won’t. Note there is no reference implementation in this process. Compilers don’t need to mimic one all-including compiler. If some feature makes the compiler worse, it won’t be implemented.

Here is my reasoning:.

  1. Since support for imaginary types is poor, old large projects would support compilers without it. Therefore there is no production code that won’t compile without it.
  2. Since support for imaginary types is poor, but they can easily be replaced by floating point types, aka the most benchmarked types in existence. In order to be faster, good code would use floating point types instead. Even if not on this compiler, on some other compiler it would make a difference.
  3. Is it worth it, to allow maybe unoptimized or poorly maintained part of the compiler, to support new bad code?

I’d like to point out that’s my view, actual reason could be different. I also don’t know the details on ISO/IEC 60559.

18

Thanks.

Not sure about compiler, but maybe this is not less readable just because it’s addition? And you miss mixed-mode operations - that’s the purpose of the imaginary type, not just arithmetic within the imaginary type. Doing this properly, you will just reinvent the complex arithmetic.

As for real life, try to address my previous example (polylog at negative floating point infinity · Issue #473 · mpmath/mpmath · GitHub) with inline version of mixed-mode arithmetic. I would expect that the source code length for polylog_continuation() will be doubled. (Or you reinvent the proper complex arithmetic.)

That’s something I would like from complex arithmetic in python: you can take a textbook formula (like \text{asinh}(z) = \log\left(z + \sqrt{1 + z^2}\right), \text{atan}(z)=i \left(\log\left(1 - i z\right) - \log\left(1 + i z\right)\right), \text{acot}(z) = \text{atan}(1/z)) and just use it in the code. No special cases if z is infinite, just above/below of the branch cut, etc. That seems to be doable with the Annex G.

But they do this, see mentioned bugs. BTW, reading the gcc thread (referenced https://sourceware.org/bugzilla/show_bug.cgi?id=15720) - this means “it indicates intent rather than completeness” (c).

IIRC, it’s also HP C compiler.

It looks that you (again) just start from the thesis you are pretended to explain)

I would agree with you, that if in violation of the standard compiler authors starting to “implement” (and add required defines to assert that) the Annex G without the imaginary type - users will have to find workarounds. Thus, in a kind of self-reinforced feedback - compilers will have less pressure to add missing types. But that doesn’t explain the beginning.

2 Likes
19

@storchaka, I prepared a PEP draft (unfortunately, that took some more time): PEP 9999: Imaginary type and IEC 60559-compatible complex arithmetic by skirpichev · Pull Request #4681 · python/peps · GitHub

Would you sponsor such PEP?

3 Likes
20

Would you sponsor such PEP?

Absolutely.

1 Like