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