root/arch/i386/math-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
  3. log2_kernel
   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,1994                                              |
   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 #include "poly.h"
  19 
  20 
  21 
  22 static void log2_kernel(FPU_REG const *arg,
  23                         Xsig *accum_result, long int *expon);
  24 
  25 
  26 /*--- poly_l2() -------------------------------------------------------------+
  27  |   Base 2 logarithm by a polynomial approximation.                         |
  28  +---------------------------------------------------------------------------*/
  29 void    poly_l2(FPU_REG const *arg, FPU_REG const *y, FPU_REG *result)
     /* [previous][next][first][last][top][bottom][index][help] */
  30 {
  31   long int             exponent, expon, expon_expon;
  32   Xsig                 accumulator, expon_accum, yaccum;
  33   char                 sign;
  34   FPU_REG              x;
  35 
  36 
  37   exponent = arg->exp - EXP_BIAS;
  38 
  39   /* From arg, make a number > sqrt(2)/2 and < sqrt(2) */
  40   if ( arg->sigh > (unsigned)0xb504f334 )
  41     {
  42       /* Treat as  sqrt(2)/2 < arg < 1 */
  43       significand(&x) = - significand(arg);
  44       x.sign = SIGN_NEG;
  45       x.tag = TW_Valid;
  46       x.exp = EXP_BIAS-1;
  47       exponent++;
  48       normalize(&x);
  49     }
  50   else
  51     {
  52       /* Treat as  1 <= arg < sqrt(2) */
  53       x.sigh = arg->sigh - 0x80000000;
  54       x.sigl = arg->sigl;
  55       x.sign = SIGN_POS;
  56       x.tag = TW_Valid;
  57       x.exp = EXP_BIAS;
  58       normalize(&x);
  59     }
  60 
  61   if ( x.tag == TW_Zero )
  62     {
  63       expon = 0;
  64       accumulator.msw = accumulator.midw = accumulator.lsw = 0;
  65     }
  66   else
  67     {
  68       log2_kernel(&x, &accumulator, &expon);
  69     }
  70 
  71   sign = exponent < 0;
  72   if ( sign ) exponent = -exponent;
  73   expon_accum.msw = exponent; expon_accum.midw = expon_accum.lsw = 0;
  74   if ( exponent )
  75     {
  76       expon_expon = 31 + norm_Xsig(&expon_accum);
  77       shr_Xsig(&accumulator, expon_expon - expon);
  78 
  79       if ( sign ^ (x.sign == SIGN_NEG) )
  80         negate_Xsig(&accumulator);
  81       add_Xsig_Xsig(&accumulator, &expon_accum);
  82     }
  83   else
  84     {
  85       expon_expon = expon;
  86       sign = x.sign;
  87     }
  88 
  89   yaccum.lsw = 0; XSIG_LL(yaccum) = significand(y);
  90   mul_Xsig_Xsig(&accumulator, &yaccum);
  91 
  92   expon_expon += round_Xsig(&accumulator);
  93 
  94   if ( accumulator.msw == 0 )
  95     {
  96       reg_move(&CONST_Z, y);
  97     }
  98   else
  99     {
 100       result->exp = expon_expon + y->exp + 1;
 101       significand(result) = XSIG_LL(accumulator);
 102       result->tag = TW_Valid; /* set the tags to Valid */
 103       result->sign = sign ^ y->sign;
 104     }
 105 
 106   return;
 107 }
 108 
 109 
 110 /*--- poly_l2p1() -----------------------------------------------------------+
 111  |   Base 2 logarithm by a polynomial approximation.                         |
 112  |   log2(x+1)                                                               |
 113  +---------------------------------------------------------------------------*/
 114 int     poly_l2p1(FPU_REG const *arg, FPU_REG const *y, FPU_REG *result)
     /* [previous][next][first][last][top][bottom][index][help] */
 115 {
 116   char                 sign;
 117   long int             exponent;
 118   Xsig                 accumulator, yaccum;
 119 
 120 
 121   sign = arg->sign;
 122 
 123   if ( arg->exp < EXP_BIAS )
 124     {
 125       log2_kernel(arg, &accumulator, &exponent);
 126 
 127       yaccum.lsw = 0;
 128       XSIG_LL(yaccum) = significand(y);
 129       mul_Xsig_Xsig(&accumulator, &yaccum);
 130 
 131       exponent += round_Xsig(&accumulator);
 132 
 133       result->exp = exponent + y->exp + 1;
 134       significand(result) = XSIG_LL(accumulator);
 135       result->tag = TW_Valid; /* set the tags to Valid */
 136       result->sign = sign ^ y->sign;
 137 
 138       return 0;
 139     }
 140   else
 141     {
 142       /* The magnitude of arg is far too large. */
 143       reg_move(y, result);
 144       if ( sign != SIGN_POS )
 145         {
 146           /* Trying to get the log of a negative number. */
 147           return 1;
 148         }
 149       else
 150         {
 151           return 0;
 152         }
 153     }
 154 
 155 }
 156 
 157 
 158 
 159 
 160 #undef HIPOWER
 161 #define HIPOWER 10
 162 static const unsigned long long logterms[HIPOWER] =
 163 {
 164   0x2a8eca5705fc2ef0LL,
 165   0xf6384ee1d01febceLL,
 166   0x093bb62877cdf642LL,
 167   0x006985d8a9ec439bLL,
 168   0x0005212c4f55a9c8LL,
 169   0x00004326a16927f0LL,
 170   0x0000038d1d80a0e7LL,
 171   0x0000003141cc80c6LL,
 172   0x00000002b1668c9fLL,
 173   0x000000002c7a46aaLL
 174 };
 175 
 176 static const unsigned long leadterm = 0xb8000000;
 177 
 178 
 179 /*--- log2_kernel() ---------------------------------------------------------+
 180  |   Base 2 logarithm by a polynomial approximation.                         |
 181  |   log2(x+1)                                                               |
 182  +---------------------------------------------------------------------------*/
 183 static void log2_kernel(FPU_REG const *arg, Xsig *accum_result,
     /* [previous][next][first][last][top][bottom][index][help] */
 184                         long int *expon)
 185 {
 186   char                 sign;
 187   long int             exponent, adj;
 188   unsigned long long   Xsq;
 189   Xsig                 accumulator, Numer, Denom, argSignif, arg_signif;
 190 
 191   sign = arg->sign;
 192 
 193   exponent = arg->exp - EXP_BIAS;
 194   Numer.lsw = Denom.lsw = 0;
 195   XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg);
 196   if ( sign == SIGN_POS )
 197     {
 198       shr_Xsig(&Denom, 2 - (1 + exponent));
 199       Denom.msw |= 0x80000000;
 200       div_Xsig(&Numer, &Denom, &argSignif);
 201     }
 202   else
 203     {
 204       shr_Xsig(&Denom, 1 - (1 + exponent));
 205       negate_Xsig(&Denom);
 206       if ( Denom.msw & 0x80000000 )
 207         {
 208           div_Xsig(&Numer, &Denom, &argSignif);
 209           exponent ++;
 210         }
 211       else
 212         {
 213           /* Denom must be 1.0 */
 214           argSignif.lsw = Numer.lsw; argSignif.midw = Numer.midw;
 215           argSignif.msw = Numer.msw;
 216         }
 217     }
 218 
 219 #ifndef PECULIAR_486
 220   /* Should check here that  |local_arg|  is within the valid range */
 221   if ( exponent >= -2 )
 222     {
 223       if ( (exponent > -2) ||
 224           (argSignif.msw > (unsigned)0xafb0ccc0) )
 225         {
 226           /* The argument is too large */
 227         }
 228     }
 229 #endif PECULIAR_486
 230 
 231   arg_signif.lsw = argSignif.lsw; XSIG_LL(arg_signif) = XSIG_LL(argSignif);
 232   adj = norm_Xsig(&argSignif);
 233   accumulator.lsw = argSignif.lsw; XSIG_LL(accumulator) = XSIG_LL(argSignif);
 234   mul_Xsig_Xsig(&accumulator, &accumulator);
 235   shr_Xsig(&accumulator, 2*(-1 - (1 + exponent + adj)));
 236   Xsq = XSIG_LL(accumulator);
 237   if ( accumulator.lsw & 0x80000000 )
 238     Xsq++;
 239 
 240   accumulator.msw = accumulator.midw = accumulator.lsw = 0;
 241   /* Do the basic fixed point polynomial evaluation */
 242   polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER-1);
 243 
 244   mul_Xsig_Xsig(&accumulator, &argSignif);
 245   shr_Xsig(&accumulator, 6 - adj);
 246 
 247   mul32_Xsig(&arg_signif, leadterm);
 248   add_two_Xsig(&accumulator, &arg_signif, &exponent);
 249 
 250   *expon = exponent + 1;
 251   accum_result->lsw = accumulator.lsw;
 252   accum_result->midw = accumulator.midw;
 253   accum_result->msw = accumulator.msw;
 254 
 255 }
/* [previous][next][first][last][top][bottom][index][help] */