Back to index

glibc  2.9
k_rem_pio2f.c
Go to the documentation of this file.
00001 /* k_rem_pio2f.c -- float version of k_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: k_rem_pio2f.c,v 1.4 1995/05/10 20:46:28 jtc Exp $";
00018 #endif
00019 
00020 #include "math.h"
00021 #include "math_private.h"
00022 
00023 /* In the float version, the input parameter x contains 8 bit
00024    integers, not 24 bit integers.  113 bit precision is not supported.  */
00025 
00026 #ifdef __STDC__
00027 static const int init_jk[] = {4,7,9}; /* initial value for jk */
00028 #else
00029 static int init_jk[] = {4,7,9}; 
00030 #endif
00031 
00032 #ifdef __STDC__
00033 static const float PIo2[] = {
00034 #else
00035 static float PIo2[] = {
00036 #endif
00037   1.5703125000e+00, /* 0x3fc90000 */
00038   4.5776367188e-04, /* 0x39f00000 */
00039   2.5987625122e-05, /* 0x37da0000 */
00040   7.5437128544e-08, /* 0x33a20000 */
00041   6.0026650317e-11, /* 0x2e840000 */
00042   7.3896444519e-13, /* 0x2b500000 */
00043   5.3845816694e-15, /* 0x27c20000 */
00044   5.6378512969e-18, /* 0x22d00000 */
00045   8.3009228831e-20, /* 0x1fc40000 */
00046   3.2756352257e-22, /* 0x1bc60000 */
00047   6.3331015649e-25, /* 0x17440000 */
00048 };
00049 
00050 #ifdef __STDC__
00051 static const float                 
00052 #else
00053 static float                
00054 #endif
00055 zero   = 0.0,
00056 one    = 1.0,
00057 two8   =  2.5600000000e+02, /* 0x43800000 */
00058 twon8  =  3.9062500000e-03; /* 0x3b800000 */
00059 
00060 #ifdef __STDC__
00061        int __kernel_rem_pio2f(float *x, float *y, int e0, int nx, int prec, const int32_t *ipio2) 
00062 #else
00063        int __kernel_rem_pio2f(x,y,e0,nx,prec,ipio2)     
00064        float x[], y[]; int e0,nx,prec; int32_t ipio2[];
00065 #endif
00066 {
00067        int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
00068        float z,fw,f[20],fq[20],q[20];
00069 
00070     /* initialize jk*/
00071        jk = init_jk[prec];
00072        jp = jk;
00073 
00074     /* determine jx,jv,q0, note that 3>q0 */
00075        jx =  nx-1;
00076        jv = (e0-3)/8; if(jv<0) jv=0;
00077        q0 =  e0-8*(jv+1);
00078 
00079     /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
00080        j = jv-jx; m = jx+jk;
00081        for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (float) ipio2[j];
00082 
00083     /* compute q[0],q[1],...q[jk] */
00084        for (i=0;i<=jk;i++) {
00085            for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
00086        }
00087 
00088        jz = jk;
00089 recompute:
00090     /* distill q[] into iq[] reversingly */
00091        for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
00092            fw    =  (float)((int32_t)(twon8* z));
00093            iq[i] =  (int32_t)(z-two8*fw);
00094            z     =  q[j-1]+fw;
00095        }
00096 
00097     /* compute n */
00098        z  = __scalbnf(z,q0);              /* actual value of z */
00099        z -= (float)8.0*__floorf(z*(float)0.125); /* trim off integer >= 8 */
00100        n  = (int32_t) z;
00101        z -= (float)n;
00102        ih = 0;
00103        if(q0>0) {    /* need iq[jz-1] to determine n */
00104            i  = (iq[jz-1]>>(8-q0)); n += i;
00105            iq[jz-1] -= i<<(8-q0);
00106            ih = iq[jz-1]>>(7-q0);
00107        } 
00108        else if(q0==0) ih = iq[jz-1]>>8;
00109        else if(z>=(float)0.5) ih=2;
00110 
00111        if(ih>0) {    /* q > 0.5 */
00112            n += 1; carry = 0;
00113            for(i=0;i<jz ;i++) {    /* compute 1-q */
00114               j = iq[i];
00115               if(carry==0) {
00116                   if(j!=0) {
00117                      carry = 1; iq[i] = 0x100- j;
00118                   }
00119               } else  iq[i] = 0xff - j;
00120            }
00121            if(q0>0) {              /* rare case: chance is 1 in 12 */
00122                switch(q0) {
00123                case 1:
00124                  iq[jz-1] &= 0x7f; break;
00125               case 2:
00126                  iq[jz-1] &= 0x3f; break;
00127                }
00128            }
00129            if(ih==2) {
00130               z = one - z;
00131               if(carry!=0) z -= __scalbnf(one,q0);
00132            }
00133        }
00134 
00135     /* check if recomputation is needed */
00136        if(z==zero) {
00137            j = 0;
00138            for (i=jz-1;i>=jk;i--) j |= iq[i];
00139            if(j==0) { /* need recomputation */
00140               for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
00141 
00142               for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
00143                   f[jx+i] = (float) ipio2[jv+i];
00144                   for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
00145                   q[i] = fw;
00146               }
00147               jz += k;
00148               goto recompute;
00149            }
00150        }
00151 
00152     /* chop off zero terms */
00153        if(z==(float)0.0) {
00154            jz -= 1; q0 -= 8;
00155            while(iq[jz]==0) { jz--; q0-=8;}
00156        } else { /* break z into 8-bit if necessary */
00157            z = __scalbnf(z,-q0);
00158            if(z>=two8) { 
00159               fw = (float)((int32_t)(twon8*z));
00160               iq[jz] = (int32_t)(z-two8*fw);
00161               jz += 1; q0 += 8;
00162               iq[jz] = (int32_t) fw;
00163            } else iq[jz] = (int32_t) z ;
00164        }
00165 
00166     /* convert integer "bit" chunk to floating-point value */
00167        fw = __scalbnf(one,q0);
00168        for(i=jz;i>=0;i--) {
00169            q[i] = fw*(float)iq[i]; fw*=twon8;
00170        }
00171 
00172     /* compute PIo2[0,...,jp]*q[jz,...,0] */
00173        for(i=jz;i>=0;i--) {
00174            for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
00175            fq[jz-i] = fw;
00176        }
00177 
00178     /* compress fq[] into y[] */
00179        switch(prec) {
00180            case 0:
00181               fw = 0.0;
00182               for (i=jz;i>=0;i--) fw += fq[i];
00183               y[0] = (ih==0)? fw: -fw; 
00184               break;
00185            case 1:
00186            case 2:
00187               fw = 0.0;
00188               for (i=jz;i>=0;i--) fw += fq[i]; 
00189               y[0] = (ih==0)? fw: -fw; 
00190               fw = fq[0]-fw;
00191               for (i=1;i<=jz;i++) fw += fq[i];
00192               y[1] = (ih==0)? fw: -fw; 
00193               break;
00194            case 3:   /* painful */
00195               for (i=jz;i>0;i--) {
00196                   fw      = fq[i-1]+fq[i]; 
00197                   fq[i]  += fq[i-1]-fw;
00198                   fq[i-1] = fw;
00199               }
00200               for (i=jz;i>1;i--) {
00201                   fw      = fq[i-1]+fq[i]; 
00202                   fq[i]  += fq[i-1]-fw;
00203                   fq[i-1] = fw;
00204               }
00205               for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; 
00206               if(ih==0) {
00207                   y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
00208               } else {
00209                   y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
00210               }
00211        }
00212        return n&7;
00213 }