Home | History | Annotate | Line # | Download | only in src 
      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 #include <sys/cdefs.h>
     14 #if defined(LIBM_SCCS) && !defined(lint)
     15 __RCSID("$NetBSD: e_atanh.c,v 1.11 2002/05/26 22:01:49 wiz Exp $");
     16 #endif
     17 
     18 /* __ieee754_atanh(x)
     19  * Method :
     20  *    1.Reduced x to positive by atanh(-x) = -atanh(x)
     21  *    2.For x>=0.5
     22  *                  1              2x                          x
     23  *	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
     24  *                  2             1 - x                      1 - x
     25  *
     26  * 	For x<0.5
     27  *	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
     28  *
     29  * Special cases:
     30  *	atanh(x) is NaN if |x| > 1 with signal;
     31  *	atanh(NaN) is that NaN with no signal;
     32  *	atanh(+-1) is +-INF with signal.
     33  *
     34  */
     35 
     36 #include "math.h"
     37 #include "math_private.h"
     38 
     39 static const double one = 1.0, huge = 1e300;
     40 
     41 static const double zero = 0.0;
     42 
     43 double
     44 __ieee754_atanh(double x)
     45 {
     46 	double t;
     47 	int32_t hx,ix;
     48 	u_int32_t lx;
     49 	EXTRACT_WORDS(hx,lx,x);
     50 	ix = hx&0x7fffffff;
     51 	if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
     52 	    return (x-x)/(x-x);
     53 	if(ix==0x3ff00000)
     54 	    return x/zero;
     55 	if(ix<0x3e300000&&(huge+x)>zero) return x;	/* x<2**-28 */
     56 	SET_HIGH_WORD(x,ix);
     57 	if(ix<0x3fe00000) {		/* x < 0.5 */
     58 	    t = x+x;
     59 	    t = 0.5*log1p(t+t*x/(one-x));
     60 	} else
     61 	    t = 0.5*log1p((x+x)/(one-x));
     62 	if(hx>=0) return t; else return -t;
     63 }
     64