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