1 /* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */
2
3 /*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14 /*
15 * See s_asinh.c for complete comments.
16 *
17 * Converted to long double by David Schultz <das@FreeBSD.ORG> and
18 * Bruce D. Evans.
19 */
20
21 #include <float.h>
22 #ifdef __i386__
23 #include <ieeefp.h>
24 #endif
25
26 #include "fpmath.h"
27 #include "math.h"
28 #include "math_private.h"
29
30 /* EXP_LARGE is the threshold above which we use asinh(x) ~= log(2x). */
31 /* EXP_TINY is the threshold below which we use asinh(x) ~= x. */
32 #if LDBL_MANT_DIG == 64
33 #define EXP_LARGE 34
34 #define EXP_TINY -34
35 #elif LDBL_MANT_DIG == 113
36 #define EXP_LARGE 58
37 #define EXP_TINY -58
38 #else
39 #error "Unsupported long double format"
40 #endif
41
42 #if LDBL_MAX_EXP != 0x4000
43 /* We also require the usual expsign encoding. */
44 #error "Unsupported long double format"
45 #endif
46
47 #define BIAS (LDBL_MAX_EXP - 1)
48
49 static const double
50 one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
51 huge= 1.00000000000000000000e+300;
52
53 #if LDBL_MANT_DIG == 64
54 static const union IEEEl2bits
55 u_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
56 #define ln2 u_ln2.e
57 #elif LDBL_MANT_DIG == 113
58 static const long double
59 ln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
60 #else
61 #error "Unsupported long double format"
62 #endif
63
64 long double
asinhl(long double x)65 asinhl(long double x)
66 {
67 long double t, w;
68 uint16_t hx, ix;
69
70 ENTERI();
71 GET_LDBL_EXPSIGN(hx, x);
72 ix = hx & 0x7fff;
73 if (ix >= 0x7fff) RETURNI(x+x); /* x is inf, NaN or misnormal */
74 if (ix < BIAS + EXP_TINY) { /* |x| < TINY, or misnormal */
75 if (huge + x > one) RETURNI(x); /* return x inexact except 0 */
76 }
77 if (ix >= BIAS + EXP_LARGE) { /* |x| >= LARGE, or misnormal */
78 w = logl(fabsl(x))+ln2;
79 } else if (ix >= 0x4000) { /* LARGE > |x| >= 2.0, or misnormal */
80 t = fabsl(x);
81 w = logl(2.0*t+one/(sqrtl(x*x+one)+t));
82 } else { /* 2.0 > |x| >= TINY, or misnormal */
83 t = x*x;
84 w =log1pl(fabsl(x)+t/(one+sqrtl(one+t)));
85 }
86 RETURNI((hx & 0x8000) == 0 ? w : -w);
87 }
88