Update.
[glibc.git] / sysdeps / libm-ieee754 / e_atanh.c
blobfa4fe675c98e7ee7fc9e4029742b2b124b4daa15
1 /* @(#)e_atanh.c 5.1 93/09/24 */
2 /*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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 * ====================================================
13 #if defined(LIBM_SCCS) && !defined(lint)
14 static char rcsid[] = "$NetBSD: e_atanh.c,v 1.8 1995/05/10 20:44:55 jtc Exp $";
15 #endif
17 /* __ieee754_atanh(x)
18 * Method :
19 * 1.Reduced x to positive by atanh(-x) = -atanh(x)
20 * 2.For x>=0.5
21 * 1 2x x
22 * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
23 * 2 1 - x 1 - x
25 * For x<0.5
26 * atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
28 * Special cases:
29 * atanh(x) is NaN if |x| > 1 with signal;
30 * atanh(NaN) is that NaN with no signal;
31 * atanh(+-1) is +-INF with signal.
35 #include "math.h"
36 #include "math_private.h"
38 #ifdef __STDC__
39 static const double one = 1.0, huge = 1e300;
40 #else
41 static double one = 1.0, huge = 1e300;
42 #endif
44 #ifdef __STDC__
45 static const double zero = 0.0;
46 #else
47 static double zero = 0.0;
48 #endif
50 #ifdef __STDC__
51 double __ieee754_atanh(double x)
52 #else
53 double __ieee754_atanh(x)
54 double x;
55 #endif
57 double t;
58 int32_t hx,ix;
59 u_int32_t lx;
60 EXTRACT_WORDS(hx,lx,x);
61 ix = hx&0x7fffffff;
62 if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
63 return (x-x)/(x-x);
64 if(ix==0x3ff00000)
65 return x/zero;
66 if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */
67 SET_HIGH_WORD(x,ix);
68 if(ix<0x3fe00000) { /* x < 0.5 */
69 t = x+x;
70 t = 0.5*__log1p(t+t*x/(one-x));
71 } else
72 t = 0.5*__log1p((x+x)/(one-x));
73 if(hx>=0) return t; else return -t;