diff options
author | Thomas Gleixner <tglx@linutronix.de> | 2007-10-11 11:16:31 +0200 |
---|---|---|
committer | Thomas Gleixner <tglx@linutronix.de> | 2007-10-11 11:16:31 +0200 |
commit | da957e111bb0c189a4a3bf8a00caaecb59ed94ca (patch) | |
tree | 6916075fdd3e28869dcd3dfa2cf160a74d1cb02e /arch/i386/math-emu/poly_sin.c | |
parent | 2ec1df4130c60d1eb49dc0fa0ed15858fede6b05 (diff) |
i386: move math-emu
Signed-off-by: Thomas Gleixner <tglx@linutronix.de>
Signed-off-by: Ingo Molnar <mingo@elte.hu>
Diffstat (limited to 'arch/i386/math-emu/poly_sin.c')
-rw-r--r-- | arch/i386/math-emu/poly_sin.c | 397 |
1 files changed, 0 insertions, 397 deletions
diff --git a/arch/i386/math-emu/poly_sin.c b/arch/i386/math-emu/poly_sin.c deleted file mode 100644 index a36313fb06f..00000000000 --- a/arch/i386/math-emu/poly_sin.c +++ /dev/null @@ -1,397 +0,0 @@ -/*---------------------------------------------------------------------------+ - | poly_sin.c | - | | - | Computation of an approximation of the sin function and the cosine | - | function by a polynomial. | - | | - | Copyright (C) 1992,1993,1994,1997,1999 | - | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia | - | E-mail billm@melbpc.org.au | - | | - | | - +---------------------------------------------------------------------------*/ - - -#include "exception.h" -#include "reg_constant.h" -#include "fpu_emu.h" -#include "fpu_system.h" -#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 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 neg_terms_h[N_COEFF_NH] = -{ - 0xaaaaaaaaaaaaaa98LL, - 0x001a01a01a019064LL, - 0x0000008f76c68a77LL, - 0x0000000000d58f5eLL -}; - - -/*--- poly_sine() -----------------------------------------------------------+ - | | - +---------------------------------------------------------------------------*/ -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; - - exponent = exponent(st0_ptr); - - 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 */ - - 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), pos_terms_l, - N_COEFF_P-1); - - 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)); - - /* Divide by four, FPU_REG compatible, etc */ - exponent = 3*exponent; - - /* The minimum exponent difference is 3 */ - shr_Xsig(&accumulator, exponent(st0_ptr) - exponent); - - negate_Xsig(&accumulator); - XSIG_LL(accumulator) += significand(st0_ptr); - - 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)) */ - - 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; - - 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); - - 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, &fixed_arg); - mul64_Xsig(&accumulator, &fixed_arg); - - shr_Xsig(&accumulator, 3); - negate_Xsig(&accumulator); - - add_Xsig_Xsig(&accumulator, &argSqrd); - - shr_Xsig(&accumulator, 1); - - 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 - */ - - /* 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) --; - - echange = round_Xsig(&accumulator); - - setexponentpos(&result, echange - 1); - } - - 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); - } -#endif /* PARANOID */ - -} - - - -/*--- poly_cos() ------------------------------------------------------------+ - | | - +---------------------------------------------------------------------------*/ -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; - -#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; - } - - 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); - } - - FPU_copy_to_reg0(&result, TAG_Valid); - -#ifdef PARANOID - if ( (exponent(&result) >= 0) - && (significand(&result) > 0x8000000000000000LL) ) - { - EXCEPTION(EX_INTERNAL|0x151); - } -#endif /* PARANOID */ - -} |