1 files changed, 222 insertions, 0 deletions

diff --git a/arch/i386/math-emu/poly_tan.c b/arch/i386/math-emu/poly_tan.c
new file mode 100644
index 000000000000..8df3e03b6e6f
--- /dev/null
+++ b/arch/i386/math-emu/poly_tan.c
@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------+
+ |  poly_tan.c                                                               |
+ |                                                                           |
+ | Compute the tan of a FPU_REG, using a polynomial approximation.           |
+ |                                                                           |
+ | Copyright (C) 1992,1993,1994,1997,1999                                    |
+ |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
+ |                       Australia.  E-mail   billm@melbpc.org.au            |
+ |                                                                           |
+ |                                                                           |
+ +---------------------------------------------------------------------------*/
+
+#include "exception.h"
+#include "reg_constant.h"
+#include "fpu_emu.h"
+#include "fpu_system.h"
+#include "control_w.h"
+#include "poly.h"
+
+
+#define HiPOWERop       3       /* odd poly, positive terms */
+static const unsigned long long oddplterm[HiPOWERop] =
+{
+  0x0000000000000000LL,
+  0x0051a1cf08fca228LL,
+  0x0000000071284ff7LL
+};
+
+#define HiPOWERon       2       /* odd poly, negative terms */
+static const unsigned long long oddnegterm[HiPOWERon] =
+{
+   0x1291a9a184244e80LL,
+   0x0000583245819c21LL
+};
+
+#define HiPOWERep       2       /* even poly, positive terms */
+static const unsigned long long evenplterm[HiPOWERep] =
+{
+  0x0e848884b539e888LL,
+  0x00003c7f18b887daLL
+};
+
+#define HiPOWERen       2       /* even poly, negative terms */
+static const unsigned long long evennegterm[HiPOWERen] =
+{
+  0xf1f0200fd51569ccLL,
+  0x003afb46105c4432LL
+};
+
+static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
+
+
+/*--- poly_tan() ------------------------------------------------------------+
+ |                                                                           |
+ +---------------------------------------------------------------------------*/
+void    poly_tan(FPU_REG *st0_ptr)
+{
+  long int              exponent;
+  int                   invert;
+  Xsig                  argSq, argSqSq, accumulatoro, accumulatore, accum,
+                        argSignif, fix_up;
+  unsigned long         adj;
+
+  exponent = exponent(st0_ptr);
+
+#ifdef PARANOID
+  if ( signnegative(st0_ptr) )  /* Can't hack a number < 0.0 */
+    { arith_invalid(0); return; }  /* Need a positive number */
+#endif /* PARANOID */
+
+  /* Split the problem into two domains, smaller and larger than pi/4 */
+  if ( (exponent == 0) || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2)) )
+    {
+      /* The argument is greater than (approx) pi/4 */
+      invert = 1;
+      accum.lsw = 0;
+      XSIG_LL(accum) = significand(st0_ptr);
+ 
+      if ( exponent == 0 )
+        {
+          /* The argument is >= 1.0 */
+          /* Put the binary point at the left. */
+          XSIG_LL(accum) <<= 1;
+        }
+      /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
+      XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
+      /* This is a special case which arises due to rounding. */
+      if ( XSIG_LL(accum) == 0xffffffffffffffffLL )
+        {
+          FPU_settag0(TAG_Valid);
+          significand(st0_ptr) = 0x8a51e04daabda360LL;
+          setexponent16(st0_ptr, (0x41 + EXTENDED_Ebias) | SIGN_Negative);
+          return;
+        }
+
+      argSignif.lsw = accum.lsw;
+      XSIG_LL(argSignif) = XSIG_LL(accum);
+      exponent = -1 + norm_Xsig(&argSignif);
+    }
+  else
+    {
+      invert = 0;
+      argSignif.lsw = 0;
+      XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr);
+ 
+      if ( exponent < -1 )
+        {
+          /* shift the argument right by the required places */
+          if ( FPU_shrx(&XSIG_LL(accum), -1-exponent) >= 0x80000000U )
+            XSIG_LL(accum) ++;  /* round up */
+        }
+    }
+
+  XSIG_LL(argSq) = XSIG_LL(accum); argSq.lsw = accum.lsw;
+  mul_Xsig_Xsig(&argSq, &argSq);
+  XSIG_LL(argSqSq) = XSIG_LL(argSq); argSqSq.lsw = argSq.lsw;
+  mul_Xsig_Xsig(&argSqSq, &argSqSq);
+
+  /* Compute the negative terms for the numerator polynomial */
+  accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
+  polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, HiPOWERon-1);
+  mul_Xsig_Xsig(&accumulatoro, &argSq);
+  negate_Xsig(&accumulatoro);
+  /* Add the positive terms */
+  polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, HiPOWERop-1);
+
+  
+  /* Compute the positive terms for the denominator polynomial */
+  accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
+  polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, HiPOWERep-1);
+  mul_Xsig_Xsig(&accumulatore, &argSq);
+  negate_Xsig(&accumulatore);
+  /* Add the negative terms */
+  polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, HiPOWERen-1);
+  /* Multiply by arg^2 */
+  mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
+  mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
+  /* de-normalize and divide by 2 */
+  shr_Xsig(&accumulatore, -2*(1+exponent) + 1);
+  negate_Xsig(&accumulatore);      /* This does 1 - accumulator */
+
+  /* Now find the ratio. */
+  if ( accumulatore.msw == 0 )
+    {
+      /* accumulatoro must contain 1.0 here, (actually, 0) but it
+         really doesn't matter what value we use because it will
+         have negligible effect in later calculations
+         */
+      XSIG_LL(accum) = 0x8000000000000000LL;
+      accum.lsw = 0;
+    }
+  else
+    {
+      div_Xsig(&accumulatoro, &accumulatore, &accum);
+    }
+
+  /* Multiply by 1/3 * arg^3 */
+  mul64_Xsig(&accum, &XSIG_LL(argSignif));
+  mul64_Xsig(&accum, &XSIG_LL(argSignif));
+  mul64_Xsig(&accum, &XSIG_LL(argSignif));
+  mul64_Xsig(&accum, &twothirds);
+  shr_Xsig(&accum, -2*(exponent+1));
+
+  /* tan(arg) = arg + accum */
+  add_two_Xsig(&accum, &argSignif, &exponent);
+
+  if ( invert )
+    {
+      /* We now have the value of tan(pi_2 - arg) where pi_2 is an
+         approximation for pi/2
+         */
+      /* The next step is to fix the answer to compensate for the
+         error due to the approximation used for pi/2
+         */
+
+      /* This is (approx) delta, the error in our approx for pi/2
+         (see above). It has an exponent of -65
+         */
+      XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
+      fix_up.lsw = 0;
+
+      if ( exponent == 0 )
+        adj = 0xffffffff;   /* We want approx 1.0 here, but
+                               this is close enough. */
+      else if ( exponent > -30 )
+        {
+          adj = accum.msw >> -(exponent+1);      /* tan */
+          adj = mul_32_32(adj, adj);             /* tan^2 */
+        }
+      else
+        adj = 0;
+      adj = mul_32_32(0x898cc517, adj);          /* delta * tan^2 */
+
+      fix_up.msw += adj;
+      if ( !(fix_up.msw & 0x80000000) )   /* did fix_up overflow ? */
+        {
+          /* Yes, we need to add an msb */
+          shr_Xsig(&fix_up, 1);
+          fix_up.msw |= 0x80000000;
+          shr_Xsig(&fix_up, 64 + exponent);
+        }
+      else
+        shr_Xsig(&fix_up, 65 + exponent);
+
+      add_two_Xsig(&accum, &fix_up, &exponent);
+
+      /* accum now contains tan(pi/2 - arg).
+         Use tan(arg) = 1.0 / tan(pi/2 - arg)
+         */
+      accumulatoro.lsw = accumulatoro.midw = 0;
+      accumulatoro.msw = 0x80000000;
+      div_Xsig(&accumulatoro, &accum, &accum);
+      exponent = - exponent - 1;
+    }
+
+  /* Transfer the result */
+  round_Xsig(&accum);
+  FPU_settag0(TAG_Valid);
+  significand(st0_ptr) = XSIG_LL(accum);
+  setexponent16(st0_ptr, exponent + EXTENDED_Ebias);  /* Result is positive. */
+
+}