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

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