| 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 |  |  | /* atanhl(x) | 
    
    | 14 |  |  |  * Method : | 
    
    | 15 |  |  |  *    1.Reduced x to positive by atanh(-x) = -atanh(x) | 
    
    | 16 |  |  |  *    2.For x>=0.5 | 
    
    | 17 |  |  |  *                   1              2x                          x | 
    
    | 18 |  |  |  *	atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) | 
    
    | 19 |  |  |  *                   2             1 - x                      1 - x | 
    
    | 20 |  |  |  * | 
    
    | 21 |  |  |  *	For x<0.5 | 
    
    | 22 |  |  |  *	atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x)) | 
    
    | 23 |  |  |  * | 
    
    | 24 |  |  |  * Special cases: | 
    
    | 25 |  |  |  *	atanhl(x) is NaN if |x| > 1 with signal; | 
    
    | 26 |  |  |  *	atanhl(NaN) is that NaN with no signal; | 
    
    | 27 |  |  |  *	atanhl(+-1) is +-INF with signal. | 
    
    | 28 |  |  |  * | 
    
    | 29 |  |  |  */ | 
    
    | 30 |  |  |  | 
    
    | 31 |  |  | #include <math.h> | 
    
    | 32 |  |  |  | 
    
    | 33 |  |  | #include "math_private.h" | 
    
    | 34 |  |  |  | 
    
    | 35 |  |  | static const long double one = 1.0, huge = 1e4900L; | 
    
    | 36 |  |  |  | 
    
    | 37 |  |  | static const long double zero = 0.0; | 
    
    | 38 |  |  |  | 
    
    | 39 |  |  | long double | 
    
    | 40 |  |  | atanhl(long double x) | 
    
    | 41 |  |  | { | 
    
    | 42 |  |  | 	long double t; | 
    
    | 43 |  |  | 	int32_t ix; | 
    
    | 44 |  |  | 	u_int32_t se,i0,i1; | 
    
    | 45 |  |  | 	GET_LDOUBLE_WORDS(se,i0,i1,x); | 
    
    | 46 |  |  | 	ix = se&0x7fff; | 
    
    | 47 |  |  | 	if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31))>0x3fff) | 
    
    | 48 |  |  | 	  /* |x|>1 */ | 
    
    | 49 |  |  | 	    return (x-x)/(x-x); | 
    
    | 50 |  |  | 	if(ix==0x3fff) | 
    
    | 51 |  |  | 	    return x/zero; | 
    
    | 52 |  |  | 	if(ix<0x3fe3&&(huge+x)>zero) return x;	/* x<2**-28 */ | 
    
    | 53 |  |  | 	SET_LDOUBLE_EXP(x,ix); | 
    
    | 54 |  |  | 	if(ix<0x3ffe) {		/* x < 0.5 */ | 
    
    | 55 |  |  | 	    t = x+x; | 
    
    | 56 |  |  | 	    t = 0.5*log1pl(t+t*x/(one-x)); | 
    
    | 57 |  |  | 	} else | 
    
    | 58 |  |  | 	    t = 0.5*log1pl((x+x)/(one-x)); | 
    
    | 59 |  |  | 	if(se<=0x7fff) return t; else return -t; | 
    
    | 60 |  |  | } |