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 }