summaryrefslogtreecommitdiffstats
path: root/arch/x86/math-emu/poly_sin.c
diff options
context:
space:
mode:
authorIngo Molnar <mingo@elte.hu>2008-01-30 13:30:11 +0100
committerIngo Molnar <mingo@elte.hu>2008-01-30 13:30:11 +0100
commit3d0d14f983b55a570b976976284df4c434af3223 (patch)
tree864f11c0ce5ee1e15acdd196018b79d0d0e2685d /arch/x86/math-emu/poly_sin.c
parenta4ec1effce83796209a0258602b0cf50026d86f2 (diff)
downloadop-kernel-dev-3d0d14f983b55a570b976976284df4c434af3223.zip
op-kernel-dev-3d0d14f983b55a570b976976284df4c434af3223.tar.gz
x86: lindent arch/i386/math-emu
lindent these files: errors lines of code errors/KLOC arch/x86/math-emu/ 2236 9424 237.2 arch/x86/math-emu/ 128 8706 14.7 no other changes. No code changed: text data bss dec hex filename 5589802 612739 3833856 10036397 9924ad vmlinux.before 5589802 612739 3833856 10036397 9924ad vmlinux.after the intent of this patch is to ease the automated tracking of kernel code quality - it's just much easier for us to maintain it if every file in arch/x86 is supposed to be clean. NOTE: it is a known problem of lindent that it causes some style damage of its own, but it's a safe tool (well, except for the gcc array range initializers extension), so we did the bulk of the changes via lindent, and did the manual fixups in a followup patch. the resulting math-emu code has been tested by Thomas Gleixner on a real 386 DX CPU as well, and it works fine. Signed-off-by: Ingo Molnar <mingo@elte.hu> Signed-off-by: Thomas Gleixner <tglx@linutronix.de>
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 a36313f..7273ae0 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 */
}
OpenPOWER on IntegriCloud