arch/i386/math-emu/poly_tan.c

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