1/*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12/*
13 Long double expansions are
14 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
15 and are incorporated herein by permission of the author. The author
16 reserves the right to distribute this material elsewhere under different
17 copying permissions. These modifications are distributed here under
18 the following terms:
19
20 This library is free software; you can redistribute it and/or
21 modify it under the terms of the GNU Lesser General Public
22 License as published by the Free Software Foundation; either
23 version 2.1 of the License, or (at your option) any later version.
24
25 This library is distributed in the hope that it will be useful,
26 but WITHOUT ANY WARRANTY; without even the implied warranty of
27 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28 Lesser General Public License for more details.
29
30 You should have received a copy of the GNU Lesser General Public
31 License along with this library; if not, see
32 <http://www.gnu.org/licenses/>. */
33
34/* __ieee754_acosl(x)
35 * Method :
36 * acos(x) = pi/2 - asin(x)
37 * acos(-x) = pi/2 + asin(x)
38 * For |x| <= 0.375
39 * acos(x) = pi/2 - asin(x)
40 * Between .375 and .5 the approximation is
41 * acos(0.4375 + x) = acos(0.4375) + x P(x) / Q(x)
42 * Between .5 and .625 the approximation is
43 * acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x)
44 * For x > 0.625,
45 * acos(x) = 2 asin(sqrt((1-x)/2))
46 * computed with an extended precision square root in the leading term.
47 * For x < -0.625
48 * acos(x) = pi - 2 asin(sqrt((1-|x|)/2))
49 *
50 * Special cases:
51 * if x is NaN, return x itself;
52 * if |x|>1, return NaN with invalid signal.
53 *
54 * Functions needed: __ieee754_sqrtl.
55 */
56
57#include <math.h>
58#include <math_private.h>
59
60static const _Float128
61 one = 1,
62 pio2_hi = L(1.5707963267948966192313216916397514420986),
63 pio2_lo = L(4.3359050650618905123985220130216759843812E-35),
64
65 /* acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x)
66 -0.0625 <= x <= 0.0625
67 peak relative error 3.3e-35 */
68
69 rS0 = L(5.619049346208901520945464704848780243887E0),
70 rS1 = L(-4.460504162777731472539175700169871920352E1),
71 rS2 = L(1.317669505315409261479577040530751477488E2),
72 rS3 = L(-1.626532582423661989632442410808596009227E2),
73 rS4 = L(3.144806644195158614904369445440583873264E1),
74 rS5 = L(9.806674443470740708765165604769099559553E1),
75 rS6 = L(-5.708468492052010816555762842394927806920E1),
76 rS7 = L(-1.396540499232262112248553357962639431922E1),
77 rS8 = L(1.126243289311910363001762058295832610344E1),
78 rS9 = L(4.956179821329901954211277873774472383512E-1),
79 rS10 = L(-3.313227657082367169241333738391762525780E-1),
80
81 sS0 = L(-4.645814742084009935700221277307007679325E0),
82 sS1 = L(3.879074822457694323970438316317961918430E1),
83 sS2 = L(-1.221986588013474694623973554726201001066E2),
84 sS3 = L(1.658821150347718105012079876756201905822E2),
85 sS4 = L(-4.804379630977558197953176474426239748977E1),
86 sS5 = L(-1.004296417397316948114344573811562952793E2),
87 sS6 = L(7.530281592861320234941101403870010111138E1),
88 sS7 = L(1.270735595411673647119592092304357226607E1),
89 sS8 = L(-1.815144839646376500705105967064792930282E1),
90 sS9 = L(-7.821597334910963922204235247786840828217E-2),
91 /* 1.000000000000000000000000000000000000000E0 */
92
93 acosr5625 = L(9.7338991014954640492751132535550279812151E-1),
94 pimacosr5625 = L(2.1682027434402468335351320579240000860757E0),
95
96 /* acos(0.4375 + x) = acos(0.4375) + x rS(x) / sS(x)
97 -0.0625 <= x <= 0.0625
98 peak relative error 2.1e-35 */
99
100 P0 = L(2.177690192235413635229046633751390484892E0),
101 P1 = L(-2.848698225706605746657192566166142909573E1),
102 P2 = L(1.040076477655245590871244795403659880304E2),
103 P3 = L(-1.400087608918906358323551402881238180553E2),
104 P4 = L(2.221047917671449176051896400503615543757E1),
105 P5 = L(9.643714856395587663736110523917499638702E1),
106 P6 = L(-5.158406639829833829027457284942389079196E1),
107 P7 = L(-1.578651828337585944715290382181219741813E1),
108 P8 = L(1.093632715903802870546857764647931045906E1),
109 P9 = L(5.448925479898460003048760932274085300103E-1),
110 P10 = L(-3.315886001095605268470690485170092986337E-1),
111 Q0 = L(-1.958219113487162405143608843774587557016E0),
112 Q1 = L(2.614577866876185080678907676023269360520E1),
113 Q2 = L(-9.990858606464150981009763389881793660938E1),
114 Q3 = L(1.443958741356995763628660823395334281596E2),
115 Q4 = L(-3.206441012484232867657763518369723873129E1),
116 Q5 = L(-1.048560885341833443564920145642588991492E2),
117 Q6 = L(6.745883931909770880159915641984874746358E1),
118 Q7 = L(1.806809656342804436118449982647641392951E1),
119 Q8 = L(-1.770150690652438294290020775359580915464E1),
120 Q9 = L(-5.659156469628629327045433069052560211164E-1),
121 /* 1.000000000000000000000000000000000000000E0 */
122
123 acosr4375 = L(1.1179797320499710475919903296900511518755E0),
124 pimacosr4375 = L(2.0236129215398221908706530535894517323217E0),
125
126 /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
127 0 <= x <= 0.5
128 peak relative error 1.9e-35 */
129 pS0 = L(-8.358099012470680544198472400254596543711E2),
130 pS1 = L(3.674973957689619490312782828051860366493E3),
131 pS2 = L(-6.730729094812979665807581609853656623219E3),
132 pS3 = L(6.643843795209060298375552684423454077633E3),
133 pS4 = L(-3.817341990928606692235481812252049415993E3),
134 pS5 = L(1.284635388402653715636722822195716476156E3),
135 pS6 = L(-2.410736125231549204856567737329112037867E2),
136 pS7 = L(2.219191969382402856557594215833622156220E1),
137 pS8 = L(-7.249056260830627156600112195061001036533E-1),
138 pS9 = L(1.055923570937755300061509030361395604448E-3),
139
140 qS0 = L(-5.014859407482408326519083440151745519205E3),
141 qS1 = L(2.430653047950480068881028451580393430537E4),
142 qS2 = L(-4.997904737193653607449250593976069726962E4),
143 qS3 = L(5.675712336110456923807959930107347511086E4),
144 qS4 = L(-3.881523118339661268482937768522572588022E4),
145 qS5 = L(1.634202194895541569749717032234510811216E4),
146 qS6 = L(-4.151452662440709301601820849901296953752E3),
147 qS7 = L(5.956050864057192019085175976175695342168E2),
148 qS8 = L(-4.175375777334867025769346564600396877176E1);
149 /* 1.000000000000000000000000000000000000000E0 */
150
151_Float128
152__ieee754_acosl (_Float128 x)
153{
154 _Float128 z, r, w, p, q, s, t, f2;
155 int32_t ix, sign;
156 ieee854_long_double_shape_type u;
157
158 u.value = x;
159 sign = u.parts32.w0;
160 ix = sign & 0x7fffffff;
161 u.parts32.w0 = ix; /* |x| */
162 if (ix >= 0x3fff0000) /* |x| >= 1 */
163 {
164 if (ix == 0x3fff0000
165 && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
166 { /* |x| == 1 */
167 if ((sign & 0x80000000) == 0)
168 return 0.0; /* acos(1) = 0 */
169 else
170 return (2.0 * pio2_hi) + (2.0 * pio2_lo); /* acos(-1)= pi */
171 }
172 return (x - x) / (x - x); /* acos(|x| > 1) is NaN */
173 }
174 else if (ix < 0x3ffe0000) /* |x| < 0.5 */
175 {
176 if (ix < 0x3f8e0000) /* |x| < 2**-113 */
177 return pio2_hi + pio2_lo;
178 if (ix < 0x3ffde000) /* |x| < .4375 */
179 {
180 /* Arcsine of x. */
181 z = x * x;
182 p = (((((((((pS9 * z
183 + pS8) * z
184 + pS7) * z
185 + pS6) * z
186 + pS5) * z
187 + pS4) * z
188 + pS3) * z
189 + pS2) * z
190 + pS1) * z
191 + pS0) * z;
192 q = (((((((( z
193 + qS8) * z
194 + qS7) * z
195 + qS6) * z
196 + qS5) * z
197 + qS4) * z
198 + qS3) * z
199 + qS2) * z
200 + qS1) * z
201 + qS0;
202 r = x + x * p / q;
203 z = pio2_hi - (r - pio2_lo);
204 return z;
205 }
206 /* .4375 <= |x| < .5 */
207 t = u.value - L(0.4375);
208 p = ((((((((((P10 * t
209 + P9) * t
210 + P8) * t
211 + P7) * t
212 + P6) * t
213 + P5) * t
214 + P4) * t
215 + P3) * t
216 + P2) * t
217 + P1) * t
218 + P0) * t;
219
220 q = (((((((((t
221 + Q9) * t
222 + Q8) * t
223 + Q7) * t
224 + Q6) * t
225 + Q5) * t
226 + Q4) * t
227 + Q3) * t
228 + Q2) * t
229 + Q1) * t
230 + Q0;
231 r = p / q;
232 if (sign & 0x80000000)
233 r = pimacosr4375 - r;
234 else
235 r = acosr4375 + r;
236 return r;
237 }
238 else if (ix < 0x3ffe4000) /* |x| < 0.625 */
239 {
240 t = u.value - L(0.5625);
241 p = ((((((((((rS10 * t
242 + rS9) * t
243 + rS8) * t
244 + rS7) * t
245 + rS6) * t
246 + rS5) * t
247 + rS4) * t
248 + rS3) * t
249 + rS2) * t
250 + rS1) * t
251 + rS0) * t;
252
253 q = (((((((((t
254 + sS9) * t
255 + sS8) * t
256 + sS7) * t
257 + sS6) * t
258 + sS5) * t
259 + sS4) * t
260 + sS3) * t
261 + sS2) * t
262 + sS1) * t
263 + sS0;
264 if (sign & 0x80000000)
265 r = pimacosr5625 - p / q;
266 else
267 r = acosr5625 + p / q;
268 return r;
269 }
270 else
271 { /* |x| >= .625 */
272 z = (one - u.value) * 0.5;
273 s = __ieee754_sqrtl (z);
274 /* Compute an extended precision square root from
275 the Newton iteration s -> 0.5 * (s + z / s).
276 The change w from s to the improved value is
277 w = 0.5 * (s + z / s) - s = (s^2 + z)/2s - s = (z - s^2)/2s.
278 Express s = f1 + f2 where f1 * f1 is exactly representable.
279 w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s .
280 s + w has extended precision. */
281 u.value = s;
282 u.parts32.w2 = 0;
283 u.parts32.w3 = 0;
284 f2 = s - u.value;
285 w = z - u.value * u.value;
286 w = w - 2.0 * u.value * f2;
287 w = w - f2 * f2;
288 w = w / (2.0 * s);
289 /* Arcsine of s. */
290 p = (((((((((pS9 * z
291 + pS8) * z
292 + pS7) * z
293 + pS6) * z
294 + pS5) * z
295 + pS4) * z
296 + pS3) * z
297 + pS2) * z
298 + pS1) * z
299 + pS0) * z;
300 q = (((((((( z
301 + qS8) * z
302 + qS7) * z
303 + qS6) * z
304 + qS5) * z
305 + qS4) * z
306 + qS3) * z
307 + qS2) * z
308 + qS1) * z
309 + qS0;
310 r = s + (w + s * p / q);
311
312 if (sign & 0x80000000)
313 w = pio2_hi + (pio2_lo - r);
314 else
315 w = r;
316 return 2.0 * w;
317 }
318}
319strong_alias (__ieee754_acosl, __acosl_finite)
320