1*05a0b428SJohn Marino /* @(#)e_acosh.c 5.1 93/09/24 */
2*05a0b428SJohn Marino /*
3*05a0b428SJohn Marino * ====================================================
4*05a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5*05a0b428SJohn Marino *
6*05a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business.
7*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this
8*05a0b428SJohn Marino * software is freely granted, provided that this notice
9*05a0b428SJohn Marino * is preserved.
10*05a0b428SJohn Marino * ====================================================
11*05a0b428SJohn Marino */
12*05a0b428SJohn Marino
13*05a0b428SJohn Marino /* acosh(x)
14*05a0b428SJohn Marino * Method :
15*05a0b428SJohn Marino * Based on
16*05a0b428SJohn Marino * acosh(x) = log [ x + sqrt(x*x-1) ]
17*05a0b428SJohn Marino * we have
18*05a0b428SJohn Marino * acosh(x) := log(x)+ln2, if x is large; else
19*05a0b428SJohn Marino * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
20*05a0b428SJohn Marino * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
21*05a0b428SJohn Marino *
22*05a0b428SJohn Marino * Special cases:
23*05a0b428SJohn Marino * acosh(x) is NaN with signal if x<1.
24*05a0b428SJohn Marino * acosh(NaN) is NaN without signal.
25*05a0b428SJohn Marino */
26*05a0b428SJohn Marino
27*05a0b428SJohn Marino #include <float.h>
28*05a0b428SJohn Marino #include <math.h>
29*05a0b428SJohn Marino
30*05a0b428SJohn Marino #include "math_private.h"
31*05a0b428SJohn Marino
32*05a0b428SJohn Marino static const double
33*05a0b428SJohn Marino one = 1.0,
34*05a0b428SJohn Marino ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
35*05a0b428SJohn Marino
36*05a0b428SJohn Marino double
acosh(double x)37*05a0b428SJohn Marino acosh(double x)
38*05a0b428SJohn Marino {
39*05a0b428SJohn Marino double t;
40*05a0b428SJohn Marino int32_t hx;
41*05a0b428SJohn Marino u_int32_t lx;
42*05a0b428SJohn Marino EXTRACT_WORDS(hx,lx,x);
43*05a0b428SJohn Marino if(hx<0x3ff00000) { /* x < 1 */
44*05a0b428SJohn Marino return (x-x)/(x-x);
45*05a0b428SJohn Marino } else if(hx >=0x41b00000) { /* x > 2**28 */
46*05a0b428SJohn Marino if(hx >=0x7ff00000) { /* x is inf of NaN */
47*05a0b428SJohn Marino return x+x;
48*05a0b428SJohn Marino } else
49*05a0b428SJohn Marino return log(x)+ln2; /* acosh(huge)=log(2x) */
50*05a0b428SJohn Marino } else if(((hx-0x3ff00000)|lx)==0) {
51*05a0b428SJohn Marino return 0.0; /* acosh(1) = 0 */
52*05a0b428SJohn Marino } else if (hx > 0x40000000) { /* 2**28 > x > 2 */
53*05a0b428SJohn Marino t=x*x;
54*05a0b428SJohn Marino return log(2.0*x-one/(x+sqrt(t-one)));
55*05a0b428SJohn Marino } else { /* 1<x<2 */
56*05a0b428SJohn Marino t = x-one;
57*05a0b428SJohn Marino return log1p(t+sqrt(2.0*t+t*t));
58*05a0b428SJohn Marino }
59*05a0b428SJohn Marino }
60*05a0b428SJohn Marino
61*05a0b428SJohn Marino #if LDBL_MANT_DIG == DBL_MANT_DIG
62*05a0b428SJohn Marino __strong_alias(acoshl, acosh);
63*05a0b428SJohn Marino #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
64