simplify complex atanh and remove singularity perturbation (#55268)

fixes #55266, and use `inv(z)`
rather than `1/z` and use `muladd` in a couple places.

---------

Co-authored-by: Mosè Giordano <giordano@users.noreply.github.com>
(cherry picked from commit b7aa5e3)

base/complex.jl

@@ -1028,24 +1028,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 *= β
@@ -1055,16 +1053,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

test/complex.jl

@@ -1214,3 +1214,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