#include "control_w.h"
#include "poly.h"
-
#define HiPOWERop 3 /* odd poly, positive terms */
-static const unsigned long long oddplterm[HiPOWERop] =
-{
- 0x0000000000000000LL,
- 0x0051a1cf08fca228LL,
- 0x0000000071284ff7LL
+static const unsigned long long oddplterm[HiPOWERop] = {
+ 0x0000000000000000LL,
+ 0x0051a1cf08fca228LL,
+ 0x0000000071284ff7LL
};
#define HiPOWERon 2 /* odd poly, negative terms */
-static const unsigned long long oddnegterm[HiPOWERon] =
-{
- 0x1291a9a184244e80LL,
- 0x0000583245819c21LL
+static const unsigned long long oddnegterm[HiPOWERon] = {
+ 0x1291a9a184244e80LL,
+ 0x0000583245819c21LL
};
#define HiPOWERep 2 /* even poly, positive terms */
-static const unsigned long long evenplterm[HiPOWERep] =
-{
- 0x0e848884b539e888LL,
- 0x00003c7f18b887daLL
+static const unsigned long long evenplterm[HiPOWERep] = {
+ 0x0e848884b539e888LL,
+ 0x00003c7f18b887daLL
};
#define HiPOWERen 2 /* even poly, negative terms */
-static const unsigned long long evennegterm[HiPOWERen] =
-{
- 0xf1f0200fd51569ccLL,
- 0x003afb46105c4432LL
+static const unsigned long long evennegterm[HiPOWERen] = {
+ 0xf1f0200fd51569ccLL,
+ 0x003afb46105c4432LL
};
static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
-
/*--- poly_tan() ------------------------------------------------------------+
| |
+---------------------------------------------------------------------------*/
-void poly_tan(FPU_REG *st0_ptr)
+void poly_tan(FPU_REG *st0_ptr)
{
- long int exponent;
- int invert;
- Xsig argSq, argSqSq, accumulatoro, accumulatore, accum,
- argSignif, fix_up;
- unsigned long adj;
+ long int exponent;
+ int invert;
+ Xsig argSq, argSqSq, accumulatoro, accumulatore, accum,
+ argSignif, fix_up;
+ unsigned long adj;
- exponent = exponent(st0_ptr);
+ exponent = exponent(st0_ptr);
#ifdef PARANOID
- if ( signnegative(st0_ptr) ) /* Can't hack a number < 0.0 */
- { arith_invalid(0); return; } /* Need a positive number */
+ if (signnegative(st0_ptr)) { /* Can't hack a number < 0.0 */
+ arith_invalid(0);
+ return;
+ } /* Need a positive number */
#endif /* PARANOID */
- /* Split the problem into two domains, smaller and larger than pi/4 */
- if ( (exponent == 0) || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2)) )
- {
- /* The argument is greater than (approx) pi/4 */
- invert = 1;
- accum.lsw = 0;
- XSIG_LL(accum) = significand(st0_ptr);
-
- if ( exponent == 0 )
- {
- /* The argument is >= 1.0 */
- /* Put the binary point at the left. */
- XSIG_LL(accum) <<= 1;
- }
- /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
- XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
- /* This is a special case which arises due to rounding. */
- if ( XSIG_LL(accum) == 0xffffffffffffffffLL )
- {
- FPU_settag0(TAG_Valid);
- significand(st0_ptr) = 0x8a51e04daabda360LL;
- setexponent16(st0_ptr, (0x41 + EXTENDED_Ebias) | SIGN_Negative);
- return;
+ /* Split the problem into two domains, smaller and larger than pi/4 */
+ if ((exponent == 0)
+ || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2))) {
+ /* The argument is greater than (approx) pi/4 */
+ invert = 1;
+ accum.lsw = 0;
+ XSIG_LL(accum) = significand(st0_ptr);
+
+ if (exponent == 0) {
+ /* The argument is >= 1.0 */
+ /* Put the binary point at the left. */
+ XSIG_LL(accum) <<= 1;
+ }
+ /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
+ XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
+ /* This is a special case which arises due to rounding. */
+ if (XSIG_LL(accum) == 0xffffffffffffffffLL) {
+ FPU_settag0(TAG_Valid);
+ significand(st0_ptr) = 0x8a51e04daabda360LL;
+ setexponent16(st0_ptr,
+ (0x41 + EXTENDED_Ebias) | SIGN_Negative);
+ return;
+ }
+
+ argSignif.lsw = accum.lsw;
+ XSIG_LL(argSignif) = XSIG_LL(accum);
+ exponent = -1 + norm_Xsig(&argSignif);
+ } else {
+ invert = 0;
+ argSignif.lsw = 0;
+ XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr);
+
+ if (exponent < -1) {
+ /* shift the argument right by the required places */
+ if (FPU_shrx(&XSIG_LL(accum), -1 - exponent) >=
+ 0x80000000U)
+ XSIG_LL(accum)++; /* round up */
+ }
}
- argSignif.lsw = accum.lsw;
- XSIG_LL(argSignif) = XSIG_LL(accum);
- exponent = -1 + norm_Xsig(&argSignif);
- }
- else
- {
- invert = 0;
- argSignif.lsw = 0;
- XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr);
-
- if ( exponent < -1 )
- {
- /* shift the argument right by the required places */
- if ( FPU_shrx(&XSIG_LL(accum), -1-exponent) >= 0x80000000U )
- XSIG_LL(accum) ++; /* round up */
- }
- }
-
- XSIG_LL(argSq) = XSIG_LL(accum); argSq.lsw = accum.lsw;
- mul_Xsig_Xsig(&argSq, &argSq);
- XSIG_LL(argSqSq) = XSIG_LL(argSq); argSqSq.lsw = argSq.lsw;
- mul_Xsig_Xsig(&argSqSq, &argSqSq);
-
- /* Compute the negative terms for the numerator polynomial */
- accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
- polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, HiPOWERon-1);
- mul_Xsig_Xsig(&accumulatoro, &argSq);
- negate_Xsig(&accumulatoro);
- /* Add the positive terms */
- polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, HiPOWERop-1);
-
-
- /* Compute the positive terms for the denominator polynomial */
- accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
- polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, HiPOWERep-1);
- mul_Xsig_Xsig(&accumulatore, &argSq);
- negate_Xsig(&accumulatore);
- /* Add the negative terms */
- polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, HiPOWERen-1);
- /* Multiply by arg^2 */
- mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
- mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
- /* de-normalize and divide by 2 */
- shr_Xsig(&accumulatore, -2*(1+exponent) + 1);
- negate_Xsig(&accumulatore); /* This does 1 - accumulator */
-
- /* Now find the ratio. */
- if ( accumulatore.msw == 0 )
- {
- /* accumulatoro must contain 1.0 here, (actually, 0) but it
- really doesn't matter what value we use because it will
- have negligible effect in later calculations
- */
- XSIG_LL(accum) = 0x8000000000000000LL;
- accum.lsw = 0;
- }
- else
- {
- div_Xsig(&accumulatoro, &accumulatore, &accum);
- }
-
- /* Multiply by 1/3 * arg^3 */
- mul64_Xsig(&accum, &XSIG_LL(argSignif));
- mul64_Xsig(&accum, &XSIG_LL(argSignif));
- mul64_Xsig(&accum, &XSIG_LL(argSignif));
- mul64_Xsig(&accum, &twothirds);
- shr_Xsig(&accum, -2*(exponent+1));
-
- /* tan(arg) = arg + accum */
- add_two_Xsig(&accum, &argSignif, &exponent);
-
- if ( invert )
- {
- /* We now have the value of tan(pi_2 - arg) where pi_2 is an
- approximation for pi/2
- */
- /* The next step is to fix the answer to compensate for the
- error due to the approximation used for pi/2
- */
-
- /* This is (approx) delta, the error in our approx for pi/2
- (see above). It has an exponent of -65
- */
- XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
- fix_up.lsw = 0;
-
- if ( exponent == 0 )
- adj = 0xffffffff; /* We want approx 1.0 here, but
- this is close enough. */
- else if ( exponent > -30 )
- {
- adj = accum.msw >> -(exponent+1); /* tan */
- adj = mul_32_32(adj, adj); /* tan^2 */
+ XSIG_LL(argSq) = XSIG_LL(accum);
+ argSq.lsw = accum.lsw;
+ mul_Xsig_Xsig(&argSq, &argSq);
+ XSIG_LL(argSqSq) = XSIG_LL(argSq);
+ argSqSq.lsw = argSq.lsw;
+ mul_Xsig_Xsig(&argSqSq, &argSqSq);
+
+ /* Compute the negative terms for the numerator polynomial */
+ accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
+ polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm,
+ HiPOWERon - 1);
+ mul_Xsig_Xsig(&accumulatoro, &argSq);
+ negate_Xsig(&accumulatoro);
+ /* Add the positive terms */
+ polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm,
+ HiPOWERop - 1);
+
+ /* Compute the positive terms for the denominator polynomial */
+ accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
+ polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm,
+ HiPOWERep - 1);
+ mul_Xsig_Xsig(&accumulatore, &argSq);
+ negate_Xsig(&accumulatore);
+ /* Add the negative terms */
+ polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm,
+ HiPOWERen - 1);
+ /* Multiply by arg^2 */
+ mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
+ mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
+ /* de-normalize and divide by 2 */
+ shr_Xsig(&accumulatore, -2 * (1 + exponent) + 1);
+ negate_Xsig(&accumulatore); /* This does 1 - accumulator */
+
+ /* Now find the ratio. */
+ if (accumulatore.msw == 0) {
+ /* accumulatoro must contain 1.0 here, (actually, 0) but it
+ really doesn't matter what value we use because it will
+ have negligible effect in later calculations
+ */
+ XSIG_LL(accum) = 0x8000000000000000LL;
+ accum.lsw = 0;
+ } else {
+ div_Xsig(&accumulatoro, &accumulatore, &accum);
}
- else
- adj = 0;
- adj = mul_32_32(0x898cc517, adj); /* delta * tan^2 */
-
- fix_up.msw += adj;
- if ( !(fix_up.msw & 0x80000000) ) /* did fix_up overflow ? */
- {
- /* Yes, we need to add an msb */
- shr_Xsig(&fix_up, 1);
- fix_up.msw |= 0x80000000;
- shr_Xsig(&fix_up, 64 + exponent);
+
+ /* Multiply by 1/3 * arg^3 */
+ mul64_Xsig(&accum, &XSIG_LL(argSignif));
+ mul64_Xsig(&accum, &XSIG_LL(argSignif));
+ mul64_Xsig(&accum, &XSIG_LL(argSignif));
+ mul64_Xsig(&accum, &twothirds);
+ shr_Xsig(&accum, -2 * (exponent + 1));
+
+ /* tan(arg) = arg + accum */
+ add_two_Xsig(&accum, &argSignif, &exponent);
+
+ if (invert) {
+ /* We now have the value of tan(pi_2 - arg) where pi_2 is an
+ approximation for pi/2
+ */
+ /* The next step is to fix the answer to compensate for the
+ error due to the approximation used for pi/2
+ */
+
+ /* This is (approx) delta, the error in our approx for pi/2
+ (see above). It has an exponent of -65
+ */
+ XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
+ fix_up.lsw = 0;
+
+ if (exponent == 0)
+ adj = 0xffffffff; /* We want approx 1.0 here, but
+ this is close enough. */
+ else if (exponent > -30) {
+ adj = accum.msw >> -(exponent + 1); /* tan */
+ adj = mul_32_32(adj, adj); /* tan^2 */
+ } else
+ adj = 0;
+ adj = mul_32_32(0x898cc517, adj); /* delta * tan^2 */
+
+ fix_up.msw += adj;
+ if (!(fix_up.msw & 0x80000000)) { /* did fix_up overflow ? */
+ /* Yes, we need to add an msb */
+ shr_Xsig(&fix_up, 1);
+ fix_up.msw |= 0x80000000;
+ shr_Xsig(&fix_up, 64 + exponent);
+ } else
+ shr_Xsig(&fix_up, 65 + exponent);
+
+ add_two_Xsig(&accum, &fix_up, &exponent);
+
+ /* accum now contains tan(pi/2 - arg).
+ Use tan(arg) = 1.0 / tan(pi/2 - arg)
+ */
+ accumulatoro.lsw = accumulatoro.midw = 0;
+ accumulatoro.msw = 0x80000000;
+ div_Xsig(&accumulatoro, &accum, &accum);
+ exponent = -exponent - 1;
}
- else
- shr_Xsig(&fix_up, 65 + exponent);
-
- add_two_Xsig(&accum, &fix_up, &exponent);
-
- /* accum now contains tan(pi/2 - arg).
- Use tan(arg) = 1.0 / tan(pi/2 - arg)
- */
- accumulatoro.lsw = accumulatoro.midw = 0;
- accumulatoro.msw = 0x80000000;
- div_Xsig(&accumulatoro, &accum, &accum);
- exponent = - exponent - 1;
- }
-
- /* Transfer the result */
- round_Xsig(&accum);
- FPU_settag0(TAG_Valid);
- significand(st0_ptr) = XSIG_LL(accum);
- setexponent16(st0_ptr, exponent + EXTENDED_Ebias); /* Result is positive. */
+
+ /* Transfer the result */
+ round_Xsig(&accum);
+ FPU_settag0(TAG_Valid);
+ significand(st0_ptr) = XSIG_LL(accum);
+ setexponent16(st0_ptr, exponent + EXTENDED_Ebias); /* Result is positive. */
}
git.openpandora.org / pandora-kernel.git / blobdiff