1 #include "libm.h" 2 3 /* tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x)) 4 * = (exp(2*x) - 1)/(exp(2*x) - 1 + 2) 5 * = (1 - exp(-2*x))/(exp(-2*x) - 1 + 2) 6 */ 7 double tanh(double x) 8 { 9 union {double f; uint64_t i;} u = {.f = x}; 10 uint32_t w; 11 int sign; 12 double_t t; 13 14 /* x = |x| */ 15 sign = u.i >> 63; 16 u.i &= (uint64_t)-1/2; 17 x = u.f; 18 w = u.i >> 32; 19 20 if (w > 0x3fe193ea) { 21 /* |x| > log(3)/2 ~= 0.5493 or nan */ 22 if (w > 0x40340000) { 23 /* |x| > 20 or nan */ 24 /* note: this branch avoids raising overflow */ 25 t = 1 - 0/x; 26 } else { 27 t = expm1(2*x); 28 t = 1 - 2/(t+2); 29 } 30 } else if (w > 0x3fd058ae) { 31 /* |x| > log(5/3)/2 ~= 0.2554 */ 32 t = expm1(2*x); 33 t = t/(t+2); 34 } else if (w >= 0x00100000) { 35 /* |x| >= 0x1p-1022, up to 2ulp error in [0.1,0.2554] */ 36 t = expm1(-2*x); 37 t = -t/(t+2); 38 } else { 39 /* |x| is subnormal */ 40 /* note: the branch above would not raise underflow in [0x1p-1023,0x1p-1022) */ 41 FORCE_EVAL((float)x); 42 t = x; 43 } 44 return sign ? -t : t; 45 } 46