Arithmetics with complex infinities is inconsistent with C/C++

[326e1fe0-4c40-4e08-a9b7-069ffdf24d01]

BPO	25453
Nosy	@malemburg, @tim-one, @mdickinson, @ericvsmith, @ezio-melotti, @stevendaprano, @asmeurer, @berkerpeksag

^{Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.}

Show more details

GitHub fields: ```python assignee = None closed_at = None created_at = labels = ['interpreter-core', 'type-feature'] title = 'Arithmetics with complex infinities is inconsistent with C/C++' updated_at = user = 'https://bugs.python.org/FrancescoBiscani' ``` bugs.python.org fields: ```python activity = actor = 'Aaron.Meurer' assignee = 'none' closed = False closed_date = None closer = None components = ['Interpreter Core'] creation = creator = 'Francesco Biscani' dependencies = [] files = [] hgrepos = [] issue_num = 25453 keywords = [] message_count = 15.0 messages = ['253283', '253284', '253286', '253289', '253290', '253293', '253295', '253296', '253297', '253301', '253304', '253305', '253306', '253308', '253337'] nosy_count = 11.0 nosy_names = ['lemburg', 'tim.peters', 'mark.dickinson', 'eric.smith', 'stutzbach', 'ezio.melotti', 'steven.daprano', 'Aaron.Meurer', 'berker.peksag', 'Francesco Biscani', 'Saksham Agrawal'] pr_nums = [] priority = 'normal' resolution = None stage = None status = 'open' superseder = None type = 'enhancement' url = 'https://bugs.python.org/issue25453' versions = ['Python 3.6'] ```

Linked PRs

• gh-124829

[326e1fe0-4c40-4e08-a9b7-069ffdf24d01]

The C++11/C99 standards define a complex infinity as a complex number in which at least one of the components is inf. Consider the Python snippet:

>>> complex(float('inf'),float('nan'))*2
(nan+nanj)

This happens because complex multiplication in Python is implemented in the most straightforward way, but the presence of a nan component "infects" both components of the result and leads to a complex nan result. See also how complex multiplication is implemented in Annex G.5.1.6 of the C99 standard.

It feels wrong that a complex infinity multiplied by a real number results in a complex nan. By comparison, the result given here by C/C++ is inf+nan*j.

Note also this:

>>> complex(float('inf'),float('nan'))+2          
(inf+nanj)

That is, addition has a different behaviour because it does not require mixing up the components of the operands.

This behaviour has unexpected consequences when one interacts with math libraries implemented in C/C++ and accessed via Python through some wrapping tool. For instance, whereas 1./(inf+nan*j) is zero in C/C++, it becomes in Python

>>> 1./complex(float('inf'),float('nan'))  
(nan+nanj)

[78e49f92-dceb-4da7-aa03-b0619e5a1fac]

Hi

Would like to work on this bug...but I am new, would like some guidance

Please help me

[ezio-melotti]

Saksham, if you would like to work on this issue you can check the devguide for more information.

[stevendaprano]

I'm not entirely satisfied that the way it is calculated by C++11/C99 is correct. (Although I can see the appeal of the C version.) Mathematically, complex multiplication (a+bj)x should be identical to (a+bj)(x+0j), but obviously in the presence of NANs that is no longer the case. So it isn't clear to me that Python is wrong to allow NANs to "infect" the real part after multiplication.

Before changing the behaviour, I'd like to hear from someone who might be able to comment on what the IEEE-754 standard may have to say about this.

[326e1fe0-4c40-4e08-a9b7-069ffdf24d01]

The best reference I could find so far is in the C99 standard:

http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1256.pdf

Annex G is titled "IEC 60559-compatible complex arithmetic". In G.3.1 it is written:

""" A complex or imaginary value with at least one infinite part is regarded as an infinity (even if its other part is a NaN). """

Later on, in G.5.1.4, it is stated:

""" The * and / operators satisfy the following infinity properties for all real, imaginary, and complex operands:

• if one operand is an infinity and the other operand is a nonzero finite number or an infinity, then the result of the * operator is an infinity;
• if the first operand is an infinity and the second operand is a finite number, then the result of the / operator is an infinity;
• if the first operand is a finite number and the second operand is an infinity, then the result of the / operator is a zero; """

So to recap, according to these definitions:

• inf + nanj is a complex infinity,
• (inf + nanj) * (2 + 0j) should give a complex infinity (i.e., a complex number in which at least one component is inf).

I am not saying that these rules are necessarily "mathematically correct". I am aware that floating point is always a quagmire to get into, and the situation here is even more complex (eh!) than usual.

But it seems to me that Python, or CPython at least, follows a pattern of copying (or relying on) the behaviour of C for floating-point operations, and it seems like in the case of complex numbers this "rule" is broken. It certainly caught me by surprise anyway :)

[78e49f92-dceb-4da7-aa03-b0619e5a1fac]

I read the first 6 chapters of the devguide.

Now how should I proceed.

[ericvsmith]

Saksham: First we need our "experts" to decide what, if any, change is needed. If we decide that a change is needed, at that point we'd look at a patch.

[ericvsmith]

Saksham: Also, it would be best to take this discussion of how to produce a patch to the python-committers mailing list: https://mail.python.org/mailman/listinfo/python-committers

[berkerpeksag]

Also, it would be best to take this discussion of how to produce a patch to the python-committers mailing list:

or the core-mentorship list: https://mail.python.org/mailman/listinfo/core-mentorship :)

[ericvsmith]

Berker Peksag added the comment:

> Also, it would be best to take this discussion of how to produce a patch to the python-committers mailing list:

or the core-mentorship list: https://mail.python.org/mailman/listinfo/core-mentorship :)

Argh! That's of course what I meant. Apologies all around.

[mdickinson]

If the proposal is to add C99 Appendix G-style handling of nan+nanj results for complex multiplication, that doesn't seem unreasonable to me, though I'm not particularly convinced that the gain in complexity is worth it. So -0 from me. Whatever happens, I'd hope that the complex multiplication operation remains simple and easily explainable, not degenerating into a mass of special cases.

I'd be opposed to adding special-case handing for float-by-complex multiplication (i.e., doing anything other than promoting the float to a complex by introducing a zero imaginary part, and then doing normal complex multiplication). But I don't think that's what Francesco is proposing.

FWIW, NumPy semantics follow those of Python, though that's hardly an independent data point.

In answer to Steven, IEEE 754 has precisely nothing to say on this subject: it simply doesn't cover complex arithmetic.

Francesco: by the way, I'd be interested to see the C source that gave you inf + nan*j for this. How did you initialise the arguments to the multiplication? C99 has a significant issue here (fixed in C2011), namely that there's no standard way to create a complex value with given real and imaginary parts. (Using real_part + I*imaginary_part doesn't count, because in the absence of the Imaginary type, which few compilers seem to implement, the real part of the "I" constant mucks up the arithmetic w.r.t. infinities, NaNs and signed zeros.)

[mdickinson]

Changing Python versions and issue type: the current behaviour is certainly deliberate, so a proposal to change it would be a feature request, which could only land in Python 3.6 or later.

[326e1fe0-4c40-4e08-a9b7-069ffdf24d01]

Hi Mark,

the original code is C++, and the inf + nanj result can be reproduced by the following snippet:

"""
#include <complex>
#include <iostream>

int main(int argc, char * argv[]) {
  std::cout << std::complex<double>(3,0) / 0. << '\n';
  return 0;
}
"""

Here is the C99 version:

"""
#include <complex.h>
#include <math.h>
#include <stdio.h>

int main(int argc, char * argv[]) {
  double complex a = 3.0 + 0.0*I;
  printf("%f + i%f\n", creal(a/0.), cimag(a/0.));
  return 0;
}
"""

This is on Linux x86_64 with GCC 4.9. Clang gives the same result. Adding the "-ffast-math" compilation flag changes the result for the C version but apparently not for the C++ one.

The original code came from an implementation of a special function f(z) which has a pole near the origin. When computing f(0), the C++ code returns inf+nan*j, which is fine (in the sense that it is a complex infinity of some kind). I then need to use this result in a larger piece of code which at one point needs to compute 1/f(0), with the expected result being 0 mathematically. This works if I implement the calculation all within C/C++, but if I switch to Python when computing 1/f(0) I get nan + nan*j as a result.

Of course I could do all sorts of stuff to improve this specific calculation and avoid the problems (possibly improving the numerical behaviour near the pole via a Taylor expansion, etc. etc. etc.).

Beware that implementing C99 behaviour will result in a noticeable overhead for complex operations. In compiled C/C++ code one usually has the option to switch on/off the -ffast-math flag to improve performance at the expense of standard-conformance, but I imagine this is not feasible in Python.

[mdickinson]

Thanks for the update.

I'm not too worried about performance: I suspect that any additional overhead would be lost in the overhead of the Python machinery. It would be worth profiling to check, of course, before any change went in. I'd expect that most performance critical code of this type would be using NumPy/SciPy anyway.

I'm still unsure about whether the code is worth changing: I do agree that the behaviour suggested by C99 Annex G is (in the abstract) an improvement on Python's current behaviour, and compatibility with C would be a plus. On the other side, introducing an incompatibility with other versions of Python and with NumPy isn't ideal. For me, the potential benefits don't really overcome the drawbacks here, but I'd like to hear other opinions.

[326e1fe0-4c40-4e08-a9b7-069ffdf24d01]

@Mark

Yes I understand that this is a thorny issue. I was kinda hoping NumPy would just forward complex arithmetic to the underlying C implementation, but apparently that's not the case (which OTOH makes sense as the use of C99 complex numbers is not that widespread).

FWIW, the quoted section in the C standard (Annex G) contains the pseudocode for standard-conforming implementations of complex multiplication and division, in case the decision is taken to change the behaviour.

[skirpichev]

This happens because complex multiplication in Python is implemented in the most straightforward way, but the presence of a nan component "infects" both components

That actually happens due to absence of mixed-mode rules for complex arithmetic in CPython. Docs says: "Python fully supports mixed arithmetic: when a binary arithmetic operator has operands of different numeric types, the operand with the “narrower” type is widened to that of the other, where integer is narrower than floating point, which is narrower than complex.":

    complex(float('inf'), float('nan'))*2 ==
    complex(float('inf'), float('nan'))*(2+0j) ==
    complex(float('inf')*2 - float('nan')*0,
            float('nan')*2 + float('inf')*0) ==
    nan+nanj

C standards (since C99 and before upcoming C2y) have special rules for complex arithmetic: when one operand is 1) a real number or 2) an imaginary number (no real part). Unfortunately, it seems no compiler (except for HP UX) have support for 2).

@mdickinson, I think we have 3 options.

a) keep things as is b) support mixed-mode arithmetic rules for 1) case c) support both 1) and 2) cases; see https://github.com/skirpichev/cpython/pull/1

b) will match behaviour of popular C/C++ compilers. This also will fix the particular case, mentioned by OP. Unfortunately, that doesn't fix neither eval(repr(x)) == x invariant for complex numbers, nor allows to use usual rectangular notation. And it's easy to find examples of analytic identities, which will be broken (e.g. asin(z) definition via log and sqrt near the branch cut) with such arithmetic rules.

The later, I think, is a major drawback CPython's implementation of complex arithmetic. I.e. we can't define cot(z) just as 1/tan(z) or acot(z) as atan(1/z).

[skirpichev]

PR with 1) variant is ready for review: https://github.com/python/cpython/pull/124829