diff --git a/base/complex.jl b/base/complex.jl index 8ac126d2c6532..095c842795d38 100644 --- a/base/complex.jl +++ b/base/complex.jl @@ -1037,24 +1037,22 @@ end function atanh(z::Complex{T}) where T z = float(z) Tf = float(T) - Ω = prevfloat(typemax(Tf)) - θ = sqrt(Ω)/4 - ρ = 1/θ x, y = reim(z) ax = abs(x) ay = abs(y) + θ = sqrt(floatmax(Tf))/4 if ax > θ || ay > θ #Prevent overflow if isnan(y) if isinf(x) return Complex(copysign(zero(x),x), y) else - return Complex(real(1/z), y) + return Complex(real(inv(z)), y) end end if isinf(y) return Complex(copysign(zero(x),x), copysign(oftype(y,pi)/2, y)) end - return Complex(real(1/z), copysign(oftype(y,pi)/2, y)) + return Complex(real(inv(z)), copysign(oftype(y,pi)/2, y)) end β = copysign(one(Tf), x) z *= β @@ -1064,16 +1062,15 @@ function atanh(z::Complex{T}) where T ξ = oftype(x, Inf) η = y else - ym = ay+ρ - ξ = log(sqrt(sqrt(4+y*y))/sqrt(ym)) - η = copysign(oftype(y,pi)/2 + atan(ym/2), y)/2 + ξ = log(sqrt(sqrt(muladd(y, y, 4)))/sqrt(ay)) + η = copysign(oftype(y,pi)/2 + atan(ay/2), y)/2 end else #Normal case - ysq = (ay+ρ)^2 + ysq = ay^2 if x == 0 ξ = x else - ξ = log1p(4x/((1-x)^2 + ysq))/4 + ξ = log1p(4x/(muladd(1-x, 1-x, ysq)))/4 end η = angle(Complex((1-x)*(1+x)-ysq, 2y))/2 end diff --git a/test/complex.jl b/test/complex.jl index d798cfe16489c..63304652ee7d8 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -1215,3 +1215,9 @@ end @test !iseven(7+0im) && isodd(7+0im) @test !iseven(6+1im) && !isodd(7+1im) end + +@testset "issue #55266" begin + for T in (Float16, Float32, Float64) + @test isapprox(atanh(1+im*floatmin(T)), Complex{T}(atanh(1+im*big(floatmin(T))))) + end +end