/* $OpenBSD: poly_l2.c,v 1.1 1996/08/27 10:32:54 downsj Exp $ */ /* * poly_l2.c * * Compute the base 2 log of a FPU_REG, using a polynomial approximation. * * * Copyright (C) 1992,1993,1994 * W. Metzenthen, 22 Parker St, Ormond, Vic 3163, * Australia. E-mail billm@vaxc.cc.monash.edu.au * All rights reserved. * * This copyright notice covers the redistribution and use of the * FPU emulator developed by W. Metzenthen. It covers only its use * in the 386BSD, FreeBSD and NetBSD operating systems. Any other * use is not permitted under this copyright. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must include information specifying * that source code for the emulator is freely available and include * either: * a) an offer to provide the source code for a nominal distribution * fee, or * b) list at least two alternative methods whereby the source * can be obtained, e.g. a publically accessible bulletin board * and an anonymous ftp site from which the software can be * downloaded. * 3. All advertising materials specifically mentioning features or use of * this emulator must acknowledge that it was developed by W. Metzenthen. * 4. The name of W. Metzenthen may not be used to endorse or promote * products derived from this software without specific prior written * permission. * * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY * AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL * W. METZENTHEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * * * The purpose of this copyright, based upon the Berkeley copyright, is to * ensure that the covered software remains freely available to everyone. * * The software (with necessary differences) is also available, but under * the terms of the GNU copyleft, for the Linux operating system and for * the djgpp ms-dos extender. * * W. Metzenthen June 1994. * * * $FreeBSD: poly_l2.c,v 1.5 1994/06/10 07:44:38 rich Exp $ * */ #include #include #include #include #define HIPOWER 9 static unsigned short 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 * arg, FPU_REG * result) { short exponent; char zero; /* flag for an Xx == 0 */ unsigned short bits, shift; 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((long long *) &(denom.sigl)); 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((long long *) &Xx.sigl, (long long *) &Xx.sigl, &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) */ *((long long *) &num.sigl) = -*((long long *) &(arg->sigl)); 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((long long *) &accum.sigl, (long long *) &Xx.sigl, (long long *) &accum.sigl); *((long long *) (&accum.sigl)) += *((long long *) (&Xx.sigl)); if (Xx.sigh > accum.sigh) { /* There was an overflow */ poly_div2((long long *) &accum.sigl); 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((long long *) &accum.sigl); accum.exp++; } /* Just negate, but throw away the sign */ *((long long *) &(accum.sigl)) = -*((long long *) &(accum.sigl)); if (exponent < 0) exponent++; else exponent--; } } shift = exponent >= 0 ? exponent : -exponent; bits = 0; if (shift) { if (accum.exp) { accum.exp++; poly_div2((long long *) &accum.sigl); } while (shift) { poly_div2((long long *) &accum.sigl); 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 * arg, FPU_REG * result) { char sign = 0; 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((long long *) &(poly_arg.sigl), (long long *) &(poly_arg.sigl), &Xsq); poly_div16(&Xsq); /* Do the basic fixed point polynomial evaluation */ polynomial((u_int *) &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; }