xref: /netbsd-src/external/gpl3/gcc.old/dist/libquadmath/math/atanhq.c (revision 627f7eb200a4419d89b531d55fccd2ee3ffdcde0)
1*627f7eb2Smrg /* s_atanhl.c -- long double version of s_atan.c.
2*627f7eb2Smrg  * Conversion to long double by Ulrich Drepper,
3*627f7eb2Smrg  * Cygnus Support, drepper@cygnus.com.
4*627f7eb2Smrg  */
5*627f7eb2Smrg 
6*627f7eb2Smrg /*
7*627f7eb2Smrg  * ====================================================
8*627f7eb2Smrg  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9*627f7eb2Smrg  *
10*627f7eb2Smrg  * Developed at SunPro, a Sun Microsystems, Inc. business.
11*627f7eb2Smrg  * Permission to use, copy, modify, and distribute this
12*627f7eb2Smrg  * software is freely granted, provided that this notice
13*627f7eb2Smrg  * is preserved.
14*627f7eb2Smrg  * ====================================================
15*627f7eb2Smrg  */
16*627f7eb2Smrg 
17*627f7eb2Smrg /* atanhq(x)
18*627f7eb2Smrg  * Method :
19*627f7eb2Smrg  *    1.Reduced x to positive by atanh(-x) = -atanh(x)
20*627f7eb2Smrg  *    2.For x>=0.5
21*627f7eb2Smrg  *                   1              2x                          x
22*627f7eb2Smrg  *	atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
23*627f7eb2Smrg  *                   2             1 - x                      1 - x
24*627f7eb2Smrg  *
25*627f7eb2Smrg  *	For x<0.5
26*627f7eb2Smrg  *	atanhl(x) = 0.5*log1pq(2x+2x*x/(1-x))
27*627f7eb2Smrg  *
28*627f7eb2Smrg  * Special cases:
29*627f7eb2Smrg  *	atanhl(x) is NaN if |x| > 1 with signal;
30*627f7eb2Smrg  *	atanhl(NaN) is that NaN with no signal;
31*627f7eb2Smrg  *	atanhl(+-1) is +-INF with signal.
32*627f7eb2Smrg  *
33*627f7eb2Smrg  */
34*627f7eb2Smrg 
35*627f7eb2Smrg #include "quadmath-imp.h"
36*627f7eb2Smrg 
37*627f7eb2Smrg static const __float128 one = 1, huge = 1e4900Q;
38*627f7eb2Smrg 
39*627f7eb2Smrg static const __float128 zero = 0;
40*627f7eb2Smrg 
41*627f7eb2Smrg __float128
atanhq(__float128 x)42*627f7eb2Smrg atanhq(__float128 x)
43*627f7eb2Smrg {
44*627f7eb2Smrg 	__float128 t;
45*627f7eb2Smrg 	uint32_t jx, ix;
46*627f7eb2Smrg 	ieee854_float128 u;
47*627f7eb2Smrg 
48*627f7eb2Smrg 	u.value = x;
49*627f7eb2Smrg 	jx = u.words32.w0;
50*627f7eb2Smrg 	ix = jx & 0x7fffffff;
51*627f7eb2Smrg 	u.words32.w0 = ix;
52*627f7eb2Smrg 	if (ix >= 0x3fff0000) /* |x| >= 1.0 or infinity or NaN */
53*627f7eb2Smrg 	  {
54*627f7eb2Smrg 	    if (u.value == one)
55*627f7eb2Smrg 	      return x/zero;
56*627f7eb2Smrg 	    else
57*627f7eb2Smrg 	      return (x-x)/(x-x);
58*627f7eb2Smrg 	  }
59*627f7eb2Smrg 	if(ix<0x3fc60000 && (huge+x)>zero)	/* x < 2^-57 */
60*627f7eb2Smrg 	  {
61*627f7eb2Smrg 	    math_check_force_underflow (x);
62*627f7eb2Smrg 	    return x;
63*627f7eb2Smrg 	  }
64*627f7eb2Smrg 
65*627f7eb2Smrg 	if(ix<0x3ffe0000) {		/* x < 0.5 */
66*627f7eb2Smrg 	    t = u.value+u.value;
67*627f7eb2Smrg 	    t = 0.5*log1pq(t+t*u.value/(one-u.value));
68*627f7eb2Smrg 	} else
69*627f7eb2Smrg 	    t = 0.5*log1pq((u.value+u.value)/(one-u.value));
70*627f7eb2Smrg 	if(jx & 0x80000000) return -t; else return t;
71*627f7eb2Smrg }
72