Make inv(::Complex) return zero when input has Inf magnitude (JuliaLa…

…ng#32050)
Liozou · Jun 5, 2019 · 67f1a47 · 67f1a47
1 parent 0bc7bb6
commit 67f1a47
Show file tree

Hide file tree

Showing 2 changed files with 47 additions and 2 deletions.
diff --git a/base/complex.jl b/base/complex.jl
@@ -259,7 +259,11 @@ julia> conj(1 + 3im)
 conj(z::Complex) = Complex(real(z),-imag(z))
 abs(z::Complex) = hypot(real(z), imag(z))
 abs2(z::Complex) = real(z)*real(z) + imag(z)*imag(z)
-inv(z::Complex) = conj(z)/abs2(z)
+function inv(z::Complex)
+ c, d = reim(z)
+ (isinf(c) | isinf(d)) && return complex(copysign(zero(c), c), flipsign(-zero(d), d))
+ complex(c, -d)/(c * c + d * d)
+end
 inv(z::Complex{<:Integer}) = inv(float(z))
 
 +(z::Complex) = Complex(+real(z), +imag(z))
@@ -346,7 +350,7 @@ function /(a::Complex{T}, b::Complex{T}) where T<:Real
 end
 
 inv(z::Complex{<:Union{Float16,Float32}}) =
- oftype(z, conj(widen(z))/abs2(widen(z)))
+ oftype(z, inv(widen(z)))
 
 /(z::Complex{T}, w::Complex{T}) where {T<:Union{Float16,Float32}} =
  oftype(z, widen(z)*inv(widen(w)))
@@ -430,6 +434,7 @@ end
 
 function inv(w::ComplexF64)
  c, d = reim(w)
+ (isinf(c) | isinf(d)) && return complex(copysign(0.0, c), flipsign(-0.0, d))
  half = 0.5
  two = 2.0
  cd = max(abs(c), abs(d))

diff --git a/test/complex.jl b/test/complex.jl
@@ -999,6 +999,46 @@ end
  end
 end
 
+@testset "division by Inf, issue#23134" begin
+ @testset "$T" for T in (Float32, Float64, BigFloat)
+ @test isequal(one(T) / complex(T(Inf)), complex(zero(T), -zero(T)))
+ @test isequal(one(T) / complex(T(Inf), one(T)), complex(zero(T), -zero(T)))
+ @test isequal(one(T) / complex(T(Inf), T(NaN)), complex(zero(T), -zero(T)))
+ @test isequal(one(T) / complex(T(Inf), T(Inf)), complex(zero(T), -zero(T)))
+
+ @test isequal(one(T) / complex(T(-Inf)), complex(-zero(T), -zero(T)))
+ @test isequal(one(T) / complex(T(-Inf), one(T)), complex(-zero(T), -zero(T)))
+ @test isequal(one(T) / complex(T(-Inf), T(NaN)), complex(-zero(T), -zero(T)))
+ @test isequal(one(T) / complex(T(-Inf), T(Inf)), complex(-zero(T), -zero(T)))
+
+ @test isequal(one(T) / complex(T(Inf),-zero(T)), complex(zero(T), zero(T)))
+ @test isequal(one(T) / complex(T(Inf),-one(T)), complex(zero(T), zero(T)))
+ @test isequal(one(T) / complex(T(Inf),T(-NaN)), complex(zero(T), zero(T)))
+ @test isequal(one(T) / complex(T(Inf),T(-Inf)), complex(zero(T), zero(T)))
+
+ @test isequal(one(T) / complex(T(-Inf),-zero(T)),complex(-zero(T), zero(T)))
+ @test isequal(one(T) / complex(T(-Inf),-one(T)), complex(-zero(T), zero(T)))
+ @test isequal(one(T) / complex(T(-Inf),T(-NaN)), complex(-zero(T), zero(T)))
+ @test isequal(one(T) / complex(T(-Inf),T(-Inf)), complex(-zero(T), zero(T)))
+
+ @test isequal(one(T) / complex(zero(T), T(Inf)), complex(zero(T), -zero(T)))
+ @test isequal(one(T) / complex(one(T), T(Inf)), complex(zero(T), -zero(T)))
+ @test isequal(one(T) / complex(T(NaN), T(Inf)), complex(zero(T), -zero(T)))
+
+ @test isequal(one(T) / complex(zero(T), T(-Inf)), complex(zero(T), zero(T)))
+ @test isequal(one(T) / complex(one(T), T(-Inf)), complex(zero(T), zero(T)))
+ @test isequal(one(T) / complex(T(NaN), T(-Inf)), complex(zero(T), zero(T)))
+
+ @test isequal(one(T) / complex(-zero(T), T(Inf)), complex(-zero(T), -zero(T)))
+ @test isequal(one(T) / complex(-one(T), T(Inf)), complex(-zero(T), -zero(T)))
+ @test isequal(one(T) / complex(T(-NaN), T(Inf)), complex(-zero(T), -zero(T)))
+
+ @test isequal(one(T) / complex(-zero(T), T(-Inf)), complex(-zero(T), zero(T)))
+ @test isequal(one(T) / complex(-one(T), T(-Inf)), complex(-zero(T), zero(T)))
+ @test isequal(one(T) / complex(T(-NaN), T(-Inf)), complex(-zero(T), zero(T)))
+ end
+end
+
 @testset "complex^real, issue #14342" begin
  for T in (Float32, Float64, BigFloat), p in (T(-21//10), -21//10)
  z = T(2)+0im