home *** CD-ROM | disk | FTP | other *** search
- /*---------------------------------------------------------------------------+
- | 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;
- }
-