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