| OLD | NEW | 
| (Empty) |  | 
 |   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-102
    3,0x1p-1022) */ | 
 |  41                 FORCE_EVAL((float)x); | 
 |  42                 t = x; | 
 |  43         } | 
 |  44         return sign ? -t : t; | 
 |  45 } | 
| OLD | NEW |