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__ 29static const __int32_t two_over_pi[] = { 30#else 31static __int32_t two_over_pi[] = { 32#endif 330xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC, 340x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62, 350x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63, 360xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A, 370x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09, 380xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29, 390xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44, 400x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41, 410x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C, 420x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8, 430x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11, 440x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF, 450x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E, 460x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5, 470xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92, 480xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08, 490x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0, 500xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3, 510x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85, 520x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80, 530x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA, 540x73, 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__ 60static const __int32_t npio2_hw[] = { 61#else 62static __int32_t npio2_hw[] = { 63#endif 640x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00, 650x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00, 660x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100, 670x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00, 680x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00, 690x4242c700, 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__ 83static const float 84#else 85static float 86#endif 87zero = 0.0000000000e+00, /* 0x00000000 */ 88half = 5.0000000000e-01, /* 0x3f000000 */ 89two8 = 2.5600000000e+02, /* 0x43800000 */ 90invpio2 = 6.3661980629e-01, /* 0x3f22f984 */ 91pio2_1 = 1.5707855225e+00, /* 0x3fc90f80 */ 92pio2_1t = 1.0804334124e-05, /* 0x37354443 */ 93pio2_2 = 1.0804273188e-05, /* 0x37354400 */ 94pio2_2t = 6.0770999344e-11, /* 0x2e85a308 */ 95pio2_3 = 6.0770943833e-11, /* 0x2e85a300 */ 96pio2_3t = 6.1232342629e-17; /* 0x248d3132 */ 97 98#ifdef __STDC__ 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