Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #110280 > unrolled thread

Unexpected NANs in complex arithmetic

Back to article view | Back to comp.lang.python

Contents

  Unexpected NANs in complex arithmetic Steven D'Aprano <steve@pearwood.info> - 2016-06-22 12:48 +1000
    Re: Unexpected NANs in complex arithmetic Chris Angelico <rosuav@gmail.com> - 2016-06-22 12:57 +1000
      Re: Unexpected NANs in complex arithmetic Dan Sommers <dan@tombstonezero.net> - 2016-06-22 03:54 +0000
        Re: Unexpected NANs in complex arithmetic Steven D'Aprano <steve+comp.lang.python@pearwood.info> - 2016-06-22 16:17 +1000
          Re: Unexpected NANs in complex arithmetic Paul Rubin <no.email@nospam.invalid> - 2016-06-21 23:25 -0700
          Re: Unexpected NANs in complex arithmetic Chris Angelico <rosuav@gmail.com> - 2016-06-22 16:30 +1000
            Re: Unexpected NANs in complex arithmetic Dan Sommers <dan@tombstonezero.net> - 2016-06-22 22:54 +0000
      Re: Unexpected NANs in complex arithmetic Steven D'Aprano <steve+comp.lang.python@pearwood.info> - 2016-06-22 16:58 +1000
        Re: Unexpected NANs in complex arithmetic Steven D'Aprano <steve+comp.lang.python@pearwood.info> - 2016-06-22 17:05 +1000
    Re: Unexpected NANs in complex arithmetic Antoon Pardon <antoon.pardon@rece.vub.ac.be> - 2016-06-22 09:14 +0200
    Re: Unexpected NANs in complex arithmetic Oscar Benjamin <oscar.j.benjamin@gmail.com> - 2016-06-22 10:18 +0100

#110280 — Unexpected NANs in complex arithmetic

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?


-- 
Steven

[toc] | [next] | [standalone]

#110281

On Wed, Jun 22, 2016 at 12:48 PM, Steven D'Aprano <steve@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

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?". 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. 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/

[toc] | [prev] | [next] | [standalone]

#110284

On Wed, 22 Jun 2016 12:57:55 +1000, Chris Angelico wrote:

> On Wed, Jun 22, 2016 at 12:48 PM, Steven D'Aprano <steve@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)

[...]

>> 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
> 
> 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?". 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. I guess that's
> a justification for it coming out as NaN.

By the time Python returns a result for inf+3j, you're already in
trouble (or perhaps Python is already in trouble).  A complex number has
a real part and an imaginary part, and inf isn't real (i.e., it's not an
element of ℝ).

[toc] | [prev] | [next] | [standalone]

#110291

On Wednesday 22 June 2016 13:54, Dan Sommers wrote:

> By the time Python returns a result for inf+3j, you're already in
> trouble (or perhaps Python is already in trouble).

I don't see why. It is possible to do perfectly sensible arithmetic on INFs.


> A complex number has
> a real part and an imaginary part, and inf isn't real (i.e., it's not an
> element of ℝ).

That's okay, because floats are not the set ℝ.


-- 
Steve

[toc] | [prev] | [next] | [standalone]

#110294

Steven D'Aprano <steve+comp.lang.python@pearwood.info> writes:
>> By the time Python returns a result for inf+3j, you're already in
>> trouble (or perhaps Python is already in trouble).
> I don't see why. It is possible to do perfectly sensible arithmetic on INFs.

We sometimes think of the real line extended by +/- inf, or the complex
plane extended by a single point at infinity.  I'm not sure how to do
sensible arithmetic on a number like inf+3j, that doesn't treat it as
the same number as inf+4j.

[toc] | [prev] | [next] | [standalone]

#110297

On Wed, Jun 22, 2016 at 4:17 PM, Steven D'Aprano
<steve+comp.lang.python@pearwood.info> wrote:
> On Wednesday 22 June 2016 13:54, Dan Sommers wrote:
>
>> By the time Python returns a result for inf+3j, you're already in
>> trouble (or perhaps Python is already in trouble).
>
> I don't see why. It is possible to do perfectly sensible arithmetic on INFs.

Technically, arithmetic on INF is a short-hand for a limit. What is
inf+1? Well, it's the limit of x+1 as x tends toward +∞, which is ∞.
So inf+1 evaluates to inf.

>>> inf = float("inf")
>>> inf+1
inf

What's 5-inf? Same again - the limit of 5-x as x tends toward +∞.
That's tending in the opposite direction, so we get infinity the other
way:

>>> 5-inf
-inf

So it's not really "doing arithmetic on infinity", as such. That said,
I'm not entirely sure why "inf+3j" means we're already in trouble -
Dan, can you elaborate? ISTM you should be able to do limit-based
arithmetic, just the same.

ChrisA

[toc] | [prev] | [next] | [standalone]

#110345

On Wed, 22 Jun 2016 16:30:49 +1000, Chris Angelico wrote:

> On Wed, Jun 22, 2016 at 4:17 PM, Steven D'Aprano
> <steve+comp.lang.python@pearwood.info> wrote:
>> On Wednesday 22 June 2016 13:54, Dan Sommers wrote:
>>
>>> By the time Python returns a result for inf+3j, you're already in
>>> trouble (or perhaps Python is already in trouble).
>>
>> I don't see why. It is possible to do perfectly sensible arithmetic on INFs.
> 
> Technically, arithmetic on INF is a short-hand for a limit. What is
> inf+1? Well, it's the limit of x+1 as x tends toward +∞, which is ∞.
> So inf+1 evaluates to inf.
> 
>>>> inf = float("inf")
>>>> inf+1
> inf
> 
> What's 5-inf? Same again - the limit of 5-x as x tends toward +∞.
> That's tending in the opposite direction, so we get infinity the other
> way:
> 
>>>> 5-inf
> -inf
> 
> So it's not really "doing arithmetic on infinity", as such. That said,
> I'm not entirely sure why "inf+3j" means we're already in trouble -
> Dan, can you elaborate? ISTM you should be able to do limit-based
> arithmetic, just the same.

What does (or can) inf+3j mean, and how can that be meaningfully
different from inf+4j, or 3+infj?  My mental model for math's ℂ and
Python's complex() is the points on a plane plus a single infinity.¹

The quantity inf+3j is not a point on the plane, nor is it that single
infinity.  To support inf+3j and 4+infj, Python's complex() would have
to be the points on a plane plus four "lines of infinity" plus something
to handle the intersections of the lines.  Maybe it's me, but that model
creates more problems than it solves.

¹ Yes, I understand that math's ℂ and Python's complex() are different.

[toc] | [prev] | [next] | [standalone]

#110298

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

> On Wed, Jun 22, 2016 at 12:48 PM, Steven D'Aprano <steve@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/

[toc] | [prev] | [next] | [standalone]

#110299

On Wednesday 22 June 2016 16:58, Steven D'Aprano wrote:

> 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.

Ah, not one second after I hit send, it struck me -- Python is casting the non-
complex factor to complex before doing the multiplication!

So:

1*(INF+2j)

which I was expanding into (1*INF) + (1*2j) = (INF+2j), gets performed by 
Python as:

(1+0j)*(INF+2j)

= 1*INF + 1*2j + 0j*INF + 0j*2j
= INF + 2j + NANj - 0
= INF + NANj



-- 
Steve

[toc] | [prev] | [next] | [standalone]

#110301

Op 22-06-16 om 04:48 schreef Steven D'Aprano:
> 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)

What I remember from complex numbers is that a multiplication
with a number that has |z| = 1, is equivallent with a rotation.

So you should be able to get the polar representation of this
"number", add in the angle of -1, being π, and convert back to
the cartesian representation.

I think seen this way, the nan part makes perfect sense.


Also the multiplication of a+bj with c+dj is (ac-bd)+(ad+bc)j
With your "numbers" this gives.

  (inf*(-1) - 3*0) + (inf*0 + 3*(-1))j

Again the nan part makes perfect sense.

-- 
Antoon

[toc] | [prev] | [next] | [standalone]

#110308

On 22 June 2016 at 08:14, Antoon Pardon <antoon.pardon@rece.vub.ac.be> wrote:
> Op 22-06-16 om 04:48 schreef Steven D'Aprano:
>> 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)
>
> Also the multiplication of a+bj with c+dj is (ac-bd)+(ad+bc)j
> With your "numbers" this gives.
>
>   (inf*(-1) - 3*0) + (inf*0 + 3*(-1))j
>
> Again the nan part makes perfect sense.

This explains presumably how complex multiplication is implemented
although it's not really a mathematical justification.

The problem is that mathematically complex numbers distribute over
addition. Complex floats (i.e. Python's complex type) for the most
part distribute approximately over addition so that the expansion
Antoon refers to will usually be correct up to floating point error:

   (1)   (x1 + j*y1) * (x2 + j*y2) ~= (x1*x2 - y1*y2) + j*(x1*y2 + x2*y1)

We can think of this as applying the distribution rule twice twice
(once for each pair of brackets). If the first complex number is real
it should have zero real part i.e. y1 = 0. This leads to

   (2)   (x1 + j*y1) * (x2 + j*y2) = x1 * (x2 + j*y2) = x1*x2 + j*x1*y2

which is I think the kind of formula Steve is expecting when being
surprised above. However we also have using the full distributed
formula (1):

    (3)   (x1 + j*y1)*(x2 + j*y2) = x1*x2 - 0*y2 + j*(x1*y2 + x2*0)

If y2 and x2 are real then the RHSs of (3) and (2) are equivalent.
However in your example x2 is infinity.

In essence the problem is that the set of numbers represented by
Python's complex type is the set of complex numbers plus a few
infinities/nans. The special values break the distributivity of
addition so that the algebra implicit in the implementation of the
complex type doesn't hold.

To see that infinities/nans break distribution over addition consider:

    inf = inf*1  = inf*(3-2) = inf*3-inf*2 = inf-inf = nan

If you did want to change these cases to match some particular idea of
the infinities as limits then special casing for pure real or
imaginary complex numbers (e.g. using (2) for complex numbers with
zero imaginary part) would probably lead to the results Steve expects.
I'm not really sure what the gain is though since we can't really do
robust algebra with infinities anyway.

--
Oscar

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.python

csiph-web