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