Alternate expression for $\operatorname{atanh}(x)$ for $\lvert x \rvert \to 1$

Question

The inverse hyperbolic tangent function $\operatorname{atanh}$ is defined as:

$$\operatorname{atanh}(x) = \frac{1}{2}\log\left(\frac{1 + x}{1 - x}\right)$$

In computers using floating point arithmetic, as $x$ approaches $+1$ or $-1$, then either $1 + x$ or $1 - x$ involve calculating differences of numbers of similar magnitude. In real terms, this causes loss of precision and inaccurate results. For example, $x$ vs. $\operatorname{atanh}(\tanh(x))$ may differ significantly for $\lvert x \rvert > 3$.

Is there an alternate expression for $\operatorname{atanh}(x)$, say for $\lvert x \rvert > 0.995$, that can reduce the error?

Answer 1

According to an old joke, you can capture an elephant by looking at it through the wrong end of a telescope, which turns the animal into a tiny elephant you can easily grab. Of course this does not work in real life, but a reversed telescoping technique does work for this problem.

Assume $x$ is close to $+1$ (for $x$ close to $-1$, simply use use the odd parity of the function to reverse the sign of $x$). Begin by applying the difference of squares factorization to render

$\dfrac{1+x}{1-x}=\dfrac{(1+x)^2}{1-x^2}.$

Then iterate:

$=\dfrac{[(1+x)^2](1+x^2)}{1-x^4}$

$=\dfrac{[(1+x)^2](1+x^2)(1+x^4)}{1-x^8}$

$=\dfrac{[(1+x)^2](1+x^2)(1+x^4)(1+x^8)}{1-x^{16}}$

$=...$

where the denominator ultimately converges to $1$, and then you have

$\tanh^{-1}(x)=\log(1+x)+\dfrac12[\log(1+x^2)+\log(1+x^4)+\log(1+x^8)+...]$

$=\log(1+x)+\dfrac12\sum_{m=1}^\infty \log(1+x^{2^m}).$

Now all the terms are well-conditioned and even for $x$ very close to $1$ the series ultimately converges rapidly. If you carry $N$ bits and thus the argument of your inverse tangent function differs from $1$ by at least $1/2^N$, then after $N+\log N/\log 2$ terms subsequent addends are sure to be too small to contribute to the computed sum. This effectively terminates the series.

To compute at least the earlier terms in the series, you may want to use this series for the natural logarithm, which is developed to converge for all positive $x$ and much more rapidly than the Taylor series when the latter converges at all:

$\log(y)=2\sum_{m=1}^\infty {\dfrac{(\frac{y-1}{y+1})^{2m-1}}{2m-1}}.$

Answer 2

You can use $$\tanh ^{-1}(1-\epsilon )=\frac{1}{2} \log \left(\frac{2-\epsilon }{\epsilon }\right)=-\frac 12 \log \left(\frac{\epsilon }{2}\right)-\sum_{n=1}^\infty\frac {\epsilon^n}{2^{n+1}\,n}$$ and use

$$\tanh ^{-1}(-1+\epsilon )=-\tanh ^{-1}(1-\epsilon )$$

The summation converges very fast and, for calculation purposes could be raplaced by some $[2p+1,2p]$ Padé approximant $P_n$ such as

$$P_1=\frac{\epsilon \left(\epsilon ^2-42 \epsilon +120\right)}{12 \left(3 \epsilon ^2-24 \epsilon +40\right)}$$ whose error is $\frac{\epsilon ^6}{76800}$

Edit

If you want to know $n$ such that $$\frac {\epsilon^n}{2^{n+1}\,n}\leq 10^{-k}$$ $$n=\left\lceil -\frac{1}{\log \left(\frac{\epsilon }{2}\right)}\,\,W\left(-\frac{10^k}{2}\,\, \log \left(\frac{\epsilon }{2}\right)\right)\right\rceil$$ where $W(.)$ is Lambert function.