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