1 files changed, 272 insertions, 0 deletions

diff --git a/arch/i386/math-emu/poly_l2.c b/arch/i386/math-emu/poly_l2.c
new file mode 100644
index 000000000000..dd00e1d5b074
--- /dev/null
+++ b/arch/i386/math-emu/poly_l2.c
@@ -0,0 +1,272 @@
+/*---------------------------------------------------------------------------+
+ |  poly_l2.c                                                                |
+ |                                                                           |
+ | Compute the base 2 log of a FPU_REG, using a polynomial approximation.    |
+ |                                                                           |
+ | Copyright (C) 1992,1993,1994,1997                                         |
+ |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
+ |                  E-mail   billm@suburbia.net                              |
+ |                                                                           |
+ |                                                                           |
+ +---------------------------------------------------------------------------*/
+
+
+#include "exception.h"
+#include "reg_constant.h"
+#include "fpu_emu.h"
+#include "fpu_system.h"
+#include "control_w.h"
+#include "poly.h"
+
+
+static void log2_kernel(FPU_REG const *arg, u_char argsign,
+                        Xsig *accum_result, long int *expon);
+
+
+/*--- poly_l2() -------------------------------------------------------------+
+ |   Base 2 logarithm by a polynomial approximation.                         |
+ +---------------------------------------------------------------------------*/
+void    poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign)
+{
+  long int             exponent, expon, expon_expon;
+  Xsig                 accumulator, expon_accum, yaccum;
+  u_char                       sign, argsign;
+  FPU_REG              x;
+  int                  tag;
+
+  exponent = exponent16(st0_ptr);
+
+  /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */
+  if ( st0_ptr->sigh > (unsigned)0xb504f334 )
+    {
+      /* Treat as  sqrt(2)/2 < st0_ptr < 1 */
+      significand(&x) = - significand(st0_ptr);
+      setexponent16(&x, -1);
+      exponent++;
+      argsign = SIGN_NEG;
+    }
+  else
+    {
+      /* Treat as  1 <= st0_ptr < sqrt(2) */
+      x.sigh = st0_ptr->sigh - 0x80000000;
+      x.sigl = st0_ptr->sigl;
+      setexponent16(&x, 0);
+      argsign = SIGN_POS;
+    }
+  tag = FPU_normalize_nuo(&x);
+
+  if ( tag == TAG_Zero )
+    {
+      expon = 0;
+      accumulator.msw = accumulator.midw = accumulator.lsw = 0;
+    }
+  else
+    {
+      log2_kernel(&x, argsign, &accumulator, &expon);
+    }
+
+  if ( exponent < 0 )
+    {
+      sign = SIGN_NEG;
+      exponent = -exponent;
+    }
+  else
+    sign = SIGN_POS;
+  expon_accum.msw = exponent; expon_accum.midw = expon_accum.lsw = 0;
+  if ( exponent )
+    {
+      expon_expon = 31 + norm_Xsig(&expon_accum);
+      shr_Xsig(&accumulator, expon_expon - expon);
+
+      if ( sign ^ argsign )
+        negate_Xsig(&accumulator);
+      add_Xsig_Xsig(&accumulator, &expon_accum);
+    }
+  else
+    {
+      expon_expon = expon;
+      sign = argsign;
+    }
+
+  yaccum.lsw = 0; XSIG_LL(yaccum) = significand(st1_ptr);
+  mul_Xsig_Xsig(&accumulator, &yaccum);
+
+  expon_expon += round_Xsig(&accumulator);
+
+  if ( accumulator.msw == 0 )
+    {
+      FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
+      return;
+    }
+
+  significand(st1_ptr) = XSIG_LL(accumulator);
+  setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1);
+
+  tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign);
+  FPU_settagi(1, tag);
+
+  set_precision_flag_up();  /* 80486 appears to always do this */
+
+  return;
+
+}
+
+
+/*--- poly_l2p1() -----------------------------------------------------------+
+ |   Base 2 logarithm by a polynomial approximation.                         |
+ |   log2(x+1)                                                               |
+ +---------------------------------------------------------------------------*/
+int     poly_l2p1(u_char sign0, u_char sign1,
+                  FPU_REG *st0_ptr, FPU_REG *st1_ptr, FPU_REG *dest)
+{
+  u_char                tag;
+  long int              exponent;
+  Xsig                  accumulator, yaccum;
+
+  if ( exponent16(st0_ptr) < 0 )
+    {
+      log2_kernel(st0_ptr, sign0, &accumulator, &exponent);
+
+      yaccum.lsw = 0;
+      XSIG_LL(yaccum) = significand(st1_ptr);
+      mul_Xsig_Xsig(&accumulator, &yaccum);
+
+      exponent += round_Xsig(&accumulator);
+
+      exponent += exponent16(st1_ptr) + 1;
+      if ( exponent < EXP_WAY_UNDER ) exponent = EXP_WAY_UNDER;
+
+      significand(dest) = XSIG_LL(accumulator);
+      setexponent16(dest, exponent);
+
+      tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1);
+      FPU_settagi(1, tag);
+
+      if ( tag == TAG_Valid )
+        set_precision_flag_up();   /* 80486 appears to always do this */
+    }
+  else
+    {
+      /* The magnitude of st0_ptr is far too large. */
+
+      if ( sign0 != SIGN_POS )
+        {
+          /* Trying to get the log of a negative number. */
+#ifdef PECULIAR_486   /* Stupid 80486 doesn't worry about log(negative). */
+          changesign(st1_ptr);
+#else
+          if ( arith_invalid(1) < 0 )
+            return 1;
+#endif /* PECULIAR_486 */
+        }
+
+      /* 80486 appears to do this */
+      if ( sign0 == SIGN_NEG )
+        set_precision_flag_down();
+      else
+        set_precision_flag_up();
+    }
+
+  if ( exponent(dest) <= EXP_UNDER )
+    EXCEPTION(EX_Underflow);
+
+  return 0;
+
+}
+
+
+
+
+#undef HIPOWER
+#define HIPOWER 10
+static const unsigned long long logterms[HIPOWER] =
+{
+  0x2a8eca5705fc2ef0LL,
+  0xf6384ee1d01febceLL,
+  0x093bb62877cdf642LL,
+  0x006985d8a9ec439bLL,
+  0x0005212c4f55a9c8LL,
+  0x00004326a16927f0LL,
+  0x0000038d1d80a0e7LL,
+  0x0000003141cc80c6LL,
+  0x00000002b1668c9fLL,
+  0x000000002c7a46aaLL
+};
+
+static const unsigned long leadterm = 0xb8000000;
+
+
+/*--- log2_kernel() ---------------------------------------------------------+
+ |   Base 2 logarithm by a polynomial approximation.                         |
+ |   log2(x+1)                                                               |
+ +---------------------------------------------------------------------------*/
+static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result,
+                        long int *expon)
+{
+  long int             exponent, adj;
+  unsigned long long   Xsq;
+  Xsig                 accumulator, Numer, Denom, argSignif, arg_signif;
+
+  exponent = exponent16(arg);
+  Numer.lsw = Denom.lsw = 0;
+  XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg);
+  if ( argsign == SIGN_POS )
+    {
+      shr_Xsig(&Denom, 2 - (1 + exponent));
+      Denom.msw |= 0x80000000;
+      div_Xsig(&Numer, &Denom, &argSignif);
+    }
+  else
+    {
+      shr_Xsig(&Denom, 1 - (1 + exponent));
+      negate_Xsig(&Denom);
+      if ( Denom.msw & 0x80000000 )
+        {
+          div_Xsig(&Numer, &Denom, &argSignif);
+          exponent ++;
+        }
+      else
+        {
+          /* Denom must be 1.0 */
+          argSignif.lsw = Numer.lsw; argSignif.midw = Numer.midw;
+          argSignif.msw = Numer.msw;
+        }
+    }
+
+#ifndef PECULIAR_486
+  /* Should check here that  |local_arg|  is within the valid range */
+  if ( exponent >= -2 )
+    {
+      if ( (exponent > -2) ||
+          (argSignif.msw > (unsigned)0xafb0ccc0) )
+        {
+          /* The argument is too large */
+        }
+    }
+#endif /* PECULIAR_486 */
+
+  arg_signif.lsw = argSignif.lsw; XSIG_LL(arg_signif) = XSIG_LL(argSignif);
+  adj = norm_Xsig(&argSignif);
+  accumulator.lsw = argSignif.lsw; XSIG_LL(accumulator) = XSIG_LL(argSignif);
+  mul_Xsig_Xsig(&accumulator, &accumulator);
+  shr_Xsig(&accumulator, 2*(-1 - (1 + exponent + adj)));
+  Xsq = XSIG_LL(accumulator);
+  if ( accumulator.lsw & 0x80000000 )
+    Xsq++;
+
+  accumulator.msw = accumulator.midw = accumulator.lsw = 0;
+  /* Do the basic fixed point polynomial evaluation */
+  polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER-1);
+
+  mul_Xsig_Xsig(&accumulator, &argSignif);
+  shr_Xsig(&accumulator, 6 - adj);
+
+  mul32_Xsig(&arg_signif, leadterm);
+  add_two_Xsig(&accumulator, &arg_signif, &exponent);
+
+  *expon = exponent + 1;
+  accum_result->lsw = accumulator.lsw;
+  accum_result->midw = accumulator.midw;
+  accum_result->msw = accumulator.msw;
+
+}