1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2020 Free Software Foundation, Inc.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <https://www.gnu.org/licenses/>.
18 */
19/******************************************************************/
20/* MODULE_NAME:uasncs.c */
21/* */
22/* FUNCTIONS: uasin */
23/* uacos */
24/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */
25/* doasin.c sincos32.c dosincos.c mpa.c */
26/* sincos.tbl asincos.tbl powtwo.tbl root.tbl */
27/* */
28/* Ultimate asin/acos routines. Given an IEEE double machine */
29/* number x, compute the correctly rounded value of */
30/* arcsin(x)or arccos(x) according to the function called. */
31/* Assumption: Machine arithmetic operations are performed in */
32/* round to nearest mode of IEEE 754 standard. */
33/* */
34/******************************************************************/
35#include "endian.h"
36#include "mydefs.h"
37#include "asincos.tbl"
38#include "root.tbl"
39#include "powtwo.tbl"
40#include "MathLib.h"
41#include "uasncs.h"
42#include <float.h>
43#include <math.h>
44#include <math_private.h>
45#include <math-underflow.h>
46#include <libm-alias-finite.h>
47
48#ifndef SECTION
49# define SECTION
50#endif
51
52void __doasin(double x, double dx, double w[]);
53void __dubsin(double x, double dx, double v[]);
54void __dubcos(double x, double dx, double v[]);
55void __docos(double x, double dx, double v[]);
56double __sin32(double x, double res, double res1);
57double __cos32(double x, double res, double res1);
58
59/***************************************************************************/
60/* An ultimate asin routine. Given an IEEE double machine number x */
61/* it computes the correctly rounded (to nearest) value of arcsin(x) */
62/***************************************************************************/
63double
64SECTION
65__ieee754_asin(double x){
66 double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2];
67 mynumber u,v;
68 int4 k,m,n;
69
70 u.x = x;
71 m = u.i[HIGH_HALF];
72 k = 0x7fffffff&m; /* no sign */
73
74 if (k < 0x3e500000)
75 {
76 math_check_force_underflow (x);
77 return x; /* for x->0 => sin(x)=x */
78 }
79 /*----------------------2^-26 <= |x| < 2^ -3 -----------------*/
80 else
81 if (k < 0x3fc00000) {
82 x2 = x*x;
83 t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
84 res = x+t; /* res=arcsin(x) according to Taylor series */
85 cor = (x-res)+t;
86 if (res == res+1.025*cor) return res;
87 else {
88 x1 = x+big;
89 xx = x*x;
90 x1 -= big;
91 x2 = x - x1;
92 p = x1*x1*x1;
93 s1 = a1.x*p;
94 s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
95 ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
96 res1 = x+s1;
97 s2 = ((x-res1)+s1)+s2;
98 res = res1+s2;
99 cor = (res1-res)+s2;
100 if (res == res+1.00014*cor) return res;
101 else {
102 __doasin(x,0,w);
103 if (w[0]==(w[0]+1.00000001*w[1])) return w[0];
104 else {
105 y=fabs(x);
106 res=fabs(w[0]);
107 res1=fabs(w[0]+1.1*w[1]);
108 return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
109 }
110 }
111 }
112 }
113 /*---------------------0.125 <= |x| < 0.5 -----------------------------*/
114 else if (k < 0x3fe00000) {
115 if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
116 else n = 11*((k&0x000fffff)>>14)+352;
117 if (m>0) xx = x - asncs.x[n];
118 else xx = -x - asncs.x[n];
119 t = asncs.x[n+1]*xx;
120 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
121 +xx*asncs.x[n+6]))))+asncs.x[n+7];
122 t+=p;
123 res =asncs.x[n+8] +t;
124 cor = (asncs.x[n+8]-res)+t;
125 if (res == res+1.05*cor) return (m>0)?res:-res;
126 else {
127 r=asncs.x[n+8]+xx*asncs.x[n+9];
128 t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
129 res = r+t;
130 cor = (r-res)+t;
131 if (res == res+1.0005*cor) return (m>0)?res:-res;
132 else {
133 res1=res+1.1*cor;
134 z=0.5*(res1-res);
135 __dubsin(res,z,w);
136 z=(w[0]-fabs(x))+w[1];
137 if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
138 else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
139 else {
140 y=fabs(x);
141 return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
142 }
143 }
144 }
145 } /* else if (k < 0x3fe00000) */
146 /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
147 else
148 if (k < 0x3fe80000) {
149 n = 1056+((k&0x000fe000)>>11)*3;
150 if (m>0) xx = x - asncs.x[n];
151 else xx = -x - asncs.x[n];
152 t = asncs.x[n+1]*xx;
153 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
154 +xx*(asncs.x[n+6]+xx*asncs.x[n+7])))))+asncs.x[n+8];
155 t+=p;
156 res =asncs.x[n+9] +t;
157 cor = (asncs.x[n+9]-res)+t;
158 if (res == res+1.01*cor) return (m>0)?res:-res;
159 else {
160 r=asncs.x[n+9]+xx*asncs.x[n+10];
161 t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
162 res = r+t;
163 cor = (r-res)+t;
164 if (res == res+1.0005*cor) return (m>0)?res:-res;
165 else {
166 res1=res+1.1*cor;
167 z=0.5*(res1-res);
168 __dubsin(res,z,w);
169 z=(w[0]-fabs(x))+w[1];
170 if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
171 else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
172 else {
173 y=fabs(x);
174 return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
175 }
176 }
177 }
178 } /* else if (k < 0x3fe80000) */
179 /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
180 else
181 if (k < 0x3fed8000) {
182 n = 992+((k&0x000fe000)>>13)*13;
183 if (m>0) xx = x - asncs.x[n];
184 else xx = -x - asncs.x[n];
185 t = asncs.x[n+1]*xx;
186 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
187 +xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+xx*asncs.x[n+8]))))))+asncs.x[n+9];
188 t+=p;
189 res =asncs.x[n+10] +t;
190 cor = (asncs.x[n+10]-res)+t;
191 if (res == res+1.01*cor) return (m>0)?res:-res;
192 else {
193 r=asncs.x[n+10]+xx*asncs.x[n+11];
194 t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
195 res = r+t;
196 cor = (r-res)+t;
197 if (res == res+1.0008*cor) return (m>0)?res:-res;
198 else {
199 res1=res+1.1*cor;
200 z=0.5*(res1-res);
201 y=hp0.x-res;
202 z=((hp0.x-y)-res)+(hp1.x-z);
203 __dubcos(y,z,w);
204 z=(w[0]-fabs(x))+w[1];
205 if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
206 else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
207 else {
208 y=fabs(x);
209 return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
210 }
211 }
212 }
213 } /* else if (k < 0x3fed8000) */
214 /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
215 else
216 if (k < 0x3fee8000) {
217 n = 884+((k&0x000fe000)>>13)*14;
218 if (m>0) xx = x - asncs.x[n];
219 else xx = -x - asncs.x[n];
220 t = asncs.x[n+1]*xx;
221 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
222 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
223 +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
224 xx*asncs.x[n+9])))))))+asncs.x[n+10];
225 t+=p;
226 res =asncs.x[n+11] +t;
227 cor = (asncs.x[n+11]-res)+t;
228 if (res == res+1.01*cor) return (m>0)?res:-res;
229 else {
230 r=asncs.x[n+11]+xx*asncs.x[n+12];
231 t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
232 res = r+t;
233 cor = (r-res)+t;
234 if (res == res+1.0007*cor) return (m>0)?res:-res;
235 else {
236 res1=res+1.1*cor;
237 z=0.5*(res1-res);
238 y=(hp0.x-res)-z;
239 z=y+hp1.x;
240 y=(y-z)+hp1.x;
241 __dubcos(z,y,w);
242 z=(w[0]-fabs(x))+w[1];
243 if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
244 else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
245 else {
246 y=fabs(x);
247 return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
248 }
249 }
250 }
251 } /* else if (k < 0x3fee8000) */
252
253 /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
254 else
255 if (k < 0x3fef0000) {
256 n = 768+((k&0x000fe000)>>13)*15;
257 if (m>0) xx = x - asncs.x[n];
258 else xx = -x - asncs.x[n];
259 t = asncs.x[n+1]*xx;
260 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
261 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
262 +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
263 xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11];
264 t+=p;
265 res =asncs.x[n+12] +t;
266 cor = (asncs.x[n+12]-res)+t;
267 if (res == res+1.01*cor) return (m>0)?res:-res;
268 else {
269 r=asncs.x[n+12]+xx*asncs.x[n+13];
270 t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
271 res = r+t;
272 cor = (r-res)+t;
273 if (res == res+1.0007*cor) return (m>0)?res:-res;
274 else {
275 res1=res+1.1*cor;
276 z=0.5*(res1-res);
277 y=(hp0.x-res)-z;
278 z=y+hp1.x;
279 y=(y-z)+hp1.x;
280 __dubcos(z,y,w);
281 z=(w[0]-fabs(x))+w[1];
282 if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
283 else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
284 else {
285 y=fabs(x);
286 return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
287 }
288 }
289 }
290 } /* else if (k < 0x3fef0000) */
291 /*--------------------0.96875 <= |x| < 1 --------------------------------*/
292 else
293 if (k<0x3ff00000) {
294 z = 0.5*((m>0)?(1.0-x):(1.0+x));
295 v.x=z;
296 k=v.i[HIGH_HALF];
297 t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
298 r=1.0-t*t*z;
299 t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
300 c=t*z;
301 t=c*(1.5-0.5*t*c);
302 y=(c+t24)-t24;
303 cc = (z-y*y)/(t+y);
304 p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
305 cor = (hp1.x - 2.0*cc)-2.0*(y+cc)*p;
306 res1 = hp0.x - 2.0*y;
307 res =res1 + cor;
308 if (res == res+1.003*((res1-res)+cor)) return (m>0)?res:-res;
309 else {
310 c=y+cc;
311 cc=(y-c)+cc;
312 __doasin(c,cc,w);
313 res1=hp0.x-2.0*w[0];
314 cor=((hp0.x-res1)-2.0*w[0])+(hp1.x-2.0*w[1]);
315 res = res1+cor;
316 cor = (res1-res)+cor;
317 if (res==(res+1.0000001*cor)) return (m>0)?res:-res;
318 else {
319 y=fabs(x);
320 res1=res+1.1*cor;
321 return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
322 }
323 }
324 } /* else if (k < 0x3ff00000) */
325 /*---------------------------- |x|>=1 -------------------------------*/
326 else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x;
327 else
328 if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x + x;
329 else {
330 u.i[HIGH_HALF]=0x7ff00000;
331 v.i[HIGH_HALF]=0x7ff00000;
332 u.i[LOW_HALF]=0;
333 v.i[LOW_HALF]=0;
334 return u.x/v.x; /* NaN */
335 }
336}
337#ifndef __ieee754_asin
338libm_alias_finite (__ieee754_asin, __asin)
339#endif
340
341/*******************************************************************/
342/* */
343/* End of arcsine, below is arccosine */
344/* */
345/*******************************************************************/
346
347double
348SECTION
349__ieee754_acos(double x)
350{
351 double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps;
352 mynumber u,v;
353 int4 k,m,n;
354 u.x = x;
355 m = u.i[HIGH_HALF];
356 k = 0x7fffffff&m;
357 /*------------------- |x|<2.77556*10^-17 ----------------------*/
358 if (k < 0x3c880000) return hp0.x;
359
360 /*----------------- 2.77556*10^-17 <= |x| < 2^-3 --------------*/
361 else
362 if (k < 0x3fc00000) {
363 x2 = x*x;
364 t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
365 r=hp0.x-x;
366 cor=(((hp0.x-r)-x)+hp1.x)-t;
367 res = r+cor;
368 cor = (r-res)+cor;
369 if (res == res+1.004*cor) return res;
370 else {
371 x1 = x+big;
372 xx = x*x;
373 x1 -= big;
374 x2 = x - x1;
375 p = x1*x1*x1;
376 s1 = a1.x*p;
377 s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
378 ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
379 res1 = x+s1;
380 s2 = ((x-res1)+s1)+s2;
381 r=hp0.x-res1;
382 cor=(((hp0.x-r)-res1)+hp1.x)-s2;
383 res = r+cor;
384 cor = (r-res)+cor;
385 if (res == res+1.00004*cor) return res;
386 else {
387 __doasin(x,0,w);
388 r=hp0.x-w[0];
389 cor=((hp0.x-r)-w[0])+(hp1.x-w[1]);
390 res=r+cor;
391 cor=(r-res)+cor;
392 if (res ==(res +1.00000001*cor)) return res;
393 else {
394 res1=res+1.1*cor;
395 return __cos32(x,res,res1);
396 }
397 }
398 }
399 } /* else if (k < 0x3fc00000) */
400 /*---------------------- 0.125 <= |x| < 0.5 --------------------*/
401 else
402 if (k < 0x3fe00000) {
403 if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
404 else n = 11*((k&0x000fffff)>>14)+352;
405 if (m>0) xx = x - asncs.x[n];
406 else xx = -x - asncs.x[n];
407 t = asncs.x[n+1]*xx;
408 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
409 xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7];
410 t+=p;
411 y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]);
412 t = (m>0)?(hp1.x-t):(hp1.x+t);
413 res = y+t;
414 if (res == res+1.02*((y-res)+t)) return res;
415 else {
416 r=asncs.x[n+8]+xx*asncs.x[n+9];
417 t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
418 if (m>0)
419 {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; }
420 else
421 {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); }
422 res = p+t;
423 cor = (p-res)+t;
424 if (res == (res+1.0002*cor)) return res;
425 else {
426 res1=res+1.1*cor;
427 z=0.5*(res1-res);
428 __docos(res,z,w);
429 z=(w[0]-x)+w[1];
430 if (z>1.0e-27) return max(res,res1);
431 else if (z<-1.0e-27) return min(res,res1);
432 else return __cos32(x,res,res1);
433 }
434 }
435 } /* else if (k < 0x3fe00000) */
436
437 /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
438 else
439 if (k < 0x3fe80000) {
440 n = 1056+((k&0x000fe000)>>11)*3;
441 if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
442 else {xx = -x - asncs.x[n]; eps=1.02; }
443 t = asncs.x[n+1]*xx;
444 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
445 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+
446 xx*asncs.x[n+7])))))+asncs.x[n+8];
447 t+=p;
448 y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]);
449 t = (m>0)?(hp1.x-t):(hp1.x+t);
450 res = y+t;
451 if (res == res+eps*((y-res)+t)) return res;
452 else {
453 r=asncs.x[n+9]+xx*asncs.x[n+10];
454 t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
455 if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0004; }
456 else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0002; }
457 res = p+t;
458 cor = (p-res)+t;
459 if (res == (res+eps*cor)) return res;
460 else {
461 res1=res+1.1*cor;
462 z=0.5*(res1-res);
463 __docos(res,z,w);
464 z=(w[0]-x)+w[1];
465 if (z>1.0e-27) return max(res,res1);
466 else if (z<-1.0e-27) return min(res,res1);
467 else return __cos32(x,res,res1);
468 }
469 }
470 } /* else if (k < 0x3fe80000) */
471
472/*------------------------- 0.75 <= |x| < 0.921875 -------------*/
473 else
474 if (k < 0x3fed8000) {
475 n = 992+((k&0x000fe000)>>13)*13;
476 if (m>0) {xx = x - asncs.x[n]; eps = 1.04; }
477 else {xx = -x - asncs.x[n]; eps = 1.01; }
478 t = asncs.x[n+1]*xx;
479 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
480 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+
481 xx*asncs.x[n+8]))))))+asncs.x[n+9];
482 t+=p;
483 y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]);
484 t = (m>0)?(hp1.x-t):(hp1.x+t);
485 res = y+t;
486 if (res == res+eps*((y-res)+t)) return res;
487 else {
488 r=asncs.x[n+10]+xx*asncs.x[n+11];
489 t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
490 if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0032; }
491 else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0008; }
492 res = p+t;
493 cor = (p-res)+t;
494 if (res == (res+eps*cor)) return res;
495 else {
496 res1=res+1.1*cor;
497 z=0.5*(res1-res);
498 __docos(res,z,w);
499 z=(w[0]-x)+w[1];
500 if (z>1.0e-27) return max(res,res1);
501 else if (z<-1.0e-27) return min(res,res1);
502 else return __cos32(x,res,res1);
503 }
504 }
505 } /* else if (k < 0x3fed8000) */
506
507/*-------------------0.921875 <= |x| < 0.953125 ------------------*/
508 else
509 if (k < 0x3fee8000) {
510 n = 884+((k&0x000fe000)>>13)*14;
511 if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
512 else {xx = -x - asncs.x[n]; eps =1.005; }
513 t = asncs.x[n+1]*xx;
514 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
515 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
516 +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
517 xx*asncs.x[n+9])))))))+asncs.x[n+10];
518 t+=p;
519 y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]);
520 t = (m>0)?(hp1.x-t):(hp1.x+t);
521 res = y+t;
522 if (res == res+eps*((y-res)+t)) return res;
523 else {
524 r=asncs.x[n+11]+xx*asncs.x[n+12];
525 t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
526 if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
527 else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
528 res = p+t;
529 cor = (p-res)+t;
530 if (res == (res+eps*cor)) return res;
531 else {
532 res1=res+1.1*cor;
533 z=0.5*(res1-res);
534 __docos(res,z,w);
535 z=(w[0]-x)+w[1];
536 if (z>1.0e-27) return max(res,res1);
537 else if (z<-1.0e-27) return min(res,res1);
538 else return __cos32(x,res,res1);
539 }
540 }
541 } /* else if (k < 0x3fee8000) */
542
543 /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
544 else
545 if (k < 0x3fef0000) {
546 n = 768+((k&0x000fe000)>>13)*15;
547 if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
548 else {xx = -x - asncs.x[n]; eps=1.005;}
549 t = asncs.x[n+1]*xx;
550 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
551 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
552 +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+xx*(asncs.x[n+9]+
553 xx*asncs.x[n+10]))))))))+asncs.x[n+11];
554 t+=p;
555 y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]);
556 t = (m>0)?(hp1.x-t):(hp1.x+t);
557 res = y+t;
558 if (res == res+eps*((y-res)+t)) return res;
559 else {
560 r=asncs.x[n+12]+xx*asncs.x[n+13];
561 t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
562 if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
563 else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
564 res = p+t;
565 cor = (p-res)+t;
566 if (res == (res+eps*cor)) return res;
567 else {
568 res1=res+1.1*cor;
569 z=0.5*(res1-res);
570 __docos(res,z,w);
571 z=(w[0]-x)+w[1];
572 if (z>1.0e-27) return max(res,res1);
573 else if (z<-1.0e-27) return min(res,res1);
574 else return __cos32(x,res,res1);
575 }
576 }
577 } /* else if (k < 0x3fef0000) */
578 /*-----------------0.96875 <= |x| < 1 ---------------------------*/
579
580 else
581 if (k<0x3ff00000) {
582 z = 0.5*((m>0)?(1.0-x):(1.0+x));
583 v.x=z;
584 k=v.i[HIGH_HALF];
585 t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
586 r=1.0-t*t*z;
587 t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
588 c=t*z;
589 t=c*(1.5-0.5*t*c);
590 y = (t27*c+c)-t27*c;
591 cc = (z-y*y)/(t+y);
592 p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
593 if (m<0) {
594 cor = (hp1.x - cc)-(y+cc)*p;
595 res1 = hp0.x - y;
596 res =res1 + cor;
597 if (res == res+1.002*((res1-res)+cor)) return (res+res);
598 else {
599 c=y+cc;
600 cc=(y-c)+cc;
601 __doasin(c,cc,w);
602 res1=hp0.x-w[0];
603 cor=((hp0.x-res1)-w[0])+(hp1.x-w[1]);
604 res = res1+cor;
605 cor = (res1-res)+cor;
606 if (res==(res+1.000001*cor)) return (res+res);
607 else {
608 res=res+res;
609 res1=res+1.2*cor;
610 return __cos32(x,res,res1);
611 }
612 }
613 }
614 else {
615 cor = cc+p*(y+cc);
616 res = y + cor;
617 if (res == res+1.03*((y-res)+cor)) return (res+res);
618 else {
619 c=y+cc;
620 cc=(y-c)+cc;
621 __doasin(c,cc,w);
622 res = w[0];
623 cor=w[1];
624 if (res==(res+1.000001*cor)) return (res+res);
625 else {
626 res=res+res;
627 res1=res+1.2*cor;
628 return __cos32(x,res,res1);
629 }
630 }
631 }
632 } /* else if (k < 0x3ff00000) */
633
634 /*---------------------------- |x|>=1 -----------------------*/
635 else
636 if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x;
637 else
638 if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x + x;
639 else {
640 u.i[HIGH_HALF]=0x7ff00000;
641 v.i[HIGH_HALF]=0x7ff00000;
642 u.i[LOW_HALF]=0;
643 v.i[LOW_HALF]=0;
644 return u.x/v.x;
645 }
646}
647#ifndef __ieee754_acos
648libm_alias_finite (__ieee754_acos, __acos)
649#endif
650