1210299Sed#include "libm.h"
2193323Sed
3193323Sed/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
4193323Sedfloat atanhf(float x) {
5193323Sed    union {
6193323Sed        float f;
7193323Sed        uint32_t i;
8193323Sed    } u = {.f = x};
9193323Sed    unsigned s = u.i >> 31;
10193323Sed    float_t y;
11193323Sed
12193323Sed    /* |x| */
13193323Sed    u.i &= 0x7fffffff;
14193323Sed    y = u.f;
15193323Sed
16193323Sed    if (u.i < 0x3f800000 - (1 << 23)) {
17210299Sed        if (u.i < 0x3f800000 - (32 << 23)) {
18193323Sed            /* handle underflow */
19193323Sed            if (u.i < (1 << 23))
20193323Sed                FORCE_EVAL((float)(y * y));
21193323Sed        } else {
22193323Sed            /* |x| < 0.5, up to 1.7ulp error */
23193323Sed            y = 0.5f * log1pf(2 * y + 2 * y * y / (1 - y));
24193323Sed        }
25193323Sed    } else {
26193323Sed        /* avoid overflow */
27193323Sed        y = 0.5f * log1pf(2 * (y / (1 - y)));
28193323Sed    }
29193323Sed    return s ? -y : y;
30193323Sed}
31193323Sed