real(atanh(1+y*im)) should be asymptotic to log(2/abs(y))/2 for real y as y approches 0.

For example, this suggests that real(atanh(1+1e-200im)) should evaluate to log(2/1e-200)/2=230.60508288968455, but it currently evaluates to 177.09910463306602, which has no correct digits.

This appears to be caused by a small positive perturbation, ρ, that is added to the imaginary part of the input in the atanh computation

so that real(atanh(1+1e-200im)) in fact computes log(2/ρ)=177.09910463306602.

This perturbation is suggested by Kahan 86, but as far as I can tell, it is not necessary in either of the places it is used, and in fact harms rather than helps accuracy when the real part of the input is 1. It may be there to avoid division by 0 in atanh(1.0+0.0im), but the Julia implementation already handles this case with an explicit branch.

For comparison, cpython computes the correct answer:

> python3
Python 3.12.3 (main, Apr  9 2024, 08:09:14) [Clang 15.0.0 (clang-1500.3.9.4)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import cmath
>>> cmath.atanh(1+1e-200j)
(230.60508288968455+0.7853981633974483j)

Version Info:

julia> versioninfo()
Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 12 × Apple M2 Pro
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)