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