Why does (inf + 0j)*1 evaluate to inf + nanj?

Question

>>> (float('inf')+0j)*1
(inf+nanj)

Why? This caused a nasty bug in my code.

Why isn't 1 the multiplicative identity, giving (inf + 0j)?

Answer 1

The 1 is converted to a complex number first, 1 + 0j, which then leads to an inf * 0 multiplication, resulting in a nan.

(inf + 0j) * 1
(inf + 0j) * (1 + 0j)
inf * 1  + inf * 0j  + 0j * 1 + 0j * 0j
#          ^ this is where it comes from
inf  + nan j  + 0j - 0
inf  + nan j

Answer 2

Mechanistically, the accepted answer is, of course, correct, but I would argue that a deeper ansswer can be given.

First, it is useful to clarify the question as @PeterCordes does in a comment: "Is there a multiplicative identity for complex numbers that does work on inf + 0j?" or in other words is what OP sees a weakness in the computer implementation of complex multiplication or is there something conceptually unsound with inf+0j

Short answer:

Using polar coordinates we can view complex multiplication as a scaling and a rotation. Rotating an infinite "arm" even by 0 degrees as in the case of multiplying by one we cannot expect to place its tip with finite precision. So indeed, there is something fundamentally not right with inf+0j, namely, that as soon as we are at infinity a finite offset becomes meaningless.

Long answer:

Background: The "big thing" around which this question revolves is the matter of extending a system of numbers (think reals or complex numbers). One reason one might want to do that is to add some concept of infinity, or to "compactify" if one happens to be a mathematician. There are other reasons, too (https://en.wikipedia.org/wiki/Galois_theory, https://en.wikipedia.org/wiki/Non-standard_analysis), but we are not interested in those here.

One point compactification

The tricky bit about such an extension is, of course, that we want these new numbers to fit into the existing arithmetic. The simplest way is to add a single element at infinity (https://en.wikipedia.org/wiki/Alexandroff_extension) and make it equal anything but zero divided by zero. This works for the reals (https://en.wikipedia.org/wiki/Projectively_extended_real_line) and the complex numbers (https://en.wikipedia.org/wiki/Riemann_sphere).

Other extensions ...

While the one point compactification is simple and mathematically sound, "richer" extensions comprising multiple infinties have been sought. The IEEE 754 standard for real floating point numbers has +inf and -inf (https://en.wikipedia.org/wiki/Extended_real_number_line). Looks natural and straightforward but already forces us to jump through hoops and invent stuff like -0 https://en.wikipedia.org/wiki/Signed_zero

... of the complex plane

What about more-than-one-inf extensions of the complex plane?

In computers, complex numbers are typically implemented by sticking two fp reals together one for the real and one for the imaginary part. That is perfectly fine as long as everything is finite. As soon, however, as infinities are considered things become tricky.

The complex plane has a natural rotational symmetry, which ties in nicely with complex arithmetic as multiplying the entire plane by e^phij is the same as a phi radian rotation around 0.

That annex G thing

Now, to keep things simple, complex fp simply uses the extensions (+/-inf, nan etc.) of the underlying real number implementation. This choice may seem so natural it isn't even perceived as a choice, but let's take a closer look at what it implies. A simple visualization of this extension of the complex plane looks like (I = infinite, f = finite, 0 = 0)

I IIIIIIIII I
             
I fffffffff I
I fffffffff I
I fffffffff I
I fffffffff I
I ffff0ffff I
I fffffffff I
I fffffffff I
I fffffffff I
I fffffffff I
             
I IIIIIIIII I

But since a true complex plane is one that respects complex multiplication, a more informative projection would be

     III    
 I         I  
    fffff    
   fffffff   
  fffffffff  
I fffffffff I
I ffff0ffff I
I fffffffff I
  fffffffff  
   fffffff   
    fffff    
 I         I 
     III    

In this projection we see the "uneven distribution" of infinities that is not only ugly but also the root of problems of the kind OP has suffered: Most infinities (those of the forms (+/-inf, finite) and (finite, +/-inf) are lumped together at the four principal directions all other directions are represented by just four infinities (+/-inf, +-inf). It shouldn't come as a surprise that extending complex multiplication to this geometry is a nightmare.

Annex G of the C99 spec tries its best to make it work, including bending the rules on how inf and nan interact (essentially inf trumps nan). OP's problem is sidestepped by not promoting reals and a proposed purely imaginary type to complex, but having the real 1 behave differently from the complex 1 doesn't strike me as a solution. Tellingly, Annex G stops short of fully specifying what the product of two infinities should be.

Can we do better?

It is tempting to try and fix these problems by choosing a better geometry of infinities. In analogy to the extended real line we could add one infinity for each direction. This construction is similar to the projective plane but doesn't lump together opposite directions. Infinities would be represented in polar coordinates inf x e^{2 omega pi i}, defining products would be straightforward. In particular, OP's problem would be solved quite naturally.

But this is where the good news ends. In a way we can be hurled back to square one by---not unreasonably---requiring that our newstyle infinities support functions that extract their real or imaginary parts. Addition is another problem; adding two nonantipodal infinities we'd have to set the angle to undefined i.e. nan (one could argue the angle must lie between the two input angles but there is no simple way of representing that "partial nan-ness")

Riemann to the rescue

In view of all this maybe the good old one point compactification is the safest thing to do. Maybe the authors of Annex G felt the same when mandating a function cproj that lumps all the infinities together.

---

Here is a related question answered by people more competent on the subject matter than I am.

Answer 3

This is an implementation detail of how complex multiplication is implemented in CPython. Unlike other languages (e.g. C or C++), CPython takes a somewhat simplistic approach:

1. ints/floats are promoted to complex numbers in multiplication
2. the simple school-formula is used, which doesn't provide desired/expected results as soon as infinite numbers are involved:

Py_complex
_Py_c_prod(Py_complex a, Py_complex b)
{
    Py_complex r;
    r.real = a.real*b.real - a.imag*b.imag;
    r.imag = a.real*b.imag + a.imag*b.real;
    return r;
}

One problematic case with the above code would be:

(0.0+1.0*j)*(inf+inf*j) = (0.0*inf-1*inf)+(0.0*inf+1.0*inf)j
                        =  nan + nan*j

However, one would like to have -inf + inf*j as result.

In this respect other languages are not far ahead: complex number multiplication was for long a time not part of the C standard, included only in C99 as appendix G, which describes how a complex multiplication should be performed - and it is not as simple as the school formula above! The C++ standard doesn't specify how complex multiplication should work, thus most compiler implementations are falling back to C-implementation, which might be C99 conforming (gcc, clang) or not (MSVC).

For the above "problematic" example, C99-compliant implementations (which are more complicated than the school formula) would give (see live) the expected result:

(0.0+1.0*j)*(inf+inf*j) = -inf + inf*j 

Even with C99 standard, an unambiguous result is not defined for all inputs and it might be different even for C99-compliant versions.

Another side effect of float not being promoted to complex in C99 is that multiplyinginf+0.0j with 1.0 or 1.0+0.0j can lead to different results (see here live):

• (inf+0.0j)*1.0 = inf+0.0j
• (inf+0.0j)*(1.0+0.0j) = inf-nanj, imaginary part being -nan and not nan (as for CPython) doesn't play a role here, because all quiet nans are equivalent (see this), even some of them have sign-bit set (and thus printed as "-", see this) and some not.

Which is at least counter-intuitive.

---

My key take-away from it is: there is nothing simple about "simple" complex number multiplication (or division) and when switching between languages or even compilers one must brace oneself for subtle bugs/differences.

Answer 4

Funny definition from Python. If we are solving this with a pen and paper I would say that expected result would be expected: (inf + 0j) as you pointed out because we know that we mean the norm of 1 so (float('inf')+0j)*1 =should= ('inf'+0j):

But that is not the case as you can see... when we run it we get:

>>> Complex( float('inf') , 0j ) * 1
result: (inf + nanj)

Python understands this *1 as a complex number and not the norm of 1 so it interprets as *(1+0j) and the error appears when we try to do inf * 0j = nanj as inf*0 can't be resolved.

What you actually want to do (assuming 1 is the norm of 1):

Recall that if z = x + iy is a complex number with real part x and imaginary part y, the complex conjugate of z is defined as z* = x − iy, and the absolute value, also called the norm of z is defined as:

Assuming 1 is the norm of 1 we should do something like:

>>> c_num = complex(float('inf'),0)
>>> value = 1
>>> realPart=(c_num.real)*value
>>> imagPart=(c_num.imag)*value
>>> complex(realPart,imagPart)
result: (inf+0j)

not very intuitive I know... but sometimes coding languages get defined in a different way from what we are used in our day to day.