Julia: incorrect (Inf + Inf*im)^2.0 and other complex^float

[stevengj]

In Julia 0.4–0.7, I get:

julia> (Inf + Inf*im)^2.0
0.0 + Inf*im

which seems wrong (it assumes that the two Inf components are the "same"). The correct answer is given for ^2:

julia> (Inf + Inf*im)^2
NaN + Inf*im

The error seems to be due to this line, which dates back to #2891 by @jiahao.

My first inclination is to remove the entire p==2 special case from this ^ method. As I commented in #24497, that seems to be a performance optimization that is probably largely superseded, since most cases that need a fast z^2 probably use (a) an integer exponent (that calls a different method) and/or (b) a literal integer exponent (for which z^2 is inlined as z*z for Complex64 and Complex128).

[stevengj]

On the other hand, we have:

julia> m = realmax(Float64)
1.7976931348623157e308

julia> (m + m*im)^2      # incorrect
NaN + Inf*im

julia> (m + m*im)^2.0    # correct
0.0 + Inf*im

I don't think that it is worth the extra checks to deal with this overflow case correctly, however. You'd have to slow down all complex multiplications to make (m + m*im) * (m + m*im) give 0.0 + Inf*im.

[stevengj]

Wow, I found some even worse bugs:

julia> (0+0im)^-3.0
0.0 + 0.0im

(!!!) and

julia> (1.0+0.0im)^1e300
ERROR: InexactError: convert(Int64, 1.0e300)
Stacktrace:
 [1] convert at ./float.jl:703 [inlined]
 [2] convert at ./int.jl:530 [inlined]
 [3] convert at ./int.jl:534 [inlined]
 [4] ^(::Complex{Float64}, ::Complex{Float64}) at ./complex.jl:657
 [5] ^(::Complex{Float64}, ::Float64) at ./promotion.jl:310

[stevengj]

I'm working at a patch that fixes this, and also makes the code 60% faster for complex^complex and 120% faster for real^complex.

[stevengj]

Another bad case?

julia> Inf^(-Inf + 0.0im)
-0.0 - 0.0im

Actually, maybe that's correct: Inf^(-Inf + 0.0im) = exp((-Inf + 0.0im) * log(Inf)) = exp((-Inf + 0.0im) * Inf) = exp(-Inf + NaN*im), and it makes sense that the latter would be zero (though the sign of zero is unclear).

Python seems to get this one wrong:

In [1]: inf = float('inf')
In [2]: inf ** (-inf + 0j)
Out[2]: (nan+nanj)

[stevengj]

Closed by #24570.