diff options
Diffstat (limited to 'arch/x86/math-emu/poly_l2.c')
| -rw-r--r-- | arch/x86/math-emu/poly_l2.c | 272 |
1 files changed, 272 insertions, 0 deletions
diff --git a/arch/x86/math-emu/poly_l2.c b/arch/x86/math-emu/poly_l2.c new file mode 100644 index 00000000000..dd00e1d5b07 --- /dev/null +++ b/arch/x86/math-emu/poly_l2.c @@ -0,0 +1,272 @@ +/*---------------------------------------------------------------------------+ + | poly_l2.c | + | | + | Compute the base 2 log of a FPU_REG, using a polynomial approximation. | + | | + | Copyright (C) 1992,1993,1994,1997 | + | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia | + | E-mail billm@suburbia.net | + | | + | | + +---------------------------------------------------------------------------*/ + + +#include "exception.h" +#include "reg_constant.h" +#include "fpu_emu.h" +#include "fpu_system.h" +#include "control_w.h" +#include "poly.h" + + +static void log2_kernel(FPU_REG const *arg, u_char argsign, + 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) +{ + 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) +{ + u_char tag; + long int exponent; + Xsig accumulator, yaccum; + + 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); + + exponent += round_Xsig(&accumulator); + + exponent += exponent16(st1_ptr) + 1; + if ( exponent < EXP_WAY_UNDER ) exponent = EXP_WAY_UNDER; + + significand(dest) = XSIG_LL(accumulator); + setexponent16(dest, exponent); + + 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 ( 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; +#endif /* PECULIAR_486 */ + } + + /* 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); + + 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 leadterm = 0xb8000000; + + +/*--- log2_kernel() ---------------------------------------------------------+ + | Base 2 logarithm by a polynomial approximation. | + | log2(x+1) | + +---------------------------------------------------------------------------*/ +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; + } + } + +#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 */ + } + } +#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; + +} |