minor speedup and accuracy improvement for atanh (JuliaLang#39267)

When `x` is not small, the numeric stability given by `log1p` and the rest of the complicated stuff actually just increases error and reduces speed. Also the small x and infinity case are just extra branches that already fall out of the remaining branches. The `<.5` branch is still both slower and has a higher max error (1.4 ULP vs .8 for the other branch), so if anyone can think of anything to improve it, I'm all ears. This PR is only a 20% speed improvement, but it's also a simplification so I think it's worthwhile.

Author: oscardssmith

base/special/hyperbolic.jl

@@ -243,14 +243,10 @@ function atanh(x::T) where T <: Union{Float32, Float64}
     # Method
     # 1.Reduced x to positive by atanh(-x) = -atanh(x)
     # 2. Find the branch and the expression to calculate and return it
-    #     a) 0 <= x < 2^-28
-    #         return x
-    #     b) 2^-28 <= x < 0.5
-    #         return 0.5*log1p(2x+2x*x/(1-x))
-    #     c) 0.5 <= x < 1
-    #         return 0.5*log1p(2x/1-x)
-    #     d) x = 1
-    #         return Inf
+    #     a) 0 <= x < 0.5
+    #         return 0.5*log1p(2x/(1-x))
+    #     b) 0.5 <= x <= 1
+    #         return 0.5*log((x+1)/(1-x))
     # Special cases:
     #    if |x| > 1 throw DomainError
     isnan(x) && return x
@@ -260,20 +256,12 @@ function atanh(x::T) where T <: Union{Float32, Float64}
     if absx > 1
         atanh_domain_error(x)
     end
-    if absx < T(2)^-28
-        # in a)
-        return x
-    end
     if absx < T(0.5)
+        # in a)
+        t = log1p(T(2)*absx/(T(1)-absx))
+    else
         # in b)
-        t = absx+absx
-        t = T(0.5)*log1p(t+t*absx/(T(1)-absx))
-    elseif absx < T(1)
-        # in c)
-        t = T(0.5)*log1p((absx + absx)/(T(1)-absx))
-    elseif absx == T(1)
-        # in d)
-        return copysign(T(Inf), x)
+        t = log((T(1)+absx)/(T(1)-absx))
     end
-    return copysign(t, x)
+    return T(0.5)*copysign(t, x)
 end