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[] = "$Id: e_rem_pio2f.c,v 1.2 1994/08/18 23:05:58 jtc Exp $"; 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<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */ 118 t = fabsf(x); 119 n = (int32_t) (t*invpio2+half); 120 fn = (float)n; 121 r = t-fn*pio2_1; 122 w = fn*pio2_1t; /* 1st round good to 40 bit */ 123 if(n<32&&(ix&0xffffff00)!=npio2_hw[n-1]) { 124 y[0] = r-w; /* quick check no cancellation */ 125 } else { 126 u_int32_t high; 127 j = ix>>23; 128 y[0] = r-w; 129 GET_FLOAT_WORD(high,y[0]); 130 i = j-((high>>23)&0xff); 131 if(i>8) { /* 2nd iteration needed, good to 57 */ 132 t = r; 133 w = fn*pio2_2; 134 r = t-w; 135 w = fn*pio2_2t-((t-r)-w); 136 y[0] = r-w; 137 GET_FLOAT_WORD(high,y[0]); 138 i = j-((high>>23)&0xff); 139 if(i>25) { /* 3rd iteration need, 74 bits acc */ 140 t = r; /* will cover all possible cases */ 141 w = fn*pio2_3; 142 r = t-w; 143 w = fn*pio2_3t-((t-r)-w); 144 y[0] = r-w; 145 } 146 } 147 } 148 y[1] = (r-y[0])-w; 149 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 150 else return n; 151 } 152 /* 153 * all other (large) arguments 154 */ 155 if(ix>=0x7f800000) { /* x is inf or NaN */ 156 y[0]=y[1]=x-x; return 0; 157 } 158 /* set z = scalbn(|x|,ilogb(x)-7) */ 159 e0 = (ix>>23)-134; /* e0 = ilogb(z)-7; */ 160 SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23))); 161 for(i=0;i<2;i++) { 162 tx[i] = (float)((int32_t)(z)); 163 z = (z-tx[i])*two8; 164 } 165 tx[2] = z; 166 nx = 3; 167 while(tx[nx-1]==zero) nx--; /* skip zero term */ 168 n = __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi); 169 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 170 return n; 171 } 172