gonzui

    1: /* @(#)e_atanh.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: /* atanh(x)
   16:  * Method :
   17:  *    1.Reduced x to positive by atanh(-x) = -atanh(x)
   18:  *    2.For x>=0.5
   19:  *                  1              2x                          x
   20:  *      atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
   21:  *                  2             1 - x                      1 - x
   22:  *      
   23:  *      For x<0.5
   24:  *      atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
   25:  *
   26:  * Special cases:
   27:  *      atanh(x) is NaN if |x| > 1 with signal;
   28:  *      atanh(NaN) is that NaN with no signal;
   29:  *      atanh(+-1) is +-INF with signal.
   30:  *
   31:  */
   32: 
   33: #include <sys/cdefs.h>
   34: #include <float.h>
   35: #include <math.h>
   36: 
   37: #include "math_private.h"
   38: 
   39: static const double one = 1.0, huge = 1e300;
   40: static const double zero = 0.0;
   41: 
   42: double
   43: atanh(double x)
   44: {
   45:         double t;
   46:         int32_t hx,ix;
   47:         u_int32_t lx;
   48:         EXTRACT_WORDS(hx,lx,x);
   49:         ix = hx&0x7fffffff;
   50:         if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
   51:             return (x-x)/(x-x);
   52:         if(ix==0x3ff00000) 
   53:             return x/zero;
   54:         if(ix<0x3e300000&&(huge+x)>zero) return x;     /* x<2**-28 */
   55:         SET_HIGH_WORD(x,ix);
   56:         if(ix<0x3fe00000) {            /* x < 0.5 */
   57:             t = x+x;
   58:             t = 0.5*log1p(t+t*x/(one-x));
   59:         } else 
   60:             t = 0.5*log1p((x+x)/(one-x));
   61:         if(hx>=0) return t; else return -t;
   62: }
   63: 
   64: #if     LDBL_MANT_DIG == 53
   65: #ifdef  lint
   66: /* PROTOLIB1 */
   67: long double atanhl(long double);
   68: #else   /* lint */
   69: __weak_alias(atanhl, atanh);
   70: #endif  /* lint */
   71: #endif  /* LDBL_MANT_DIG == 53 */