root/drivers/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   billm@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 const     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 const *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   unsigned 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 
  75   /* shift num left, lose the ms bit */
  76   num.sigh <<= 1;
  77   if ( num.sigl & 0x80000000 ) num.sigh |= 1;
  78   num.sigl <<= 1;
  79 
  80   denom.sigl = num.sigl;
  81   denom.sigh = num.sigh;
  82   poly_div4(&significand(&denom));
  83   denom.sigh += 0x80000000;                     /* set the msb */
  84   Xx.exp = EXP_BIAS;  /* needed to prevent errors in div routine */
  85   reg_u_div(&num, &denom, &Xx, FULL_PRECISION);
  86 
  87   zero = !(Xx.sigh | Xx.sigl);
  88   
  89   mul64(&significand(&Xx), &significand(&Xx), &Xsq);
  90   poly_div16(&Xsq);
  91 
  92   accum.exp = -1;               /* exponent of accum */
  93 
  94   /* Do the basic fixed point polynomial evaluation */
  95   polynomial((unsigned *)&accum.sigl, (unsigned *)&Xsq, lterms, HIPOWER-1);
  96 
  97   if ( !exponent )
  98     {
  99       /* If the exponent is zero, then we would lose precision by
 100          sticking to fixed point computation here */
 101       /* We need to re-compute Xx because of loss of precision. */
 102       FPU_REG   lXx;
 103       char      sign;
 104       
 105       sign = accum.sign;
 106       accum.sign = 0;
 107 
 108       /* make accum compatible and normalize */
 109       accum.exp = EXP_BIAS + accum.exp;
 110       normalize(&accum);
 111 
 112       if ( zero )
 113         {
 114           reg_move(&CONST_Z, result);
 115         }
 116       else
 117         {
 118           /* we need to re-compute lXx to better accuracy */
 119           num.tag = TW_Valid;           /* set the tags to Vaild */
 120           num.sign = 0;         /* set the sign to positive */
 121           num.exp = EXP_BIAS - 1;
 122           if ( sign )
 123             {
 124               /* The argument is of the form 1-x */
 125               /* Use  1-1/(1-x) = x/(1-x) */
 126               significand(&num) = - significand(arg);
 127               normalize(&num);
 128               reg_div(&num, arg, &num, FULL_PRECISION);
 129             }
 130           else
 131             {
 132               normalize(&num);
 133             }
 134 
 135           denom.tag = TW_Valid; /* set the tags to Valid */
 136           denom.sign = SIGN_POS;        /* set the sign to positive */
 137           denom.exp = EXP_BIAS;
 138           
 139           reg_div(&num, &denom, &lXx, FULL_PRECISION);
 140 
 141           reg_u_mul(&lXx, &accum, &accum, FULL_PRECISION);
 142 
 143           reg_u_add(&lXx, &accum, result, FULL_PRECISION);
 144           
 145           normalize(result);
 146         }
 147 
 148       result->sign = sign;
 149       return;
 150     }
 151 
 152   mul64(&significand(&accum),
 153         &significand(&Xx), &significand(&accum));
 154 
 155   significand(&accum) += significand(&Xx);
 156 
 157   if ( Xx.sigh > accum.sigh )
 158     {
 159       /* There was an overflow */
 160 
 161       poly_div2(&significand(&accum));
 162       accum.sigh |= 0x80000000;
 163       accum.exp++;
 164     }
 165 
 166   /* When we add the exponent to the accum result later, we will
 167      require that their signs are the same. Here we ensure that
 168      this is so. */
 169   if ( exponent && ((exponent < 0) ^ (accum.sign)) )
 170     {
 171       /* signs are different */
 172 
 173       accum.sign = !accum.sign;
 174 
 175       /* An exceptional case is when accum is zero */
 176       if ( accum.sigl | accum.sigh )
 177         {
 178           /* find 1-accum */
 179           /* Shift to get exponent == 0 */
 180           if ( accum.exp < 0 )
 181             {
 182               poly_div2(&significand(&accum));
 183               accum.exp++;
 184             }
 185           /* Just negate, but throw away the sign */
 186           significand(&accum) = - significand(&accum);
 187           if ( exponent < 0 )
 188             exponent++;
 189           else
 190             exponent--;
 191         }
 192     }
 193 
 194   shift = exponent >= 0 ? exponent : -exponent ;
 195   bits = 0;
 196   if ( shift )
 197     {
 198       if ( accum.exp )
 199         {
 200           accum.exp++;
 201           poly_div2(&significand(&accum));
 202         }
 203       while ( shift )
 204         {
 205           poly_div2(&significand(&accum));
 206           if ( shift & 1)
 207             accum.sigh |= 0x80000000;
 208           shift >>= 1;
 209           bits++;
 210         }
 211     }
 212 
 213   /* Convert to 64 bit signed-compatible */
 214   accum.exp += bits + EXP_BIAS - 1;
 215 
 216   reg_move(&accum, result);
 217   normalize(result);
 218 
 219   return;
 220 }
 221 
 222 
 223 /*--- poly_l2p1() -----------------------------------------------------------+
 224  |   Base 2 logarithm by a polynomial approximation.                         |
 225  |   log2(x+1)                                                               |
 226  +---------------------------------------------------------------------------*/
 227 int     poly_l2p1(FPU_REG const *arg, FPU_REG *result)
     /* [previous][next][first][last][top][bottom][index][help] */
 228 {
 229   char          sign = 0;
 230   unsigned long long     Xsq;
 231   FPU_REG       arg_pl1, denom, accum, local_arg, poly_arg;
 232 
 233 
 234   sign = arg->sign;
 235 
 236   reg_add(arg, &CONST_1, &arg_pl1, FULL_PRECISION);
 237 
 238   if ( (arg_pl1.sign) | (arg_pl1.tag) )
 239     {                   /* We need a valid positive number! */
 240       return 1;
 241     }
 242 
 243   reg_add(&CONST_1, &arg_pl1, &denom, FULL_PRECISION);
 244   reg_div(arg, &denom, &local_arg, FULL_PRECISION);
 245   local_arg.sign = 0;   /* Make the sign positive */
 246 
 247   /* Now we need to check that  |local_arg| is less than
 248      3-2*sqrt(2) = 0.17157.. = .0xafb0ccc0 * 2^-2 */
 249 
 250   if ( local_arg.exp >= EXP_BIAS - 3 )
 251     {
 252       if ( (local_arg.exp > EXP_BIAS - 3) ||
 253           (local_arg.sigh > (unsigned)0xafb0ccc0) )
 254         {
 255           /* The argument is large */
 256           poly_l2(&arg_pl1, result); return 0;
 257         }
 258     }
 259 
 260   /* Make a copy of local_arg */
 261   reg_move(&local_arg, &poly_arg);
 262 
 263   /* Get poly_arg bits aligned as required */
 264   shrx((unsigned *)&(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));
 265 
 266   mul64(&significand(&poly_arg), &significand(&poly_arg), &Xsq);
 267   poly_div16(&Xsq);
 268 
 269   /* Do the basic fixed point polynomial evaluation */
 270   polynomial(&(accum.sigl), (unsigned *)&Xsq, lterms, HIPOWER-1);
 271 
 272   accum.tag = TW_Valid; /* set the tags to Valid */
 273   accum.sign = SIGN_POS;        /* and make accum positive */
 274 
 275   /* make accum compatible and normalize */
 276   accum.exp = EXP_BIAS - 1;
 277   normalize(&accum);
 278 
 279   reg_u_mul(&local_arg, &accum, &accum, FULL_PRECISION);
 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] */