[issue25453] Arithmetics with complex infinities is inconsistent with C/C++
Francesco Biscani added the comment: @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. -- ___ Python tracker <rep...@bugs.python.org> <http://bugs.python.org/issue25453> ___ ___ Python-bugs-list mailing list Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue25453] Arithmetics with complex infinities is inconsistent with C/C++
Francesco Biscani added the comment: 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 :) -- ___ Python tracker <rep...@bugs.python.org> <http://bugs.python.org/issue25453> ___ ___ Python-bugs-list mailing list Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue25453] Arithmetics with complex infinities is inconsistent with C/C++
Francesco Biscani added the comment:
Hi Mark,
the original code is C++, and the inf + nanj result can be reproduced by the
following snippet:
"""
#include
#include
int main(int argc, char * argv[]) {
std::cout << std::complex(3,0) / 0. << '\n';
return 0;
}
"""
Here is the C99 version:
"""
#include
#include
#include
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.
--
___
Python tracker <rep...@bugs.python.org>
<http://bugs.python.org/issue25453>
___
___
Python-bugs-list mailing list
Unsubscribe:
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue25453] Arithmetics with complex infinities is inconsistent with C/C++
New submission from Francesco Biscani:
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)
------
messages: 253283
nosy: Francesco Biscani
priority: normal
severity: normal
status: open
title: Arithmetics with complex infinities is inconsistent with C/C++
type: behavior
versions: Python 2.7, Python 3.3
___
Python tracker <rep...@bugs.python.org>
<http://bugs.python.org/issue25453>
___
___
Python-bugs-list mailing list
Unsubscribe:
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
