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