• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /* ef_rem_pio2.c -- float version of e_rem_pio2.c
2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3  */
4 
5 /*
6  * ====================================================
7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8  *
9  * Developed at SunPro, a Sun Microsystems, Inc. business.
10  * Permission to use, copy, modify, and distribute this
11  * software is freely granted, provided that this notice
12  * is preserved.
13  * ====================================================
14  *
15  */
16 
17 /* __ieee754_rem_pio2f(x,y)
18  *
19  * return the remainder of x rem pi/2 in y[0]+y[1]
20  * use __kernel_rem_pio2f()
21  */
22 
23 #include "fdlibm.h"
24 
25 /*
26  * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
27  */
28 #ifdef __STDC__
29 static const __int32_t two_over_pi[] = {
30 #else
31 static __int32_t two_over_pi[] = {
32 #endif
33 0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
34 0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
35 0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
36 0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
37 0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
38 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
39 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
40 0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
41 0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
42 0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
43 0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
44 0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
45 0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
46 0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
47 0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
48 0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
49 0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
50 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
51 0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
52 0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
53 0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
54 0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
55 };
56 
57 /* This array is like the one in e_rem_pio2.c, but the numbers are
58    single precision and the last 8 bits are forced to 0.  */
59 #ifdef __STDC__
60 static const __int32_t npio2_hw[] = {
61 #else
62 static __int32_t npio2_hw[] = {
63 #endif
64 0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
65 0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
66 0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
67 0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
68 0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
69 0x4242c700, 0x42490f00
70 };
71 
72 /*
73  * invpio2:  24 bits of 2/pi
74  * pio2_1:   first  17 bit of pi/2
75  * pio2_1t:  pi/2 - pio2_1
76  * pio2_2:   second 17 bit of pi/2
77  * pio2_2t:  pi/2 - (pio2_1+pio2_2)
78  * pio2_3:   third  17 bit of pi/2
79  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
80  */
81 
82 #ifdef __STDC__
83 static const float
84 #else
85 static float
86 #endif
87 zero =  0.0000000000e+00, /* 0x00000000 */
88 half =  5.0000000000e-01, /* 0x3f000000 */
89 two8 =  2.5600000000e+02, /* 0x43800000 */
90 invpio2 =  6.3661980629e-01, /* 0x3f22f984 */
91 pio2_1  =  1.5707855225e+00, /* 0x3fc90f80 */
92 pio2_1t =  1.0804334124e-05, /* 0x37354443 */
93 pio2_2  =  1.0804273188e-05, /* 0x37354400 */
94 pio2_2t =  6.0770999344e-11, /* 0x2e85a308 */
95 pio2_3  =  6.0770943833e-11, /* 0x2e85a300 */
96 pio2_3t =  6.1232342629e-17; /* 0x248d3132 */
97 
98 #ifdef __STDC__
__ieee754_rem_pio2f(float x,float * y)99 	__int32_t __ieee754_rem_pio2f(float x, float *y)
100 #else
101 	__int32_t __ieee754_rem_pio2f(x,y)
102 	float x,y[];
103 #endif
104 {
105 	float z,w,t,r,fn;
106 	float tx[3];
107 	__int32_t i,j,n,ix,hx;
108 	int e0,nx;
109 
110 	GET_FLOAT_WORD(hx,x);
111 	ix = hx&0x7fffffff;
112 	if(ix<=0x3f490fd8)   /* |x| ~<= pi/4 , no need for reduction */
113 	    {y[0] = x; y[1] = 0; return 0;}
114 	if(ix<0x4016cbe4) {  /* |x| < 3pi/4, special case with n=+-1 */
115 	    if(hx>0) {
116 		z = x - pio2_1;
117 		if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
118 		    y[0] = z - pio2_1t;
119 		    y[1] = (z-y[0])-pio2_1t;
120 		} else {		/* near pi/2, use 24+24+24 bit pi */
121 		    z -= pio2_2;
122 		    y[0] = z - pio2_2t;
123 		    y[1] = (z-y[0])-pio2_2t;
124 		}
125 		return 1;
126 	    } else {	/* negative x */
127 		z = x + pio2_1;
128 		if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
129 		    y[0] = z + pio2_1t;
130 		    y[1] = (z-y[0])+pio2_1t;
131 		} else {		/* near pi/2, use 24+24+24 bit pi */
132 		    z += pio2_2;
133 		    y[0] = z + pio2_2t;
134 		    y[1] = (z-y[0])+pio2_2t;
135 		}
136 		return -1;
137 	    }
138 	}
139 	if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
140 	    t  = fabsf(x);
141 	    n  = (__int32_t) (t*invpio2+half);
142 	    fn = (float)n;
143 	    r  = t-fn*pio2_1;
144 	    w  = fn*pio2_1t;	/* 1st round good to 40 bit */
145 	    if(n<32&&(ix&0xffffff00)!=npio2_hw[n-1]) {
146 		y[0] = r-w;	/* quick check no cancellation */
147 	    } else {
148 	        __uint32_t high;
149 	        j  = ix>>23;
150 	        y[0] = r-w;
151 		GET_FLOAT_WORD(high,y[0]);
152 	        i = j-((high>>23)&0xff);
153 	        if(i>8) {  /* 2nd iteration needed, good to 57 */
154 		    t  = r;
155 		    w  = fn*pio2_2;
156 		    r  = t-w;
157 		    w  = fn*pio2_2t-((t-r)-w);
158 		    y[0] = r-w;
159 		    GET_FLOAT_WORD(high,y[0]);
160 		    i = j-((high>>23)&0xff);
161 		    if(i>25)  {	/* 3rd iteration need, 74 bits acc */
162 		    	t  = r;	/* will cover all possible cases */
163 		    	w  = fn*pio2_3;
164 		    	r  = t-w;
165 		    	w  = fn*pio2_3t-((t-r)-w);
166 		    	y[0] = r-w;
167 		    }
168 		}
169 	    }
170 	    y[1] = (r-y[0])-w;
171 	    if(hx<0) 	{y[0] = -y[0]; y[1] = -y[1]; return -n;}
172 	    else	 return n;
173 	}
174     /*
175      * all other (large) arguments
176      */
177 	if(!FLT_UWORD_IS_FINITE(ix)) {
178 	    y[0]=y[1]=x-x; return 0;
179 	}
180     /* set z = scalbn(|x|,ilogb(x)-7) */
181 	e0 	= (int)((ix>>23)-134);	/* e0 = ilogb(z)-7; */
182 	SET_FLOAT_WORD(z, ix - ((__int32_t)e0<<23));
183 	for(i=0;i<2;i++) {
184 		tx[i] = (float)((__int32_t)(z));
185 		z     = (z-tx[i])*two8;
186 	}
187 	tx[2] = z;
188 	nx = 3;
189 	while(tx[nx-1]==zero) nx--;	/* skip zero term */
190 	n  =  __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi);
191 	if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
192 	return n;
193 }
194