RTAI API: base/math/e_rem_pio2.c Source File

00001 /* @(#)e_rem_pio2.c 5.1 93/09/24 */ 00002 /* 00003 * ==================================================== 00004 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 00005 * 00006 * Developed at SunPro, a Sun Microsystems, Inc. business. 00007 * Permission to use, copy, modify, and distribute this 00008 * software is freely granted, provided that this notice 00009 * is preserved. 00010 * ==================================================== 00011 */ 00012 00013 #if defined(LIBM_SCCS) && !defined(lint) 00014 static char rcsid[] = "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $"; 00015 #endif 00016 00017 /* __ieee754_rem_pio2(x,y) 00018 * 00019 * return the remainder of x rem pi/2 in y[0]+y[1] 00020 * use __kernel_rem_pio2() 00021 */ 00022 00023 #include "math.h" 00024 #include "mathP.h" 00025 00026 /* 00027 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 00028 */ 00029 #ifdef __STDC__ 00030 static const int32_t two_over_pi[] = { 00031 #else 00032 static int32_t two_over_pi[] = { 00033 #endif 00034 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 00035 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 00036 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 00037 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 00038 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 00039 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 00040 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 00041 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 00042 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 00043 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 00044 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 00045 }; 00046 00047 #ifdef __STDC__ 00048 static const int32_t npio2_hw[] = { 00049 #else 00050 static int32_t npio2_hw[] = { 00051 #endif 00052 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 00053 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 00054 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 00055 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C, 00056 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, 00057 0x404858EB, 0x404921FB, 00058 }; 00059 00060 /* 00061 * invpio2: 53 bits of 2/pi 00062 * pio2_1: first 33 bit of pi/2 00063 * pio2_1t: pi/2 - pio2_1 00064 * pio2_2: second 33 bit of pi/2 00065 * pio2_2t: pi/2 - (pio2_1+pio2_2) 00066 * pio2_3: third 33 bit of pi/2 00067 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) 00068 */ 00069 00070 #ifdef __STDC__ 00071 static const double 00072 #else 00073 static double 00074 #endif 00075 zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 00076 half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 00077 two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ 00078 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ 00079 pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ 00080 pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ 00081 pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ 00082 pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ 00083 pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ 00084 pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ 00085 00086 #ifdef __STDC__ 00087 int32_t __ieee754_rem_pio2(double x, double *y) 00088 #else 00089 int32_t __ieee754_rem_pio2(x,y) 00090 double x,y[]; 00091 #endif 00092 { 00093 double z,w,t,r,fn; 00094 double tx[3]; 00095 int32_t e0,i,j,nx,n,ix,hx; 00096 u_int32_t low; 00097 00098 /* XXX z is not properly initialized without this. Is this a 00099 * bug? --ds */ 00100 z = 0; 00101 00102 GET_HIGH_WORD(hx,x); /* high word of x */ 00103 ix = hx&0x7fffffff; 00104 if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ 00105 {y[0] = x; y[1] = 0; return 0;} 00106 if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ 00107 if(hx>0) { 00108 z = x - pio2_1; 00109 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 00110 y[0] = z - pio2_1t; 00111 y[1] = (z-y[0])-pio2_1t; 00112 } else { /* near pi/2, use 33+33+53 bit pi */ 00113 z -= pio2_2; 00114 y[0] = z - pio2_2t; 00115 y[1] = (z-y[0])-pio2_2t; 00116 } 00117 return 1; 00118 } else { /* negative x */ 00119 z = x + pio2_1; 00120 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 00121 y[0] = z + pio2_1t; 00122 y[1] = (z-y[0])+pio2_1t; 00123 } else { /* near pi/2, use 33+33+53 bit pi */ 00124 z += pio2_2; 00125 y[0] = z + pio2_2t; 00126 y[1] = (z-y[0])+pio2_2t; 00127 } 00128 return -1; 00129 } 00130 } 00131 if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 00132 t = fabs(x); 00133 n = (int32_t) (t*invpio2+half); 00134 fn = (double)n; 00135 r = t-fn*pio2_1; 00136 w = fn*pio2_1t; /* 1st round good to 85 bit */ 00137 if(n<32&&ix!=npio2_hw[n-1]) { 00138 y[0] = r-w; /* quick check no cancellation */ 00139 } else { 00140 u_int32_t high; 00141 j = ix>>20; 00142 y[0] = r-w; 00143 GET_HIGH_WORD(high,y[0]); 00144 i = j-((high>>20)&0x7ff); 00145 if(i>16) { /* 2nd iteration needed, good to 118 */ 00146 t = r; 00147 w = fn*pio2_2; 00148 r = t-w; 00149 w = fn*pio2_2t-((t-r)-w); 00150 y[0] = r-w; 00151 GET_HIGH_WORD(high,y[0]); 00152 i = j-((high>>20)&0x7ff); 00153 if(i>49) { /* 3rd iteration need, 151 bits acc */ 00154 t = r; /* will cover all possible cases */ 00155 w = fn*pio2_3; 00156 r = t-w; 00157 w = fn*pio2_3t-((t-r)-w); 00158 y[0] = r-w; 00159 } 00160 } 00161 } 00162 y[1] = (r-y[0])-w; 00163 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 00164 else return n; 00165 } 00166 /* 00167 * all other (large) arguments 00168 */ 00169 if(ix>=0x7ff00000) { /* x is inf or NaN */ 00170 y[0]=y[1]=x-x; return 0; 00171 } 00172 /* set z = scalbn(|x|,ilogb(x)-23) */ 00173 GET_LOW_WORD(low,x); 00174 SET_LOW_WORD(z,low); 00175 e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ 00176 SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20))); 00177 for(i=0;i<2;i++) { 00178 tx[i] = (double)((int32_t)(z)); 00179 z = (z-tx[i])*two24; 00180 } 00181 tx[2] = z; 00182 nx = 3; 00183 while(tx[nx-1]==zero) nx--; /* skip zero term */ 00184 n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi); 00185 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 00186 return n; 00187 }