Unexpected NANs in complex arithmetic

[Steven D'Aprano]

On Wednesday 22 June 2016 12:57, Chris Angelico wrote:

> On Wed, Jun 22, 2016 at 12:48 PM, Steven D'Aprano <steve at pearwood.info>
> wrote:
>> I'm doing some arithmetic on complex numbers involving INFs, and getting
>> unexpected NANs.
>>
>> py> INF = float('inf')
>> py> z = INF + 3j
>> py> z
>> (inf+3j)
>> py> -z
>> (-inf-3j)
>>
>> So far, nothing unexpected has occurred. But:
>>
>> py> -1*z  # should be the same as -z
>> (-inf+nanj)
>>
>>
>> And even more strange:
>>
>> py> 1*z
>> (inf+nanj)
>>
>>
>>
>> Is this the right behaviour? If so, what's the justification for it?
> 
> I've no idea, so I Googled StackExchange [1] and found this:
> 
> http://math.stackexchange.com/questions/585766/what-is-infinity-in-complex-
plane-and-what-are-operation-with-infinity-extended

Not actually very helpful, since complex is not modelled on the extended 
complex numbers. The complex data type has a large but finite number of 
infinities: ±INF+0.1j, ±INF+0.2j, ±INF+0.3j, etc, and a large but finite number 
of NANs. The extended complex numbers has a single infinity, and no NANs.


> Notably this:
> 
> """
> To some extent, +∞+∞ and −∞−∞ also play this role on the real axis:
> they are not a destination, they are road signs that tell us to go in
> a certain direction and never stop.
> """
> 
> So when your real part is float("inf"), what you're really saying is
> "Go as far as you possibly can in the positive direction, then keep
> going (because you haven't run out of numbers yet), and tell me, what
> is 1*z tending towards?".

In Real ℝ and Complex ℂ arithmetic, 1*z must equal z for all z's, or else 
arithmetic is inconsistent and anything is possible. That also applies to the 
extended ℝ and ℂ, both of which have a single infinity. Since neither ℝ nor ℂ 
have anything like NANs, there's no problem there.

In extended ℂ, we cannot distinguish between different infinities: the entire 
complex plane is bent around in a sphere so that all the "infinities" are at a 
single point. That is *not* what the complex() data type models:

py> INF+2j == INF+3j  # two distinct infinities (there are many more)
False

IEEE-754 floats are not ℝ, and INF does double-duty as "a really big but finite 
number which has overflowed", and "an actual infinity, whatever that means". So 
long as you remember that INF can mean either of those, you can usually reason 
about the behaviour of IEEE-754 INF.

The complex() type should, I expect, be treated more or less as a two-tuple of 
z = (x, y) where x and y are both floats. So:

1*(x+yj) should equal (1*x) + (1*y)j, which is well-defined. If y is finite, 
1*y should remain finite even if x is INF.

I would understand if 1*(NAN+2j) returned NAN+NANj. I'd disagree, but 
understand. But it makes no sense to me to introduce a NAN into a calculation 
just because you multiply by 1, even if it includes an INF. And multiplying by 
-1 should be identical to negating.


> Infinity isn't a number, and the imaginary
> part of your complex number isn't really a factor in figuring out
> where you're likely to end up with your multiplication.

If you consider ∞+2j, that's equivalent to:

- starting at the origin, move up two spaces, then to the right forever;

which is clearly not the same as ∞-2j:

- starting at the origin, move down two spaces, then to the right forever.

That's because the complex() plane tries to model ℂ with-INFs-and-NANs (note 
plurals) rather than the extended Complex plane, a.k.a. the Riemann sphere.

If complex() were modelled after the extended complex plane, then all complex 
numbers with ±INF as either the real or imaginary part should be considered 
equal:

INF+2j == INF-2j == INF+INFj == -INF-INFj


Putting a NAN in their certainly doesn't accomplish that.


> I guess that's
> a justification for it coming out as NaN.
> 
> ChrisA
> 
> [1]
> [http://www.theallium.com/engineering/computer-programming-to-be-officially-
renamed-googling-stackoverflow/