2 /* @(#)e_rem_pio2.c 1.4 95/01/18 */
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
7 * Developed at SunSoft, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
11 * ====================================================
13 * Optimized by Bruce D. Evans.
16 //#include <sys/cdefs.h>
17 //__FBSDID("$FreeBSD$");
19 /* __ieee754_rem_pio2(x,y)
21 * return the remainder of x rem pi/2 in y[0]+y[1]
22 * use __kernel_rem_pio2()
27 #include "math_private.h"
30 * invpio2: 53 bits of 2/pi
31 * pio2_1: first 33 bit of pi/2
32 * pio2_1t: pi/2 - pio2_1
33 * pio2_2: second 33 bit of pi/2
34 * pio2_2t: pi/2 - (pio2_1+pio2_2)
35 * pio2_3: third 33 bit of pi/2
36 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
40 zero
= 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
41 two24
= 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
42 invpio2
= 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
43 pio2_1
= 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
44 pio2_1t
= 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
45 pio2_2
= 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
46 pio2_2t
= 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
47 pio2_3
= 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
48 pio2_3t
= 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
50 #ifdef INLINE_REM_PIO2
54 __ieee754_rem_pio2(double x
, double *y
)
58 int32_t e0
,i
,j
,nx
,n
,ix
,hx
;
61 GET_HIGH_WORD(hx
,x
); /* high word of x */
63 #if 0 /* Must be handled in caller. */
64 if(ix
<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
65 {y
[0] = x
; y
[1] = 0; return 0;}
67 if (ix
<= 0x400f6a7a) { /* |x| ~<= 5pi/4 */
68 if ((ix
& 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */
69 goto medium
; /* cancellation -- use medium case */
70 if (ix
<= 0x4002d97c) { /* |x| ~<= 3pi/4 */
72 z
= x
- pio2_1
; /* one round good to 85 bits */
74 y
[1] = (z
-y
[0])-pio2_1t
;
79 y
[1] = (z
-y
[0])+pio2_1t
;
86 y
[1] = (z
-y
[0])-2*pio2_1t
;
91 y
[1] = (z
-y
[0])+2*pio2_1t
;
96 if (ix
<= 0x401c463b) { /* |x| ~<= 9pi/4 */
97 if (ix
<= 0x4015fdbc) { /* |x| ~<= 7pi/4 */
98 if (ix
== 0x4012d97c) /* |x| ~= 3pi/2 */
102 y
[0] = z
- 3*pio2_1t
;
103 y
[1] = (z
-y
[0])-3*pio2_1t
;
107 y
[0] = z
+ 3*pio2_1t
;
108 y
[1] = (z
-y
[0])+3*pio2_1t
;
112 if (ix
== 0x401921fb) /* |x| ~= 4pi/2 */
116 y
[0] = z
- 4*pio2_1t
;
117 y
[1] = (z
-y
[0])-4*pio2_1t
;
121 y
[0] = z
+ 4*pio2_1t
;
122 y
[1] = (z
-y
[0])+4*pio2_1t
;
127 if(ix
<0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */
129 fn
= rnint((double_t
)x
*invpio2
);
132 w
= fn
*pio2_1t
; /* 1st round good to 85 bit */
137 GET_HIGH_WORD(high
,y
[0]);
138 i
= j
-((high
>>20)&0x7ff);
139 if(i
>16) { /* 2nd iteration needed, good to 118 */
143 w
= fn
*pio2_2t
-((t
-r
)-w
);
145 GET_HIGH_WORD(high
,y
[0]);
146 i
= j
-((high
>>20)&0x7ff);
147 if(i
>49) { /* 3rd iteration need, 151 bits acc */
148 t
= r
; /* will cover all possible cases */
151 w
= fn
*pio2_3t
-((t
-r
)-w
);
160 * all other (large) arguments
162 if(ix
>=0x7ff00000) { /* x is inf or NaN */
163 y
[0]=y
[1]=x
-x
; return 0;
165 /* set z = scalbn(|x|,ilogb(x)-23) */
167 e0
= (ix
>>20)-1046; /* e0 = ilogb(z)-23; */
168 INSERT_WORDS(z
, ix
- ((int32_t)(e0
<<20)), low
);
170 tx
[i
] = (double)((int32_t)(z
));
175 while(tx
[nx
-1]==zero
) nx
--; /* skip zero term */
176 n
= __kernel_rem_pio2(tx
,ty
,e0
,nx
,1);
177 if(hx
<0) {y
[0] = -ty
[0]; y
[1] = -ty
[1]; return -n
;}
178 y
[0] = ty
[0]; y
[1] = ty
[1]; return n
;