root/kernel/FPU-emu/poly_l2.c

/* [previous][next][first][last][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. poly_l2
  2. poly_l2p1

   1 /*---------------------------------------------------------------------------+
   2  |  poly_l2.c                                                                |
   3  |                                                                           |
   4  | Compute the base 2 log of a FPU_REG, using a polynomial approximation.    |
   5  |                                                                           |
   6  | Copyright (C) 1992    W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
   7  |                       Australia.  E-mail apm233m@vaxc.cc.monash.edu.au    |
   8  |                                                                           |
   9  |                                                                           |
  10  +---------------------------------------------------------------------------*/
  11 
  12 
  13 #include "exception.h"
  14 #include "reg_constant.h"
  15 #include "fpu_emu.h"
  16 
  17 
  18 
  19 #define HIPOWER 9
  20 static unsigned short   lterms[HIPOWER][4] =
  21         {
  22         /* Ideal computation with these coeffs gives about
  23            64.6 bit rel accuracy. */
  24         { 0xe177, 0xb82f, 0x7652, 0x7154 },
  25         { 0xee0f, 0xe80f, 0x2770, 0x7b1c },
  26         { 0x0fc0, 0xbe87, 0xb143, 0x49dd },
  27         { 0x78b9, 0xdadd, 0xec54, 0x34c2 },
  28         { 0x003a, 0x5de9, 0x628b, 0x2909 },
  29         { 0x5588, 0xed16, 0x4abf, 0x2193 },
  30         { 0xb461, 0x85f7, 0x347a, 0x1c6a },
  31         { 0x0975, 0x87b3, 0xd5bf, 0x1876 },
  32         { 0xe85c, 0xcec9, 0x84e7, 0x187d }
  33         };
  34 
  35 
  36 
  37 
  38 /*--- poly_l2() -------------------------------------------------------------+
  39  |   Base 2 logarithm by a polynomial approximation.                         |
  40  +---------------------------------------------------------------------------*/
  41 void    poly_l2(FPU_REG *arg, FPU_REG *result)
     /* [previous][next][first][last][top][bottom][index][help] */
  42 {
  43   short           exponent;
  44   char            zero;         /* flag for an Xx == 0 */
  45   unsigned short  bits, shift;
  46   long long       Xsq;
  47   FPU_REG         accum, denom, num, Xx;
  48 
  49 
  50   exponent = arg->exp - EXP_BIAS;
  51 
  52   accum.tag = TW_Valid; /* set the tags to Valid */
  53 
  54   if ( arg->sigh > (unsigned)0xb504f334 )
  55     {
  56       /* This is good enough for the computation of the polynomial
  57          sum, but actually results in a loss of precision for
  58          the computation of Xx. This will matter only if exponent
  59          becomes zero. */
  60       exponent++;
  61       accum.sign = 1;   /* sign to negative */
  62       num.exp = EXP_BIAS;  /* needed to prevent errors in div routine */
  63       reg_u_div((long long *)&(CONST_1.sigl), (long long *)&(arg->sigl), &num);
  64     }
  65   else
  66     {
  67       accum.sign = 0;   /* set the sign to positive */
  68       num.sigl = arg->sigl;             /* copy the mantissa */
  69       num.sigh = arg->sigh;
  70     }
  71 
  72   /* shift num left, lose the ms bit */
  73   num.sigh <<= 1;
  74   if ( num.sigl & 0x80000000 ) num.sigh |= 1;
  75   num.sigl <<= 1;
  76 
  77   denom.sigl = num.sigl;
  78   denom.sigh = num.sigh;
  79   poly_div4((long long *)&(denom.sigl));
  80   denom.sigh += 0x80000000;                     /* set the msb */
  81   Xx.exp = EXP_BIAS;  /* needed to prevent errors in div routine */
  82   reg_u_div((long long *)&num.sigl, (long long *)&(denom.sigl), &Xx);
  83 
  84   zero = !(Xx.sigh | Xx.sigl);
  85   
  86   mul64((long long *)&Xx.sigl, (long long *)&Xx.sigl, &Xsq);
  87   poly_div16(&Xsq);
  88 
  89   accum.exp = -1;               /* exponent of accum */
  90 
  91   /* Do the basic fixed point polynomial evaluation */
  92   polynomial((unsigned *)&accum.sigl, (unsigned *)&Xsq, lterms, HIPOWER-1);
  93 
  94   if ( !exponent )
  95     {
  96       /* If the exponent is zero, then we would lose precision by
  97          sticking to fixed point computation here */
  98       /* We need to re-compute Xx because of loss of precision. */
  99       FPU_REG   lXx;
 100       char      sign;
 101       
 102       sign = accum.sign;
 103       accum.sign = 0;
 104 
 105       /* make accum compatible and normalize */
 106       accum.exp = EXP_BIAS + accum.exp;
 107       normalize(&accum);
 108 
 109       if ( zero )
 110         {
 111           reg_move(&CONST_Z, result);
 112         }
 113       else
 114         {
 115           /* we need to re-compute lXx to better accuracy */
 116           num.tag = TW_Valid;           /* set the tags to Vaild */
 117           num.sign = 0;         /* set the sign to positive */
 118           num.exp = EXP_BIAS - 1;
 119           if ( sign )
 120             {
 121               /* The argument is of the form 1-x */
 122               /* Use  1-1/(1-x) = x/(1-x) */
 123               *((long long *)&num.sigl) = - *((long long *)&(arg->sigl));
 124               normalize(&num);
 125               reg_div(&num, arg, &num);
 126             }
 127           else
 128             {
 129               normalize(&num);
 130             }
 131           
 132           denom.tag = TW_Valid; /* set the tags to Valid */
 133           denom.sign = SIGN_POS;        /* set the sign to positive */
 134           denom.exp = EXP_BIAS;
 135           
 136           reg_div(&num, &denom, &lXx);
 137 
 138           reg_u_mul(&lXx, &accum, &accum);
 139           accum.exp += - EXP_BIAS + 1;
 140 
 141           reg_u_add(&lXx, &accum, result);
 142           
 143           normalize(result);
 144         }
 145       result->sign = sign;
 146       return;
 147     }
 148 
 149   mul64((long long *)&accum.sigl,
 150         (long long *)&Xx.sigl, (long long *)&accum.sigl);
 151 
 152   *((long long *)(&accum.sigl)) += *((long long *)(&Xx.sigl));
 153 
 154   if ( Xx.sigh > accum.sigh )
 155     {
 156       /* There was an overflow */
 157 
 158       poly_div2((long long *)&accum.sigl);
 159       accum.sigh |= 0x80000000;
 160       accum.exp++;
 161     }
 162 
 163   /* When we add the exponent to the accum result later, we will
 164      require that their signs are the same. Here we ensure that
 165      this is so. */
 166   if ( exponent && ((exponent < 0) ^ (accum.sign)) )
 167     {
 168       /* signs are different */
 169 
 170       accum.sign = !accum.sign;
 171 
 172       /* An exceptional case is when accum is zero */
 173       if ( accum.sigl | accum.sigh )
 174         {
 175           /* find 1-accum */
 176           /* Shift to get exponent == 0 */
 177           if ( accum.exp < 0 )
 178             {
 179               poly_div2((long long *)&accum.sigl);
 180               accum.exp++;
 181             }
 182           /* Just negate, but throw away the sign */
 183           *((long long *)&(accum.sigl)) = - *((long long *)&(accum.sigl));
 184           if ( exponent < 0 )
 185             exponent++;
 186           else
 187             exponent--;
 188         }
 189     }
 190 
 191   shift = exponent >= 0 ? exponent : -exponent ;
 192   bits = 0;
 193   if ( shift )
 194     {
 195       if ( accum.exp )
 196         {
 197           accum.exp++;
 198           poly_div2((long long *)&accum.sigl);
 199         }
 200       while ( shift )
 201         {
 202           poly_div2((long long *)&accum.sigl);
 203           if ( shift & 1)
 204             accum.sigh |= 0x80000000;
 205           shift >>= 1;
 206           bits++;
 207         }
 208     }
 209 
 210   /* Convert to 64 bit signed-compatible */
 211   accum.exp += bits + EXP_BIAS - 1;
 212 
 213   reg_move(&accum, result);
 214   normalize(result);
 215 
 216   return;
 217 }
 218 
 219 
 220 /*--- poly_l2p1() -----------------------------------------------------------+
 221  |   Base 2 logarithm by a polynomial approximation.                         |
 222  |   log2(x+1)                                                               |
 223  +---------------------------------------------------------------------------*/
 224 int     poly_l2p1(FPU_REG *arg, FPU_REG *result)
     /* [previous][next][first][last][top][bottom][index][help] */
 225 {
 226   char          sign = 0;
 227   long long     Xsq;
 228   FPU_REG       arg_pl1, denom, accum, local_arg, poly_arg;
 229 
 230 
 231   sign = arg->sign;
 232 
 233   reg_add(arg, &CONST_1, &arg_pl1);
 234 
 235   if ( (arg_pl1.sign) | (arg_pl1.tag) )
 236     {                   /* We need a valid positive number! */
 237       return 1;
 238     }
 239 
 240   reg_add(&CONST_1, &arg_pl1, &denom);
 241   reg_div(arg, &denom, &local_arg);
 242   local_arg.sign = 0;   /* Make the sign positive */
 243 
 244   /* Now we need to check that  |local_arg| is less than
 245      3-2*sqrt(2) = 0.17157.. = .0xafb0ccc0 * 2^-2 */
 246 
 247   if ( local_arg.exp >= EXP_BIAS - 3 )
 248     {
 249       if ( (local_arg.exp > EXP_BIAS - 3) ||
 250           (local_arg.sigh > (unsigned)0xafb0ccc0) )
 251         {
 252           /* The argument is large */
 253           poly_l2(&arg_pl1, result); return 0;
 254         }
 255     }
 256 
 257   /* Make a copy of local_arg */
 258   reg_move(&local_arg, &poly_arg);
 259 
 260   /* Get poly_arg bits aligned as required */
 261   shrx((unsigned *)&(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));
 262 
 263   mul64((long long *)&(poly_arg.sigl), (long long *)&(poly_arg.sigl), &Xsq);
 264   poly_div16(&Xsq);
 265 
 266   /* Do the basic fixed point polynomial evaluation */
 267   polynomial(&(accum.sigl), (unsigned *)&Xsq, lterms, HIPOWER-1);
 268 
 269   accum.tag = TW_Valid; /* set the tags to Valid */
 270   accum.sign = SIGN_POS;        /* and make accum positive */
 271 
 272   /* make accum compatible and normalize */
 273   accum.exp = EXP_BIAS - 1;
 274   normalize(&accum);
 275 
 276   reg_u_mul(&local_arg, &accum, &accum);
 277   accum.exp -= EXP_BIAS - 1;
 278 
 279   reg_u_add(&local_arg, &accum, result);
 280 
 281   /* Multiply the result by 2 */
 282   result->exp++;
 283 
 284   result->sign = sign;
 285   
 286   return 0;
 287 }

/* [previous][next][first][last][top][bottom][index][help] */