289 lines
7.8 KiB
C
289 lines
7.8 KiB
C
/*---------------------------------------------------------------------------+
|
|
| poly_l2.c |
|
|
| |
|
|
| Compute the base 2 log of a FPU_REG, using a polynomial approximation. |
|
|
| |
|
|
| Copyright (C) 1992,1993 |
|
|
| W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
|
|
| Australia. E-mail billm@vaxc.cc.monash.edu.au |
|
|
| |
|
|
| |
|
|
+---------------------------------------------------------------------------*/
|
|
|
|
|
|
#include "exception.h"
|
|
#include "reg_constant.h"
|
|
#include "fpu_emu.h"
|
|
#include "control_w.h"
|
|
|
|
|
|
|
|
#define HIPOWER 9
|
|
static unsigned short const lterms[HIPOWER][4] =
|
|
{
|
|
/* Ideal computation with these coeffs gives about
|
|
64.6 bit rel accuracy. */
|
|
{ 0xe177, 0xb82f, 0x7652, 0x7154 },
|
|
{ 0xee0f, 0xe80f, 0x2770, 0x7b1c },
|
|
{ 0x0fc0, 0xbe87, 0xb143, 0x49dd },
|
|
{ 0x78b9, 0xdadd, 0xec54, 0x34c2 },
|
|
{ 0x003a, 0x5de9, 0x628b, 0x2909 },
|
|
{ 0x5588, 0xed16, 0x4abf, 0x2193 },
|
|
{ 0xb461, 0x85f7, 0x347a, 0x1c6a },
|
|
{ 0x0975, 0x87b3, 0xd5bf, 0x1876 },
|
|
{ 0xe85c, 0xcec9, 0x84e7, 0x187d }
|
|
};
|
|
|
|
|
|
|
|
|
|
/*--- poly_l2() -------------------------------------------------------------+
|
|
| Base 2 logarithm by a polynomial approximation. |
|
|
+---------------------------------------------------------------------------*/
|
|
void poly_l2(FPU_REG const *arg, FPU_REG *result)
|
|
{
|
|
short exponent;
|
|
char zero; /* flag for an Xx == 0 */
|
|
unsigned short bits, shift;
|
|
unsigned long long Xsq;
|
|
FPU_REG accum, denom, num, Xx;
|
|
|
|
|
|
exponent = arg->exp - EXP_BIAS;
|
|
|
|
accum.tag = TW_Valid; /* set the tags to Valid */
|
|
|
|
if ( arg->sigh > (unsigned)0xb504f334 )
|
|
{
|
|
/* This is good enough for the computation of the polynomial
|
|
sum, but actually results in a loss of precision for
|
|
the computation of Xx. This will matter only if exponent
|
|
becomes zero. */
|
|
exponent++;
|
|
accum.sign = 1; /* sign to negative */
|
|
num.exp = EXP_BIAS; /* needed to prevent errors in div routine */
|
|
reg_u_div(&CONST_1, arg, &num, FULL_PRECISION);
|
|
}
|
|
else
|
|
{
|
|
accum.sign = 0; /* set the sign to positive */
|
|
num.sigl = arg->sigl; /* copy the mantissa */
|
|
num.sigh = arg->sigh;
|
|
}
|
|
|
|
|
|
/* shift num left, lose the ms bit */
|
|
num.sigh <<= 1;
|
|
if ( num.sigl & 0x80000000 ) num.sigh |= 1;
|
|
num.sigl <<= 1;
|
|
|
|
denom.sigl = num.sigl;
|
|
denom.sigh = num.sigh;
|
|
poly_div4(&significand(&denom));
|
|
denom.sigh += 0x80000000; /* set the msb */
|
|
Xx.exp = EXP_BIAS; /* needed to prevent errors in div routine */
|
|
reg_u_div(&num, &denom, &Xx, FULL_PRECISION);
|
|
|
|
zero = !(Xx.sigh | Xx.sigl);
|
|
|
|
mul64(&significand(&Xx), &significand(&Xx), &Xsq);
|
|
poly_div16(&Xsq);
|
|
|
|
accum.exp = -1; /* exponent of accum */
|
|
|
|
/* Do the basic fixed point polynomial evaluation */
|
|
polynomial((unsigned *)&accum.sigl, (unsigned *)&Xsq, lterms, HIPOWER-1);
|
|
|
|
if ( !exponent )
|
|
{
|
|
/* If the exponent is zero, then we would lose precision by
|
|
sticking to fixed point computation here */
|
|
/* We need to re-compute Xx because of loss of precision. */
|
|
FPU_REG lXx;
|
|
char sign;
|
|
|
|
sign = accum.sign;
|
|
accum.sign = 0;
|
|
|
|
/* make accum compatible and normalize */
|
|
accum.exp = EXP_BIAS + accum.exp;
|
|
normalize(&accum);
|
|
|
|
if ( zero )
|
|
{
|
|
reg_move(&CONST_Z, result);
|
|
}
|
|
else
|
|
{
|
|
/* we need to re-compute lXx to better accuracy */
|
|
num.tag = TW_Valid; /* set the tags to Vaild */
|
|
num.sign = 0; /* set the sign to positive */
|
|
num.exp = EXP_BIAS - 1;
|
|
if ( sign )
|
|
{
|
|
/* The argument is of the form 1-x */
|
|
/* Use 1-1/(1-x) = x/(1-x) */
|
|
significand(&num) = - significand(arg);
|
|
normalize(&num);
|
|
reg_div(&num, arg, &num, FULL_PRECISION);
|
|
}
|
|
else
|
|
{
|
|
normalize(&num);
|
|
}
|
|
|
|
denom.tag = TW_Valid; /* set the tags to Valid */
|
|
denom.sign = SIGN_POS; /* set the sign to positive */
|
|
denom.exp = EXP_BIAS;
|
|
|
|
reg_div(&num, &denom, &lXx, FULL_PRECISION);
|
|
|
|
reg_u_mul(&lXx, &accum, &accum, FULL_PRECISION);
|
|
|
|
reg_u_add(&lXx, &accum, result, FULL_PRECISION);
|
|
|
|
normalize(result);
|
|
}
|
|
|
|
result->sign = sign;
|
|
return;
|
|
}
|
|
|
|
mul64(&significand(&accum),
|
|
&significand(&Xx), &significand(&accum));
|
|
|
|
significand(&accum) += significand(&Xx);
|
|
|
|
if ( Xx.sigh > accum.sigh )
|
|
{
|
|
/* There was an overflow */
|
|
|
|
poly_div2(&significand(&accum));
|
|
accum.sigh |= 0x80000000;
|
|
accum.exp++;
|
|
}
|
|
|
|
/* When we add the exponent to the accum result later, we will
|
|
require that their signs are the same. Here we ensure that
|
|
this is so. */
|
|
if ( exponent && ((exponent < 0) ^ (accum.sign)) )
|
|
{
|
|
/* signs are different */
|
|
|
|
accum.sign = !accum.sign;
|
|
|
|
/* An exceptional case is when accum is zero */
|
|
if ( accum.sigl | accum.sigh )
|
|
{
|
|
/* find 1-accum */
|
|
/* Shift to get exponent == 0 */
|
|
if ( accum.exp < 0 )
|
|
{
|
|
poly_div2(&significand(&accum));
|
|
accum.exp++;
|
|
}
|
|
/* Just negate, but throw away the sign */
|
|
significand(&accum) = - significand(&accum);
|
|
if ( exponent < 0 )
|
|
exponent++;
|
|
else
|
|
exponent--;
|
|
}
|
|
}
|
|
|
|
shift = exponent >= 0 ? exponent : -exponent ;
|
|
bits = 0;
|
|
if ( shift )
|
|
{
|
|
if ( accum.exp )
|
|
{
|
|
accum.exp++;
|
|
poly_div2(&significand(&accum));
|
|
}
|
|
while ( shift )
|
|
{
|
|
poly_div2(&significand(&accum));
|
|
if ( shift & 1)
|
|
accum.sigh |= 0x80000000;
|
|
shift >>= 1;
|
|
bits++;
|
|
}
|
|
}
|
|
|
|
/* Convert to 64 bit signed-compatible */
|
|
accum.exp += bits + EXP_BIAS - 1;
|
|
|
|
reg_move(&accum, result);
|
|
normalize(result);
|
|
|
|
return;
|
|
}
|
|
|
|
|
|
/*--- poly_l2p1() -----------------------------------------------------------+
|
|
| Base 2 logarithm by a polynomial approximation. |
|
|
| log2(x+1) |
|
|
+---------------------------------------------------------------------------*/
|
|
int poly_l2p1(FPU_REG const *arg, FPU_REG *result)
|
|
{
|
|
char sign = 0;
|
|
unsigned long long Xsq;
|
|
FPU_REG arg_pl1, denom, accum, local_arg, poly_arg;
|
|
|
|
|
|
sign = arg->sign;
|
|
|
|
reg_add(arg, &CONST_1, &arg_pl1, FULL_PRECISION);
|
|
|
|
if ( (arg_pl1.sign) | (arg_pl1.tag) )
|
|
{ /* We need a valid positive number! */
|
|
return 1;
|
|
}
|
|
|
|
reg_add(&CONST_1, &arg_pl1, &denom, FULL_PRECISION);
|
|
reg_div(arg, &denom, &local_arg, FULL_PRECISION);
|
|
local_arg.sign = 0; /* Make the sign positive */
|
|
|
|
/* Now we need to check that |local_arg| is less than
|
|
3-2*sqrt(2) = 0.17157.. = .0xafb0ccc0 * 2^-2 */
|
|
|
|
if ( local_arg.exp >= EXP_BIAS - 3 )
|
|
{
|
|
if ( (local_arg.exp > EXP_BIAS - 3) ||
|
|
(local_arg.sigh > (unsigned)0xafb0ccc0) )
|
|
{
|
|
/* The argument is large */
|
|
poly_l2(&arg_pl1, result); return 0;
|
|
}
|
|
}
|
|
|
|
/* Make a copy of local_arg */
|
|
reg_move(&local_arg, &poly_arg);
|
|
|
|
/* Get poly_arg bits aligned as required */
|
|
shrx((unsigned *)&(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));
|
|
|
|
mul64(&significand(&poly_arg), &significand(&poly_arg), &Xsq);
|
|
poly_div16(&Xsq);
|
|
|
|
/* Do the basic fixed point polynomial evaluation */
|
|
polynomial(&(accum.sigl), (unsigned *)&Xsq, lterms, HIPOWER-1);
|
|
|
|
accum.tag = TW_Valid; /* set the tags to Valid */
|
|
accum.sign = SIGN_POS; /* and make accum positive */
|
|
|
|
/* make accum compatible and normalize */
|
|
accum.exp = EXP_BIAS - 1;
|
|
normalize(&accum);
|
|
|
|
reg_u_mul(&local_arg, &accum, &accum, FULL_PRECISION);
|
|
|
|
reg_u_add(&local_arg, &accum, result, FULL_PRECISION);
|
|
|
|
/* Multiply the result by 2 */
|
|
result->exp++;
|
|
|
|
result->sign = sign;
|
|
|
|
return 0;
|
|
}
|