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