Division of finite values by infinite complex values should always return zero

[LordZorgoth]

Since Issue 9790 was closed and it was mentioned that corner cases should be separate issues, I am reposting one particular element of that issue where Julia is currently both wrong and, according to a comment in the thread, inconsistent with C++.

Division of any finite value by any complex infinite value should return 0. Currently, we have the following behavior in Julia 0.6.0:

julia> 1/(Inf+im*Inf)
NaN + NaN*im

julia> 1/(Inf+im*NaN)
NaN + NaN*im

julia> 1/(Inf+0*im)
0.0 - 0.0im

The third case is correct, but both of the other two cases are mathematically wrong. The first case is wrong. The second case is also wrong because if NaN is replaced by some real value, or by Inf, or by -Inf, it is clear that the result should be zero in each case. Further, it is inconsistent with the behavior of C++.

[p-zubieta]

It was just recently discussed that we actually follow the ISO standard (I believe the cases you posted above also fit the standard), see this #22983 (comment). It is in fact the C and C++ compilers (at least gcc and g++) that do not follow it.

[martinholters]

Wasn't that about division by (complex) zero?

[LordZorgoth]

The discussion in that thread, as far as I can tell, was indeed about division by zero, not infinity. I can't access the complex division standard because of the paywall, so I don't know what, if anything, it says about division by infinity. It would be quite odd if it is expected to return NaN when it is clearly zero, though I suppose the reasoning would be to avoid an operation to check for infinity in the denominator when performing complex division.

Incidentally, dividing by complex zero should certainly return NaN, even from a mathematician's perspective, unless you have one ComplexInfinity, which we have already decided we aren't doing. You don't know what direction the limit to zero is coming from, so you don't know what the argument is, which for us means that we don't know whether the real and imaginary parts should be Inf, -Inf, or NaN, thus making NaN the only logical choice. That is a very distinct case from division of a finite value by a complex infinity, which should certainly return 0 from any mathematical perspective.

[LordZorgoth]

Oh, I see now; The specific algorithm for complex division is defined in the standard posted on the old issue, and it doesn't cover the special case where the denominator is infinite, which would result in NaNs being returned. Presumably it is defined that way so that the behavior of the fastest languages is consistent with that of slow languages.

[p-zubieta]

Wasn't that about division by (complex) zero?

Sorry, yes, you're right.

I can't access the complex division standard because of the paywall, so I don't know what, if anything, it says about division by infinity.

http://standards.iso.org/ittf/PubliclyAvailableStandards/. Anyway, the relevant information about complex division is already in the other issue I linked.

The specific algorithm for complex division is defined, and it doesn't cover the special case where the denominator is infinite, which would result in NaNs being returned.

It defines a fallback for (x + y*im) / (z + w*im) as ( (x*z + y*w) + (y*z + z*w)*im ) / (z*z + w*w) which will give the results on the OP if I'm not mistaken.

[p-zubieta]

Though, I'd like to point that having 1/(Inf+im*Inf) == 0.0 seems sensible. As for 1/(Inf+im*NaN) == 0.0 I'm not so sure.

[LordZorgoth]

From a mathematical perspective, you definitely do want 1/(Inf+im*NaN) to be zero. Inf+im*NaN is still an infinite number. You can't add an imaginary anything to real infinity and get something that isn't infinite. An example of where Inf+im*NaN comes from is (1+0*im)/0. That's definitely infinite. The reason it has an imaginary part of NaN is simple because we don't know the relative sizes of the zeroes in the imaginary part if we are taking some kind of limit... but that doesn't matter when we're dividing a finite number by that quantity, because we certainly do know that the real part is infinite.

Further, from a computational perspective, if you are already checking for 1/(Inf+im*Inf) to make it return zero, it is actually more efficient to make 1/(Inf+im*NaN) return zero, as you only have to find one infinity, not two.

There's certainly an argument for returning NaN for all divisions by complex infinity for the sake of speed and consistency with the standard, but if 1/(Inf+im*Inf)=0, 1/(Inf+im*NaN) should definitely be zero as well.

[simonbyrne]

For the record, both Octave and R give the equivalent of 1/Complex(Inf, Inf) == 1/Complex(Inf, NaN) == 0. It would else help deal with some edge cases of #29840.

[stevengj]

Python gives:

>>> 1/complex(inf,inf)
(nan+nanj)

ISO/IEC 10967-3:2006(E) seems to be silent on this issue: section 5.2.6.3 (Complex division) only defines the result for inputs in c(F), which is defined as complex numbers over "non-special" floating-point values (i.e. not including NaN, ±Inf, and -0), unless I am missing some discussion elsewhere in the standard.

[stevengj]

gcc 6.3 on GNU/Linux and clang on MacOS both give 0-0i from 1 / (INFINITY + INFINITY*I):

#include <complex.h>
#include <stdio.h>
#include <math.h>
int main(void) {
     double complex z = 1 / (INFINITY + INFINITY*I);
     printf("%g%+gi\n", creal(z), cimag(z));
     return 0;
}

Similarly, g++ gives 0-0i for 1.0 / std::complex<double>(INFINITY, INFINITY) and 0+0i for 1.0 / std::complex<double>(INFINITY, NAN).

[stevengj]

Overall, I can't find any reason not to give something that == 0 from division of finite value by ±Inf ± Inf*im.

[giordano]

The code above with icc:

$ icc -std=c99 complex.c && ./a.out 
complex.c(5): warning #1189: complex floating-point operation result is out of range
       double complex z = 1 / (INFINITY + INFINITY*I);
                                                  ^

-nan-nani

and this C++ code

#include <complex>
#include <iostream>

int main(void) {
  std::cout << 1.0 / std::complex<double>(INFINITY, INFINITY) << std::endl;
  return 0;
}

with icc gives:

$ icc complex.cpp && ./a.out 
(-nan,-nan)

In Mathematica, dividing 1 by ComplexInfinity gives 0:

In[1]:= 1 / ComplexInfinity

Out[1]= 0