git.gag.com Git - fw/altos/blob - src/math/ef_rem_pio2.c

   1 /* ef_rem_pio2.c -- float version of e_rem_pio2.c
   2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   3  */
   4
   5 /*
   6  * ====================================================
   7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   8  *
   9  * Developed at SunPro, a Sun Microsystems, Inc. business.
  10  * Permission to use, copy, modify, and distribute this
  11  * software is freely granted, provided that this notice
  12  * is preserved.
  13  * ====================================================
  14  *
  15  */
  16
  17 /* __ieee754_rem_pio2f(x,y)
  18  *
  19  * return the remainder of x rem pi/2 in y[0]+y[1]
  20  * use __kernel_rem_pio2f()
  21  */
  22
  23 #include "fdlibm.h"
  24
  25 /*
  26  * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
  27  */
  28 #ifdef __STDC__
  29 static const __int32_t two_over_pi[] = {
  30 #else
  31 static __int32_t two_over_pi[] = {
  32 #endif
  33 0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
  34 0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
  35 0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
  36 0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
  37 0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
  38 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
  39 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
  40 0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
  41 0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
  42 0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
  43 0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
  44 0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
  45 0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
  46 0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
  47 0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
  48 0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
  49 0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
  50 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
  51 0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
  52 0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
  53 0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
  54 0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
  55 };
  56
  57 /* This array is like the one in e_rem_pio2.c, but the numbers are
  58    single precision and the last 8 bits are forced to 0.  */
  59 #ifdef __STDC__
  60 static const __int32_t npio2_hw[] = {
  61 #else
  62 static __int32_t npio2_hw[] = {
  63 #endif
  64 0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
  65 0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
  66 0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
  67 0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
  68 0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
  69 0x4242c700, 0x42490f00
  70 };
  71
  72 /*
  73  * invpio2:  24 bits of 2/pi
  74  * pio2_1:   first  17 bit of pi/2
  75  * pio2_1t:  pi/2 - pio2_1
  76  * pio2_2:   second 17 bit of pi/2
  77  * pio2_2t:  pi/2 - (pio2_1+pio2_2)
  78  * pio2_3:   third  17 bit of pi/2
  79  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
  80  */
  81
  82 #ifdef __STDC__
  83 static const float
  84 #else
  85 static float
  86 #endif
  87 zero =  0.0000000000e+00, /* 0x00000000 */
  88 half =  5.0000000000e-01, /* 0x3f000000 */
  89 two8 =  2.5600000000e+02, /* 0x43800000 */
  90 invpio2 =  6.3661980629e-01, /* 0x3f22f984 */
  91 pio2_1  =  1.5707855225e+00, /* 0x3fc90f80 */
  92 pio2_1t =  1.0804334124e-05, /* 0x37354443 */
  93 pio2_2  =  1.0804273188e-05, /* 0x37354400 */
  94 pio2_2t =  6.0770999344e-11, /* 0x2e85a308 */
  95 pio2_3  =  6.0770943833e-11, /* 0x2e85a300 */
  96 pio2_3t =  6.1232342629e-17; /* 0x248d3132 */
  97
  98 #ifdef __STDC__
  99         __int32_t __ieee754_rem_pio2f(float x, float *y)
 100 #else
 101         __int32_t __ieee754_rem_pio2f(x,y)
 102         float x,y[];
 103 #endif
 104 {
 105         float z,w,t,r,fn;
 106         float tx[3];
 107         __int32_t i,j,n,ix,hx;
 108         int e0,nx;
 109
 110         GET_FLOAT_WORD(hx,x);
 111         ix = hx&0x7fffffff;
 112         if(ix<=0x3f490fd8)   /* |x| ~<= pi/4 , no need for reduction */
 113             {y[0] = x; y[1] = 0; return 0;}
 114         if(ix<0x4016cbe4) {  /* |x| < 3pi/4, special case with n=+-1 */
 115             if(hx>0) {
 116                 z = x - pio2_1;
 117                 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
 118                     y[0] = z - pio2_1t;
 119                     y[1] = (z-y[0])-pio2_1t;
 120                 } else {                /* near pi/2, use 24+24+24 bit pi */
 121                     z -= pio2_2;
 122                     y[0] = z - pio2_2t;
 123                     y[1] = (z-y[0])-pio2_2t;
 124                 }
 125                 return 1;
 126             } else {    /* negative x */
 127                 z = x + pio2_1;
 128                 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
 129                     y[0] = z + pio2_1t;
 130                     y[1] = (z-y[0])+pio2_1t;
 131                 } else {                /* near pi/2, use 24+24+24 bit pi */
 132                     z += pio2_2;
 133                     y[0] = z + pio2_2t;
 134                     y[1] = (z-y[0])+pio2_2t;
 135                 }
 136                 return -1;
 137             }
 138         }
 139         if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
 140             t  = fabsf(x);
 141             n  = (__int32_t) (t*invpio2+half);
 142             fn = (float)n;
 143             r  = t-fn*pio2_1;
 144             w  = fn*pio2_1t;    /* 1st round good to 40 bit */
 145             if(n<32&&(ix&(__int32_t)0xffffff00)!=npio2_hw[n-1]) {
 146                 y[0] = r-w;     /* quick check no cancellation */
 147             } else {
 148                 __uint32_t high;
 149                 j  = ix>>23;
 150                 y[0] = r-w;
 151                 GET_FLOAT_WORD(high,y[0]);
 152                 i = j-((high>>23)&0xff);
 153                 if(i>8) {  /* 2nd iteration needed, good to 57 */
 154                     t  = r;
 155                     w  = fn*pio2_2;
 156                     r  = t-w;
 157                     w  = fn*pio2_2t-((t-r)-w);
 158                     y[0] = r-w;
 159                     GET_FLOAT_WORD(high,y[0]);
 160                     i = j-((high>>23)&0xff);
 161                     if(i>25)  { /* 3rd iteration need, 74 bits acc */
 162                         t  = r; /* will cover all possible cases */
 163                         w  = fn*pio2_3;
 164                         r  = t-w;
 165                         w  = fn*pio2_3t-((t-r)-w);
 166                         y[0] = r-w;
 167                     }
 168                 }
 169             }
 170             y[1] = (r-y[0])-w;
 171             if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}
 172             else         return n;
 173         }
 174     /*
 175      * all other (large) arguments
 176      */
 177         if(!FLT_UWORD_IS_FINITE(ix)) {
 178             y[0]=y[1]=x-x; return 0;
 179         }
 180     /* set z = scalbn(|x|,ilogb(x)-7) */
 181         e0      = (int)((ix>>23)-134);  /* e0 = ilogb(z)-7; */
 182         SET_FLOAT_WORD(z, ix - ((__int32_t)e0<<23));
 183         for(i=0;i<2;i++) {
 184                 tx[i] = (float)((__int32_t)(z));
 185                 z     = (z-tx[i])*two8;
 186         }
 187         tx[2] = z;
 188         nx = 3;
 189         while(tx[nx-1]==zero) nx--;     /* skip zero term */
 190         n  =  __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi);
 191         if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
 192         return n;
 193 }