This issue tracker has been migrated to GitHub, and is currently read-only.
For more information, see the GitHub FAQs in the Python's Developer Guide.

classification
Title: Arithmetics with complex infinities is inconsistent with C/C++
Type: enhancement Stage:
Components: Interpreter Core Versions: Python 3.6
process
Status: open Resolution:
Dependencies: Superseder:
Assigned To: Nosy List: Aaron.Meurer, Francesco Biscani, Saksham Agrawal, berker.peksag, eric.smith, ezio.melotti, lemburg, mark.dickinson, steven.daprano, stutzbach, tim.peters
Priority: normal Keywords:

Created on 2015-10-21 12:20 by Francesco Biscani, last changed 2022-04-11 14:58 by admin.

Messages (15)
msg253283 - (view) Author: Francesco Biscani (Francesco Biscani) Date: 2015-10-21 12:20
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)
msg253284 - (view) Author: Saksham Agrawal (Saksham Agrawal) Date: 2015-10-21 12:28
Hi

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

Please help me
msg253286 - (view) Author: Ezio Melotti (ezio.melotti) * (Python committer) Date: 2015-10-21 13:03
Saksham, if you would like to work on this issue you can check the devguide for more information.
msg253289 - (view) Author: Steven D'Aprano (steven.daprano) * (Python committer) Date: 2015-10-21 14:40
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.
msg253290 - (view) Author: Francesco Biscani (Francesco Biscani) Date: 2015-10-21 15:19
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 :)
msg253293 - (view) Author: Saksham Agrawal (Saksham Agrawal) Date: 2015-10-21 15:29
I read the first 6 chapters of the devguide.

Now how should I proceed.
msg253295 - (view) Author: Eric V. Smith (eric.smith) * (Python committer) Date: 2015-10-21 15:45
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.
msg253296 - (view) Author: Eric V. Smith (eric.smith) * (Python committer) Date: 2015-10-21 15:50
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
msg253297 - (view) Author: Berker Peksag (berker.peksag) * (Python committer) Date: 2015-10-21 17:06
> 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 :)
msg253301 - (view) Author: Eric V. Smith (eric.smith) * (Python committer) Date: 2015-10-21 17:23
> 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.
msg253304 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2015-10-21 18:04
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.)
msg253305 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2015-10-21 18:06
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.
msg253306 - (view) Author: Francesco Biscani (Francesco Biscani) Date: 2015-10-21 18:46
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.
msg253308 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2015-10-21 19:02
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.
msg253337 - (view) Author: Francesco Biscani (Francesco Biscani) Date: 2015-10-22 16:29
@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.
History
Date User Action Args
2022-04-11 14:58:23adminsetgithub: 69639
2016-05-16 21:29:52Aaron.Meurersetnosy: + Aaron.Meurer
2015-10-22 16:29:32Francesco Biscanisetmessages: + msg253337
2015-10-21 19:02:17mark.dickinsonsetmessages: + msg253308
2015-10-21 18:46:03Francesco Biscanisetmessages: + msg253306
2015-10-21 18:22:46mark.dickinsonsetnosy: + tim.peters
2015-10-21 18:06:46mark.dickinsonsettype: behavior -> enhancement
messages: + msg253305
versions: - Python 2.7, Python 3.4, Python 3.5
2015-10-21 18:04:53mark.dickinsonsetmessages: + msg253304
2015-10-21 17:23:14eric.smithsetmessages: + msg253301
2015-10-21 17:06:37berker.peksagsetnosy: + berker.peksag
messages: + msg253297
2015-10-21 15:50:03eric.smithsetmessages: + msg253296
2015-10-21 15:45:43eric.smithsetmessages: + msg253295
2015-10-21 15:29:52Saksham Agrawalsetmessages: + msg253293
2015-10-21 15:19:15Francesco Biscanisetmessages: + msg253290
2015-10-21 14:40:40steven.dapranosetnosy: + steven.daprano
messages: + msg253289
2015-10-21 13:03:50ezio.melottisetversions: + Python 3.4, Python 3.5, Python 3.6, - Python 3.3
nosy: + ezio.melotti, stutzbach, lemburg, eric.smith, mark.dickinson

messages: + msg253286

components: + Interpreter Core
2015-10-21 12:28:24Saksham Agrawalsetnosy: + Saksham Agrawal
messages: + msg253284
2015-10-21 12:20:43Francesco Biscanicreate