1131722Sobrien#include "libm.h"
2218822Sdim
3131722Sobrien/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
4131722Sobriendouble atanh(double x) {
5131722Sobrien    union {
6131722Sobrien        double f;
7131722Sobrien        uint64_t i;
8131722Sobrien    } u = {.f = x};
9131722Sobrien    unsigned e = u.i >> 52 & 0x7ff;
10131722Sobrien    unsigned s = u.i >> 63;
11131722Sobrien    double_t y;
12131722Sobrien
13131722Sobrien    /* |x| */
14131722Sobrien    u.i &= (uint64_t)-1 / 2;
15131722Sobrien    y = u.f;
16131722Sobrien
17131722Sobrien    if (e < 0x3ff - 1) {
18218822Sdim        if (e < 0x3ff - 32) {
19131722Sobrien            /* handle underflow */
20218822Sdim            if (e == 0)
21218822Sdim                FORCE_EVAL((float)y);
22218822Sdim        } else {
23218822Sdim            /* |x| < 0.5, up to 1.7ulp error */
24218822Sdim            y = 0.5 * log1p(2 * y + 2 * y * y / (1 - y));
25218822Sdim        }
26218822Sdim    } else {
27218822Sdim        /* avoid overflow */
28218822Sdim        y = 0.5 * log1p(2 * (y / (1 - y)));
29131722Sobrien    }
30131722Sobrien    return s ? -y : y;
31218822Sdim}
32