1*05a0b428SJohn Marino /* @(#)e_atanh.c 5.1 93/09/24 */
2*05a0b428SJohn Marino /*
3*05a0b428SJohn Marino * ====================================================
4*05a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5*05a0b428SJohn Marino *
6*05a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business.
7*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this
8*05a0b428SJohn Marino * software is freely granted, provided that this notice
9*05a0b428SJohn Marino * is preserved.
10*05a0b428SJohn Marino * ====================================================
11*05a0b428SJohn Marino */
12*05a0b428SJohn Marino
13*05a0b428SJohn Marino /* atanh(x)
14*05a0b428SJohn Marino * Method :
15*05a0b428SJohn Marino * 1.Reduced x to positive by atanh(-x) = -atanh(x)
16*05a0b428SJohn Marino * 2.For x>=0.5
17*05a0b428SJohn Marino * 1 2x x
18*05a0b428SJohn Marino * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
19*05a0b428SJohn Marino * 2 1 - x 1 - x
20*05a0b428SJohn Marino *
21*05a0b428SJohn Marino * For x<0.5
22*05a0b428SJohn Marino * atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
23*05a0b428SJohn Marino *
24*05a0b428SJohn Marino * Special cases:
25*05a0b428SJohn Marino * atanh(x) is NaN if |x| > 1 with signal;
26*05a0b428SJohn Marino * atanh(NaN) is that NaN with no signal;
27*05a0b428SJohn Marino * atanh(+-1) is +-INF with signal.
28*05a0b428SJohn Marino *
29*05a0b428SJohn Marino */
30*05a0b428SJohn Marino
31*05a0b428SJohn Marino #include <float.h>
32*05a0b428SJohn Marino #include <math.h>
33*05a0b428SJohn Marino
34*05a0b428SJohn Marino #include "math_private.h"
35*05a0b428SJohn Marino
36*05a0b428SJohn Marino static const double one = 1.0, huge = 1e300;
37*05a0b428SJohn Marino static const double zero = 0.0;
38*05a0b428SJohn Marino
39*05a0b428SJohn Marino double
atanh(double x)40*05a0b428SJohn Marino atanh(double x)
41*05a0b428SJohn Marino {
42*05a0b428SJohn Marino double t;
43*05a0b428SJohn Marino int32_t hx,ix;
44*05a0b428SJohn Marino u_int32_t lx;
45*05a0b428SJohn Marino EXTRACT_WORDS(hx,lx,x);
46*05a0b428SJohn Marino ix = hx&0x7fffffff;
47*05a0b428SJohn Marino if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
48*05a0b428SJohn Marino return (x-x)/(x-x);
49*05a0b428SJohn Marino if(ix==0x3ff00000)
50*05a0b428SJohn Marino return x/zero;
51*05a0b428SJohn Marino if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */
52*05a0b428SJohn Marino SET_HIGH_WORD(x,ix);
53*05a0b428SJohn Marino if(ix<0x3fe00000) { /* x < 0.5 */
54*05a0b428SJohn Marino t = x+x;
55*05a0b428SJohn Marino t = 0.5*log1p(t+t*x/(one-x));
56*05a0b428SJohn Marino } else
57*05a0b428SJohn Marino t = 0.5*log1p((x+x)/(one-x));
58*05a0b428SJohn Marino if(hx>=0) return t; else return -t;
59*05a0b428SJohn Marino }
60*05a0b428SJohn Marino
61*05a0b428SJohn Marino #if LDBL_MANT_DIG == DBL_MANT_DIG
62*05a0b428SJohn Marino __strong_alias(atanhl, atanh);
63*05a0b428SJohn Marino #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
64