/*---------------------------------------------------------------------------+
| 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;
}
