2 * IBM Accurate Mathematical Library
3 * Copyright (c) International Business Machines Corp., 2001
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU Lesser General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
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"
43 void __doasin(double x
, double dx
, double w
[]);
44 void __dubsin(double x
, double dx
, double v
[]);
45 void __dubcos(double x
, double dx
, double v
[]);
46 void __docos(double x
, double dx
, double v
[]);
47 double __sin32(double x
, double res
, double res1
);
48 double __cos32(double x
, double res
, double res1
);
50 /***************************************************************************/
51 /* An ultimate asin routine. Given an IEEE double machine number x */
52 /* it computes the correctly rounded (to nearest) value of arcsin(x) */
53 /***************************************************************************/
54 double __ieee754_asin(double x
){
55 double x1
,x2
,xx
,s1
,s2
,res1
,p
,t
,res
,r
,cor
,cc
,y
,c
,z
,w
[2];
64 k
= 0x7fffffff&m
; /* no sign */
66 if (k
< 0x3e500000) return x
; /* for x->0 => sin(x)=x */
67 /*----------------------2^-26 <= |x| < 2^ -3 -----------------*/
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 */
74 if (res
== res
+1.025*cor
) return res
;
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
;
85 s2
= ((x
-res1
)+s1
)+s2
;
88 if (res
== res
+1.00014*cor
) return res
;
91 if (w
[0]==(w
[0]+1.00000001*w
[1])) return w
[0];
95 res1
=ABS(w
[0]+1.1*w
[1]);
96 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
101 /*---------------------0.125 <= |x| < 0.5 -----------------------------*/
102 else if (k
< 0x3fe00000) {
103 if (k
<0x3fd00000) n
= 11*((k
&0x000fffff)>>15);
104 else n
= 11*((k
&0x000fffff)>>14)+352;
105 if (m
>0) xx
= x
- asncs
.x
[n
];
106 else xx
= -x
- asncs
.x
[n
];
108 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+xx
*(asncs
.x
[n
+5]
109 +xx
*asncs
.x
[n
+6]))))+asncs
.x
[n
+7];
111 res
=asncs
.x
[n
+8] +t
;
112 cor
= (asncs
.x
[n
+8]-res
)+t
;
113 if (res
== res
+1.05*cor
) return (m
>0)?res
:-res
;
115 r
=asncs
.x
[n
+8]+xx
*asncs
.x
[n
+9];
116 t
=((asncs
.x
[n
+8]-r
)+xx
*asncs
.x
[n
+9])+(p
+xx
*asncs
.x
[n
+10]);
119 if (res
== res
+1.0005*cor
) return (m
>0)?res
:-res
;
124 z
=(w
[0]-ABS(x
))+w
[1];
125 if (z
>1.0e-27) return (m
>0)?min(res
,res1
):-min(res
,res1
);
126 else if (z
<-1.0e-27) return (m
>0)?max(res
,res1
):-max(res
,res1
);
129 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
133 } /* else if (k < 0x3fe00000) */
134 /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
136 if (k
< 0x3fe80000) {
137 n
= 1056+((k
&0x000fe000)>>11)*3;
138 if (m
>0) xx
= x
- asncs
.x
[n
];
139 else xx
= -x
- asncs
.x
[n
];
141 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+xx
*(asncs
.x
[n
+5]
142 +xx
*(asncs
.x
[n
+6]+xx
*asncs
.x
[n
+7])))))+asncs
.x
[n
+8];
144 res
=asncs
.x
[n
+9] +t
;
145 cor
= (asncs
.x
[n
+9]-res
)+t
;
146 if (res
== res
+1.01*cor
) return (m
>0)?res
:-res
;
148 r
=asncs
.x
[n
+9]+xx
*asncs
.x
[n
+10];
149 t
=((asncs
.x
[n
+9]-r
)+xx
*asncs
.x
[n
+10])+(p
+xx
*asncs
.x
[n
+11]);
152 if (res
== res
+1.0005*cor
) return (m
>0)?res
:-res
;
157 z
=(w
[0]-ABS(x
))+w
[1];
158 if (z
>1.0e-27) return (m
>0)?min(res
,res1
):-min(res
,res1
);
159 else if (z
<-1.0e-27) return (m
>0)?max(res
,res1
):-max(res
,res1
);
162 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
166 } /* else if (k < 0x3fe80000) */
167 /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
169 if (k
< 0x3fed8000) {
170 n
= 992+((k
&0x000fe000)>>13)*13;
171 if (m
>0) xx
= x
- asncs
.x
[n
];
172 else xx
= -x
- asncs
.x
[n
];
174 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+xx
*(asncs
.x
[n
+5]
175 +xx
*(asncs
.x
[n
+6]+xx
*(asncs
.x
[n
+7]+xx
*asncs
.x
[n
+8]))))))+asncs
.x
[n
+9];
177 res
=asncs
.x
[n
+10] +t
;
178 cor
= (asncs
.x
[n
+10]-res
)+t
;
179 if (res
== res
+1.01*cor
) return (m
>0)?res
:-res
;
181 r
=asncs
.x
[n
+10]+xx
*asncs
.x
[n
+11];
182 t
=((asncs
.x
[n
+10]-r
)+xx
*asncs
.x
[n
+11])+(p
+xx
*asncs
.x
[n
+12]);
185 if (res
== res
+1.0008*cor
) return (m
>0)?res
:-res
;
190 z
=((hp0
.x
-y
)-res
)+(hp1
.x
-z
);
192 z
=(w
[0]-ABS(x
))+w
[1];
193 if (z
>1.0e-27) return (m
>0)?min(res
,res1
):-min(res
,res1
);
194 else if (z
<-1.0e-27) return (m
>0)?max(res
,res1
):-max(res
,res1
);
197 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
201 } /* else if (k < 0x3fed8000) */
202 /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
204 if (k
< 0x3fee8000) {
205 n
= 884+((k
&0x000fe000)>>13)*14;
206 if (m
>0) xx
= x
- asncs
.x
[n
];
207 else xx
= -x
- asncs
.x
[n
];
209 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
210 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
211 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+
212 xx
*asncs
.x
[n
+9])))))))+asncs
.x
[n
+10];
214 res
=asncs
.x
[n
+11] +t
;
215 cor
= (asncs
.x
[n
+11]-res
)+t
;
216 if (res
== res
+1.01*cor
) return (m
>0)?res
:-res
;
218 r
=asncs
.x
[n
+11]+xx
*asncs
.x
[n
+12];
219 t
=((asncs
.x
[n
+11]-r
)+xx
*asncs
.x
[n
+12])+(p
+xx
*asncs
.x
[n
+13]);
222 if (res
== res
+1.0007*cor
) return (m
>0)?res
:-res
;
230 z
=(w
[0]-ABS(x
))+w
[1];
231 if (z
>1.0e-27) return (m
>0)?min(res
,res1
):-min(res
,res1
);
232 else if (z
<-1.0e-27) return (m
>0)?max(res
,res1
):-max(res
,res1
);
235 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
239 } /* else if (k < 0x3fee8000) */
241 /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
243 if (k
< 0x3fef0000) {
244 n
= 768+((k
&0x000fe000)>>13)*15;
245 if (m
>0) xx
= x
- asncs
.x
[n
];
246 else xx
= -x
- asncs
.x
[n
];
248 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
249 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
250 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+
251 xx
*(asncs
.x
[n
+9]+xx
*asncs
.x
[n
+10]))))))))+asncs
.x
[n
+11];
253 res
=asncs
.x
[n
+12] +t
;
254 cor
= (asncs
.x
[n
+12]-res
)+t
;
255 if (res
== res
+1.01*cor
) return (m
>0)?res
:-res
;
257 r
=asncs
.x
[n
+12]+xx
*asncs
.x
[n
+13];
258 t
=((asncs
.x
[n
+12]-r
)+xx
*asncs
.x
[n
+13])+(p
+xx
*asncs
.x
[n
+14]);
261 if (res
== res
+1.0007*cor
) return (m
>0)?res
:-res
;
269 z
=(w
[0]-ABS(x
))+w
[1];
270 if (z
>1.0e-27) return (m
>0)?min(res
,res1
):-min(res
,res1
);
271 else if (z
<-1.0e-27) return (m
>0)?max(res
,res1
):-max(res
,res1
);
274 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
278 } /* else if (k < 0x3fef0000) */
279 /*--------------------0.96875 <= |x| < 1 --------------------------------*/
282 z
= 0.5*((m
>0)?(1.0-x
):(1.0+x
));
285 t
=inroot
[(k
&0x001fffff)>>14]*powtwo
[511-(k
>>21)];
287 t
= t
*(rt0
+r
*(rt1
+r
*(rt2
+r
*rt3
)));
292 p
=(((((f6
*z
+f5
)*z
+f4
)*z
+f3
)*z
+f2
)*z
+f1
)*z
;
293 cor
= (hp1
.x
- 2.0*cc
)-2.0*(y
+cc
)*p
;
294 res1
= hp0
.x
- 2.0*y
;
296 if (res
== res
+1.003*((res1
-res
)+cor
)) return (m
>0)?res
:-res
;
302 cor
=((hp0
.x
-res1
)-2.0*w
[0])+(hp1
.x
-2.0*w
[1]);
304 cor
= (res1
-res
)+cor
;
305 if (res
==(res
+1.0000001*cor
)) return (m
>0)?res
:-res
;
309 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
312 } /* else if (k < 0x3ff00000) */
313 /*---------------------------- |x|>=1 -------------------------------*/
314 else if (k
==0x3ff00000 && u
.i
[LOW_HALF
]==0) return (m
>0)?hp0
.x
:-hp0
.x
;
316 if (k
>0x7ff00000 || (k
== 0x7ff00000 && u
.i
[LOW_HALF
] != 0)) return x
;
318 u
.i
[HIGH_HALF
]=0x7ff00000;
319 v
.i
[HIGH_HALF
]=0x7ff00000;
322 return u
.x
/v
.x
; /* NaN */
326 /*******************************************************************/
328 /* End of arcsine, below is arccosine */
330 /*******************************************************************/
332 double __ieee754_acos(double x
)
334 double x1
,x2
,xx
,s1
,s2
,res1
,p
,t
,res
,r
,cor
,cc
,y
,c
,z
,w
[2],eps
;
346 /*------------------- |x|<2.77556*10^-17 ----------------------*/
347 if (k
< 0x3c880000) return hp0
.x
;
349 /*----------------- 2.77556*10^-17 <= |x| < 2^-3 --------------*/
351 if (k
< 0x3fc00000) {
353 t
= (((((f6
*x2
+ f5
)*x2
+ f4
)*x2
+ f3
)*x2
+ f2
)*x2
+ f1
)*(x2
*x
);
355 cor
=(((hp0
.x
-r
)-x
)+hp1
.x
)-t
;
358 if (res
== res
+1.004*cor
) return res
;
366 s2
= ((((((c7
*xx
+ c6
)*xx
+ c5
)*xx
+ c4
)*xx
+ c3
)*xx
+ c2
)*xx
*xx
*x
+
367 ((a1
.x
+a2
.x
)*x2
*x2
+ 0.5*x1
*x
)*x2
) + a2
.x
*p
;
369 s2
= ((x
-res1
)+s1
)+s2
;
371 cor
=(((hp0
.x
-r
)-res1
)+hp1
.x
)-s2
;
374 if (res
== res
+1.00004*cor
) return res
;
378 cor
=((hp0
.x
-r
)-w
[0])+(hp1
.x
-w
[1]);
381 if (res
==(res
+1.00000001*cor
)) return res
;
384 return __cos32(x
,res
,res1
);
388 } /* else if (k < 0x3fc00000) */
389 /*---------------------- 0.125 <= |x| < 0.5 --------------------*/
391 if (k
< 0x3fe00000) {
392 if (k
<0x3fd00000) n
= 11*((k
&0x000fffff)>>15);
393 else n
= 11*((k
&0x000fffff)>>14)+352;
394 if (m
>0) xx
= x
- asncs
.x
[n
];
395 else xx
= -x
- asncs
.x
[n
];
397 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
398 xx
*(asncs
.x
[n
+5]+xx
*asncs
.x
[n
+6]))))+asncs
.x
[n
+7];
400 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+8]):(hp0
.x
+asncs
.x
[n
+8]);
401 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
403 if (res
== res
+1.02*((y
-res
)+t
)) return res
;
405 r
=asncs
.x
[n
+8]+xx
*asncs
.x
[n
+9];
406 t
=((asncs
.x
[n
+8]-r
)+xx
*asncs
.x
[n
+9])+(p
+xx
*asncs
.x
[n
+10]);
408 {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; }
410 {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); }
413 if (res
== (res
+1.0002*cor
)) return res
;
419 if (z
>1.0e-27) return max(res
,res1
);
420 else if (z
<-1.0e-27) return min(res
,res1
);
421 else return __cos32(x
,res
,res1
);
424 } /* else if (k < 0x3fe00000) */
426 /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
428 if (k
< 0x3fe80000) {
429 n
= 1056+((k
&0x000fe000)>>11)*3;
430 if (m
>0) {xx
= x
- asncs
.x
[n
]; eps
=1.04; }
431 else {xx
= -x
- asncs
.x
[n
]; eps
=1.02; }
433 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
434 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]+
435 xx
*asncs
.x
[n
+7])))))+asncs
.x
[n
+8];
437 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+9]):(hp0
.x
+asncs
.x
[n
+9]);
438 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
440 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
442 r
=asncs
.x
[n
+9]+xx
*asncs
.x
[n
+10];
443 t
=((asncs
.x
[n
+9]-r
)+xx
*asncs
.x
[n
+10])+(p
+xx
*asncs
.x
[n
+11]);
444 if (m
>0) {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; eps
=1.0004; }
445 else {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); eps
=1.0002; }
448 if (res
== (res
+eps
*cor
)) return res
;
454 if (z
>1.0e-27) return max(res
,res1
);
455 else if (z
<-1.0e-27) return min(res
,res1
);
456 else return __cos32(x
,res
,res1
);
459 } /* else if (k < 0x3fe80000) */
461 /*------------------------- 0.75 <= |x| < 0.921875 -------------*/
463 if (k
< 0x3fed8000) {
464 n
= 992+((k
&0x000fe000)>>13)*13;
465 if (m
>0) {xx
= x
- asncs
.x
[n
]; eps
= 1.04; }
466 else {xx
= -x
- asncs
.x
[n
]; eps
= 1.01; }
468 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
469 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]+xx
*(asncs
.x
[n
+7]+
470 xx
*asncs
.x
[n
+8]))))))+asncs
.x
[n
+9];
472 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+10]):(hp0
.x
+asncs
.x
[n
+10]);
473 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
475 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
477 r
=asncs
.x
[n
+10]+xx
*asncs
.x
[n
+11];
478 t
=((asncs
.x
[n
+10]-r
)+xx
*asncs
.x
[n
+11])+(p
+xx
*asncs
.x
[n
+12]);
479 if (m
>0) {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; eps
=1.0032; }
480 else {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); eps
=1.0008; }
483 if (res
== (res
+eps
*cor
)) return res
;
489 if (z
>1.0e-27) return max(res
,res1
);
490 else if (z
<-1.0e-27) return min(res
,res1
);
491 else return __cos32(x
,res
,res1
);
494 } /* else if (k < 0x3fed8000) */
496 /*-------------------0.921875 <= |x| < 0.953125 ------------------*/
498 if (k
< 0x3fee8000) {
499 n
= 884+((k
&0x000fe000)>>13)*14;
500 if (m
>0) {xx
= x
- asncs
.x
[n
]; eps
=1.04; }
501 else {xx
= -x
- asncs
.x
[n
]; eps
=1.005; }
503 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
504 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
505 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+
506 xx
*asncs
.x
[n
+9])))))))+asncs
.x
[n
+10];
508 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+11]):(hp0
.x
+asncs
.x
[n
+11]);
509 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
511 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
513 r
=asncs
.x
[n
+11]+xx
*asncs
.x
[n
+12];
514 t
=((asncs
.x
[n
+11]-r
)+xx
*asncs
.x
[n
+12])+(p
+xx
*asncs
.x
[n
+13]);
515 if (m
>0) {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; eps
=1.0030; }
516 else {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); eps
=1.0005; }
519 if (res
== (res
+eps
*cor
)) return res
;
525 if (z
>1.0e-27) return max(res
,res1
);
526 else if (z
<-1.0e-27) return min(res
,res1
);
527 else return __cos32(x
,res
,res1
);
530 } /* else if (k < 0x3fee8000) */
532 /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
534 if (k
< 0x3fef0000) {
535 n
= 768+((k
&0x000fe000)>>13)*15;
536 if (m
>0) {xx
= x
- asncs
.x
[n
]; eps
=1.04; }
537 else {xx
= -x
- asncs
.x
[n
]; eps
=1.005;}
539 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
540 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
541 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+xx
*(asncs
.x
[n
+9]+
542 xx
*asncs
.x
[n
+10]))))))))+asncs
.x
[n
+11];
544 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+12]):(hp0
.x
+asncs
.x
[n
+12]);
545 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
547 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
549 r
=asncs
.x
[n
+12]+xx
*asncs
.x
[n
+13];
550 t
=((asncs
.x
[n
+12]-r
)+xx
*asncs
.x
[n
+13])+(p
+xx
*asncs
.x
[n
+14]);
551 if (m
>0) {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; eps
=1.0030; }
552 else {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); eps
=1.0005; }
555 if (res
== (res
+eps
*cor
)) return res
;
561 if (z
>1.0e-27) return max(res
,res1
);
562 else if (z
<-1.0e-27) return min(res
,res1
);
563 else return __cos32(x
,res
,res1
);
566 } /* else if (k < 0x3fef0000) */
567 /*-----------------0.96875 <= |x| < 1 ---------------------------*/
571 z
= 0.5*((m
>0)?(1.0-x
):(1.0+x
));
574 t
=inroot
[(k
&0x001fffff)>>14]*powtwo
[511-(k
>>21)];
576 t
= t
*(rt0
+r
*(rt1
+r
*(rt2
+r
*rt3
)));
581 p
=(((((f6
*z
+f5
)*z
+f4
)*z
+f3
)*z
+f2
)*z
+f1
)*z
;
583 cor
= (hp1
.x
- cc
)-(y
+cc
)*p
;
586 if (res
== res
+1.002*((res1
-res
)+cor
)) return (res
+res
);
592 cor
=((hp0
.x
-res1
)-w
[0])+(hp1
.x
-w
[1]);
594 cor
= (res1
-res
)+cor
;
595 if (res
==(res
+1.000001*cor
)) return (res
+res
);
599 return __cos32(x
,res
,res1
);
606 if (res
== res
+1.03*((y
-res
)+cor
)) return (res
+res
);
613 if (res
==(res
+1.000001*cor
)) return (res
+res
);
617 return __cos32(x
,res
,res1
);
621 } /* else if (k < 0x3ff00000) */
623 /*---------------------------- |x|>=1 -----------------------*/
625 if (k
==0x3ff00000 && u
.i
[LOW_HALF
]==0) return (m
>0)?0:2.0*hp0
.x
;
627 if (k
>0x7ff00000 || (k
== 0x7ff00000 && u
.i
[LOW_HALF
] != 0)) return x
;
629 u
.i
[HIGH_HALF
]=0x7ff00000;
630 v
.i
[HIGH_HALF
]=0x7ff00000;