nsz Git - libm/blob - src/math/tanhl.c

   1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_tanhl.c */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12 /* tanhl(x)
  13  * Return the Hyperbolic Tangent of x
  14  *
  15  * Method :
  16  *                                      x    -x
  17  *                                     e  - e
  18  *      0. tanhl(x) is defined to be -----------
  19  *                                      x    -x
  20  *                                     e  + e
  21  *      1. reduce x to non-negative by tanhl(-x) = -tanhl(x).
  22  *      2.  0      <= x <= 2**-55 : tanhl(x) := x*(one+x)
  23  *                                               -t
  24  *          2**-55 <  x <=  1     : tanhl(x) := -----; t = expm1l(-2x)
  25  *                                              t + 2
  26  *                                                    2
  27  *          1      <= x <=  23.0  : tanhl(x) := 1-  ----- ; t=expm1l(2x)
  28  *                                                  t + 2
  29  *          23.0   <  x <= INF    : tanhl(x) := 1.
  30  *
  31  * Special cases:
  32  *      tanhl(NaN) is NaN;
  33  *      only tanhl(0)=0 is exact for finite argument.
  34  */
  35
  36 #include "libm.h"
  37
  38 #if LD64
  39 long double tanhl(long double x)
  40 {
  41         return tanh(x);
  42 }
  43 #elif LD80
  44 static const long double one=1.0, two=2.0, tiny = 1.0e-4900L;
  45
  46 long double tanhl(long double x)
  47 {
  48         long double t,z;
  49         int32_t se;
  50         uint32_t jj0,jj1,ix;
  51
  52         /* High word of |x|. */
  53         GET_LDOUBLE_WORDS(se, jj0, jj1, x);
  54         ix = se & 0x7fff;
  55
  56         /* x is INF or NaN */
  57         if (ix == 0x7fff) {
  58                 /* for NaN it's not important which branch: tanhl(NaN) = NaN */
  59                 if (se & 0x8000)
  60                         return one/x-one;  /* tanhl(-inf)= -1; */
  61                 return one/x+one;          /* tanhl(+inf)= +1 */
  62         }
  63
  64         /* |x| < 23 */
  65         if (ix < 0x4003 || (ix == 0x4003 && jj0 < 0xb8000000u)) {
  66                 if ((ix|jj0|jj1) == 0) /* x == +- 0 */
  67                         return x;
  68                 if (ix < 0x3fc8)       /* |x| < 2**-55 */
  69                         return x*(one+tiny);  /* tanh(small) = small */
  70                 if (ix >= 0x3fff) {    /* |x| >= 1  */
  71                         t = expm1l(two*fabsl(x));
  72                         z = one - two/(t+two);
  73                 } else {
  74                         t = expm1l(-two*fabsl(x));
  75                         z = -t/(t+two);
  76                 }
  77         /* |x| > 23, return +-1 */
  78         } else {
  79                 z = one - tiny;  /* raise inexact flag */
  80         }
  81         return se & 0x8000 ? -z : z;
  82 }
  83 #endif