x86/math-emu/poly_tan.c

0001 // SPDX-License-Identifier: GPL-2.0
0002 /*---------------------------------------------------------------------------+
0003  |  poly_tan.c                                                               |
0004  |                                                                           |
0005  | Compute the tan of a FPU_REG, using a polynomial approximation.           |
0006  |                                                                           |
0007  | Copyright (C) 1992,1993,1994,1997,1999                                    |
0008  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
0009  |                       Australia.  E-mail   billm@melbpc.org.au            |
0010  |                                                                           |
0011  |                                                                           |
0012  +---------------------------------------------------------------------------*/
0013
0014 #include "exception.h"
0015 #include "reg_constant.h"
0016 #include "fpu_emu.h"
0017 #include "fpu_system.h"
0018 #include "control_w.h"
0019 #include "poly.h"
0020
0021 #define HiPOWERop   3   /* odd poly, positive terms */
0022 static const unsigned long long oddplterm[HiPOWERop] = {
0023     0x0000000000000000LL,
0024     0x0051a1cf08fca228LL,
0025     0x0000000071284ff7LL
0026 };
0027
0028 #define HiPOWERon   2   /* odd poly, negative terms */
0029 static const unsigned long long oddnegterm[HiPOWERon] = {
0030     0x1291a9a184244e80LL,
0031     0x0000583245819c21LL
0032 };
0033
0034 #define HiPOWERep   2   /* even poly, positive terms */
0035 static const unsigned long long evenplterm[HiPOWERep] = {
0036     0x0e848884b539e888LL,
0037     0x00003c7f18b887daLL
0038 };
0039
0040 #define HiPOWERen   2   /* even poly, negative terms */
0041 static const unsigned long long evennegterm[HiPOWERen] = {
0042     0xf1f0200fd51569ccLL,
0043     0x003afb46105c4432LL
0044 };
0045
0046 static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
0047
0048 /*--- poly_tan() ------------------------------------------------------------+
0049  |                                                                           |
0050  +---------------------------------------------------------------------------*/
0051 void poly_tan(FPU_REG *st0_ptr)
0052 {
0053     long int exponent;
0054     int invert;
0055     Xsig argSq, argSqSq, accumulatoro, accumulatore, accum,
0056         argSignif, fix_up;
0057     unsigned long adj;
0058
0059     exponent = exponent(st0_ptr);
0060
0061 #ifdef PARANOID
0062     if (signnegative(st0_ptr)) {    /* Can't hack a number < 0.0 */
0063         arith_invalid(0);
0064         return;
0065     }           /* Need a positive number */
0066 #endif /* PARANOID */
0067
0068     /* Split the problem into two domains, smaller and larger than pi/4 */
0069     if ((exponent == 0)
0070         || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2))) {
0071         /* The argument is greater than (approx) pi/4 */
0072         invert = 1;
0073         accum.lsw = 0;
0074         XSIG_LL(accum) = significand(st0_ptr);
0075
0076         if (exponent == 0) {
0077             /* The argument is >= 1.0 */
0078             /* Put the binary point at the left. */
0079             XSIG_LL(accum) <<= 1;
0080         }
0081         /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
0082         XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
0083         /* This is a special case which arises due to rounding. */
0084         if (XSIG_LL(accum) == 0xffffffffffffffffLL) {
0085             FPU_settag0(TAG_Valid);
0086             significand(st0_ptr) = 0x8a51e04daabda360LL;
0087             setexponent16(st0_ptr,
0088                       (0x41 + EXTENDED_Ebias) | SIGN_Negative);
0089             return;
0090         }
0091
0092         argSignif.lsw = accum.lsw;
0093         XSIG_LL(argSignif) = XSIG_LL(accum);
0094         exponent = -1 + norm_Xsig(&argSignif);
0095     } else {
0096         invert = 0;
0097         argSignif.lsw = 0;
0098         XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr);
0099
0100         if (exponent < -1) {
0101             /* shift the argument right by the required places */
0102             if (FPU_shrx(&XSIG_LL(accum), -1 - exponent) >=
0103                 0x80000000U)
0104                 XSIG_LL(accum)++;   /* round up */
0105         }
0106     }
0107
0108     XSIG_LL(argSq) = XSIG_LL(accum);
0109     argSq.lsw = accum.lsw;
0110     mul_Xsig_Xsig(&argSq, &argSq);
0111     XSIG_LL(argSqSq) = XSIG_LL(argSq);
0112     argSqSq.lsw = argSq.lsw;
0113     mul_Xsig_Xsig(&argSqSq, &argSqSq);
0114
0115     /* Compute the negative terms for the numerator polynomial */
0116     accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
0117     polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm,
0118             HiPOWERon - 1);
0119     mul_Xsig_Xsig(&accumulatoro, &argSq);
0120     negate_Xsig(&accumulatoro);
0121     /* Add the positive terms */
0122     polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm,
0123             HiPOWERop - 1);
0124
0125     /* Compute the positive terms for the denominator polynomial */
0126     accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
0127     polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm,
0128             HiPOWERep - 1);
0129     mul_Xsig_Xsig(&accumulatore, &argSq);
0130     negate_Xsig(&accumulatore);
0131     /* Add the negative terms */
0132     polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm,
0133             HiPOWERen - 1);
0134     /* Multiply by arg^2 */
0135     mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
0136     mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
0137     /* de-normalize and divide by 2 */
0138     shr_Xsig(&accumulatore, -2 * (1 + exponent) + 1);
0139     negate_Xsig(&accumulatore); /* This does 1 - accumulator */
0140
0141     /* Now find the ratio. */
0142     if (accumulatore.msw == 0) {
0143         /* accumulatoro must contain 1.0 here, (actually, 0) but it
0144            really doesn't matter what value we use because it will
0145            have negligible effect in later calculations
0146          */
0147         XSIG_LL(accum) = 0x8000000000000000LL;
0148         accum.lsw = 0;
0149     } else {
0150         div_Xsig(&accumulatoro, &accumulatore, &accum);
0151     }
0152
0153     /* Multiply by 1/3 * arg^3 */
0154     mul64_Xsig(&accum, &XSIG_LL(argSignif));
0155     mul64_Xsig(&accum, &XSIG_LL(argSignif));
0156     mul64_Xsig(&accum, &XSIG_LL(argSignif));
0157     mul64_Xsig(&accum, &twothirds);
0158     shr_Xsig(&accum, -2 * (exponent + 1));
0159
0160     /* tan(arg) = arg + accum */
0161     add_two_Xsig(&accum, &argSignif, &exponent);
0162
0163     if (invert) {
0164         /* We now have the value of tan(pi_2 - arg) where pi_2 is an
0165            approximation for pi/2
0166          */
0167         /* The next step is to fix the answer to compensate for the
0168            error due to the approximation used for pi/2
0169          */
0170
0171         /* This is (approx) delta, the error in our approx for pi/2
0172            (see above). It has an exponent of -65
0173          */
0174         XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
0175         fix_up.lsw = 0;
0176
0177         if (exponent == 0)
0178             adj = 0xffffffff;   /* We want approx 1.0 here, but
0179                            this is close enough. */
0180         else if (exponent > -30) {
0181             adj = accum.msw >> -(exponent + 1); /* tan */
0182             adj = mul_32_32(adj, adj);  /* tan^2 */
0183         } else
0184             adj = 0;
0185         adj = mul_32_32(0x898cc517, adj);   /* delta * tan^2 */
0186
0187         fix_up.msw += adj;
0188         if (!(fix_up.msw & 0x80000000)) {   /* did fix_up overflow ? */
0189             /* Yes, we need to add an msb */
0190             shr_Xsig(&fix_up, 1);
0191             fix_up.msw |= 0x80000000;
0192             shr_Xsig(&fix_up, 64 + exponent);
0193         } else
0194             shr_Xsig(&fix_up, 65 + exponent);
0195
0196         add_two_Xsig(&accum, &fix_up, &exponent);
0197
0198         /* accum now contains tan(pi/2 - arg).
0199            Use tan(arg) = 1.0 / tan(pi/2 - arg)
0200          */
0201         accumulatoro.lsw = accumulatoro.midw = 0;
0202         accumulatoro.msw = 0x80000000;
0203         div_Xsig(&accumulatoro, &accum, &accum);
0204         exponent = -exponent - 1;
0205     }
0206
0207     /* Transfer the result */
0208     round_Xsig(&accum);
0209     FPU_settag0(TAG_Valid);
0210     significand(st0_ptr) = XSIG_LL(accum);
0211     setexponent16(st0_ptr, exponent + EXTENDED_Ebias);  /* Result is positive. */
0212
0213 }