e_rem_pio2f.c revision 50476
1117397Skan/* e_rem_pio2f.c -- float version of e_rem_pio2.c 2117397Skan * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 3117397Skan */ 4117397Skan 5117397Skan/* 6117397Skan * ==================================================== 7117397Skan * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8117397Skan * 9117397Skan * Developed at SunPro, a Sun Microsystems, Inc. business. 10117397Skan * Permission to use, copy, modify, and distribute this 11117397Skan * software is freely granted, provided that this notice 12117397Skan * is preserved. 13117397Skan * ==================================================== 14117397Skan */ 15117397Skan 16117397Skan#ifndef lint 17117397Skanstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/e_rem_pio2f.c 50476 1999-08-28 00:22:10Z peter $"; 18169691Skan#endif 19117397Skan 20117397Skan/* __ieee754_rem_pio2f(x,y) 21117397Skan * 22117397Skan * return the remainder of x rem pi/2 in y[0]+y[1] 23117397Skan * use __kernel_rem_pio2f() 24117397Skan */ 25117397Skan 26117397Skan#include "math.h" 27117397Skan#include "math_private.h" 28117397Skan 29117397Skan/* 30117397Skan * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 31117397Skan */ 32117397Skan#ifdef __STDC__ 33117397Skanstatic const int32_t two_over_pi[] = { 34117397Skan#else 35117397Skanstatic int32_t two_over_pi[] = { 36117397Skan#endif 37117397Skan0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC, 38117397Skan0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62, 39169691Skan0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63, 40169691Skan0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A, 41169691Skan0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09, 42117397Skan0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29, 43117397Skan0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44, 44117397Skan0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41, 45117397Skan0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C, 46117397Skan0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8, 47117397Skan0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11, 48117397Skan0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF, 49117397Skan0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E, 50117397Skan0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5, 51117397Skan0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92, 52117397Skan0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08, 53117397Skan0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0, 54117397Skan0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3, 55117397Skan0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85, 56117397Skan0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80, 57117397Skan0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA, 58117397Skan0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B, 59117397Skan}; 60117397Skan 61117397Skan/* This array is like the one in e_rem_pio2.c, but the numbers are 62169691Skan single precision and the last 8 bits are forced to 0. */ 63169691Skan#ifdef __STDC__ 64static const int32_t npio2_hw[] = { 65#else 66static int32_t npio2_hw[] = { 67#endif 680x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00, 690x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00, 700x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100, 710x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00, 720x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00, 730x4242c700, 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__ 87static const float 88#else 89static float 90#endif 91zero = 0.0000000000e+00, /* 0x00000000 */ 92half = 5.0000000000e-01, /* 0x3f000000 */ 93two8 = 2.5600000000e+02, /* 0x43800000 */ 94invpio2 = 6.3661980629e-01, /* 0x3f22f984 */ 95pio2_1 = 1.5707855225e+00, /* 0x3fc90f80 */ 96pio2_1t = 1.0804334124e-05, /* 0x37354443 */ 97pio2_2 = 1.0804273188e-05, /* 0x37354400 */ 98pio2_2t = 6.0770999344e-11, /* 0x2e85a308 */ 99pio2_3 = 6.0770943833e-11, /* 0x2e85a300 */ 100pio2_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