Subversion Repositories or1k

/*---------------------------------------------------------------------------+

 |  poly_tan.c                                                               |

 |                                                                           |

 | Compute the tan of a FPU_REG, using a polynomial approximation.           |

 |                                                                           |

 | Copyright (C) 1992,1993,1994                                              |

 |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |

 |                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au    |

 |                                                                           |

 |                                                                           |

 +---------------------------------------------------------------------------*/

#include "exception.h"

#include "reg_constant.h"

#include "fpu_emu.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 const *arg, FPU_REG *result)

{

  long int              exponent;

  int                   invert;

  Xsig                  argSq, argSqSq, accumulatoro, accumulatore, accum,

                        argSignif, fix_up;

  unsigned long         adj;

  exponent = arg->exp - EXP_BIAS;

#ifdef PARANOID

  if ( arg->sign != 0 )  /* Can't hack a number < 0.0 */

    { arith_invalid(result); return; }  /* Need a positive number */

#endif PARANOID

  /* Split the problem into two domains, smaller and larger than pi/4 */

  if ( (exponent == 0) || ((exponent == -1) && (arg->sigh > 0xc90fdaa2)) )

    {

      /* The argument is greater than (approx) pi/4 */

      invert = 1;

      accum.lsw = 0;

      XSIG_LL(accum) = significand(arg);

      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);

      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(arg);

      if ( exponent < -1 )

        {

          /* shift the argument right by the required places */

          if ( 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 */

          mul_32_32(adj, adj, &adj);           /* tan^2 */

        }

      else

        adj = 0;

      mul_32_32(0x898cc517, adj, &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);

  *(short *)&(result->sign) = 0;

  significand(result) = XSIG_LL(accum);

  result->exp = EXP_BIAS + exponent;

}