Quelle  poly_atan.c   Sprache: C

// SPDX-License-Identifier: GPL-2.0
/*---------------------------------------------------------------------------+
 |  poly_atan.c                                                              |
 |                                                                           |
 | Compute the arctan of a FPU_REG, using a polynomial approximation.        |
 |                                                                           |
 | Copyright (C) 1992,1993,1994,1997                                         |
 |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
 |                  E-mail   billm@suburbia.net                              |
 |                                                                           |
 |                                                                           |
 +---------------------------------------------------------------------------*/

#include "exception.h"
#include "reg_constant.h"
#include "fpu_emu.h"
#include "fpu_system.h"
#include "status_w.h"
#include "control_w.h"
#include "poly.h"

#define HIPOWERon 6 /* odd poly, negative terms */
static const unsigned long long oddnegterms[HIPOWERon] = {
 0x0000000000000000LL, /* Dummy (not for - 1.0) */
 0x015328437f756467LL,
 0x0005dda27b73dec6LL,
 0x0000226bf2bfb91aLL,
 0x000000ccc439c5f7LL,
 0x0000000355438407LL
};

#define HIPOWERop 6 /* odd poly, positive terms */
static const unsigned long long oddplterms[HIPOWERop] = {
/*  0xaaaaaaaaaaaaaaabLL,  transferred to fixedpterm[] */
 0x0db55a71875c9ac2LL,
 0x0029fce2d67880b0LL,
 0x0000dfd3908b4596LL,
 0x00000550fd61dab4LL,
 0x0000001c9422b3f9LL,
 0x000000003e3301e1LL
};

static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;

static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);

static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);

/*--- poly_atan() -----------------------------------------------------------+
 |                                                                           |
 +---------------------------------------------------------------------------*/
void poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
        FPU_REG *st1_ptr, u_char st1_tag)
{
 u_char transformed, inverted, sign1, sign2;
 int exponent;
 long int dummy_exp;
 Xsig accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq;
 u_char tag;

 sign1 = getsign(st0_ptr);
 sign2 = getsign(st1_ptr);
 if (st0_tag == TAG_Valid) {
  exponent = exponent(st0_ptr);
 } else {
  /* This gives non-compatible stack contents... */
  FPU_to_exp16(st0_ptr, st0_ptr);
  exponent = exponent16(st0_ptr);
 }
 if (st1_tag == TAG_Valid) {
  exponent -= exponent(st1_ptr);
 } else {
  /* This gives non-compatible stack contents... */
  FPU_to_exp16(st1_ptr, st1_ptr);
  exponent -= exponent16(st1_ptr);
 }

 if ((exponent < 0) || ((exponent == 0) &&
          ((st0_ptr->sigh < st1_ptr->sigh) ||
    ((st0_ptr->sigh == st1_ptr->sigh) &&
     (st0_ptr->sigl < st1_ptr->sigl))))) {
  inverted = 1;
  Numer.lsw = Denom.lsw = 0;
  XSIG_LL(Numer) = significand(st0_ptr);
  XSIG_LL(Denom) = significand(st1_ptr);
 } else {
  inverted = 0;
  exponent = -exponent;
  Numer.lsw = Denom.lsw = 0;
  XSIG_LL(Numer) = significand(st1_ptr);
  XSIG_LL(Denom) = significand(st0_ptr);
 }
 div_Xsig(&Numer, &Denom, &argSignif);
 exponent += norm_Xsig(&argSignif);

 if ((exponent >= -1)
     || ((exponent == -2) && (argSignif.msw > 0xd413ccd0))) {
  /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
  /* Convert the argument by an identity for atan */
  transformed = 1;

  if (exponent >= 0) {
#ifdef PARANOID
   if (!((exponent == 0) &&
         (argSignif.lsw == 0) && (argSignif.midw == 0) &&
         (argSignif.msw == 0x80000000))) {
    EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic error */
    return;
   }
#endif /* PARANOID */
   argSignif.msw = 0; /* Make the transformed arg -> 0.0 */
  } else {
   Numer.lsw = Denom.lsw = argSignif.lsw;
   XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);

   if (exponent < -1)
    shr_Xsig(&Numer, -1 - exponent);
   negate_Xsig(&Numer);

   shr_Xsig(&Denom, -exponent);
   Denom.msw |= 0x80000000;

   div_Xsig(&Numer, &Denom, &argSignif);

   exponent = -1 + norm_Xsig(&argSignif);
  }
 } else {
  transformed = 0;
 }

 argSq.lsw = argSignif.lsw;
 argSq.midw = argSignif.midw;
 argSq.msw = argSignif.msw;
 mul_Xsig_Xsig(&argSq, &argSq);

 argSqSq.lsw = argSq.lsw;
 argSqSq.midw = argSq.midw;
 argSqSq.msw = argSq.msw;
 mul_Xsig_Xsig(&argSqSq, &argSqSq);

 accumulatore.lsw = argSq.lsw;
 XSIG_LL(accumulatore) = XSIG_LL(argSq);

 shr_Xsig(&argSq, 2 * (-1 - exponent - 1));
 shr_Xsig(&argSqSq, 4 * (-1 - exponent - 1));

 /* Now have argSq etc with binary point at the left
   .1xxxxxxxx */

 /* Do the basic fixed point polynomial evaluation */
 accumulator.msw = accumulator.midw = accumulator.lsw = 0;
 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
   oddplterms, HIPOWERop - 1);
 mul64_Xsig(&accumulator, &XSIG_LL(argSq));
 negate_Xsig(&accumulator);
 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms,
   HIPOWERon - 1);
 negate_Xsig(&accumulator);
 add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);

 mul64_Xsig(&accumulatore, &denomterm);
 shr_Xsig(&accumulatore, 1 + 2 * (-1 - exponent));
 accumulatore.msw |= 0x80000000;

 div_Xsig(&accumulator, &accumulatore, &accumulator);

 mul_Xsig_Xsig(&accumulator, &argSignif);
 mul_Xsig_Xsig(&accumulator, &argSq);

 shr_Xsig(&accumulator, 3);
 negate_Xsig(&accumulator);
 add_Xsig_Xsig(&accumulator, &argSignif);

 if (transformed) {
  /* compute pi/4 - accumulator */
  shr_Xsig(&accumulator, -1 - exponent);
  negate_Xsig(&accumulator);
  add_Xsig_Xsig(&accumulator, &pi_signif);
  exponent = -1;
 }

 if (inverted) {
  /* compute pi/2 - accumulator */
  shr_Xsig(&accumulator, -exponent);
  negate_Xsig(&accumulator);
  add_Xsig_Xsig(&accumulator, &pi_signif);
  exponent = 0;
 }

 if (sign1) {
  /* compute pi - accumulator */
  shr_Xsig(&accumulator, 1 - exponent);
  negate_Xsig(&accumulator);
  add_Xsig_Xsig(&accumulator, &pi_signif);
  exponent = 1;
 }

 exponent += round_Xsig(&accumulator);

 significand(st1_ptr) = XSIG_LL(accumulator);
 setexponent16(st1_ptr, exponent);

 tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
 FPU_settagi(1, tag);

 set_precision_flag_up(); /* We do not really know if up or down,
   use this as the default. */

}