2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001 Free Software Foundation
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, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 /******************************************************************/
21 /* MODULE_NAME:uasncs.c */
23 /* FUNCTIONS: uasin */
25 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */
26 /* doasin.c sincos32.c dosincos.c mpa.c */
27 /* sincos.tbl asincos.tbl powtwo.tbl root.tbl */
29 /* Ultimate asin/acos routines. Given an IEEE double machine */
30 /* number x, compute the correctly rounded value of */
31 /* arcsin(x)or arccos(x) according to the function called. */
32 /* Assumption: Machine arithmetic operations are performed in */
33 /* round to nearest mode of IEEE 754 standard. */
35 /******************************************************************/
38 #include "asincos.tbl"
43 #include "math_private.h"
45 void __doasin(double x
, double dx
, double w
[]);
46 void __dubsin(double x
, double dx
, double v
[]);
47 void __dubcos(double x
, double dx
, double v
[]);
48 void __docos(double x
, double dx
, double v
[]);
49 double __sin32(double x
, double res
, double res1
);
50 double __cos32(double x
, double res
, double res1
);
52 /***************************************************************************/
53 /* An ultimate asin routine. Given an IEEE double machine number x */
54 /* it computes the correctly rounded (to nearest) value of arcsin(x) */
55 /***************************************************************************/
56 double __ieee754_asin(double x
){
57 double x1
,x2
,xx
,s1
,s2
,res1
,p
,t
,res
,r
,cor
,cc
,y
,c
,z
,w
[2];
66 k
= 0x7fffffff&m
; /* no sign */
68 if (k
< 0x3e500000) return x
; /* for x->0 => sin(x)=x */
69 /*----------------------2^-26 <= |x| < 2^ -3 -----------------*/
73 t
= (((((f6
*x2
+ f5
)*x2
+ f4
)*x2
+ f3
)*x2
+ f2
)*x2
+ f1
)*(x2
*x
);
74 res
= x
+t
; /* res=arcsin(x) according to Taylor series */
76 if (res
== res
+1.025*cor
) return res
;
84 s2
= ((((((c7
*xx
+ c6
)*xx
+ c5
)*xx
+ c4
)*xx
+ c3
)*xx
+ c2
)*xx
*xx
*x
+
85 ((a1
.x
+a2
.x
)*x2
*x2
+ 0.5*x1
*x
)*x2
) + a2
.x
*p
;
87 s2
= ((x
-res1
)+s1
)+s2
;
90 if (res
== res
+1.00014*cor
) return res
;
93 if (w
[0]==(w
[0]+1.00000001*w
[1])) return w
[0];
97 res1
=ABS(w
[0]+1.1*w
[1]);
98 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
103 /*---------------------0.125 <= |x| < 0.5 -----------------------------*/
104 else if (k
< 0x3fe00000) {
105 if (k
<0x3fd00000) n
= 11*((k
&0x000fffff)>>15);
106 else n
= 11*((k
&0x000fffff)>>14)+352;
107 if (m
>0) xx
= x
- asncs
.x
[n
];
108 else xx
= -x
- asncs
.x
[n
];
110 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+xx
*(asncs
.x
[n
+5]
111 +xx
*asncs
.x
[n
+6]))))+asncs
.x
[n
+7];
113 res
=asncs
.x
[n
+8] +t
;
114 cor
= (asncs
.x
[n
+8]-res
)+t
;
115 if (res
== res
+1.05*cor
) return (m
>0)?res
:-res
;
117 r
=asncs
.x
[n
+8]+xx
*asncs
.x
[n
+9];
118 t
=((asncs
.x
[n
+8]-r
)+xx
*asncs
.x
[n
+9])+(p
+xx
*asncs
.x
[n
+10]);
121 if (res
== res
+1.0005*cor
) return (m
>0)?res
:-res
;
126 z
=(w
[0]-ABS(x
))+w
[1];
127 if (z
>1.0e-27) return (m
>0)?min(res
,res1
):-min(res
,res1
);
128 else if (z
<-1.0e-27) return (m
>0)?max(res
,res1
):-max(res
,res1
);
131 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
135 } /* else if (k < 0x3fe00000) */
136 /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
138 if (k
< 0x3fe80000) {
139 n
= 1056+((k
&0x000fe000)>>11)*3;
140 if (m
>0) xx
= x
- asncs
.x
[n
];
141 else xx
= -x
- asncs
.x
[n
];
143 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+xx
*(asncs
.x
[n
+5]
144 +xx
*(asncs
.x
[n
+6]+xx
*asncs
.x
[n
+7])))))+asncs
.x
[n
+8];
146 res
=asncs
.x
[n
+9] +t
;
147 cor
= (asncs
.x
[n
+9]-res
)+t
;
148 if (res
== res
+1.01*cor
) return (m
>0)?res
:-res
;
150 r
=asncs
.x
[n
+9]+xx
*asncs
.x
[n
+10];
151 t
=((asncs
.x
[n
+9]-r
)+xx
*asncs
.x
[n
+10])+(p
+xx
*asncs
.x
[n
+11]);
154 if (res
== res
+1.0005*cor
) return (m
>0)?res
:-res
;
159 z
=(w
[0]-ABS(x
))+w
[1];
160 if (z
>1.0e-27) return (m
>0)?min(res
,res1
):-min(res
,res1
);
161 else if (z
<-1.0e-27) return (m
>0)?max(res
,res1
):-max(res
,res1
);
164 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
168 } /* else if (k < 0x3fe80000) */
169 /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
171 if (k
< 0x3fed8000) {
172 n
= 992+((k
&0x000fe000)>>13)*13;
173 if (m
>0) xx
= x
- asncs
.x
[n
];
174 else xx
= -x
- asncs
.x
[n
];
176 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+xx
*(asncs
.x
[n
+5]
177 +xx
*(asncs
.x
[n
+6]+xx
*(asncs
.x
[n
+7]+xx
*asncs
.x
[n
+8]))))))+asncs
.x
[n
+9];
179 res
=asncs
.x
[n
+10] +t
;
180 cor
= (asncs
.x
[n
+10]-res
)+t
;
181 if (res
== res
+1.01*cor
) return (m
>0)?res
:-res
;
183 r
=asncs
.x
[n
+10]+xx
*asncs
.x
[n
+11];
184 t
=((asncs
.x
[n
+10]-r
)+xx
*asncs
.x
[n
+11])+(p
+xx
*asncs
.x
[n
+12]);
187 if (res
== res
+1.0008*cor
) return (m
>0)?res
:-res
;
192 z
=((hp0
.x
-y
)-res
)+(hp1
.x
-z
);
194 z
=(w
[0]-ABS(x
))+w
[1];
195 if (z
>1.0e-27) return (m
>0)?min(res
,res1
):-min(res
,res1
);
196 else if (z
<-1.0e-27) return (m
>0)?max(res
,res1
):-max(res
,res1
);
199 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
203 } /* else if (k < 0x3fed8000) */
204 /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
206 if (k
< 0x3fee8000) {
207 n
= 884+((k
&0x000fe000)>>13)*14;
208 if (m
>0) xx
= x
- asncs
.x
[n
];
209 else xx
= -x
- asncs
.x
[n
];
211 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
212 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
213 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+
214 xx
*asncs
.x
[n
+9])))))))+asncs
.x
[n
+10];
216 res
=asncs
.x
[n
+11] +t
;
217 cor
= (asncs
.x
[n
+11]-res
)+t
;
218 if (res
== res
+1.01*cor
) return (m
>0)?res
:-res
;
220 r
=asncs
.x
[n
+11]+xx
*asncs
.x
[n
+12];
221 t
=((asncs
.x
[n
+11]-r
)+xx
*asncs
.x
[n
+12])+(p
+xx
*asncs
.x
[n
+13]);
224 if (res
== res
+1.0007*cor
) return (m
>0)?res
:-res
;
232 z
=(w
[0]-ABS(x
))+w
[1];
233 if (z
>1.0e-27) return (m
>0)?min(res
,res1
):-min(res
,res1
);
234 else if (z
<-1.0e-27) return (m
>0)?max(res
,res1
):-max(res
,res1
);
237 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
241 } /* else if (k < 0x3fee8000) */
243 /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
245 if (k
< 0x3fef0000) {
246 n
= 768+((k
&0x000fe000)>>13)*15;
247 if (m
>0) xx
= x
- asncs
.x
[n
];
248 else xx
= -x
- asncs
.x
[n
];
250 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
251 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
252 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+
253 xx
*(asncs
.x
[n
+9]+xx
*asncs
.x
[n
+10]))))))))+asncs
.x
[n
+11];
255 res
=asncs
.x
[n
+12] +t
;
256 cor
= (asncs
.x
[n
+12]-res
)+t
;
257 if (res
== res
+1.01*cor
) return (m
>0)?res
:-res
;
259 r
=asncs
.x
[n
+12]+xx
*asncs
.x
[n
+13];
260 t
=((asncs
.x
[n
+12]-r
)+xx
*asncs
.x
[n
+13])+(p
+xx
*asncs
.x
[n
+14]);
263 if (res
== res
+1.0007*cor
) return (m
>0)?res
:-res
;
271 z
=(w
[0]-ABS(x
))+w
[1];
272 if (z
>1.0e-27) return (m
>0)?min(res
,res1
):-min(res
,res1
);
273 else if (z
<-1.0e-27) return (m
>0)?max(res
,res1
):-max(res
,res1
);
276 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
280 } /* else if (k < 0x3fef0000) */
281 /*--------------------0.96875 <= |x| < 1 --------------------------------*/
284 z
= 0.5*((m
>0)?(1.0-x
):(1.0+x
));
287 t
=inroot
[(k
&0x001fffff)>>14]*powtwo
[511-(k
>>21)];
289 t
= t
*(rt0
+r
*(rt1
+r
*(rt2
+r
*rt3
)));
294 p
=(((((f6
*z
+f5
)*z
+f4
)*z
+f3
)*z
+f2
)*z
+f1
)*z
;
295 cor
= (hp1
.x
- 2.0*cc
)-2.0*(y
+cc
)*p
;
296 res1
= hp0
.x
- 2.0*y
;
298 if (res
== res
+1.003*((res1
-res
)+cor
)) return (m
>0)?res
:-res
;
304 cor
=((hp0
.x
-res1
)-2.0*w
[0])+(hp1
.x
-2.0*w
[1]);
306 cor
= (res1
-res
)+cor
;
307 if (res
==(res
+1.0000001*cor
)) return (m
>0)?res
:-res
;
311 return (m
>0)?__sin32(y
,res
,res1
):-__sin32(y
,res
,res1
);
314 } /* else if (k < 0x3ff00000) */
315 /*---------------------------- |x|>=1 -------------------------------*/
316 else if (k
==0x3ff00000 && u
.i
[LOW_HALF
]==0) return (m
>0)?hp0
.x
:-hp0
.x
;
318 if (k
>0x7ff00000 || (k
== 0x7ff00000 && u
.i
[LOW_HALF
] != 0)) return x
;
320 u
.i
[HIGH_HALF
]=0x7ff00000;
321 v
.i
[HIGH_HALF
]=0x7ff00000;
324 return u
.x
/v
.x
; /* NaN */
328 /*******************************************************************/
330 /* End of arcsine, below is arccosine */
332 /*******************************************************************/
334 double __ieee754_acos(double x
)
336 double x1
,x2
,xx
,s1
,s2
,res1
,p
,t
,res
,r
,cor
,cc
,y
,c
,z
,w
[2],eps
;
348 /*------------------- |x|<2.77556*10^-17 ----------------------*/
349 if (k
< 0x3c880000) return hp0
.x
;
351 /*----------------- 2.77556*10^-17 <= |x| < 2^-3 --------------*/
353 if (k
< 0x3fc00000) {
355 t
= (((((f6
*x2
+ f5
)*x2
+ f4
)*x2
+ f3
)*x2
+ f2
)*x2
+ f1
)*(x2
*x
);
357 cor
=(((hp0
.x
-r
)-x
)+hp1
.x
)-t
;
360 if (res
== res
+1.004*cor
) return res
;
368 s2
= ((((((c7
*xx
+ c6
)*xx
+ c5
)*xx
+ c4
)*xx
+ c3
)*xx
+ c2
)*xx
*xx
*x
+
369 ((a1
.x
+a2
.x
)*x2
*x2
+ 0.5*x1
*x
)*x2
) + a2
.x
*p
;
371 s2
= ((x
-res1
)+s1
)+s2
;
373 cor
=(((hp0
.x
-r
)-res1
)+hp1
.x
)-s2
;
376 if (res
== res
+1.00004*cor
) return res
;
380 cor
=((hp0
.x
-r
)-w
[0])+(hp1
.x
-w
[1]);
383 if (res
==(res
+1.00000001*cor
)) return res
;
386 return __cos32(x
,res
,res1
);
390 } /* else if (k < 0x3fc00000) */
391 /*---------------------- 0.125 <= |x| < 0.5 --------------------*/
393 if (k
< 0x3fe00000) {
394 if (k
<0x3fd00000) n
= 11*((k
&0x000fffff)>>15);
395 else n
= 11*((k
&0x000fffff)>>14)+352;
396 if (m
>0) xx
= x
- asncs
.x
[n
];
397 else xx
= -x
- asncs
.x
[n
];
399 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
400 xx
*(asncs
.x
[n
+5]+xx
*asncs
.x
[n
+6]))))+asncs
.x
[n
+7];
402 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+8]):(hp0
.x
+asncs
.x
[n
+8]);
403 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
405 if (res
== res
+1.02*((y
-res
)+t
)) return res
;
407 r
=asncs
.x
[n
+8]+xx
*asncs
.x
[n
+9];
408 t
=((asncs
.x
[n
+8]-r
)+xx
*asncs
.x
[n
+9])+(p
+xx
*asncs
.x
[n
+10]);
410 {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; }
412 {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); }
415 if (res
== (res
+1.0002*cor
)) return res
;
421 if (z
>1.0e-27) return max(res
,res1
);
422 else if (z
<-1.0e-27) return min(res
,res1
);
423 else return __cos32(x
,res
,res1
);
426 } /* else if (k < 0x3fe00000) */
428 /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
430 if (k
< 0x3fe80000) {
431 n
= 1056+((k
&0x000fe000)>>11)*3;
432 if (m
>0) {xx
= x
- asncs
.x
[n
]; eps
=1.04; }
433 else {xx
= -x
- asncs
.x
[n
]; eps
=1.02; }
435 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
436 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]+
437 xx
*asncs
.x
[n
+7])))))+asncs
.x
[n
+8];
439 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+9]):(hp0
.x
+asncs
.x
[n
+9]);
440 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
442 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
444 r
=asncs
.x
[n
+9]+xx
*asncs
.x
[n
+10];
445 t
=((asncs
.x
[n
+9]-r
)+xx
*asncs
.x
[n
+10])+(p
+xx
*asncs
.x
[n
+11]);
446 if (m
>0) {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; eps
=1.0004; }
447 else {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); eps
=1.0002; }
450 if (res
== (res
+eps
*cor
)) return res
;
456 if (z
>1.0e-27) return max(res
,res1
);
457 else if (z
<-1.0e-27) return min(res
,res1
);
458 else return __cos32(x
,res
,res1
);
461 } /* else if (k < 0x3fe80000) */
463 /*------------------------- 0.75 <= |x| < 0.921875 -------------*/
465 if (k
< 0x3fed8000) {
466 n
= 992+((k
&0x000fe000)>>13)*13;
467 if (m
>0) {xx
= x
- asncs
.x
[n
]; eps
= 1.04; }
468 else {xx
= -x
- asncs
.x
[n
]; eps
= 1.01; }
470 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
471 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]+xx
*(asncs
.x
[n
+7]+
472 xx
*asncs
.x
[n
+8]))))))+asncs
.x
[n
+9];
474 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+10]):(hp0
.x
+asncs
.x
[n
+10]);
475 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
477 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
479 r
=asncs
.x
[n
+10]+xx
*asncs
.x
[n
+11];
480 t
=((asncs
.x
[n
+10]-r
)+xx
*asncs
.x
[n
+11])+(p
+xx
*asncs
.x
[n
+12]);
481 if (m
>0) {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; eps
=1.0032; }
482 else {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); eps
=1.0008; }
485 if (res
== (res
+eps
*cor
)) return res
;
491 if (z
>1.0e-27) return max(res
,res1
);
492 else if (z
<-1.0e-27) return min(res
,res1
);
493 else return __cos32(x
,res
,res1
);
496 } /* else if (k < 0x3fed8000) */
498 /*-------------------0.921875 <= |x| < 0.953125 ------------------*/
500 if (k
< 0x3fee8000) {
501 n
= 884+((k
&0x000fe000)>>13)*14;
502 if (m
>0) {xx
= x
- asncs
.x
[n
]; eps
=1.04; }
503 else {xx
= -x
- asncs
.x
[n
]; eps
=1.005; }
505 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
506 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
507 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+
508 xx
*asncs
.x
[n
+9])))))))+asncs
.x
[n
+10];
510 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+11]):(hp0
.x
+asncs
.x
[n
+11]);
511 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
513 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
515 r
=asncs
.x
[n
+11]+xx
*asncs
.x
[n
+12];
516 t
=((asncs
.x
[n
+11]-r
)+xx
*asncs
.x
[n
+12])+(p
+xx
*asncs
.x
[n
+13]);
517 if (m
>0) {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; eps
=1.0030; }
518 else {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); eps
=1.0005; }
521 if (res
== (res
+eps
*cor
)) return res
;
527 if (z
>1.0e-27) return max(res
,res1
);
528 else if (z
<-1.0e-27) return min(res
,res1
);
529 else return __cos32(x
,res
,res1
);
532 } /* else if (k < 0x3fee8000) */
534 /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
536 if (k
< 0x3fef0000) {
537 n
= 768+((k
&0x000fe000)>>13)*15;
538 if (m
>0) {xx
= x
- asncs
.x
[n
]; eps
=1.04; }
539 else {xx
= -x
- asncs
.x
[n
]; eps
=1.005;}
541 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
542 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
543 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+xx
*(asncs
.x
[n
+9]+
544 xx
*asncs
.x
[n
+10]))))))))+asncs
.x
[n
+11];
546 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+12]):(hp0
.x
+asncs
.x
[n
+12]);
547 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
549 if (res
== res
+eps
*((y
-res
)+t
)) return res
;
551 r
=asncs
.x
[n
+12]+xx
*asncs
.x
[n
+13];
552 t
=((asncs
.x
[n
+12]-r
)+xx
*asncs
.x
[n
+13])+(p
+xx
*asncs
.x
[n
+14]);
553 if (m
>0) {p
= hp0
.x
-r
; t
= (((hp0
.x
-p
)-r
)-t
)+hp1
.x
; eps
=1.0030; }
554 else {p
= hp0
.x
+r
; t
= ((hp0
.x
-p
)+r
)+(hp1
.x
+t
); eps
=1.0005; }
557 if (res
== (res
+eps
*cor
)) return res
;
563 if (z
>1.0e-27) return max(res
,res1
);
564 else if (z
<-1.0e-27) return min(res
,res1
);
565 else return __cos32(x
,res
,res1
);
568 } /* else if (k < 0x3fef0000) */
569 /*-----------------0.96875 <= |x| < 1 ---------------------------*/
573 z
= 0.5*((m
>0)?(1.0-x
):(1.0+x
));
576 t
=inroot
[(k
&0x001fffff)>>14]*powtwo
[511-(k
>>21)];
578 t
= t
*(rt0
+r
*(rt1
+r
*(rt2
+r
*rt3
)));
583 p
=(((((f6
*z
+f5
)*z
+f4
)*z
+f3
)*z
+f2
)*z
+f1
)*z
;
585 cor
= (hp1
.x
- cc
)-(y
+cc
)*p
;
588 if (res
== res
+1.002*((res1
-res
)+cor
)) return (res
+res
);
594 cor
=((hp0
.x
-res1
)-w
[0])+(hp1
.x
-w
[1]);
596 cor
= (res1
-res
)+cor
;
597 if (res
==(res
+1.000001*cor
)) return (res
+res
);
601 return __cos32(x
,res
,res1
);
608 if (res
== res
+1.03*((y
-res
)+cor
)) return (res
+res
);
615 if (res
==(res
+1.000001*cor
)) return (res
+res
);
619 return __cos32(x
,res
,res1
);
623 } /* else if (k < 0x3ff00000) */
625 /*---------------------------- |x|>=1 -----------------------*/
627 if (k
==0x3ff00000 && u
.i
[LOW_HALF
]==0) return (m
>0)?0:2.0*hp0
.x
;
629 if (k
>0x7ff00000 || (k
== 0x7ff00000 && u
.i
[LOW_HALF
] != 0)) return x
;
631 u
.i
[HIGH_HALF
]=0x7ff00000;
632 v
.i
[HIGH_HALF
]=0x7ff00000;