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) */ 4float atanhf(float x) 5{ 6 union {float f; uint32_t i;} u = {.f = x}; 7 unsigned s = u.i >> 31; 8 float_t y; 9 10 /* |x| */ 11 u.i &= 0x7fffffff; 12 y = u.f; 13 14 if (u.i < 0x3f800000 - (1<<23)) { 15 if (u.i < 0x3f800000 - (32<<23)) { 16 /* handle underflow */ 17 if (u.i < (1<<23)) 18 FORCE_EVAL((float)(y*y)); 19 } else { 20 /* |x| < 0.5, up to 1.7ulp error */ 21 y = 0.5f*log1pf(2*y + 2*y*y/(1-y)); 22 } 23 } else { 24 /* avoid overflow */ 25 y = 0.5f*log1pf(2*(y/(1-y))); 26 } 27 return s ? -y : y; 28} 29