#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
-} ;
+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] =
-{
+static const unsigned long long oddplterms[HIPOWERop] = {
/* 0xaaaaaaaaaaaaaaabLL, transferred to fixedpterm[] */
- 0x0db55a71875c9ac2LL,
- 0x0029fce2d67880b0LL,
- 0x0000dfd3908b4596LL,
- 0x00000550fd61dab4LL,
- 0x0000001c9422b3f9LL,
- 0x000000003e3301e1LL
+ 0x0db55a71875c9ac2LL,
+ 0x0029fce2d67880b0LL,
+ 0x0000dfd3908b4596LL,
+ 0x00000550fd61dab4LL,
+ 0x0000001c9422b3f9LL,
+ 0x000000003e3301e1LL
};
static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
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)
+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 )
- {
+ 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;
- }
+ 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 */
+ 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;
}
- 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);
+
+ if (sign1) {
+ /* compute pi - accumulator */
+ shr_Xsig(&accumulator, 1 - exponent);
+ negate_Xsig(&accumulator);
+ add_Xsig_Xsig(&accumulator, &pi_signif);
+ exponent = 1;
}
- }
- 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. */
+
+ 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. */
}
git.openpandora.org / pandora-kernel.git / blobdiff