arch/i386/math-emu/poly_l2.c

   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,1997                                         |
   7  |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
   8  |                  E-mail   billm@suburbia.net                              |
   9  |                                                                           |
  10  |                                                                           |
  11  +---------------------------------------------------------------------------*/
  12
  13
  14 #include "exception.h"
  15 #include "reg_constant.h"
  16 #include "fpu_emu.h"
  17 #include "fpu_system.h"
  18 #include "control_w.h"
  19 #include "poly.h"
  20
  21
  22 static void log2_kernel(FPU_REG const *arg, u_char argsign,
  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 *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign)
  30 {
  31   long int             exponent, expon, expon_expon;
  32   Xsig                 accumulator, expon_accum, yaccum;
  33   u_char                       sign, argsign;
  34   FPU_REG              x;
  35   int                  tag;
  36
  37   exponent = exponent16(st0_ptr);
  38
  39   /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */
  40   if ( st0_ptr->sigh > (unsigned)0xb504f334 )
  41     {
  42       /* Treat as  sqrt(2)/2 < st0_ptr < 1 */
  43       significand(&x) = - significand(st0_ptr);
  44       setexponent16(&x, -1);
  45       exponent++;
  46       argsign = SIGN_NEG;
  47     }
  48   else
  49     {
  50       /* Treat as  1 <= st0_ptr < sqrt(2) */
  51       x.sigh = st0_ptr->sigh - 0x80000000;
  52       x.sigl = st0_ptr->sigl;
  53       setexponent16(&x, 0);
  54       argsign = SIGN_POS;
  55     }
  56   tag = FPU_normalize_nuo(&x);
  57
  58   if ( tag == TAG_Zero )
  59     {
  60       expon = 0;
  61       accumulator.msw = accumulator.midw = accumulator.lsw = 0;
  62     }
  63   else
  64     {
  65       log2_kernel(&x, argsign, &accumulator, &expon);
  66     }
  67
  68   if ( exponent < 0 )
  69     {
  70       sign = SIGN_NEG;
  71       exponent = -exponent;
  72     }
  73   else
  74     sign = SIGN_POS;
  75   expon_accum.msw = exponent; expon_accum.midw = expon_accum.lsw = 0;
  76   if ( exponent )
  77     {
  78       expon_expon = 31 + norm_Xsig(&expon_accum);
  79       shr_Xsig(&accumulator, expon_expon - expon);
  80
  81       if ( sign ^ argsign )
  82         negate_Xsig(&accumulator);
  83       add_Xsig_Xsig(&accumulator, &expon_accum);
  84     }
  85   else
  86     {
  87       expon_expon = expon;
  88       sign = argsign;
  89     }
  90
  91   yaccum.lsw = 0; XSIG_LL(yaccum) = significand(st1_ptr);
  92   mul_Xsig_Xsig(&accumulator, &yaccum);
  93
  94   expon_expon += round_Xsig(&accumulator);
  95
  96   if ( accumulator.msw == 0 )
  97     {
  98       FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
  99       return;
 100     }
 101
 102   significand(st1_ptr) = XSIG_LL(accumulator);
 103   setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1);
 104
 105   tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign);
 106   FPU_settagi(1, tag);
 107
 108   set_precision_flag_up();  /* 80486 appears to always do this */
 109
 110   return;
 111
 112 }
 113
 114
 115 /*--- poly_l2p1() -----------------------------------------------------------+
 116  |   Base 2 logarithm by a polynomial approximation.                         |
 117  |   log2(x+1)                                                               |
 118  +---------------------------------------------------------------------------*/
 119 int     poly_l2p1(u_char sign0, u_char sign1,
 120                   FPU_REG *st0_ptr, FPU_REG *st1_ptr, FPU_REG *dest)
 121 {
 122   u_char                tag;
 123   long int              exponent;
 124   Xsig                  accumulator, yaccum;
 125
 126   if ( exponent16(st0_ptr) < 0 )
 127     {
 128       log2_kernel(st0_ptr, sign0, &accumulator, &exponent);
 129
 130       yaccum.lsw = 0;
 131       XSIG_LL(yaccum) = significand(st1_ptr);
 132       mul_Xsig_Xsig(&accumulator, &yaccum);
 133
 134       exponent += round_Xsig(&accumulator);
 135
 136       exponent += exponent16(st1_ptr) + 1;
 137       if ( exponent < EXP_WAY_UNDER ) exponent = EXP_WAY_UNDER;
 138
 139       significand(dest) = XSIG_LL(accumulator);
 140       setexponent16(dest, exponent);
 141
 142       tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1);
 143       FPU_settagi(1, tag);
 144
 145       if ( tag == TAG_Valid )
 146         set_precision_flag_up();   /* 80486 appears to always do this */
 147     }
 148   else
 149     {
 150       /* The magnitude of st0_ptr is far too large. */
 151
 152       if ( sign0 != SIGN_POS )
 153         {
 154           /* Trying to get the log of a negative number. */
 155 #ifdef PECULIAR_486   /* Stupid 80486 doesn't worry about log(negative). */
 156           changesign(st1_ptr);
 157 #else
 158           if ( arith_invalid(1) < 0 )
 159             return 1;
 160 #endif /* PECULIAR_486 */
 161         }
 162
 163       /* 80486 appears to do this */
 164       if ( sign0 == SIGN_NEG )
 165         set_precision_flag_down();
 166       else
 167         set_precision_flag_up();
 168     }
 169
 170   if ( exponent(dest) <= EXP_UNDER )
 171     EXCEPTION(EX_Underflow);
 172
 173   return 0;
 174
 175 }
 176
 177
 178
 179
 180 #undef HIPOWER
 181 #define HIPOWER 10
 182 static const unsigned long long logterms[HIPOWER] =
 183 {
 184   0x2a8eca5705fc2ef0LL,
 185   0xf6384ee1d01febceLL,
 186   0x093bb62877cdf642LL,
 187   0x006985d8a9ec439bLL,
 188   0x0005212c4f55a9c8LL,
 189   0x00004326a16927f0LL,
 190   0x0000038d1d80a0e7LL,
 191   0x0000003141cc80c6LL,
 192   0x00000002b1668c9fLL,
 193   0x000000002c7a46aaLL
 194 };
 195
 196 static const unsigned long leadterm = 0xb8000000;
 197
 198
 199 /*--- log2_kernel() ---------------------------------------------------------+
 200  |   Base 2 logarithm by a polynomial approximation.                         |
 201  |   log2(x+1)                                                               |
 202  +---------------------------------------------------------------------------*/
 203 static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result,
 204                         long int *expon)
 205 {
 206   long int             exponent, adj;
 207   unsigned long long   Xsq;
 208   Xsig                 accumulator, Numer, Denom, argSignif, arg_signif;
 209
 210   exponent = exponent16(arg);
 211   Numer.lsw = Denom.lsw = 0;
 212   XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg);
 213   if ( argsign == SIGN_POS )
 214     {
 215       shr_Xsig(&Denom, 2 - (1 + exponent));
 216       Denom.msw |= 0x80000000;
 217       div_Xsig(&Numer, &Denom, &argSignif);
 218     }
 219   else
 220     {
 221       shr_Xsig(&Denom, 1 - (1 + exponent));
 222       negate_Xsig(&Denom);
 223       if ( Denom.msw & 0x80000000 )
 224         {
 225           div_Xsig(&Numer, &Denom, &argSignif);
 226           exponent ++;
 227         }
 228       else
 229         {
 230           /* Denom must be 1.0 */
 231           argSignif.lsw = Numer.lsw; argSignif.midw = Numer.midw;
 232           argSignif.msw = Numer.msw;
 233         }
 234     }
 235
 236 #ifndef PECULIAR_486
 237   /* Should check here that  |local_arg|  is within the valid range */
 238   if ( exponent >= -2 )
 239     {
 240       if ( (exponent > -2) ||
 241           (argSignif.msw > (unsigned)0xafb0ccc0) )
 242         {
 243           /* The argument is too large */
 244         }
 245     }
 246 #endif /* PECULIAR_486 */
 247
 248   arg_signif.lsw = argSignif.lsw; XSIG_LL(arg_signif) = XSIG_LL(argSignif);
 249   adj = norm_Xsig(&argSignif);
 250   accumulator.lsw = argSignif.lsw; XSIG_LL(accumulator) = XSIG_LL(argSignif);
 251   mul_Xsig_Xsig(&accumulator, &accumulator);
 252   shr_Xsig(&accumulator, 2*(-1 - (1 + exponent + adj)));
 253   Xsq = XSIG_LL(accumulator);
 254   if ( accumulator.lsw & 0x80000000 )
 255     Xsq++;
 256
 257   accumulator.msw = accumulator.midw = accumulator.lsw = 0;
 258   /* Do the basic fixed point polynomial evaluation */
 259   polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER-1);
 260
 261   mul_Xsig_Xsig(&accumulator, &argSignif);
 262   shr_Xsig(&accumulator, 6 - adj);
 263
 264   mul32_Xsig(&arg_signif, leadterm);
 265   add_two_Xsig(&accumulator, &arg_signif, &exponent);
 266
 267   *expon = exponent + 1;
 268   accum_result->lsw = accumulator.lsw;
 269   accum_result->midw = accumulator.midw;
 270   accum_result->msw = accumulator.msw;
 271
 272 }