1
2 /* @(#)e_acosh.c 5.1 93/09/24 */
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 /* acosh(x)
16 * Method :
17 * Based on
18 * acosh(x) = log [ x + sqrt(x*x-1) ]
19 * we have
20 * acosh(x) := log(x)+ln2, if x is large; else
21 * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
22 * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
23 *
24 * Special cases:
25 * acosh(x) is NaN with signal if x<1.
26 * acosh(NaN) is NaN without signal.
27 */
28
29 #include "fdlibm.h"
30
31 #ifdef _NEED_FLOAT64
32
33 static const __float64
34 one = _F_64(1.0),
35 ln2 = _F_64(6.93147180559945286227e-01); /* 0x3FE62E42, 0xFEFA39EF */
36
37 __float64
acosh64(__float64 x)38 acosh64(__float64 x)
39 {
40 __float64 t;
41 __int32_t hx;
42 __uint32_t lx;
43 EXTRACT_WORDS(hx, lx, x);
44 if (hx < 0x3ff00000) { /* x < 1 */
45 return __math_invalid(x);
46 } else if (hx >= 0x41b00000) { /* x > 2**28 */
47 if (hx >= 0x7ff00000) { /* x is inf of NaN */
48 return x + x;
49 } else
50 return log64(x) + ln2; /* acosh(huge)=log(2x) */
51 } else if (((hx - 0x3ff00000) | lx) == 0) {
52 return _F_64(0.0); /* acosh(1) = 0 */
53 } else if (hx > 0x40000000) { /* 2**28 > x > 2 */
54 t = x * x;
55 return log64(_F_64(2.0) * x - one / (x + sqrt64(t - one)));
56 } else { /* 1<x<2 */
57 t = x - one;
58 return log1p64(t + sqrt64(_F_64(2.0) * t + t * t));
59 }
60 }
61
62 _MATH_ALIAS_d_d(acosh)
63
64 #endif /* _NEED_FLOAT64 */
65
