3 /* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
6 union {double f
; uint64_t i
;} u
= {.f
= x
};
7 unsigned e
= u
.i
>> 52 & 0x7ff;
8 unsigned s
= u
.i
>> 63;
12 u
.i
&= (uint64_t)-1/2;
17 /* handle underflow */
21 /* |x| < 0.5, up to 1.7ulp error */
22 y
= 0.5*log1p(2*y
+ 2*y
*y
/(1-y
));
26 y
= 0.5*log1p(2*(y
/(1-y
)));