xref: /openbsd-src/lib/libm/src/ld128/e_atanhl.c (revision 49393c004c040ee201e6408db68882c3fe4cb110)
1*49393c00Smartynas /* @(#)e_atanh.c 5.1 93/09/24 */
2*49393c00Smartynas /*
3*49393c00Smartynas  * ====================================================
4*49393c00Smartynas  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5*49393c00Smartynas  *
6*49393c00Smartynas  * Developed at SunPro, a Sun Microsystems, Inc. business.
7*49393c00Smartynas  * Permission to use, copy, modify, and distribute this
8*49393c00Smartynas  * software is freely granted, provided that this notice
9*49393c00Smartynas  * is preserved.
10*49393c00Smartynas  * ====================================================
11*49393c00Smartynas  */
12*49393c00Smartynas 
13*49393c00Smartynas /* atanhl(x)
14*49393c00Smartynas  * Method :
15*49393c00Smartynas  *    1.Reduced x to positive by atanh(-x) = -atanh(x)
16*49393c00Smartynas  *    2.For x>=0.5
17*49393c00Smartynas  *                   1              2x                          x
18*49393c00Smartynas  *	atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
19*49393c00Smartynas  *                   2             1 - x                      1 - x
20*49393c00Smartynas  *
21*49393c00Smartynas  * 	For x<0.5
22*49393c00Smartynas  *	atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x))
23*49393c00Smartynas  *
24*49393c00Smartynas  * Special cases:
25*49393c00Smartynas  *	atanhl(x) is NaN if |x| > 1 with signal;
26*49393c00Smartynas  *	atanhl(NaN) is that NaN with no signal;
27*49393c00Smartynas  *	atanhl(+-1) is +-INF with signal.
28*49393c00Smartynas  *
29*49393c00Smartynas  */
30*49393c00Smartynas 
31*49393c00Smartynas #include <math.h>
32*49393c00Smartynas 
33*49393c00Smartynas #include "math_private.h"
34*49393c00Smartynas 
35*49393c00Smartynas static const long double one = 1.0L, huge = 1e4900L;
36*49393c00Smartynas 
37*49393c00Smartynas static const long double zero = 0.0L;
38*49393c00Smartynas 
39*49393c00Smartynas long double
atanhl(long double x)40*49393c00Smartynas atanhl(long double x)
41*49393c00Smartynas {
42*49393c00Smartynas 	long double t;
43*49393c00Smartynas 	u_int32_t jx, ix;
44*49393c00Smartynas 	ieee_quad_shape_type u;
45*49393c00Smartynas 
46*49393c00Smartynas 	u.value = x;
47*49393c00Smartynas 	jx = u.parts32.mswhi;
48*49393c00Smartynas 	ix = jx & 0x7fffffff;
49*49393c00Smartynas 	u.parts32.mswhi = ix;
50*49393c00Smartynas 	if (ix >= 0x3fff0000) /* |x| >= 1.0 or infinity or NaN */
51*49393c00Smartynas 	  {
52*49393c00Smartynas 	    if (u.value == one)
53*49393c00Smartynas 	      return x/zero;
54*49393c00Smartynas 	    else
55*49393c00Smartynas 	      return (x-x)/(x-x);
56*49393c00Smartynas 	  }
57*49393c00Smartynas 	if(ix<0x3fc60000 && (huge+x)>zero) return x;	/* x < 2^-57 */
58*49393c00Smartynas 
59*49393c00Smartynas 	if(ix<0x3ffe0000) {		/* x < 0.5 */
60*49393c00Smartynas 	    t = u.value+u.value;
61*49393c00Smartynas 	    t = 0.5*log1pl(t+t*u.value/(one-u.value));
62*49393c00Smartynas 	} else
63*49393c00Smartynas 	    t = 0.5*log1pl((u.value+u.value)/(one-u.value));
64*49393c00Smartynas 	if(jx & 0x80000000) return -t; else return t;
65*49393c00Smartynas }
66