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;
  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                         /* raise inexact if x!=+-inf and handle nan */
  26                         t = 1 + 0/(x + 0x1p-120f);
  27                 } else {
  28                         t = expm1(2*x);
  29                         t = 1 - 2/(t+2);
  30                 }
  31         } else if (w > 0x3fd058ae) {
  32                 /* |x| > log(5/3)/2 ~= 0.2554 */
  33                 t = expm1(2*x);
  34                 t = t/(t+2);
  35         } else {
  36                 /* |x| is small, up to 2ulp error in [0.1,0.2554] */
  37                 t = expm1(-2*x);
  38                 t = -t/(t+2);
  39         }
  40         return sign ? -t : t;
  41 }