nsz Git - musl/blob - src/math/tanh.c

   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 }