Back to index

glibc  2.9
e_rem_pio2f.c
Go to the documentation of this file.
00001 /* e_rem_pio2f.c -- float version of e_rem_pio2.c
00002  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
00003  */
00004 
00005 /*
00006  * ====================================================
00007  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00008  *
00009  * Developed at SunPro, a Sun Microsystems, Inc. business.
00010  * Permission to use, copy, modify, and distribute this
00011  * software is freely granted, provided that this notice
00012  * is preserved.
00013  * ====================================================
00014  */
00015 
00016 #if defined(LIBM_SCCS) && !defined(lint)
00017 static char rcsid[] = "$NetBSD: e_rem_pio2f.c,v 1.5 1995/05/10 20:46:03 jtc Exp $";
00018 #endif
00019 
00020 /* __ieee754_rem_pio2f(x,y)
00021  *
00022  * return the remainder of x rem pi/2 in y[0]+y[1]
00023  * use __kernel_rem_pio2f()
00024  */
00025 
00026 #include "math.h"
00027 #include "math_private.h"
00028 
00029 /*
00030  * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
00031  */
00032 #ifdef __STDC__
00033 static const int32_t two_over_pi[] = {
00034 #else
00035 static int32_t two_over_pi[] = {
00036 #endif
00037 0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
00038 0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
00039 0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
00040 0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
00041 0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
00042 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
00043 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
00044 0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
00045 0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
00046 0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
00047 0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
00048 0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
00049 0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
00050 0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
00051 0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
00052 0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
00053 0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
00054 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
00055 0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
00056 0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
00057 0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
00058 0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
00059 };
00060 
00061 /* This array is like the one in e_rem_pio2.c, but the numbers are
00062    single precision and the last 8 bits are forced to 0.  */
00063 #ifdef __STDC__
00064 static const int32_t npio2_hw[] = {
00065 #else
00066 static int32_t npio2_hw[] = {
00067 #endif
00068 0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
00069 0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
00070 0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
00071 0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
00072 0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
00073 0x4242c700, 0x42490f00
00074 };
00075 
00076 /*
00077  * invpio2:  24 bits of 2/pi
00078  * pio2_1:   first  17 bit of pi/2
00079  * pio2_1t:  pi/2 - pio2_1
00080  * pio2_2:   second 17 bit of pi/2
00081  * pio2_2t:  pi/2 - (pio2_1+pio2_2)
00082  * pio2_3:   third  17 bit of pi/2
00083  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
00084  */
00085 
00086 #ifdef __STDC__
00087 static const float
00088 #else
00089 static float
00090 #endif
00091 zero =  0.0000000000e+00, /* 0x00000000 */
00092 half =  5.0000000000e-01, /* 0x3f000000 */
00093 two8 =  2.5600000000e+02, /* 0x43800000 */
00094 invpio2 =  6.3661980629e-01, /* 0x3f22f984 */
00095 pio2_1  =  1.5707855225e+00, /* 0x3fc90f80 */
00096 pio2_1t =  1.0804334124e-05, /* 0x37354443 */
00097 pio2_2  =  1.0804273188e-05, /* 0x37354400 */
00098 pio2_2t =  6.0770999344e-11, /* 0x2e85a308 */
00099 pio2_3  =  6.0770943833e-11, /* 0x2e85a300 */
00100 pio2_3t =  6.1232342629e-17; /* 0x248d3132 */
00101 
00102 #ifdef __STDC__
00103        int32_t __ieee754_rem_pio2f(float x, float *y)
00104 #else
00105        int32_t __ieee754_rem_pio2f(x,y)
00106        float x,y[];
00107 #endif
00108 {
00109        float z,w,t,r,fn;
00110        float tx[3];
00111        int32_t e0,i,j,nx,n,ix,hx;
00112 
00113        GET_FLOAT_WORD(hx,x);
00114        ix = hx&0x7fffffff;
00115        if(ix<=0x3f490fd8)   /* |x| ~<= pi/4 , no need for reduction */
00116            {y[0] = x; y[1] = 0; return 0;}
00117        if(ix<0x4016cbe4) {  /* |x| < 3pi/4, special case with n=+-1 */
00118            if(hx>0) {
00119               z = x - pio2_1;
00120               if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
00121                   y[0] = z - pio2_1t;
00122                   y[1] = (z-y[0])-pio2_1t;
00123               } else {             /* near pi/2, use 24+24+24 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            } else {  /* negative x */
00130               z = x + pio2_1;
00131               if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
00132                   y[0] = z + pio2_1t;
00133                   y[1] = (z-y[0])+pio2_1t;
00134               } else {             /* near pi/2, use 24+24+24 bit pi */
00135                   z += pio2_2;
00136                   y[0] = z + pio2_2t;
00137                   y[1] = (z-y[0])+pio2_2t;
00138               }
00139               return -1;
00140            }
00141        }
00142        if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
00143            t  = fabsf(x);
00144            n  = (int32_t) (t*invpio2+half);
00145            fn = (float)n;
00146            r  = t-fn*pio2_1;
00147            w  = fn*pio2_1t; /* 1st round good to 40 bit */
00148            if(n<32&&(int32_t)(ix&0xffffff00)!=npio2_hw[n-1]) {
00149               y[0] = r-w;   /* quick check no cancellation */
00150            } else {
00151                u_int32_t high;
00152                j  = ix>>23;
00153                y[0] = r-w;
00154               GET_FLOAT_WORD(high,y[0]);
00155                i = j-((high>>23)&0xff);
00156                if(i>8) {  /* 2nd iteration needed, good to 57 */
00157                   t  = r;
00158                   w  = fn*pio2_2;
00159                   r  = t-w;
00160                   w  = fn*pio2_2t-((t-r)-w);
00161                   y[0] = r-w;
00162                   GET_FLOAT_WORD(high,y[0]);
00163                   i = j-((high>>23)&0xff);
00164                   if(i>25)  {      /* 3rd iteration need, 74 bits acc */
00165                      t  = r;       /* will cover all possible cases */
00166                      w  = fn*pio2_3;
00167                      r  = t-w;
00168                      w  = fn*pio2_3t-((t-r)-w);
00169                      y[0] = r-w;
00170                   }
00171               }
00172            }
00173            y[1] = (r-y[0])-w;
00174            if(hx<0)  {y[0] = -y[0]; y[1] = -y[1]; return -n;}
00175            else       return n;
00176        }
00177     /*
00178      * all other (large) arguments
00179      */
00180        if(ix>=0x7f800000) {        /* x is inf or NaN */
00181            y[0]=y[1]=x-x; return 0;
00182        }
00183     /* set z = scalbn(|x|,ilogb(x)-7) */
00184        e0     = (ix>>23)-134;             /* e0 = ilogb(z)-7; */
00185        SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23)));
00186        for(i=0;i<2;i++) {
00187               tx[i] = (float)((int32_t)(z));
00188               z     = (z-tx[i])*two8;
00189        }
00190        tx[2] = z;
00191        nx = 3;
00192        while(tx[nx-1]==zero) nx--; /* skip zero term */
00193        n  =  __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi);
00194        if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
00195        return n;
00196 }