arch/x86/math-emu/poly_atan.c

   1 // SPDX-License-Identifier: GPL-2.0
   2 /*---------------------------------------------------------------------------+
   3  |  poly_atan.c                                                              |
   4  |                                                                           |
   5  | Compute the arctan of a FPU_REG, using a polynomial approximation.        |
   6  |                                                                           |
   7  | Copyright (C) 1992,1993,1994,1997                                         |
   8  |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
   9  |                  E-mail   billm@suburbia.net                              |
  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 "status_w.h"
  19 #include "control_w.h"
  20 #include "poly.h"
  21
  22 #define HIPOWERon       6       /* odd poly, negative terms */
  23 static const unsigned long long oddnegterms[HIPOWERon] = {
  24         0x0000000000000000LL,   /* Dummy (not for - 1.0) */
  25         0x015328437f756467LL,
  26         0x0005dda27b73dec6LL,
  27         0x0000226bf2bfb91aLL,
  28         0x000000ccc439c5f7LL,
  29         0x0000000355438407LL
  30 };
  31
  32 #define HIPOWERop       6       /* odd poly, positive terms */
  33 static const unsigned long long oddplterms[HIPOWERop] = {
  34 /*  0xaaaaaaaaaaaaaaabLL,  transferred to fixedpterm[] */
  35         0x0db55a71875c9ac2LL,
  36         0x0029fce2d67880b0LL,
  37         0x0000dfd3908b4596LL,
  38         0x00000550fd61dab4LL,
  39         0x0000001c9422b3f9LL,
  40         0x000000003e3301e1LL
  41 };
  42
  43 static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
  44
  45 static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
  46
  47 static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
  48
  49 /*--- poly_atan() -----------------------------------------------------------+
  50  |                                                                           |
  51  +---------------------------------------------------------------------------*/
  52 void poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
  53                FPU_REG *st1_ptr, u_char st1_tag)
  54 {
  55         u_char transformed, inverted, sign1, sign2;
  56         int exponent;
  57         long int dummy_exp;
  58         Xsig accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq;
  59         u_char tag;
  60
  61         sign1 = getsign(st0_ptr);
  62         sign2 = getsign(st1_ptr);
  63         if (st0_tag == TAG_Valid) {
  64                 exponent = exponent(st0_ptr);
  65         } else {
  66                 /* This gives non-compatible stack contents... */
  67                 FPU_to_exp16(st0_ptr, st0_ptr);
  68                 exponent = exponent16(st0_ptr);
  69         }
  70         if (st1_tag == TAG_Valid) {
  71                 exponent -= exponent(st1_ptr);
  72         } else {
  73                 /* This gives non-compatible stack contents... */
  74                 FPU_to_exp16(st1_ptr, st1_ptr);
  75                 exponent -= exponent16(st1_ptr);
  76         }
  77
  78         if ((exponent < 0) || ((exponent == 0) &&
  79                                ((st0_ptr->sigh < st1_ptr->sigh) ||
  80                                 ((st0_ptr->sigh == st1_ptr->sigh) &&
  81                                  (st0_ptr->sigl < st1_ptr->sigl))))) {
  82                 inverted = 1;
  83                 Numer.lsw = Denom.lsw = 0;
  84                 XSIG_LL(Numer) = significand(st0_ptr);
  85                 XSIG_LL(Denom) = significand(st1_ptr);
  86         } else {
  87                 inverted = 0;
  88                 exponent = -exponent;
  89                 Numer.lsw = Denom.lsw = 0;
  90                 XSIG_LL(Numer) = significand(st1_ptr);
  91                 XSIG_LL(Denom) = significand(st0_ptr);
  92         }
  93         div_Xsig(&Numer, &Denom, &argSignif);
  94         exponent += norm_Xsig(&argSignif);
  95
  96         if ((exponent >= -1)
  97             || ((exponent == -2) && (argSignif.msw > 0xd413ccd0))) {
  98                 /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
  99                 /* Convert the argument by an identity for atan */
 100                 transformed = 1;
 101
 102                 if (exponent >= 0) {
 103 #ifdef PARANOID
 104                         if (!((exponent == 0) &&
 105                               (argSignif.lsw == 0) && (argSignif.midw == 0) &&
 106                               (argSignif.msw == 0x80000000))) {
 107                                 EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic error */
 108                                 return;
 109                         }
 110 #endif /* PARANOID */
 111                         argSignif.msw = 0;      /* Make the transformed arg -> 0.0 */
 112                 } else {
 113                         Numer.lsw = Denom.lsw = argSignif.lsw;
 114                         XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
 115
 116                         if (exponent < -1)
 117                                 shr_Xsig(&Numer, -1 - exponent);
 118                         negate_Xsig(&Numer);
 119
 120                         shr_Xsig(&Denom, -exponent);
 121                         Denom.msw |= 0x80000000;
 122
 123                         div_Xsig(&Numer, &Denom, &argSignif);
 124
 125                         exponent = -1 + norm_Xsig(&argSignif);
 126                 }
 127         } else {
 128                 transformed = 0;
 129         }
 130
 131         argSq.lsw = argSignif.lsw;
 132         argSq.midw = argSignif.midw;
 133         argSq.msw = argSignif.msw;
 134         mul_Xsig_Xsig(&argSq, &argSq);
 135
 136         argSqSq.lsw = argSq.lsw;
 137         argSqSq.midw = argSq.midw;
 138         argSqSq.msw = argSq.msw;
 139         mul_Xsig_Xsig(&argSqSq, &argSqSq);
 140
 141         accumulatore.lsw = argSq.lsw;
 142         XSIG_LL(accumulatore) = XSIG_LL(argSq);
 143
 144         shr_Xsig(&argSq, 2 * (-1 - exponent - 1));
 145         shr_Xsig(&argSqSq, 4 * (-1 - exponent - 1));
 146
 147         /* Now have argSq etc with binary point at the left
 148            .1xxxxxxxx */
 149
 150         /* Do the basic fixed point polynomial evaluation */
 151         accumulator.msw = accumulator.midw = accumulator.lsw = 0;
 152         polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
 153                         oddplterms, HIPOWERop - 1);
 154         mul64_Xsig(&accumulator, &XSIG_LL(argSq));
 155         negate_Xsig(&accumulator);
 156         polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms,
 157                         HIPOWERon - 1);
 158         negate_Xsig(&accumulator);
 159         add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
 160
 161         mul64_Xsig(&accumulatore, &denomterm);
 162         shr_Xsig(&accumulatore, 1 + 2 * (-1 - exponent));
 163         accumulatore.msw |= 0x80000000;
 164
 165         div_Xsig(&accumulator, &accumulatore, &accumulator);
 166
 167         mul_Xsig_Xsig(&accumulator, &argSignif);
 168         mul_Xsig_Xsig(&accumulator, &argSq);
 169
 170         shr_Xsig(&accumulator, 3);
 171         negate_Xsig(&accumulator);
 172         add_Xsig_Xsig(&accumulator, &argSignif);
 173
 174         if (transformed) {
 175                 /* compute pi/4 - accumulator */
 176                 shr_Xsig(&accumulator, -1 - exponent);
 177                 negate_Xsig(&accumulator);
 178                 add_Xsig_Xsig(&accumulator, &pi_signif);
 179                 exponent = -1;
 180         }
 181
 182         if (inverted) {
 183                 /* compute pi/2 - accumulator */
 184                 shr_Xsig(&accumulator, -exponent);
 185                 negate_Xsig(&accumulator);
 186                 add_Xsig_Xsig(&accumulator, &pi_signif);
 187                 exponent = 0;
 188         }
 189
 190         if (sign1) {
 191                 /* compute pi - accumulator */
 192                 shr_Xsig(&accumulator, 1 - exponent);
 193                 negate_Xsig(&accumulator);
 194                 add_Xsig_Xsig(&accumulator, &pi_signif);
 195                 exponent = 1;
 196         }
 197
 198         exponent += round_Xsig(&accumulator);
 199
 200         significand(st1_ptr) = XSIG_LL(accumulator);
 201         setexponent16(st1_ptr, exponent);
 202
 203         tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
 204         FPU_settagi(1, tag);
 205
 206         set_precision_flag_up();        /* We do not really know if up or down,
 207                                            use this as the default. */
 208
 209 }