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 logarithm of a 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(REG *arg, 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   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       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 /*        normalize(&accum); ********/
 141 
 142           reg_u_add(&lXx, &accum, result);
 143           
 144           normalize(result);
 145         }
 146       result->sign = sign;
 147       return;
 148     }
 149 
 150   mul64((long long *)&accum.sigl,
 151         (long long *)&Xx.sigl, (long long *)&accum.sigl);
 152 
 153   *((long long *)(&accum.sigl)) += *((long long *)(&Xx.sigl));
 154 
 155   if ( Xx.sigh > accum.sigh )
 156     {
 157       /* There was an overflow */
 158 
 159       poly_div2((long long *)&accum.sigl);
 160       accum.sigh |= 0x80000000;
 161       accum.exp++;
 162     }
 163 
 164   /* When we add the exponent to the accum result later, we will
 165      require that their signs are the same. Here we ensure that
 166      this is so. */
 167   if ( exponent && ((exponent < 0) ^ (accum.sign)) )
 168     {
 169       /* signs are different */
 170 
 171       accum.sign = !accum.sign;
 172 
 173       /* An exceptional case is when accum is zero */
 174       if ( accum.sigl | accum.sigh )
 175         {
 176           /* find 1-accum */
 177           /* Shift to get exponent == 0 */
 178           if ( accum.exp < 0 )
 179             {
 180               poly_div2((long long *)&accum.sigl);
 181               accum.exp++;
 182             }
 183           /* Just negate, but throw away the sign */
 184           *((long long *)&(accum.sigl)) = - *((long long *)&(accum.sigl));
 185           if ( exponent < 0 )
 186             exponent++;
 187           else
 188             exponent--;
 189         }
 190     }
 191 
 192   shift = exponent >= 0 ? exponent : -exponent ;
 193   bits = 0;
 194   if ( shift )
 195     {
 196       if ( accum.exp )
 197         {
 198           accum.exp++;
 199           poly_div2((long long *)&accum.sigl);
 200         }
 201       while ( shift )
 202         {
 203           poly_div2((long long *)&accum.sigl);
 204           if ( shift & 1)
 205             accum.sigh |= 0x80000000;
 206           shift >>= 1;
 207           bits++;
 208         }
 209     }
 210 
 211   /* Convert to 64 bit signed-compatible */
 212   accum.exp += bits + EXP_BIAS - 1;
 213 
 214   reg_move(&accum, result);
 215   normalize(result);
 216 
 217   return;
 218 }
 219 
 220 
 221 /*--- poly_l2p1() -----------------------------------------------------------+
 222  |   Base 2 logarithm by a polynomial approximation.                         |
 223  |   log2(x+1)                                                               |
 224  +---------------------------------------------------------------------------*/
 225 int     poly_l2p1(REG *arg, REG *result)
     /* [previous][next][first][last][top][bottom][index][help] */
 226 {
 227   char          sign = 0;
 228   long long     Xsq;
 229   REG           arg_pl1, denom, accum, local_arg, poly_arg;
 230 
 231 
 232   sign = arg->sign;
 233 
 234   reg_add(arg, &CONST_1, &arg_pl1);
 235 
 236   if ( (arg_pl1.sign) | (arg_pl1.tag) )
 237     {                   /* We need a valid positive number! */
 238       return 1;
 239     }
 240 
 241   reg_add(&CONST_1, &arg_pl1, &denom);
 242   reg_div(arg, &denom, &local_arg);
 243   local_arg.sign = 0;   /* Make the sign positive */
 244 
 245   /* Now we need to check that  |local_arg| is less than
 246      3-2*sqrt(2) = 0.17157.. = .0xafb0ccc0 * 2^-2 */
 247 
 248   if ( local_arg.exp >= EXP_BIAS - 3 )
 249     {
 250       if ( (local_arg.exp > EXP_BIAS - 3) ||
 251           (local_arg.sigh > (unsigned)0xafb0ccc0) )
 252         {
 253           /* The argument is large */
 254           poly_l2(&arg_pl1, result); return 0;
 255         }
 256     }
 257 
 258   /* Make a copy of local_arg */
 259   reg_move(&local_arg, &poly_arg);
 260 
 261   /* Get poly_arg bits aligned as required */
 262   shrx((unsigned *)&(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));
 263 
 264   mul64((long long *)&(poly_arg.sigl), (long long *)&(poly_arg.sigl), &Xsq);
 265   poly_div16(&Xsq);
 266 
 267   /* Do the basic fixed point polynomial evaluation */
 268   polynomial(&(accum.sigl), (unsigned *)&Xsq, lterms, HIPOWER-1);
 269 
 270   accum.tag = TW_Valid; /* set the tags to Valid */
 271   accum.sign = SIGN_POS;        /* and make accum positive */
 272 
 273   /* make accum compatible and normalize */
 274   accum.exp = EXP_BIAS - 1;
 275   normalize(&accum);
 276 
 277   reg_u_mul(&local_arg, &accum, &accum);
 278   accum.exp -= EXP_BIAS - 1;
 279 
 280   reg_u_add(&local_arg, &accum, result);
 281 
 282   /* Multiply the result by 2 */
 283   result->exp++;
 284 
 285   result->sign = sign;
 286   
 287   return 0;
 288 }

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