gonzui

    1: /* @(#)s_tanh.c 5.1 93/09/24 */
    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: 
   13: /* LINTLIBRARY */
   14: 
   15: /* Tanh(x)
   16:  * Return the Hyperbolic Tangent of x
   17:  *
   18:  * Method :
   19:  *                                     x    -x
   20:  *                                    e  - e
   21:  *      0. tanh(x) is defined to be -----------
   22:  *                                     x    -x
   23:  *                                    e  + e
   24:  *      1. reduce x to non-negative by tanh(-x) = -tanh(x).
   25:  *      2.  0      <= x <= 2**-55 : tanh(x) := x*(one+x)
   26:  *                                              -t
   27:  *          2**-55 <  x <=  1     : tanh(x) := -----; t = expm1(-2x)
   28:  *                                             t + 2
   29:  *                                                   2
   30:  *          1      <= x <=  22.0  : tanh(x) := 1-  ----- ; t=expm1(2x)
   31:  *                                                 t + 2
   32:  *          22.0   <  x <= INF    : tanh(x) := 1.
   33:  *
   34:  * Special cases:
   35:  *      tanh(NaN) is NaN;
   36:  *      only tanh(0)=0 is exact for finite argument.
   37:  */
   38: 
   39: #include <sys/cdefs.h>
   40: #include <float.h>
   41: #include <math.h>
   42: 
   43: #include "math_private.h"
   44: 
   45: static const double one=1.0, two=2.0, tiny = 1.0e-300;
   46: 
   47: double
   48: tanh(double x)
   49: {
   50:         double t,z;
   51:         int32_t jx,ix,lx;
   52: 
   53:     /* High word of |x|. */
   54:         EXTRACT_WORDS(jx,lx,x);
   55:         ix = jx&0x7fffffff;
   56: 
   57:     /* x is INF or NaN */
   58:         if(ix>=0x7ff00000) { 
   59:             if (jx>=0) return one/x+one;    /* tanh(+-inf)=+-1 */
   60:             else       return one/x-one;    /* tanh(NaN) = NaN */
   61:         }
   62: 
   63:     /* |x| < 22 */
   64:         if (ix < 0x40360000) {         /* |x|<22 */
   65:             if ((ix | lx) == 0)
   66:                 return x;             /* x == +-0 */
   67:             if (ix<0x3c800000)                 /* |x|<2**-55 */
   68:                 return x*(one+x);     /* tanh(small) = small */
   69:             if (ix>=0x3ff00000) {      /* |x|>=1  */
   70:                 t = expm1(two*fabs(x));
   71:                 z = one - two/(t+two);
   72:             } else {
   73:                 t = expm1(-two*fabs(x));
   74:                 z= -t/(t+two);
   75:             }
   76:     /* |x| > 22, return +-1 */
   77:         } else {
   78:             z = one - tiny;            /* raised inexact flag */
   79:         }
   80:         return (jx>=0)? z: -z;
   81: }
   82: 
   83: #if     LDBL_MANT_DIG == 53
   84: #ifdef  lint
   85: /* PROTOLIB1 */
   86: long double tanhl(long double);
   87: #else   /* lint */
   88: __weak_alias(tanhl, tanh);
   89: #endif  /* lint */
   90: #endif  /* LDBL_MANT_DIG == 53 */