Why `inv(Complex(0.0, 0.0))` returns `Complex(NaN, NaN)`?

In Julia, inv(Complex(0.0, 0.0)) is evaluated as

julia> c = inv(Complex(0.0, 0.0))
NaN + NaN*im

julia> isnan(c)
true

julia> isinf(c)
false

julia> versioninfo()
Julia Version 1.8.4
Commit 00177ebc4fc (2022-12-23 21:32 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × AMD Ryzen 7 2700X Eight-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, znver1)
  Threads: 1 on 16 virtual cores

But I think it should be

julia> c = inv(Complex(0.0, 0.0))
Inf - Inf*im

julia> isnan(c)
false

julia> isinf(c)
true

Because

\begin{aligned} \frac{1}{a+bi} &= \frac{a}{a^2+b^2} + i \frac{-b}{a^2+b^2} \\ \lim_{a \to +0\\b \to +0}\frac{a}{a^2+b^2} &= +\infty \\ \lim_{a \to +0\\b \to +0}\frac{-b}{a^2+b^2} &= -\infty \end{aligned}

Is there any reason for the NaN value? If so, it would be nice to have more documentation about this.

Hello!

Your expressions are the result of evaluating limits. However, even your limits are not well defined, since both a and b could be approximated from the other side (negative), being both real numbers.

In computing, to deal with these situations in an uniform way, the IEEE 754 standard defines the inverse of 0.0 as NaN (not-a-number).

julia> inv(0.0)
Inf

julia> 1/0.0
Inf

Do you mean that the standard defines it for complex zeros?

1 Like

If you look at

\frac{a}{a^2+b^2}

If you let b go to zero first, then a, the limit is \infty. But if you let a go to zero first, then the expression is always zero (0/b). So the limit isn’t well defined, but for a different reason.

As for whether one approaches the limit from positive or negative, you have this:

julia> inv(0.0)
Inf

julia> inv(-0.0)
-Inf
1 Like

That is not correct. The multiplicative inverse of 0.0 with IEEE_754 is Inf.

julia> inv(0.0)
Inf

julia> inv(-0.0)
-Inf


julia> 0.0 / 0.0
NaN

There are more than two sides here in the complex plane. In the figure below, x is the real component, y is the imaginary component, and z is the angle of the complex number.

julia> using Makie, CairoMakie

julia> begin
           fig = Figure();
           ax3d = Axis3(fig[1, 1]; aspect = (1, 1, 1),
                  perspectiveness = 0.5, azimuth = 2.19, elevation = 0.9)
           x = -1:0.1:1
           y = -1:0.1:1
           z = angle.(1 ./ (x .+ y'*im))
           p = surface!(x, y, z)
           ax3d.xlabel = "Real"
           ax3d.ylabel = "Imag"
           ax3d.zlabel = "Angle of Multiplicative Inverse"
           save("multiplicative_inverse.png", fig)
       end

I think the limit argument is not very convincing here.

1 Like

If you let b go to zero first, then a, the limit is . But if you let a go to zero first, then the expression is always zero ( 0/b ).

Ahh, sorry I was stupid.

What I was trying to say is, it is better to define inv function that is commutative with abs function.

julia> c1 = Complex(3,4)
3 + 4im

julia> inv(abs(c1))
0.2

julia> abs(inv(c1))
0.2

julia> c2 = Complex(0,0)
0 + 0im

julia> inv(abs(c2))
Inf

julia> abs(inv(c2))  # This should be Inf
NaN

In the theory of complex functions, we sometimes equate all infinite points as \infty = 1/0 (See Riemannian sphere), so returning NaN is not useful here.

The choice for the return value of inv(Complex(0.0, 0.0)) will be the following:

  • Complex(NaN, NaN)
    • The current output.
    • As @mkitti mentioned, the limit of arguments does not exist. So returning NaN makes a little sense (similar to 0.0 / 0.0).
  • Complex(Inf, 0.0) or Complex(0.0, Inf)
    • If we need inv∘abs == abs∘inv (or ), these will be acceptable.
    • But choosing the position of 0.0 is not straightforward, so this is not a good way.
  • Complex(Inf, Inf)
    • This is okay.
  • Complex(Inf, -Inf)
    • This is better because this is consistent with inv(Complex(Inf,Inf)).
julia> inv(Complex(Inf,Inf))
0.0 - 0.0im

julia> inv(Complex(-Inf,Inf))
-0.0 - 0.0im

julia> inv(Complex(Inf,-Inf))
0.0 + 0.0im

julia> inv(Complex(-Inf,-Inf))
-0.0 + 0.0im

julia> (inv∘inv)(Complex(Inf,Inf))  # I also hope `inv∘inv ≈ id`
NaN + NaN*im
1 Like

OMG, true. Sorry, 7 AM…

My sense is that we are somewhat limited by our Cartesian description of complex values. If we moved to polar description, I think we might have a few more options. That is if we described the complex values in terms of abs and angle such that

θ = angle(C)
R = abs(C)
C = exp(θ*1im)

Then we could have R == ∞ for some θ dependent on the original angle.

There was a detailed discussion of this topic, including a comparison with several languages/compilers and the ISO standard in: should complex/0 division match complex/(0+0im)? · Issue #22983 · JuliaLang/julia · GitHub

TLDR: Julia’s behavior seems to be consistent with the ISO standard, and with the viewpoint that 0.0 + 0.0im represents a limit of x + iy as x,y \to 0^+ but in which the order of the limits is undefined (hence giving NaN from division in which the resultt cannot be determined). It’s also consistent with gfortran, but other languages make different choices.

10 Likes

Now that we’ve established that this already has been investigated and that we should not change this, perhaps we could discuss how to make a custom Complex type with the desired behavior?

Probably a lot easier just to define your own division operator /′, or check for zero denominators on a case-by-case basis since the desired behavior is probably application-dependent. Depends on what your practical use-case (if any) is.

2 Likes

Thank you all for the conversation!
I made a package RiemannSphereOperations.jl for complex number operations on the Riemann sphere.

julia> using RiemannSphereOperations

julia> inv′(complex(0,0))
Inf - Inf*im

julia> inv′(complex(0//1, 0//1))
1//0 - 1//0*im

julia> complex(Inf,-Inf) ==′ complex(-Inf,0)
true

julia> complex(Inf,-Inf) ==′ complex(1,0)
false
1 Like