root/arch/i386/math-emu/poly_tan.c

/* */

DEFINITIONS

1. poly_tan

   1 /*---------------------------------------------------------------------------+
   2  |  poly_tan.c                                                               |
   3  |                                                                           |
   4  | Compute the tan 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 #include "exception.h"
  14 #include "reg_constant.h"
  15 #include "fpu_emu.h"
  16 #include "control_w.h"
  17 #include "poly.h"
  18 
  19 
  20 #define HiPOWERop       3       /* odd poly, positive terms */
  21 static const unsigned long long oddplterm[HiPOWERop] =
  22 {
  23   0x0000000000000000LL,
  24   0x0051a1cf08fca228LL,
  25   0x0000000071284ff7LL
  26 };
  27 
  28 #define HiPOWERon       2       /* odd poly, negative terms */
  29 static const unsigned long long oddnegterm[HiPOWERon] =
  30 {
  31    0x1291a9a184244e80LL,
  32    0x0000583245819c21LL
  33 };
  34 
  35 #define HiPOWERep       2       /* even poly, positive terms */
  36 static const unsigned long long evenplterm[HiPOWERep] =
  37 {
  38   0x0e848884b539e888LL,
  39   0x00003c7f18b887daLL
  40 };
  41 
  42 #define HiPOWERen       2       /* even poly, negative terms */
  43 static const unsigned long long evennegterm[HiPOWERen] =
  44 {
  45   0xf1f0200fd51569ccLL,
  46   0x003afb46105c4432LL
  47 };
  48 
  49 static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
  50 
  51 
  52 /*--- poly_tan() ------------------------------------------------------------+
  53  |                                                                           |
  54  +---------------------------------------------------------------------------*/
  55 void    poly_tan(FPU_REG const *arg, FPU_REG *result)
     /*  */
  56 {
  57   long int              exponent;
  58   int                   invert;
  59   Xsig                  argSq, argSqSq, accumulatoro, accumulatore, accum,
  60                         argSignif, fix_up;
  61   unsigned long         adj;
  62 
  63   exponent = arg->exp - EXP_BIAS;
  64 
  65 #ifdef PARANOID
  66   if ( arg->sign != 0 ) /* Can't hack a number < 0.0 */
  67     { arith_invalid(result); return; }  /* Need a positive number */
  68 #endif PARANOID
  69 
  70   /* Split the problem into two domains, smaller and larger than pi/4 */
  71   if ( (exponent == 0) || ((exponent == -1) && (arg->sigh > 0xc90fdaa2)) )
  72     {
  73       /* The argument is greater than (approx) pi/4 */
  74       invert = 1;
  75       accum.lsw = 0;
  76       XSIG_LL(accum) = significand(arg);
  77  
  78       if ( exponent == 0 )
  79         {
  80           /* The argument is >= 1.0 */
  81           /* Put the binary point at the left. */
  82           XSIG_LL(accum) <<= 1;
  83         }
  84       /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
  85       XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
  86 
  87       argSignif.lsw = accum.lsw;
  88       XSIG_LL(argSignif) = XSIG_LL(accum);
  89       exponent = -1 + norm_Xsig(&argSignif);
  90     }
  91   else
  92     {
  93       invert = 0;
  94       argSignif.lsw = 0;
  95       XSIG_LL(accum) = XSIG_LL(argSignif) = significand(arg);
  96  
  97       if ( exponent < -1 )
  98         {
  99           /* shift the argument right by the required places */
 100           if ( shrx(&XSIG_LL(accum), -1-exponent) >= 0x80000000U )
 101             XSIG_LL(accum) ++;  /* round up */
 102         }
 103     }
 104 
 105   XSIG_LL(argSq) = XSIG_LL(accum); argSq.lsw = accum.lsw;
 106   mul_Xsig_Xsig(&argSq, &argSq);
 107   XSIG_LL(argSqSq) = XSIG_LL(argSq); argSqSq.lsw = argSq.lsw;
 108   mul_Xsig_Xsig(&argSqSq, &argSqSq);
 109 
 110   /* Compute the negative terms for the numerator polynomial */
 111   accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
 112   polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, HiPOWERon-1);
 113   mul_Xsig_Xsig(&accumulatoro, &argSq);
 114   negate_Xsig(&accumulatoro);
 115   /* Add the positive terms */
 116   polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, HiPOWERop-1);
 117 
 118   
 119   /* Compute the positive terms for the denominator polynomial */
 120   accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
 121   polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, HiPOWERep-1);
 122   mul_Xsig_Xsig(&accumulatore, &argSq);
 123   negate_Xsig(&accumulatore);
 124   /* Add the negative terms */
 125   polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, HiPOWERen-1);
 126   /* Multiply by arg^2 */
 127   mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
 128   mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
 129   /* de-normalize and divide by 2 */
 130   shr_Xsig(&accumulatore, -2*(1+exponent) + 1);
 131   negate_Xsig(&accumulatore);      /* This does 1 - accumulator */
 132 
 133   /* Now find the ratio. */
 134   if ( accumulatore.msw == 0 )
 135     {
 136       /* accumulatoro must contain 1.0 here, (actually, 0) but it
 137          really doesn't matter what value we use because it will
 138          have negligible effect in later calculations
 139          */
 140       XSIG_LL(accum) = 0x8000000000000000LL;
 141       accum.lsw = 0;
 142     }
 143   else
 144     {
 145       div_Xsig(&accumulatoro, &accumulatore, &accum);
 146     }
 147 
 148   /* Multiply by 1/3 * arg^3 */
 149   mul64_Xsig(&accum, &XSIG_LL(argSignif));
 150   mul64_Xsig(&accum, &XSIG_LL(argSignif));
 151   mul64_Xsig(&accum, &XSIG_LL(argSignif));
 152   mul64_Xsig(&accum, &twothirds);
 153   shr_Xsig(&accum, -2*(exponent+1));
 154 
 155   /* tan(arg) = arg + accum */
 156   add_two_Xsig(&accum, &argSignif, &exponent);
 157 
 158   if ( invert )
 159     {
 160       /* We now have the value of tan(pi_2 - arg) where pi_2 is an
 161          approximation for pi/2
 162          */
 163       /* The next step is to fix the answer to compensate for the
 164          error due to the approximation used for pi/2
 165          */
 166 
 167       /* This is (approx) delta, the error in our approx for pi/2
 168          (see above). It has an exponent of -65
 169          */
 170       XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
 171       fix_up.lsw = 0;
 172 
 173       if ( exponent == 0 )
 174         adj = 0xffffffff;   /* We want approx 1.0 here, but
 175                                this is close enough. */
 176       else if ( exponent > -30 )
 177         {
 178           adj = accum.msw >> -(exponent+1);      /* tan */
 179           mul_32_32(adj, adj, &adj);           /* tan^2 */
 180         }
 181       else
 182         adj = 0;
 183       mul_32_32(0x898cc517, adj, &adj);        /* delta * tan^2 */
 184 
 185       fix_up.msw += adj;
 186       if ( !(fix_up.msw & 0x80000000) )   /* did fix_up overflow ? */
 187         {
 188           /* Yes, we need to add an msb */
 189           shr_Xsig(&fix_up, 1);
 190           fix_up.msw |= 0x80000000;
 191           shr_Xsig(&fix_up, 64 + exponent);
 192         }
 193       else
 194         shr_Xsig(&fix_up, 65 + exponent);
 195 
 196       add_two_Xsig(&accum, &fix_up, &exponent);
 197 
 198       /* accum now contains tan(pi/2 - arg).
 199          Use tan(arg) = 1.0 / tan(pi/2 - arg)
 200          */
 201       accumulatoro.lsw = accumulatoro.midw = 0;
 202       accumulatoro.msw = 0x80000000;
 203       div_Xsig(&accumulatoro, &accum, &accum);
 204       exponent = - exponent - 1;
 205     }
 206 
 207   /* Transfer the result */
 208   round_Xsig(&accum);
 209   *(short *)&(result->sign) = 0;
 210   significand(result) = XSIG_LL(accum);
 211   result->exp = EXP_BIAS + exponent;
 212 
 213 }

/* */