arch/x86/math-emu/poly_tan.c

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