/*---------------------------------------------------------------------------+ | poly_l2.c | | $Id: poly_l2.c,v 1.4 2001-10-06 03:53:46 bdenney Exp $ | | | 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, s32 *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) { s32 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, 0); 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; s32 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 u64 logterms[HIPOWER] = { BX_CONST64(0x2a8eca5705fc2ef0), BX_CONST64(0xf6384ee1d01febce), BX_CONST64(0x093bb62877cdf642), BX_CONST64(0x006985d8a9ec439b), BX_CONST64(0x0005212c4f55a9c8), BX_CONST64(0x00004326a16927f0), BX_CONST64(0x0000038d1d80a0e7), BX_CONST64(0x0000003141cc80c6), BX_CONST64(0x00000002b1668c9f), BX_CONST64(0x000000002c7a46aa) }; static const u32 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, s32 *expon) { s32 exponent, adj; u64 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; }