]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/e_asin.c
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2017 Free Software Foundation, Inc.
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.
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.
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/>.
19 /******************************************************************/
20 /* MODULE_NAME:uasncs.c */
22 /* FUNCTIONS: uasin */
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 */
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. */
34 /******************************************************************/
37 #include "asincos.tbl"
44 #include <math_private.h>
50 void __doasin(double x
, double dx
, double w
[]);
51 void __dubsin(double x
, double dx
, double v
[]);
52 void __dubcos(double x
, double dx
, double v
[]);
53 void __docos(double x
, double dx
, double v
[]);
54 double __sin32(double x
, double res
, double res1
);
55 double __cos32(double x
, double res
, double res1
);
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 /***************************************************************************/
63 __ieee754_asin(double x
){
64 double x1
,x2
,xx
,s1
,s2
,res1
,p
,t
,res
,r
,cor
,cc
,y
,c
,z
,w
[2];
70 k
= 0x7fffffff&m
; /* no sign */
74 math_check_force_underflow (x
);
75 return x
; /* for x->0 => sin(x)=x */
77 /*----------------------2^-26 <= |x| < 2^ -3 -----------------*/
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 */
84 if (res
== res
+1.025*cor
) return res
;
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
;
95 s2
= ((x
-res1
)+s1
)+s2
;
98 if (res
== res
+1.00014*cor
) return res
;
101 if (w
[0]==(w
[0]+1.00000001*w
[1])) return w
[0];
105 res1
=fabs(w
[0]+1.1*w
[1]);
106 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
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
];
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];
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
;
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]);
129 if (res
== res
+1.0005*cor
) return (m
>0)?res
:-res
;
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
);
139 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
143 } /* else if (k < 0x3fe00000) */
144 /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
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
];
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];
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
;
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]);
162 if (res
== res
+1.0005*cor
) return (m
>0)?res
:-res
;
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
);
172 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
176 } /* else if (k < 0x3fe80000) */
177 /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
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
];
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];
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
;
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]);
195 if (res
== res
+1.0008*cor
) return (m
>0)?res
:-res
;
200 z
=((hp0
.x
-y
)-res
)+(hp1
.x
-z
);
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
);
207 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
211 } /* else if (k < 0x3fed8000) */
212 /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
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
];
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];
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
;
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]);
232 if (res
== res
+1.0007*cor
) return (m
>0)?res
:-res
;
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
);
245 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
249 } /* else if (k < 0x3fee8000) */
251 /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
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
];
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];
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
;
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]);
271 if (res
== res
+1.0007*cor
) return (m
>0)?res
:-res
;
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
);
284 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
288 } /* else if (k < 0x3fef0000) */
289 /*--------------------0.96875 <= |x| < 1 --------------------------------*/
292 z
= 0.5*((m
>0)?(1.0-x
):(1.0+x
));
295 t
=inroot
[(k
&0x001fffff)>>14]*powtwo
[511-(k
>>21)];
297 t
= t
*(rt0
+r
*(rt1
+r
*(rt2
+r
*rt3
)));
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
;
306 if (res
== res
+1.003*((res1
-res
)+cor
)) return (m
>0)?res
:-res
;
312 cor
=((hp0
.x
-res1
)-2.0*w
[0])+(hp1
.x
-2.0*w
[1]);
314 cor
= (res1
-res
)+cor
;
315 if (res
==(res
+1.0000001*cor
)) return (m
>0)?res
:-res
;
319 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
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
;
326 if (k
>0x7ff00000 || (k
== 0x7ff00000 && u
.i
[LOW_HALF
] != 0)) return x
+ x
;
328 u
.i
[HIGH_HALF
]=0x7ff00000;
329 v
.i
[HIGH_HALF
]=0x7ff00000;
332 return u
.x
/v
.x
; /* NaN */
335 #ifndef __ieee754_asin
336 strong_alias (__ieee754_asin
, __asin_finite
)
339 /*******************************************************************/
341 /* End of arcsine, below is arccosine */
343 /*******************************************************************/
347 __ieee754_acos(double x
)
349 double x1
,x2
,xx
,s1
,s2
,res1
,p
,t
,res
,r
,cor
,cc
,y
,c
,z
,w
[2],eps
;
355 /*------------------- |x|<2.77556*10^-17 ----------------------*/
356 if (k
< 0x3c880000) return hp0
.x
;
358 /*----------------- 2.77556*10^-17 <= |x| < 2^-3 --------------*/
360 if (k
< 0x3fc00000) {
362 t
= (((((f6
*x2
+ f5
)*x2
+ f4
)*x2
+ f3
)*x2
+ f2
)*x2
+ f1
)*(x2
*x
);
364 cor
=(((hp0
.x
-r
)-x
)+hp1
.x
)-t
;
367 if (res
== res
+1.004*cor
) return res
;
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
;
378 s2
= ((x
-res1
)+s1
)+s2
;
380 cor
=(((hp0
.x
-r
)-res1
)+hp1
.x
)-s2
;
383 if (res
== res
+1.00004*cor
) return res
;
387 cor
=((hp0
.x
-r
)-w
[0])+(hp1
.x
-w
[1]);
390 if (res
==(res
+1.00000001*cor
)) return res
;
393 return __cos32(x
,res
,res1
);
397 } /* else if (k < 0x3fc00000) */
398 /*---------------------- 0.125 <= |x| < 0.5 --------------------*/
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
];
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];
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
);
412 if (res
== res
+1.02*((y
-res
)+t
)) return res
;
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]);
417 {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; }
419 {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); }
422 if (res
== (res
+1.0002*cor
)) return res
;
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
);
433 } /* else if (k < 0x3fe00000) */
435 /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
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; }
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];
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
);
449 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
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; }
457 if (res
== (res
+eps
*cor
)) return res
;
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
);
468 } /* else if (k < 0x3fe80000) */
470 /*------------------------- 0.75 <= |x| < 0.921875 -------------*/
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; }
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];
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
);
484 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
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; }
492 if (res
== (res
+eps
*cor
)) return res
;
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
);
503 } /* else if (k < 0x3fed8000) */
505 /*-------------------0.921875 <= |x| < 0.953125 ------------------*/
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; }
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];
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
);
520 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
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; }
528 if (res
== (res
+eps
*cor
)) return res
;
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
);
539 } /* else if (k < 0x3fee8000) */
541 /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
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;}
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];
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
);
556 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
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; }
564 if (res
== (res
+eps
*cor
)) return res
;
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
);
575 } /* else if (k < 0x3fef0000) */
576 /*-----------------0.96875 <= |x| < 1 ---------------------------*/
580 z
= 0.5*((m
>0)?(1.0-x
):(1.0+x
));
583 t
=inroot
[(k
&0x001fffff)>>14]*powtwo
[511-(k
>>21)];
585 t
= t
*(rt0
+r
*(rt1
+r
*(rt2
+r
*rt3
)));
590 p
=(((((f6
*z
+f5
)*z
+f4
)*z
+f3
)*z
+f2
)*z
+f1
)*z
;
592 cor
= (hp1
.x
- cc
)-(y
+cc
)*p
;
595 if (res
== res
+1.002*((res1
-res
)+cor
)) return (res
+res
);
601 cor
=((hp0
.x
-res1
)-w
[0])+(hp1
.x
-w
[1]);
603 cor
= (res1
-res
)+cor
;
604 if (res
==(res
+1.000001*cor
)) return (res
+res
);
608 return __cos32(x
,res
,res1
);
615 if (res
== res
+1.03*((y
-res
)+cor
)) return (res
+res
);
622 if (res
==(res
+1.000001*cor
)) return (res
+res
);
626 return __cos32(x
,res
,res1
);
630 } /* else if (k < 0x3ff00000) */
632 /*---------------------------- |x|>=1 -----------------------*/
634 if (k
==0x3ff00000 && u
.i
[LOW_HALF
]==0) return (m
>0)?0:2.0*hp0
.x
;
636 if (k
>0x7ff00000 || (k
== 0x7ff00000 && u
.i
[LOW_HALF
] != 0)) return x
+ x
;
638 u
.i
[HIGH_HALF
]=0x7ff00000;
639 v
.i
[HIGH_HALF
]=0x7ff00000;
645 #ifndef __ieee754_acos
646 strong_alias (__ieee754_acos
, __acos_finite
)