1*998b640bSDavid Schultz /* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */
2*998b640bSDavid Schultz
3*998b640bSDavid Schultz /*
4*998b640bSDavid Schultz * ====================================================
5*998b640bSDavid Schultz * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6*998b640bSDavid Schultz *
7*998b640bSDavid Schultz * Developed at SunSoft, a Sun Microsystems, Inc. business.
8*998b640bSDavid Schultz * Permission to use, copy, modify, and distribute this
9*998b640bSDavid Schultz * software is freely granted, provided that this notice
10*998b640bSDavid Schultz * is preserved.
11*998b640bSDavid Schultz * ====================================================
12*998b640bSDavid Schultz *
13*998b640bSDavid Schultz */
14*998b640bSDavid Schultz
15*998b640bSDavid Schultz /*
16*998b640bSDavid Schultz * See e_atanh.c for complete comments.
17*998b640bSDavid Schultz *
18*998b640bSDavid Schultz * Converted to long double by David Schultz <das@FreeBSD.ORG> and
19*998b640bSDavid Schultz * Bruce D. Evans.
20*998b640bSDavid Schultz */
21*998b640bSDavid Schultz
22*998b640bSDavid Schultz #include <float.h>
23*998b640bSDavid Schultz #ifdef __i386__
24*998b640bSDavid Schultz #include <ieeefp.h>
25*998b640bSDavid Schultz #endif
26*998b640bSDavid Schultz
27*998b640bSDavid Schultz #include "fpmath.h"
28*998b640bSDavid Schultz #include "math.h"
29*998b640bSDavid Schultz #include "math_private.h"
30*998b640bSDavid Schultz
31*998b640bSDavid Schultz /* EXP_TINY is the threshold below which we use atanh(x) ~= x. */
32*998b640bSDavid Schultz #if LDBL_MANT_DIG == 64
33*998b640bSDavid Schultz #define EXP_TINY -34
34*998b640bSDavid Schultz #elif LDBL_MANT_DIG == 113
35*998b640bSDavid Schultz #define EXP_TINY -58
36*998b640bSDavid Schultz #else
37*998b640bSDavid Schultz #error "Unsupported long double format"
38*998b640bSDavid Schultz #endif
39*998b640bSDavid Schultz
40*998b640bSDavid Schultz #if LDBL_MAX_EXP != 0x4000
41*998b640bSDavid Schultz /* We also require the usual expsign encoding. */
42*998b640bSDavid Schultz #error "Unsupported long double format"
43*998b640bSDavid Schultz #endif
44*998b640bSDavid Schultz
45*998b640bSDavid Schultz #define BIAS (LDBL_MAX_EXP - 1)
46*998b640bSDavid Schultz
47*998b640bSDavid Schultz static const double one = 1.0, huge = 1e300;
48*998b640bSDavid Schultz static const double zero = 0.0;
49*998b640bSDavid Schultz
50*998b640bSDavid Schultz long double
atanhl(long double x)51*998b640bSDavid Schultz atanhl(long double x)
52*998b640bSDavid Schultz {
53*998b640bSDavid Schultz long double t;
54*998b640bSDavid Schultz uint16_t hx, ix;
55*998b640bSDavid Schultz
56*998b640bSDavid Schultz ENTERI();
57*998b640bSDavid Schultz GET_LDBL_EXPSIGN(hx, x);
58*998b640bSDavid Schultz ix = hx & 0x7fff;
59*998b640bSDavid Schultz if (ix >= 0x3fff) /* |x| >= 1, or NaN or misnormal */
60*998b640bSDavid Schultz RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x));
61*998b640bSDavid Schultz if (ix < BIAS + EXP_TINY && (huge + x) > zero)
62*998b640bSDavid Schultz RETURNI(x); /* x is tiny */
63*998b640bSDavid Schultz SET_LDBL_EXPSIGN(x, ix);
64*998b640bSDavid Schultz if (ix < 0x3ffe) { /* |x| < 0.5, or misnormal */
65*998b640bSDavid Schultz t = x+x;
66*998b640bSDavid Schultz t = 0.5*log1pl(t+t*x/(one-x));
67*998b640bSDavid Schultz } else
68*998b640bSDavid Schultz t = 0.5*log1pl((x+x)/(one-x));
69*998b640bSDavid Schultz RETURNI((hx & 0x8000) == 0 ? t : -t);
70*998b640bSDavid Schultz }
71