initial cmath code and minor libm.h update
[libm] / src / math / atanhl.c
1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_atanh.c */
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 /* atanhl(x)
13  * Method :
14  *    1.Reduced x to positive by atanh(-x) = -atanh(x)
15  *    2.For x>=0.5
16  *                   1              2x                          x
17  *      atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
18  *                   2             1 - x                      1 - x
19  *
20  *      For x<0.5
21  *      atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x))
22  *
23  * Special cases:
24  *      atanhl(x) is NaN if |x| > 1 with signal;
25  *      atanhl(NaN) is that NaN with no signal;
26  *      atanhl(+-1) is +-INF with signal.
27  */
28
29 #include "libm.h"
30
31 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
32 long double atanhl(long double x)
33 {
34         return atanh(x);
35 }
36 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
37 static const long double zero = 0.0, one = 1.0, huge = 1e4900L;
38
39 long double atanhl(long double x)
40 {
41         long double t;
42         int32_t ix;
43         uint32_t se,i0,i1;
44
45         GET_LDOUBLE_WORDS(se, i0, i1, x);
46         ix = se & 0x7fff;
47         if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31)) > 0x3fff)
48                 /* |x| > 1 */
49                 return (x-x)/(x-x);
50         if (ix == 0x3fff)
51                 return x/zero;
52         if (ix < 0x3fe3 && huge+x > zero)  /* x < 2**-28 */
53                 return x;
54         SET_LDOUBLE_EXP(x, ix);
55         if (ix < 0x3ffe) {  /* x < 0.5 */
56                 t = x + x;
57                 t = 0.5*log1pl(t + t*x/(one-x));
58         } else
59                 t = 0.5*log1pl((x + x)/(one - x));
60         if (se <= 0x7fff)
61                 return t;
62         return -t;
63 }
64 #endif