diff options
Diffstat (limited to 'arch/x86/math-emu/poly_l2.c')
| -rw-r--r-- | arch/x86/math-emu/poly_l2.c | 376 |
1 files changed, 174 insertions, 202 deletions
diff --git a/arch/x86/math-emu/poly_l2.c b/arch/x86/math-emu/poly_l2.c index dd00e1d5b074..8e2ff4b28a0a 100644 --- a/arch/x86/math-emu/poly_l2.c +++ b/arch/x86/math-emu/poly_l2.c @@ -10,7 +10,6 @@ | | +---------------------------------------------------------------------------*/ - #include "exception.h" #include "reg_constant.h" #include "fpu_emu.h" @@ -18,184 +17,163 @@ #include "control_w.h" #include "poly.h" - static void log2_kernel(FPU_REG const *arg, u_char argsign, - Xsig *accum_result, long int *expon); - + Xsig * accum_result, long int *expon); /*--- poly_l2() -------------------------------------------------------------+ | Base 2 logarithm by a polynomial approximation. | +---------------------------------------------------------------------------*/ -void poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign) +void poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign) { - long int exponent, expon, expon_expon; - Xsig accumulator, expon_accum, yaccum; - u_char sign, argsign; - FPU_REG x; - int tag; - - exponent = exponent16(st0_ptr); - - /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */ - if ( st0_ptr->sigh > (unsigned)0xb504f334 ) - { - /* Treat as sqrt(2)/2 < st0_ptr < 1 */ - significand(&x) = - significand(st0_ptr); - setexponent16(&x, -1); - exponent++; - argsign = SIGN_NEG; - } - else - { - /* Treat as 1 <= st0_ptr < sqrt(2) */ - x.sigh = st0_ptr->sigh - 0x80000000; - x.sigl = st0_ptr->sigl; - setexponent16(&x, 0); - argsign = SIGN_POS; - } - tag = FPU_normalize_nuo(&x); - - if ( tag == TAG_Zero ) - { - expon = 0; - accumulator.msw = accumulator.midw = accumulator.lsw = 0; - } - else - { - log2_kernel(&x, argsign, &accumulator, &expon); - } - - if ( exponent < 0 ) - { - sign = SIGN_NEG; - exponent = -exponent; - } - else - sign = SIGN_POS; - expon_accum.msw = exponent; expon_accum.midw = expon_accum.lsw = 0; - if ( exponent ) - { - expon_expon = 31 + norm_Xsig(&expon_accum); - shr_Xsig(&accumulator, expon_expon - expon); - - if ( sign ^ argsign ) - negate_Xsig(&accumulator); - add_Xsig_Xsig(&accumulator, &expon_accum); - } - else - { - expon_expon = expon; - sign = argsign; - } - - yaccum.lsw = 0; XSIG_LL(yaccum) = significand(st1_ptr); - mul_Xsig_Xsig(&accumulator, &yaccum); - - expon_expon += round_Xsig(&accumulator); - - if ( accumulator.msw == 0 ) - { - FPU_copy_to_reg1(&CONST_Z, TAG_Zero); - return; - } - - significand(st1_ptr) = XSIG_LL(accumulator); - setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1); - - tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign); - FPU_settagi(1, tag); - - set_precision_flag_up(); /* 80486 appears to always do this */ - - return; + long int exponent, expon, expon_expon; + Xsig accumulator, expon_accum, yaccum; + u_char sign, argsign; + FPU_REG x; + int tag; + + exponent = exponent16(st0_ptr); + + /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */ + if (st0_ptr->sigh > (unsigned)0xb504f334) { + /* Treat as sqrt(2)/2 < st0_ptr < 1 */ + significand(&x) = -significand(st0_ptr); + setexponent16(&x, -1); + exponent++; + argsign = SIGN_NEG; + } else { + /* Treat as 1 <= st0_ptr < sqrt(2) */ + x.sigh = st0_ptr->sigh - 0x80000000; + x.sigl = st0_ptr->sigl; + setexponent16(&x, 0); + argsign = SIGN_POS; + } + tag = FPU_normalize_nuo(&x); -} + if (tag == TAG_Zero) { + expon = 0; + accumulator.msw = accumulator.midw = accumulator.lsw = 0; + } else { + log2_kernel(&x, argsign, &accumulator, &expon); + } + + if (exponent < 0) { + sign = SIGN_NEG; + exponent = -exponent; + } else + sign = SIGN_POS; + expon_accum.msw = exponent; + expon_accum.midw = expon_accum.lsw = 0; + if (exponent) { + expon_expon = 31 + norm_Xsig(&expon_accum); + shr_Xsig(&accumulator, expon_expon - expon); + + if (sign ^ argsign) + negate_Xsig(&accumulator); + add_Xsig_Xsig(&accumulator, &expon_accum); + } else { + expon_expon = expon; + sign = argsign; + } + + yaccum.lsw = 0; + XSIG_LL(yaccum) = significand(st1_ptr); + mul_Xsig_Xsig(&accumulator, &yaccum); + + expon_expon += round_Xsig(&accumulator); + + if (accumulator.msw == 0) { + FPU_copy_to_reg1(&CONST_Z, TAG_Zero); + return; + } + + significand(st1_ptr) = XSIG_LL(accumulator); + setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1); + tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign); + FPU_settagi(1, tag); + + set_precision_flag_up(); /* 80486 appears to always do this */ + + return; + +} /*--- poly_l2p1() -----------------------------------------------------------+ | Base 2 logarithm by a polynomial approximation. | | log2(x+1) | +---------------------------------------------------------------------------*/ -int poly_l2p1(u_char sign0, u_char sign1, - FPU_REG *st0_ptr, FPU_REG *st1_ptr, FPU_REG *dest) +int poly_l2p1(u_char sign0, u_char sign1, + FPU_REG * st0_ptr, FPU_REG * st1_ptr, FPU_REG * dest) { - u_char tag; - long int exponent; - Xsig accumulator, yaccum; + u_char tag; + long int exponent; + Xsig accumulator, yaccum; - if ( exponent16(st0_ptr) < 0 ) - { - log2_kernel(st0_ptr, sign0, &accumulator, &exponent); + if (exponent16(st0_ptr) < 0) { + log2_kernel(st0_ptr, sign0, &accumulator, &exponent); - yaccum.lsw = 0; - XSIG_LL(yaccum) = significand(st1_ptr); - mul_Xsig_Xsig(&accumulator, &yaccum); + yaccum.lsw = 0; + XSIG_LL(yaccum) = significand(st1_ptr); + mul_Xsig_Xsig(&accumulator, &yaccum); - exponent += round_Xsig(&accumulator); + exponent += round_Xsig(&accumulator); - exponent += exponent16(st1_ptr) + 1; - if ( exponent < EXP_WAY_UNDER ) exponent = EXP_WAY_UNDER; + exponent += exponent16(st1_ptr) + 1; + if (exponent < EXP_WAY_UNDER) + exponent = EXP_WAY_UNDER; - significand(dest) = XSIG_LL(accumulator); - setexponent16(dest, exponent); + significand(dest) = XSIG_LL(accumulator); + setexponent16(dest, exponent); - tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1); - FPU_settagi(1, tag); + tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1); + FPU_settagi(1, tag); - if ( tag == TAG_Valid ) - set_precision_flag_up(); /* 80486 appears to always do this */ - } - else - { - /* The magnitude of st0_ptr is far too large. */ + if (tag == TAG_Valid) + set_precision_flag_up(); /* 80486 appears to always do this */ + } else { + /* The magnitude of st0_ptr is far too large. */ - if ( sign0 != SIGN_POS ) - { - /* Trying to get the log of a negative number. */ -#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */ - changesign(st1_ptr); + if (sign0 != SIGN_POS) { + /* Trying to get the log of a negative number. */ +#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */ + changesign(st1_ptr); #else - if ( arith_invalid(1) < 0 ) - return 1; + if (arith_invalid(1) < 0) + return 1; #endif /* PECULIAR_486 */ - } + } - /* 80486 appears to do this */ - if ( sign0 == SIGN_NEG ) - set_precision_flag_down(); - else - set_precision_flag_up(); - } + /* 80486 appears to do this */ + if (sign0 == SIGN_NEG) + set_precision_flag_down(); + else + set_precision_flag_up(); + } - if ( exponent(dest) <= EXP_UNDER ) - EXCEPTION(EX_Underflow); + if (exponent(dest) <= EXP_UNDER) + EXCEPTION(EX_Underflow); - return 0; + return 0; } - - - #undef HIPOWER #define HIPOWER 10 -static const unsigned long long logterms[HIPOWER] = -{ - 0x2a8eca5705fc2ef0LL, - 0xf6384ee1d01febceLL, - 0x093bb62877cdf642LL, - 0x006985d8a9ec439bLL, - 0x0005212c4f55a9c8LL, - 0x00004326a16927f0LL, - 0x0000038d1d80a0e7LL, - 0x0000003141cc80c6LL, - 0x00000002b1668c9fLL, - 0x000000002c7a46aaLL +static const unsigned long long logterms[HIPOWER] = { + 0x2a8eca5705fc2ef0LL, + 0xf6384ee1d01febceLL, + 0x093bb62877cdf642LL, + 0x006985d8a9ec439bLL, + 0x0005212c4f55a9c8LL, + 0x00004326a16927f0LL, + 0x0000038d1d80a0e7LL, + 0x0000003141cc80c6LL, + 0x00000002b1668c9fLL, + 0x000000002c7a46aaLL }; static const unsigned long leadterm = 0xb8000000; - /*--- log2_kernel() ---------------------------------------------------------+ | Base 2 logarithm by a polynomial approximation. | | log2(x+1) | @@ -203,70 +181,64 @@ static const unsigned long leadterm = 0xb8000000; static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result, long int *expon) { - long int exponent, adj; - unsigned long long Xsq; - Xsig accumulator, Numer, Denom, argSignif, arg_signif; - - exponent = exponent16(arg); - Numer.lsw = Denom.lsw = 0; - XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg); - if ( argsign == SIGN_POS ) - { - shr_Xsig(&Denom, 2 - (1 + exponent)); - Denom.msw |= 0x80000000; - div_Xsig(&Numer, &Denom, &argSignif); - } - else - { - shr_Xsig(&Denom, 1 - (1 + exponent)); - negate_Xsig(&Denom); - if ( Denom.msw & 0x80000000 ) - { - div_Xsig(&Numer, &Denom, &argSignif); - exponent ++; - } - else - { - /* Denom must be 1.0 */ - argSignif.lsw = Numer.lsw; argSignif.midw = Numer.midw; - argSignif.msw = Numer.msw; + long int exponent, adj; + unsigned long long Xsq; + Xsig accumulator, Numer, Denom, argSignif, arg_signif; + + exponent = exponent16(arg); + Numer.lsw = Denom.lsw = 0; + XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg); + if (argsign == SIGN_POS) { + shr_Xsig(&Denom, 2 - (1 + exponent)); + Denom.msw |= 0x80000000; + div_Xsig(&Numer, &Denom, &argSignif); + } else { + shr_Xsig(&Denom, 1 - (1 + exponent)); + negate_Xsig(&Denom); + if (Denom.msw & 0x80000000) { + div_Xsig(&Numer, &Denom, &argSignif); + exponent++; + } else { + /* Denom must be 1.0 */ + argSignif.lsw = Numer.lsw; + argSignif.midw = Numer.midw; + argSignif.msw = Numer.msw; + } } - } #ifndef PECULIAR_486 - /* Should check here that |local_arg| is within the valid range */ - if ( exponent >= -2 ) - { - if ( (exponent > -2) || - (argSignif.msw > (unsigned)0xafb0ccc0) ) - { - /* The argument is too large */ + /* Should check here that |local_arg| is within the valid range */ + if (exponent >= -2) { + if ((exponent > -2) || (argSignif.msw > (unsigned)0xafb0ccc0)) { + /* The argument is too large */ + } } - } #endif /* PECULIAR_486 */ - arg_signif.lsw = argSignif.lsw; XSIG_LL(arg_signif) = XSIG_LL(argSignif); - adj = norm_Xsig(&argSignif); - accumulator.lsw = argSignif.lsw; XSIG_LL(accumulator) = XSIG_LL(argSignif); - mul_Xsig_Xsig(&accumulator, &accumulator); - shr_Xsig(&accumulator, 2*(-1 - (1 + exponent + adj))); - Xsq = XSIG_LL(accumulator); - if ( accumulator.lsw & 0x80000000 ) - Xsq++; - - accumulator.msw = accumulator.midw = accumulator.lsw = 0; - /* Do the basic fixed point polynomial evaluation */ - polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER-1); - - mul_Xsig_Xsig(&accumulator, &argSignif); - shr_Xsig(&accumulator, 6 - adj); - - mul32_Xsig(&arg_signif, leadterm); - add_two_Xsig(&accumulator, &arg_signif, &exponent); - - *expon = exponent + 1; - accum_result->lsw = accumulator.lsw; - accum_result->midw = accumulator.midw; - accum_result->msw = accumulator.msw; + arg_signif.lsw = argSignif.lsw; + XSIG_LL(arg_signif) = XSIG_LL(argSignif); + adj = norm_Xsig(&argSignif); + accumulator.lsw = argSignif.lsw; + XSIG_LL(accumulator) = XSIG_LL(argSignif); + mul_Xsig_Xsig(&accumulator, &accumulator); + shr_Xsig(&accumulator, 2 * (-1 - (1 + exponent + adj))); + Xsq = XSIG_LL(accumulator); + if (accumulator.lsw & 0x80000000) + Xsq++; + + accumulator.msw = accumulator.midw = accumulator.lsw = 0; + /* Do the basic fixed point polynomial evaluation */ + polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER - 1); + + mul_Xsig_Xsig(&accumulator, &argSignif); + shr_Xsig(&accumulator, 6 - adj); + + mul32_Xsig(&arg_signif, leadterm); + add_two_Xsig(&accumulator, &arg_signif, &exponent); + + *expon = exponent + 1; + accum_result->lsw = accumulator.lsw; + accum_result->midw = accumulator.midw; + accum_result->msw = accumulator.msw; } |