aboutsummaryrefslogtreecommitdiff
path: root/arch/x86/math-emu/poly_sin.c
diff options
context:
space:
mode:
Diffstat (limited to 'arch/x86/math-emu/poly_sin.c')
-rw-r--r--arch/x86/math-emu/poly_sin.c599
1 files changed, 290 insertions, 309 deletions
diff --git a/arch/x86/math-emu/poly_sin.c b/arch/x86/math-emu/poly_sin.c
index a36313fb06f..b862039c728 100644
--- a/arch/x86/math-emu/poly_sin.c
+++ b/arch/x86/math-emu/poly_sin.c
@@ -11,7 +11,6 @@
| |
+---------------------------------------------------------------------------*/
-
#include "exception.h"
#include "reg_constant.h"
#include "fpu_emu.h"
@@ -19,379 +18,361 @@
#include "control_w.h"
#include "poly.h"
-
#define N_COEFF_P 4
#define N_COEFF_N 4
-static const unsigned long long pos_terms_l[N_COEFF_P] =
-{
- 0xaaaaaaaaaaaaaaabLL,
- 0x00d00d00d00cf906LL,
- 0x000006b99159a8bbLL,
- 0x000000000d7392e6LL
+static const unsigned long long pos_terms_l[N_COEFF_P] = {
+ 0xaaaaaaaaaaaaaaabLL,
+ 0x00d00d00d00cf906LL,
+ 0x000006b99159a8bbLL,
+ 0x000000000d7392e6LL
};
-static const unsigned long long neg_terms_l[N_COEFF_N] =
-{
- 0x2222222222222167LL,
- 0x0002e3bc74aab624LL,
- 0x0000000b09229062LL,
- 0x00000000000c7973LL
+static const unsigned long long neg_terms_l[N_COEFF_N] = {
+ 0x2222222222222167LL,
+ 0x0002e3bc74aab624LL,
+ 0x0000000b09229062LL,
+ 0x00000000000c7973LL
};
-
-
#define N_COEFF_PH 4
#define N_COEFF_NH 4
-static const unsigned long long pos_terms_h[N_COEFF_PH] =
-{
- 0x0000000000000000LL,
- 0x05b05b05b05b0406LL,
- 0x000049f93edd91a9LL,
- 0x00000000c9c9ed62LL
+static const unsigned long long pos_terms_h[N_COEFF_PH] = {
+ 0x0000000000000000LL,
+ 0x05b05b05b05b0406LL,
+ 0x000049f93edd91a9LL,
+ 0x00000000c9c9ed62LL
};
-static const unsigned long long neg_terms_h[N_COEFF_NH] =
-{
- 0xaaaaaaaaaaaaaa98LL,
- 0x001a01a01a019064LL,
- 0x0000008f76c68a77LL,
- 0x0000000000d58f5eLL
+static const unsigned long long neg_terms_h[N_COEFF_NH] = {
+ 0xaaaaaaaaaaaaaa98LL,
+ 0x001a01a01a019064LL,
+ 0x0000008f76c68a77LL,
+ 0x0000000000d58f5eLL
};
-
/*--- poly_sine() -----------------------------------------------------------+
| |
+---------------------------------------------------------------------------*/
-void poly_sine(FPU_REG *st0_ptr)
+void poly_sine(FPU_REG *st0_ptr)
{
- int exponent, echange;
- Xsig accumulator, argSqrd, argTo4;
- unsigned long fix_up, adj;
- unsigned long long fixed_arg;
- FPU_REG result;
+ int exponent, echange;
+ Xsig accumulator, argSqrd, argTo4;
+ unsigned long fix_up, adj;
+ unsigned long long fixed_arg;
+ FPU_REG result;
- exponent = exponent(st0_ptr);
+ exponent = exponent(st0_ptr);
- accumulator.lsw = accumulator.midw = accumulator.msw = 0;
+ accumulator.lsw = accumulator.midw = accumulator.msw = 0;
- /* Split into two ranges, for arguments below and above 1.0 */
- /* The boundary between upper and lower is approx 0.88309101259 */
- if ( (exponent < -1) || ((exponent == -1) && (st0_ptr->sigh <= 0xe21240aa)) )
- {
- /* The argument is <= 0.88309101259 */
+ /* Split into two ranges, for arguments below and above 1.0 */
+ /* The boundary between upper and lower is approx 0.88309101259 */
+ if ((exponent < -1)
+ || ((exponent == -1) && (st0_ptr->sigh <= 0xe21240aa))) {
+ /* The argument is <= 0.88309101259 */
+
+ argSqrd.msw = st0_ptr->sigh;
+ argSqrd.midw = st0_ptr->sigl;
+ argSqrd.lsw = 0;
+ mul64_Xsig(&argSqrd, &significand(st0_ptr));
+ shr_Xsig(&argSqrd, 2 * (-1 - exponent));
+ argTo4.msw = argSqrd.msw;
+ argTo4.midw = argSqrd.midw;
+ argTo4.lsw = argSqrd.lsw;
+ mul_Xsig_Xsig(&argTo4, &argTo4);
- argSqrd.msw = st0_ptr->sigh; argSqrd.midw = st0_ptr->sigl; argSqrd.lsw = 0;
- mul64_Xsig(&argSqrd, &significand(st0_ptr));
- shr_Xsig(&argSqrd, 2*(-1-exponent));
- argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
- argTo4.lsw = argSqrd.lsw;
- mul_Xsig_Xsig(&argTo4, &argTo4);
+ polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
+ N_COEFF_N - 1);
+ mul_Xsig_Xsig(&accumulator, &argSqrd);
+ negate_Xsig(&accumulator);
- polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
- N_COEFF_N-1);
- mul_Xsig_Xsig(&accumulator, &argSqrd);
- negate_Xsig(&accumulator);
+ polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
+ N_COEFF_P - 1);
- polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
- N_COEFF_P-1);
+ shr_Xsig(&accumulator, 2); /* Divide by four */
+ accumulator.msw |= 0x80000000; /* Add 1.0 */
- shr_Xsig(&accumulator, 2); /* Divide by four */
- accumulator.msw |= 0x80000000; /* Add 1.0 */
+ mul64_Xsig(&accumulator, &significand(st0_ptr));
+ mul64_Xsig(&accumulator, &significand(st0_ptr));
+ mul64_Xsig(&accumulator, &significand(st0_ptr));
- mul64_Xsig(&accumulator, &significand(st0_ptr));
- mul64_Xsig(&accumulator, &significand(st0_ptr));
- mul64_Xsig(&accumulator, &significand(st0_ptr));
+ /* Divide by four, FPU_REG compatible, etc */
+ exponent = 3 * exponent;
- /* Divide by four, FPU_REG compatible, etc */
- exponent = 3*exponent;
+ /* The minimum exponent difference is 3 */
+ shr_Xsig(&accumulator, exponent(st0_ptr) - exponent);
- /* The minimum exponent difference is 3 */
- shr_Xsig(&accumulator, exponent(st0_ptr) - exponent);
+ negate_Xsig(&accumulator);
+ XSIG_LL(accumulator) += significand(st0_ptr);
- negate_Xsig(&accumulator);
- XSIG_LL(accumulator) += significand(st0_ptr);
+ echange = round_Xsig(&accumulator);
- echange = round_Xsig(&accumulator);
+ setexponentpos(&result, exponent(st0_ptr) + echange);
+ } else {
+ /* The argument is > 0.88309101259 */
+ /* We use sin(st(0)) = cos(pi/2-st(0)) */
- setexponentpos(&result, exponent(st0_ptr) + echange);
- }
- else
- {
- /* The argument is > 0.88309101259 */
- /* We use sin(st(0)) = cos(pi/2-st(0)) */
+ fixed_arg = significand(st0_ptr);
- fixed_arg = significand(st0_ptr);
+ if (exponent == 0) {
+ /* The argument is >= 1.0 */
- if ( exponent == 0 )
- {
- /* The argument is >= 1.0 */
+ /* Put the binary point at the left. */
+ fixed_arg <<= 1;
+ }
+ /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
+ fixed_arg = 0x921fb54442d18469LL - fixed_arg;
+ /* There is a special case which arises due to rounding, to fix here. */
+ if (fixed_arg == 0xffffffffffffffffLL)
+ fixed_arg = 0;
- /* Put the binary point at the left. */
- fixed_arg <<= 1;
- }
- /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
- fixed_arg = 0x921fb54442d18469LL - fixed_arg;
- /* There is a special case which arises due to rounding, to fix here. */
- if ( fixed_arg == 0xffffffffffffffffLL )
- fixed_arg = 0;
+ XSIG_LL(argSqrd) = fixed_arg;
+ argSqrd.lsw = 0;
+ mul64_Xsig(&argSqrd, &fixed_arg);
- XSIG_LL(argSqrd) = fixed_arg; argSqrd.lsw = 0;
- mul64_Xsig(&argSqrd, &fixed_arg);
+ XSIG_LL(argTo4) = XSIG_LL(argSqrd);
+ argTo4.lsw = argSqrd.lsw;
+ mul_Xsig_Xsig(&argTo4, &argTo4);
- XSIG_LL(argTo4) = XSIG_LL(argSqrd); argTo4.lsw = argSqrd.lsw;
- mul_Xsig_Xsig(&argTo4, &argTo4);
+ polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
+ N_COEFF_NH - 1);
+ mul_Xsig_Xsig(&accumulator, &argSqrd);
+ negate_Xsig(&accumulator);
- polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
- N_COEFF_NH-1);
- mul_Xsig_Xsig(&accumulator, &argSqrd);
- negate_Xsig(&accumulator);
+ polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
+ N_COEFF_PH - 1);
+ negate_Xsig(&accumulator);
- polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
- N_COEFF_PH-1);
- negate_Xsig(&accumulator);
+ mul64_Xsig(&accumulator, &fixed_arg);
+ mul64_Xsig(&accumulator, &fixed_arg);
- mul64_Xsig(&accumulator, &fixed_arg);
- mul64_Xsig(&accumulator, &fixed_arg);
+ shr_Xsig(&accumulator, 3);
+ negate_Xsig(&accumulator);
- shr_Xsig(&accumulator, 3);
- negate_Xsig(&accumulator);
+ add_Xsig_Xsig(&accumulator, &argSqrd);
- add_Xsig_Xsig(&accumulator, &argSqrd);
+ shr_Xsig(&accumulator, 1);
- shr_Xsig(&accumulator, 1);
+ accumulator.lsw |= 1; /* A zero accumulator here would cause problems */
+ negate_Xsig(&accumulator);
- accumulator.lsw |= 1; /* A zero accumulator here would cause problems */
- negate_Xsig(&accumulator);
+ /* The basic computation is complete. Now fix the answer to
+ compensate for the error due to the approximation used for
+ pi/2
+ */
- /* The basic computation is complete. Now fix the answer to
- compensate for the error due to the approximation used for
- pi/2
- */
+ /* This has an exponent of -65 */
+ fix_up = 0x898cc517;
+ /* The fix-up needs to be improved for larger args */
+ if (argSqrd.msw & 0xffc00000) {
+ /* Get about 32 bit precision in these: */
+ fix_up -= mul_32_32(0x898cc517, argSqrd.msw) / 6;
+ }
+ fix_up = mul_32_32(fix_up, LL_MSW(fixed_arg));
- /* This has an exponent of -65 */
- fix_up = 0x898cc517;
- /* The fix-up needs to be improved for larger args */
- if ( argSqrd.msw & 0xffc00000 )
- {
- /* Get about 32 bit precision in these: */
- fix_up -= mul_32_32(0x898cc517, argSqrd.msw) / 6;
- }
- fix_up = mul_32_32(fix_up, LL_MSW(fixed_arg));
+ adj = accumulator.lsw; /* temp save */
+ accumulator.lsw -= fix_up;
+ if (accumulator.lsw > adj)
+ XSIG_LL(accumulator)--;
- adj = accumulator.lsw; /* temp save */
- accumulator.lsw -= fix_up;
- if ( accumulator.lsw > adj )
- XSIG_LL(accumulator) --;
+ echange = round_Xsig(&accumulator);
- echange = round_Xsig(&accumulator);
-
- setexponentpos(&result, echange - 1);
- }
+ setexponentpos(&result, echange - 1);
+ }
- significand(&result) = XSIG_LL(accumulator);
- setsign(&result, getsign(st0_ptr));
- FPU_copy_to_reg0(&result, TAG_Valid);
+ significand(&result) = XSIG_LL(accumulator);
+ setsign(&result, getsign(st0_ptr));
+ FPU_copy_to_reg0(&result, TAG_Valid);
#ifdef PARANOID
- if ( (exponent(&result) >= 0)
- && (significand(&result) > 0x8000000000000000LL) )
- {
- EXCEPTION(EX_INTERNAL|0x150);
- }
+ if ((exponent(&result) >= 0)
+ && (significand(&result) > 0x8000000000000000LL)) {
+ EXCEPTION(EX_INTERNAL | 0x150);
+ }
#endif /* PARANOID */
}
-
-
/*--- poly_cos() ------------------------------------------------------------+
| |
+---------------------------------------------------------------------------*/
-void poly_cos(FPU_REG *st0_ptr)
+void poly_cos(FPU_REG *st0_ptr)
{
- FPU_REG result;
- long int exponent, exp2, echange;
- Xsig accumulator, argSqrd, fix_up, argTo4;
- unsigned long long fixed_arg;
+ FPU_REG result;
+ long int exponent, exp2, echange;
+ Xsig accumulator, argSqrd, fix_up, argTo4;
+ unsigned long long fixed_arg;
#ifdef PARANOID
- if ( (exponent(st0_ptr) > 0)
- || ((exponent(st0_ptr) == 0)
- && (significand(st0_ptr) > 0xc90fdaa22168c234LL)) )
- {
- EXCEPTION(EX_Invalid);
- FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
- return;
- }
-#endif /* PARANOID */
-
- exponent = exponent(st0_ptr);
-
- accumulator.lsw = accumulator.midw = accumulator.msw = 0;
-
- if ( (exponent < -1) || ((exponent == -1) && (st0_ptr->sigh <= 0xb00d6f54)) )
- {
- /* arg is < 0.687705 */
-
- argSqrd.msw = st0_ptr->sigh; argSqrd.midw = st0_ptr->sigl;
- argSqrd.lsw = 0;
- mul64_Xsig(&argSqrd, &significand(st0_ptr));
-
- if ( exponent < -1 )
- {
- /* shift the argument right by the required places */
- shr_Xsig(&argSqrd, 2*(-1-exponent));
- }
-
- argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
- argTo4.lsw = argSqrd.lsw;
- mul_Xsig_Xsig(&argTo4, &argTo4);
-
- polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
- N_COEFF_NH-1);
- mul_Xsig_Xsig(&accumulator, &argSqrd);
- negate_Xsig(&accumulator);
-
- polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
- N_COEFF_PH-1);
- negate_Xsig(&accumulator);
-
- mul64_Xsig(&accumulator, &significand(st0_ptr));
- mul64_Xsig(&accumulator, &significand(st0_ptr));
- shr_Xsig(&accumulator, -2*(1+exponent));
-
- shr_Xsig(&accumulator, 3);
- negate_Xsig(&accumulator);
-
- add_Xsig_Xsig(&accumulator, &argSqrd);
-
- shr_Xsig(&accumulator, 1);
-
- /* It doesn't matter if accumulator is all zero here, the
- following code will work ok */
- negate_Xsig(&accumulator);
-
- if ( accumulator.lsw & 0x80000000 )
- XSIG_LL(accumulator) ++;
- if ( accumulator.msw == 0 )
- {
- /* The result is 1.0 */
- FPU_copy_to_reg0(&CONST_1, TAG_Valid);
- return;
- }
- else
- {
- significand(&result) = XSIG_LL(accumulator);
-
- /* will be a valid positive nr with expon = -1 */
- setexponentpos(&result, -1);
- }
- }
- else
- {
- fixed_arg = significand(st0_ptr);
-
- if ( exponent == 0 )
- {
- /* The argument is >= 1.0 */
-
- /* Put the binary point at the left. */
- fixed_arg <<= 1;
- }
- /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
- fixed_arg = 0x921fb54442d18469LL - fixed_arg;
- /* There is a special case which arises due to rounding, to fix here. */
- if ( fixed_arg == 0xffffffffffffffffLL )
- fixed_arg = 0;
-
- exponent = -1;
- exp2 = -1;
-
- /* A shift is needed here only for a narrow range of arguments,
- i.e. for fixed_arg approx 2^-32, but we pick up more... */
- if ( !(LL_MSW(fixed_arg) & 0xffff0000) )
- {
- fixed_arg <<= 16;
- exponent -= 16;
- exp2 -= 16;
+ if ((exponent(st0_ptr) > 0)
+ || ((exponent(st0_ptr) == 0)
+ && (significand(st0_ptr) > 0xc90fdaa22168c234LL))) {
+ EXCEPTION(EX_Invalid);
+ FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
+ return;
}
+#endif /* PARANOID */
- XSIG_LL(argSqrd) = fixed_arg; argSqrd.lsw = 0;
- mul64_Xsig(&argSqrd, &fixed_arg);
-
- if ( exponent < -1 )
- {
- /* shift the argument right by the required places */
- shr_Xsig(&argSqrd, 2*(-1-exponent));
- }
-
- argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
- argTo4.lsw = argSqrd.lsw;
- mul_Xsig_Xsig(&argTo4, &argTo4);
-
- polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
- N_COEFF_N-1);
- mul_Xsig_Xsig(&accumulator, &argSqrd);
- negate_Xsig(&accumulator);
-
- polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
- N_COEFF_P-1);
-
- shr_Xsig(&accumulator, 2); /* Divide by four */
- accumulator.msw |= 0x80000000; /* Add 1.0 */
-
- mul64_Xsig(&accumulator, &fixed_arg);
- mul64_Xsig(&accumulator, &fixed_arg);
- mul64_Xsig(&accumulator, &fixed_arg);
-
- /* Divide by four, FPU_REG compatible, etc */
- exponent = 3*exponent;
-
- /* The minimum exponent difference is 3 */
- shr_Xsig(&accumulator, exp2 - exponent);
-
- negate_Xsig(&accumulator);
- XSIG_LL(accumulator) += fixed_arg;
-
- /* The basic computation is complete. Now fix the answer to
- compensate for the error due to the approximation used for
- pi/2
- */
-
- /* This has an exponent of -65 */
- XSIG_LL(fix_up) = 0x898cc51701b839a2ll;
- fix_up.lsw = 0;
-
- /* The fix-up needs to be improved for larger args */
- if ( argSqrd.msw & 0xffc00000 )
- {
- /* Get about 32 bit precision in these: */
- fix_up.msw -= mul_32_32(0x898cc517, argSqrd.msw) / 2;
- fix_up.msw += mul_32_32(0x898cc517, argTo4.msw) / 24;
+ exponent = exponent(st0_ptr);
+
+ accumulator.lsw = accumulator.midw = accumulator.msw = 0;
+
+ if ((exponent < -1)
+ || ((exponent == -1) && (st0_ptr->sigh <= 0xb00d6f54))) {
+ /* arg is < 0.687705 */
+
+ argSqrd.msw = st0_ptr->sigh;
+ argSqrd.midw = st0_ptr->sigl;
+ argSqrd.lsw = 0;
+ mul64_Xsig(&argSqrd, &significand(st0_ptr));
+
+ if (exponent < -1) {
+ /* shift the argument right by the required places */
+ shr_Xsig(&argSqrd, 2 * (-1 - exponent));
+ }
+
+ argTo4.msw = argSqrd.msw;
+ argTo4.midw = argSqrd.midw;
+ argTo4.lsw = argSqrd.lsw;
+ mul_Xsig_Xsig(&argTo4, &argTo4);
+
+ polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
+ N_COEFF_NH - 1);
+ mul_Xsig_Xsig(&accumulator, &argSqrd);
+ negate_Xsig(&accumulator);
+
+ polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
+ N_COEFF_PH - 1);
+ negate_Xsig(&accumulator);
+
+ mul64_Xsig(&accumulator, &significand(st0_ptr));
+ mul64_Xsig(&accumulator, &significand(st0_ptr));
+ shr_Xsig(&accumulator, -2 * (1 + exponent));
+
+ shr_Xsig(&accumulator, 3);
+ negate_Xsig(&accumulator);
+
+ add_Xsig_Xsig(&accumulator, &argSqrd);
+
+ shr_Xsig(&accumulator, 1);
+
+ /* It doesn't matter if accumulator is all zero here, the
+ following code will work ok */
+ negate_Xsig(&accumulator);
+
+ if (accumulator.lsw & 0x80000000)
+ XSIG_LL(accumulator)++;
+ if (accumulator.msw == 0) {
+ /* The result is 1.0 */
+ FPU_copy_to_reg0(&CONST_1, TAG_Valid);
+ return;
+ } else {
+ significand(&result) = XSIG_LL(accumulator);
+
+ /* will be a valid positive nr with expon = -1 */
+ setexponentpos(&result, -1);
+ }
+ } else {
+ fixed_arg = significand(st0_ptr);
+
+ if (exponent == 0) {
+ /* The argument is >= 1.0 */
+
+ /* Put the binary point at the left. */
+ fixed_arg <<= 1;
+ }
+ /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
+ fixed_arg = 0x921fb54442d18469LL - fixed_arg;
+ /* There is a special case which arises due to rounding, to fix here. */
+ if (fixed_arg == 0xffffffffffffffffLL)
+ fixed_arg = 0;
+
+ exponent = -1;
+ exp2 = -1;
+
+ /* A shift is needed here only for a narrow range of arguments,
+ i.e. for fixed_arg approx 2^-32, but we pick up more... */
+ if (!(LL_MSW(fixed_arg) & 0xffff0000)) {
+ fixed_arg <<= 16;
+ exponent -= 16;
+ exp2 -= 16;
+ }
+
+ XSIG_LL(argSqrd) = fixed_arg;
+ argSqrd.lsw = 0;
+ mul64_Xsig(&argSqrd, &fixed_arg);
+
+ if (exponent < -1) {
+ /* shift the argument right by the required places */
+ shr_Xsig(&argSqrd, 2 * (-1 - exponent));
+ }
+
+ argTo4.msw = argSqrd.msw;
+ argTo4.midw = argSqrd.midw;
+ argTo4.lsw = argSqrd.lsw;
+ mul_Xsig_Xsig(&argTo4, &argTo4);
+
+ polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
+ N_COEFF_N - 1);
+ mul_Xsig_Xsig(&accumulator, &argSqrd);
+ negate_Xsig(&accumulator);
+
+ polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
+ N_COEFF_P - 1);
+
+ shr_Xsig(&accumulator, 2); /* Divide by four */
+ accumulator.msw |= 0x80000000; /* Add 1.0 */
+
+ mul64_Xsig(&accumulator, &fixed_arg);
+ mul64_Xsig(&accumulator, &fixed_arg);
+ mul64_Xsig(&accumulator, &fixed_arg);
+
+ /* Divide by four, FPU_REG compatible, etc */
+ exponent = 3 * exponent;
+
+ /* The minimum exponent difference is 3 */
+ shr_Xsig(&accumulator, exp2 - exponent);
+
+ negate_Xsig(&accumulator);
+ XSIG_LL(accumulator) += fixed_arg;
+
+ /* The basic computation is complete. Now fix the answer to
+ compensate for the error due to the approximation used for
+ pi/2
+ */
+
+ /* This has an exponent of -65 */
+ XSIG_LL(fix_up) = 0x898cc51701b839a2ll;
+ fix_up.lsw = 0;
+
+ /* The fix-up needs to be improved for larger args */
+ if (argSqrd.msw & 0xffc00000) {
+ /* Get about 32 bit precision in these: */
+ fix_up.msw -= mul_32_32(0x898cc517, argSqrd.msw) / 2;
+ fix_up.msw += mul_32_32(0x898cc517, argTo4.msw) / 24;
+ }
+
+ exp2 += norm_Xsig(&accumulator);
+ shr_Xsig(&accumulator, 1); /* Prevent overflow */
+ exp2++;
+ shr_Xsig(&fix_up, 65 + exp2);
+
+ add_Xsig_Xsig(&accumulator, &fix_up);
+
+ echange = round_Xsig(&accumulator);
+
+ setexponentpos(&result, exp2 + echange);
+ significand(&result) = XSIG_LL(accumulator);
}
- exp2 += norm_Xsig(&accumulator);
- shr_Xsig(&accumulator, 1); /* Prevent overflow */
- exp2++;
- shr_Xsig(&fix_up, 65 + exp2);
-
- add_Xsig_Xsig(&accumulator, &fix_up);
-
- echange = round_Xsig(&accumulator);
-
- setexponentpos(&result, exp2 + echange);
- significand(&result) = XSIG_LL(accumulator);
- }
-
- FPU_copy_to_reg0(&result, TAG_Valid);
+ FPU_copy_to_reg0(&result, TAG_Valid);
#ifdef PARANOID
- if ( (exponent(&result) >= 0)
- && (significand(&result) > 0x8000000000000000LL) )
- {
- EXCEPTION(EX_INTERNAL|0x151);
- }
+ if ((exponent(&result) >= 0)
+ && (significand(&result) > 0x8000000000000000LL)) {
+ EXCEPTION(EX_INTERNAL | 0x151);
+ }
#endif /* PARANOID */
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-24 23:40:16 +0000