1 /* md_acosh.c
2 *
3 * Inverse hyperbolic cosine
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, md_acosh();
10 *
11 * y = md_acosh( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns inverse hyperbolic cosine of argument.
18 *
19 * If 1 <= x < 1.5, a rational approximation
20 *
21 * sqrt(z) * P(z)/Q(z)
22 *
23 * where z = x-1, is used. Otherwise,
24 *
25 * md_acosh(x) = md_log( x + sqrt( (x-1)(x+1) ).
26 *
27 *
28 *
29 * ACCURACY:
30 *
31 * Relative error:
32 * arithmetic domain # trials peak rms
33 * DEC 1,3 30000 4.2e-17 1.1e-17
34 * IEEE 1,3 30000 4.6e-16 8.7e-17
35 *
36 *
37 * ERROR MESSAGES:
38 *
39 * message condition value returned
40 * md_acosh domain |x| < 1 NAN
41 *
42 */
43
44 /* md_acosh.c */
45
46 /*
47 Cephes Math Library Release 2.8: June, 2000
48 Copyright 1984, 1995, 2000 by Stephen L. Moshier
49 */
50
51
52 /* md_acosh(z) = sqrt(x) * R(x), z = x + 1, interval 0 < x < 0.5 */
53
54 #include "mconf.h"
55
56 #ifdef UNK
57 static double P[] = {
58 1.18801130533544501356E2,
59 3.94726656571334401102E3,
60 3.43989375926195455866E4,
61 1.08102874834699867335E5,
62 1.10855947270161294369E5
63 };
64 static double Q[] = {
65 /* 1.00000000000000000000E0,*/
66 1.86145380837903397292E2,
67 4.15352677227719831579E3,
68 2.97683430363289370382E4,
69 8.29725251988426222434E4,
70 7.83869920495893927727E4
71 };
72 #endif
73
74 #ifdef DEC
75 static unsigned short P[] = {
76 0041755,0115055,0144002,0146444,
77 0043166,0132103,0155150,0150302,
78 0044006,0057360,0003021,0162753,
79 0044323,0021557,0175225,0056253,
80 0044330,0101771,0040046,0006636
81 };
82 static unsigned short Q[] = {
83 /*0040200,0000000,0000000,0000000,*/
84 0042072,0022467,0126670,0041232,
85 0043201,0146066,0152142,0034015,
86 0043750,0110257,0121165,0026100,
87 0044242,0007103,0034667,0033173,
88 0044231,0014576,0175573,0017472
89 };
90 #endif
91
92 #ifdef IBMPC
93 static unsigned short P[] = {
94 0x59a4,0xb900,0xb345,0x405d,
95 0x1a18,0x7b4d,0xd688,0x40ae,
96 0x3cbd,0x00c2,0xcbde,0x40e0,
97 0xab95,0xff52,0x646d,0x40fa,
98 0xc1b4,0x2804,0x107f,0x40fb
99 };
100 static unsigned short Q[] = {
101 /*0x0000,0x0000,0x0000,0x3ff0,*/
102 0x0853,0xf5b7,0x44a6,0x4067,
103 0x4702,0xda8c,0x3986,0x40b0,
104 0xa588,0xf44e,0x1215,0x40dd,
105 0xe6cf,0x6736,0x41c8,0x40f4,
106 0x63e7,0xdf6f,0x232f,0x40f3
107 };
108 #endif
109
110 #ifdef MIEEE
111 static unsigned short P[] = {
112 0x405d,0xb345,0xb900,0x59a4,
113 0x40ae,0xd688,0x7b4d,0x1a18,
114 0x40e0,0xcbde,0x00c2,0x3cbd,
115 0x40fa,0x646d,0xff52,0xab95,
116 0x40fb,0x107f,0x2804,0xc1b4
117 };
118 static unsigned short Q[] = {
119 0x4067,0x44a6,0xf5b7,0x0853,
120 0x40b0,0x3986,0xda8c,0x4702,
121 0x40dd,0x1215,0xf44e,0xa588,
122 0x40f4,0x41c8,0x6736,0xe6cf,
123 0x40f3,0x232f,0xdf6f,0x63e7,
124 };
125 #endif
126
127 #ifdef ANSIPROT
128 extern double polevl ( double, void *, int );
129 extern double p1evl ( double, void *, int );
130 extern double md_log ( double );
131 extern double sqrt ( double );
132 #else
133 double md_log(), sqrt(), polevl(), p1evl();
134 #endif
135 extern double LOGE2, INFINITY, NAN;
136
md_acosh(x)137 double md_acosh(x)
138 double x;
139 {
140 double a, z;
141
142 if( x < 1.0 )
143 {
144 mtherr( "md_acosh", DOMAIN );
145 return(NAN);
146 }
147
148 if( x > 1.0e8 )
149 {
150 #ifdef INFINITIES
151 if( x == INFINITY )
152 return( INFINITY );
153 #endif
154 return( md_log(x) + LOGE2 );
155 }
156
157 z = x - 1.0;
158
159 if( z < 0.5 )
160 {
161 a = sqrt(z) * (polevl(z, P, 4) / p1evl(z, Q, 5) );
162 return( a );
163 }
164
165 a = sqrt( z*(x+1.0) );
166 return( md_log(x + a) );
167 }
168