I need to calculate the inverse hyperbolic tangent with great precision. The equation in question is -atanh(1/ (1 + 10**-x))
where x
should be on the order of 240 and return the approximate value of -276.65.
I tried some libraries that python provides, like numpy, math and mpmath, however, the largest "x" value I could use was 17. For "x" larger than 17, I get -inf or a domain error.
All you need is in https://en.wikipedia.org/wiki/Hyperbolic_functions.
You have tanh(x) = sinh(x)/cosh(x)
but coth(x) = cosh(x)/sinh(x)
. So instead of artanh(1/(1+ε))
you can write arcoth(1+ε)
.
Then you can express that using the logarithm. On Wikipedia you find arcoth(x) = ½ ln((x+1)/(x-1))
. So arcoth(1+ε)
becomes ½ ln((2+ε)/ε)
.
That is something you can feed into Python. Remembering that you had a negative sign in front of the whole expression you can write it as
eps = 10**(-x) # e.g. eps = 1e-240
res = -0.5 * math.log((2 + eps)/eps)
If eps is as small as you indicate, the numerator will round to 2
but the denominator should have enough precision to make the whole thing work. It does for your example: I got -276.6567847495655
where the correct result would be more like -276.65678474956545473686759
(computed using RealIntervalField(1024)
in Sage to get guarantees on numeric error) so all the digits we got from double precision arithmetic are pretty much correct.