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 
32 
33 static const long double one = 1.0L, huge = 1e4900L;
34 
35 static const long double zero = 0.0L;
36 
37 long double
atanhl(long double x)38 atanhl(long double x)
39 {
40 	long double t;
41 	u_int32_t jx, ix;
42 	ieee_quad_shape_type u;
43 
44 	u.value = x;
45 	jx = u.parts32.mswhi;
46 	ix = jx & 0x7fffffff;
47 	u.parts32.mswhi = ix;
48 	if (ix >= 0x3fff0000) /* |x| >= 1.0 or infinity or NaN */
49 	  {
50 	    if (u.value == one)
51                 return __math_divzerol(ix != jx);
52 	    else
53                 return __math_invalidl(x);
54 	  }
55 	if(ix<0x3fc60000 && (huge+x)>zero) return x;	/* x < 2^-57 */
56 
57 	if(ix<0x3ffe0000) {		/* x < 0.5 */
58 	    t = u.value+u.value;
59 	    t = 0.5L*log1pl(t+t*u.value/(one-u.value));
60 	} else
61 	    t = 0.5L*log1pl((u.value+u.value)/(one-u.value));
62 	if(jx & 0x80000000) return -t; else return t;
63 }
64