x86/math-emu/poly_atan.c

0001 // SPDX-License-Identifier: GPL-2.0
0002 /*---------------------------------------------------------------------------+
0003  |  poly_atan.c                                                              |
0004  |                                                                           |
0005  | Compute the arctan of a FPU_REG, using a polynomial approximation.        |
0006  |                                                                           |
0007  | Copyright (C) 1992,1993,1994,1997                                         |
0008  |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
0009  |                  E-mail   billm@suburbia.net                              |
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 "status_w.h"
0019 #include "control_w.h"
0020 #include "poly.h"
0021
0022 #define HIPOWERon   6   /* odd poly, negative terms */
0023 static const unsigned long long oddnegterms[HIPOWERon] = {
0024     0x0000000000000000LL,   /* Dummy (not for - 1.0) */
0025     0x015328437f756467LL,
0026     0x0005dda27b73dec6LL,
0027     0x0000226bf2bfb91aLL,
0028     0x000000ccc439c5f7LL,
0029     0x0000000355438407LL
0030 };
0031
0032 #define HIPOWERop   6   /* odd poly, positive terms */
0033 static const unsigned long long oddplterms[HIPOWERop] = {
0034 /*  0xaaaaaaaaaaaaaaabLL,  transferred to fixedpterm[] */
0035     0x0db55a71875c9ac2LL,
0036     0x0029fce2d67880b0LL,
0037     0x0000dfd3908b4596LL,
0038     0x00000550fd61dab4LL,
0039     0x0000001c9422b3f9LL,
0040     0x000000003e3301e1LL
0041 };
0042
0043 static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
0044
0045 static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
0046
0047 static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
0048
0049 /*--- poly_atan() -----------------------------------------------------------+
0050  |                                                                           |
0051  +---------------------------------------------------------------------------*/
0052 void poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
0053            FPU_REG *st1_ptr, u_char st1_tag)
0054 {
0055     u_char transformed, inverted, sign1, sign2;
0056     int exponent;
0057     long int dummy_exp;
0058     Xsig accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq;
0059     u_char tag;
0060
0061     sign1 = getsign(st0_ptr);
0062     sign2 = getsign(st1_ptr);
0063     if (st0_tag == TAG_Valid) {
0064         exponent = exponent(st0_ptr);
0065     } else {
0066         /* This gives non-compatible stack contents... */
0067         FPU_to_exp16(st0_ptr, st0_ptr);
0068         exponent = exponent16(st0_ptr);
0069     }
0070     if (st1_tag == TAG_Valid) {
0071         exponent -= exponent(st1_ptr);
0072     } else {
0073         /* This gives non-compatible stack contents... */
0074         FPU_to_exp16(st1_ptr, st1_ptr);
0075         exponent -= exponent16(st1_ptr);
0076     }
0077
0078     if ((exponent < 0) || ((exponent == 0) &&
0079                    ((st0_ptr->sigh < st1_ptr->sigh) ||
0080                 ((st0_ptr->sigh == st1_ptr->sigh) &&
0081                  (st0_ptr->sigl < st1_ptr->sigl))))) {
0082         inverted = 1;
0083         Numer.lsw = Denom.lsw = 0;
0084         XSIG_LL(Numer) = significand(st0_ptr);
0085         XSIG_LL(Denom) = significand(st1_ptr);
0086     } else {
0087         inverted = 0;
0088         exponent = -exponent;
0089         Numer.lsw = Denom.lsw = 0;
0090         XSIG_LL(Numer) = significand(st1_ptr);
0091         XSIG_LL(Denom) = significand(st0_ptr);
0092     }
0093     div_Xsig(&Numer, &Denom, &argSignif);
0094     exponent += norm_Xsig(&argSignif);
0095
0096     if ((exponent >= -1)
0097         || ((exponent == -2) && (argSignif.msw > 0xd413ccd0))) {
0098         /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
0099         /* Convert the argument by an identity for atan */
0100         transformed = 1;
0101
0102         if (exponent >= 0) {
0103 #ifdef PARANOID
0104             if (!((exponent == 0) &&
0105                   (argSignif.lsw == 0) && (argSignif.midw == 0) &&
0106                   (argSignif.msw == 0x80000000))) {
0107                 EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic error */
0108                 return;
0109             }
0110 #endif /* PARANOID */
0111             argSignif.msw = 0;  /* Make the transformed arg -> 0.0 */
0112         } else {
0113             Numer.lsw = Denom.lsw = argSignif.lsw;
0114             XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
0115
0116             if (exponent < -1)
0117                 shr_Xsig(&Numer, -1 - exponent);
0118             negate_Xsig(&Numer);
0119
0120             shr_Xsig(&Denom, -exponent);
0121             Denom.msw |= 0x80000000;
0122
0123             div_Xsig(&Numer, &Denom, &argSignif);
0124
0125             exponent = -1 + norm_Xsig(&argSignif);
0126         }
0127     } else {
0128         transformed = 0;
0129     }
0130
0131     argSq.lsw = argSignif.lsw;
0132     argSq.midw = argSignif.midw;
0133     argSq.msw = argSignif.msw;
0134     mul_Xsig_Xsig(&argSq, &argSq);
0135
0136     argSqSq.lsw = argSq.lsw;
0137     argSqSq.midw = argSq.midw;
0138     argSqSq.msw = argSq.msw;
0139     mul_Xsig_Xsig(&argSqSq, &argSqSq);
0140
0141     accumulatore.lsw = argSq.lsw;
0142     XSIG_LL(accumulatore) = XSIG_LL(argSq);
0143
0144     shr_Xsig(&argSq, 2 * (-1 - exponent - 1));
0145     shr_Xsig(&argSqSq, 4 * (-1 - exponent - 1));
0146
0147     /* Now have argSq etc with binary point at the left
0148        .1xxxxxxxx */
0149
0150     /* Do the basic fixed point polynomial evaluation */
0151     accumulator.msw = accumulator.midw = accumulator.lsw = 0;
0152     polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
0153             oddplterms, HIPOWERop - 1);
0154     mul64_Xsig(&accumulator, &XSIG_LL(argSq));
0155     negate_Xsig(&accumulator);
0156     polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms,
0157             HIPOWERon - 1);
0158     negate_Xsig(&accumulator);
0159     add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
0160
0161     mul64_Xsig(&accumulatore, &denomterm);
0162     shr_Xsig(&accumulatore, 1 + 2 * (-1 - exponent));
0163     accumulatore.msw |= 0x80000000;
0164
0165     div_Xsig(&accumulator, &accumulatore, &accumulator);
0166
0167     mul_Xsig_Xsig(&accumulator, &argSignif);
0168     mul_Xsig_Xsig(&accumulator, &argSq);
0169
0170     shr_Xsig(&accumulator, 3);
0171     negate_Xsig(&accumulator);
0172     add_Xsig_Xsig(&accumulator, &argSignif);
0173
0174     if (transformed) {
0175         /* compute pi/4 - accumulator */
0176         shr_Xsig(&accumulator, -1 - exponent);
0177         negate_Xsig(&accumulator);
0178         add_Xsig_Xsig(&accumulator, &pi_signif);
0179         exponent = -1;
0180     }
0181
0182     if (inverted) {
0183         /* compute pi/2 - accumulator */
0184         shr_Xsig(&accumulator, -exponent);
0185         negate_Xsig(&accumulator);
0186         add_Xsig_Xsig(&accumulator, &pi_signif);
0187         exponent = 0;
0188     }
0189
0190     if (sign1) {
0191         /* compute pi - accumulator */
0192         shr_Xsig(&accumulator, 1 - exponent);
0193         negate_Xsig(&accumulator);
0194         add_Xsig_Xsig(&accumulator, &pi_signif);
0195         exponent = 1;
0196     }
0197
0198     exponent += round_Xsig(&accumulator);
0199
0200     significand(st1_ptr) = XSIG_LL(accumulator);
0201     setexponent16(st1_ptr, exponent);
0202
0203     tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
0204     FPU_settagi(1, tag);
0205
0206     set_precision_flag_up();    /* We do not really know if up or down,
0207                        use this as the default. */
0208
0209 }