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