Back to index

glibc  2.9
e_rem_pio2.c
Go to the documentation of this 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 "math_private.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        GET_HIGH_WORD(hx,x);        /* high word of x */
00099        ix = hx&0x7fffffff;
00100        if(ix<=0x3fe921fb)   /* |x| ~<= pi/4 , no need for reduction */
00101            {y[0] = x; y[1] = 0; return 0;}
00102        if(ix<0x4002d97c) {  /* |x| < 3pi/4, special case with n=+-1 */
00103            if(hx>0) { 
00104               z = x - pio2_1;
00105               if(ix!=0x3ff921fb) {        /* 33+53 bit pi is good enough */
00106                   y[0] = z - pio2_1t;
00107                   y[1] = (z-y[0])-pio2_1t;
00108               } else {             /* near pi/2, use 33+33+53 bit pi */
00109                   z -= pio2_2;
00110                   y[0] = z - pio2_2t;
00111                   y[1] = (z-y[0])-pio2_2t;
00112               }
00113               return 1;
00114            } else {  /* negative x */
00115               z = x + pio2_1;
00116               if(ix!=0x3ff921fb) {        /* 33+53 bit pi is good enough */
00117                   y[0] = z + pio2_1t;
00118                   y[1] = (z-y[0])+pio2_1t;
00119               } else {             /* near pi/2, use 33+33+53 bit pi */
00120                   z += pio2_2;
00121                   y[0] = z + pio2_2t;
00122                   y[1] = (z-y[0])+pio2_2t;
00123               }
00124               return -1;
00125            }
00126        }
00127        if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
00128            t  = fabs(x);
00129            n  = (int32_t) (t*invpio2+half);
00130            fn = (double)n;
00131            r  = t-fn*pio2_1;
00132            w  = fn*pio2_1t; /* 1st round good to 85 bit */
00133            if(n<32&&ix!=npio2_hw[n-1]) {  
00134               y[0] = r-w;   /* quick check no cancellation */
00135            } else {
00136                u_int32_t high;
00137                j  = ix>>20;
00138                y[0] = r-w; 
00139               GET_HIGH_WORD(high,y[0]);
00140                i = j-((high>>20)&0x7ff);
00141                if(i>16) {  /* 2nd iteration needed, good to 118 */
00142                   t  = r;
00143                   w  = fn*pio2_2;  
00144                   r  = t-w;
00145                   w  = fn*pio2_2t-((t-r)-w);     
00146                   y[0] = r-w;
00147                   GET_HIGH_WORD(high,y[0]);
00148                   i = j-((high>>20)&0x7ff);
00149                   if(i>49)  {      /* 3rd iteration need, 151 bits acc */
00150                      t  = r;       /* will cover all possible cases */
00151                      w  = fn*pio2_3;      
00152                      r  = t-w;
00153                      w  = fn*pio2_3t-((t-r)-w);  
00154                      y[0] = r-w;
00155                   }
00156               }
00157            }
00158            y[1] = (r-y[0])-w;
00159            if(hx<0)  {y[0] = -y[0]; y[1] = -y[1]; return -n;}
00160            else       return n;
00161        }
00162     /* 
00163      * all other (large) arguments
00164      */
00165        if(ix>=0x7ff00000) {        /* x is inf or NaN */
00166            y[0]=y[1]=x-x; return 0;
00167        }
00168     /* set z = scalbn(|x|,ilogb(x)-23) */
00169        GET_LOW_WORD(low,x);
00170        SET_LOW_WORD(z,low);
00171        e0     = (ix>>20)-1046;     /* e0 = ilogb(z)-23; */
00172        SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20)));
00173        for(i=0;i<2;i++) {
00174               tx[i] = (double)((int32_t)(z));
00175               z     = (z-tx[i])*two24;
00176        }
00177        tx[2] = z;
00178        nx = 3;
00179        while(tx[nx-1]==zero) nx--; /* skip zero term */
00180        n  =  __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
00181        if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
00182        return n;
00183 }