1#include "libm.h"
2
3/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
4double atanh(double x) {
5    union {
6        double f;
7        uint64_t i;
8    } u = {.f = x};
9    unsigned e = u.i >> 52 & 0x7ff;
10    unsigned s = u.i >> 63;
11    double_t y;
12
13    /* |x| */
14    u.i &= (uint64_t)-1 / 2;
15    y = u.f;
16
17    if (e < 0x3ff - 1) {
18        if (e < 0x3ff - 32) {
19            /* handle underflow */
20            if (e == 0)
21                FORCE_EVAL((float)y);
22        } else {
23            /* |x| < 0.5, up to 1.7ulp error */
24            y = 0.5 * log1p(2 * y + 2 * y * y / (1 - y));
25        }
26    } else {
27        /* avoid overflow */
28        y = 0.5 * log1p(2 * (y / (1 - y)));
29    }
30    return s ? -y : y;
31}
32