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 SunSoft, 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 #include <sys/cdefs.h>
16
17 #include "namespace.h"
18
19 #include <float.h>
20 #include <machine/ieee.h>
21
22 #include "math.h"
23 #include "math_private.h"
24
25 __weak_alias(acoshl, _acoshl)
26
27 #ifdef __HAVE_LONG_DOUBLE
28
29 /*
30 * See e_acosh.c for complete comments.
31 *
32 * Converted to long double by David Schultz <das@FreeBSD.ORG> and
33 * Bruce D. Evans.
34 */
35
36 /* EXP_LARGE is the threshold above which we use acosh(x) ~= log(2x). */
37 #if LDBL_MANT_DIG == 64
38 #define EXP_LARGE 34
39 #elif LDBL_MANT_DIG == 113
40 #define EXP_LARGE 58
41 #else
42 #error "Unsupported long double format"
43 #endif
44
45 #if LDBL_MAX_EXP != 0x4000
46 /* We also require the usual expsign encoding. */
47 #error "Unsupported long double format"
48 #endif
49
50 #define BIAS (LDBL_MAX_EXP - 1)
51
52 static const double
53 one = 1.0;
54
55 #if LDBL_MANT_DIG == 64
56 static const union ieee_ext_u
57 u_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
58 #define ln2 u_ln2.extu_ld
59 #elif LDBL_MANT_DIG == 113
60 static const long double
61 ln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
62 #else
63 #error "Unsupported long double format"
64 #endif
65
66 long double
acoshl(long double x)67 acoshl(long double x)
68 {
69 long double t;
70 int16_t hx;
71
72 ENTERI();
73 GET_LDBL_EXPSIGN(hx, x);
74 if (hx < 0x3fff) { /* x < 1, or misnormal */
75 RETURNI((x-x)/(x-x));
76 } else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */
77 if (hx >= 0x7fff) { /* x is inf, NaN or misnormal */
78 RETURNI(x+x);
79 } else
80 RETURNI(logl(x)+ln2); /* acosh(huge)=log(2x), or misnormal */
81 } else if (hx == 0x3fff && x == 1) {
82 RETURNI(0.0); /* acosh(1) = 0 */
83 } else if (hx >= 0x4000) { /* LARGE > x >= 2, or misnormal */
84 t=x*x;
85 RETURNI(logl(2.0*x-one/(x+sqrtl(t-one))));
86 } else { /* 1<x<2 */
87 t = x-one;
88 RETURNI(log1pl(t+sqrtl(2.0*t+t*t)));
89 }
90 }
91 #else
92 long double
93 acoshl(long double x)
94 {
95 return acosh(x);
96 }
97 #endif
98