4 * Compute the base 2 log of a FPU_REG, using a polynomial approximation.
7 * Copyright (C) 1992,1993,1994
8 * W. Metzenthen, 22 Parker St, Ormond, Vic 3163,
9 * Australia. E-mail billm@vaxc.cc.monash.edu.au
10 * All rights reserved.
12 * This copyright notice covers the redistribution and use of the
13 * FPU emulator developed by W. Metzenthen. It covers only its use
14 * in the 386BSD, FreeBSD and NetBSD operating systems. Any other
15 * use is not permitted under this copyright.
17 * Redistribution and use in source and binary forms, with or without
18 * modification, are permitted provided that the following conditions
20 * 1. Redistributions of source code must retain the above copyright
21 * notice, this list of conditions and the following disclaimer.
22 * 2. Redistributions in binary form must include information specifying
23 * that source code for the emulator is freely available and include
25 * a) an offer to provide the source code for a nominal distribution
27 * b) list at least two alternative methods whereby the source
28 * can be obtained, e.g. a publically accessible bulletin board
29 * and an anonymous ftp site from which the software can be
31 * 3. All advertising materials specifically mentioning features or use of
32 * this emulator must acknowledge that it was developed by W. Metzenthen.
33 * 4. The name of W. Metzenthen may not be used to endorse or promote
34 * products derived from this software without specific prior written
37 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
38 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
39 * AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
40 * W. METZENTHEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
41 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
42 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
43 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
44 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
45 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
46 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
49 * The purpose of this copyright, based upon the Berkeley copyright, is to
50 * ensure that the covered software remains freely available to everyone.
52 * The software (with necessary differences) is also available, but under
53 * the terms of the GNU copyleft, for the Linux operating system and for
54 * the djgpp ms-dos extender.
56 * W. Metzenthen June 1994.
59 * $FreeBSD: src/sys/gnu/i386/fpemul/poly_l2.c,v 1.10 1999/08/28 00:42:53 peter Exp $
60 * $DragonFly: src/sys/i386/gnu/fpemul/Attic/poly_l2.c,v 1.3 2003/08/07 21:17:20 dillon Exp $
65 #include "reg_constant.h"
66 #include "control_w.h"
71 static unsigned short lterms[HIPOWER][4] =
73 /* Ideal computation with these coeffs gives about 64.6 bit rel
75 {0xe177, 0xb82f, 0x7652, 0x7154},
76 {0xee0f, 0xe80f, 0x2770, 0x7b1c},
77 {0x0fc0, 0xbe87, 0xb143, 0x49dd},
78 {0x78b9, 0xdadd, 0xec54, 0x34c2},
79 {0x003a, 0x5de9, 0x628b, 0x2909},
80 {0x5588, 0xed16, 0x4abf, 0x2193},
81 {0xb461, 0x85f7, 0x347a, 0x1c6a},
82 {0x0975, 0x87b3, 0xd5bf, 0x1876},
83 {0xe85c, 0xcec9, 0x84e7, 0x187d}
89 /*--- poly_l2() -------------------------------------------------------------+
90 | Base 2 logarithm by a polynomial approximation. |
91 +---------------------------------------------------------------------------*/
93 poly_l2(FPU_REG * arg, FPU_REG * result)
96 char zero; /* flag for an Xx == 0 */
97 unsigned short bits, shift;
99 FPU_REG accum, denom, num, Xx;
102 exponent = arg->exp - EXP_BIAS;
104 accum.tag = TW_Valid; /* set the tags to Valid */
106 if (arg->sigh > (unsigned) 0xb504f334) {
107 /* This is good enough for the computation of the polynomial
108 * sum, but actually results in a loss of precision for the
109 * computation of Xx. This will matter only if exponent
112 accum.sign = 1; /* sign to negative */
113 num.exp = EXP_BIAS; /* needed to prevent errors in div
115 reg_u_div(&CONST_1, arg, &num, FULL_PRECISION);
117 accum.sign = 0; /* set the sign to positive */
118 num.sigl = arg->sigl; /* copy the mantissa */
119 num.sigh = arg->sigh;
123 /* shift num left, lose the ms bit */
125 if (num.sigl & 0x80000000)
129 denom.sigl = num.sigl;
130 denom.sigh = num.sigh;
131 poly_div4((long long *) &(denom.sigl));
132 denom.sigh += 0x80000000; /* set the msb */
133 Xx.exp = EXP_BIAS; /* needed to prevent errors in div routine */
134 reg_u_div(&num, &denom, &Xx, FULL_PRECISION);
136 zero = !(Xx.sigh | Xx.sigl);
138 mul64((long long *) &Xx.sigl, (long long *) &Xx.sigl, &Xsq);
141 accum.exp = -1; /* exponent of accum */
143 /* Do the basic fixed point polynomial evaluation */
144 polynomial((unsigned *) &accum.sigl, (unsigned *) &Xsq, lterms, HIPOWER - 1);
147 /* If the exponent is zero, then we would lose precision by
148 * sticking to fixed point computation here */
149 /* We need to re-compute Xx because of loss of precision. */
156 /* make accum compatible and normalize */
157 accum.exp = EXP_BIAS + accum.exp;
161 reg_move(&CONST_Z, result);
163 /* we need to re-compute lXx to better accuracy */
164 num.tag = TW_Valid; /* set the tags to Vaild */
165 num.sign = 0; /* set the sign to positive */
166 num.exp = EXP_BIAS - 1;
168 /* The argument is of the form 1-x */
169 /* Use 1-1/(1-x) = x/(1-x) */
170 *((long long *) &num.sigl) = -*((long long *) &(arg->sigl));
172 reg_div(&num, arg, &num, FULL_PRECISION);
177 denom.tag = TW_Valid; /* set the tags to Valid */
178 denom.sign = SIGN_POS; /* set the sign to positive */
179 denom.exp = EXP_BIAS;
181 reg_div(&num, &denom, &lXx, FULL_PRECISION);
183 reg_u_mul(&lXx, &accum, &accum, FULL_PRECISION);
185 reg_u_add(&lXx, &accum, result, FULL_PRECISION);
193 mul64((long long *) &accum.sigl,
194 (long long *) &Xx.sigl, (long long *) &accum.sigl);
196 *((long long *) (&accum.sigl)) += *((long long *) (&Xx.sigl));
198 if (Xx.sigh > accum.sigh) {
199 /* There was an overflow */
201 poly_div2((long long *) &accum.sigl);
202 accum.sigh |= 0x80000000;
205 /* When we add the exponent to the accum result later, we will require
206 * that their signs are the same. Here we ensure that this is so. */
207 if (exponent && ((exponent < 0) ^ (accum.sign))) {
208 /* signs are different */
210 accum.sign = !accum.sign;
212 /* An exceptional case is when accum is zero */
213 if (accum.sigl | accum.sigh) {
215 /* Shift to get exponent == 0 */
217 poly_div2((long long *) &accum.sigl);
220 /* Just negate, but throw away the sign */
221 *((long long *) &(accum.sigl)) = -*((long long *) &(accum.sigl));
228 shift = exponent >= 0 ? exponent : -exponent;
233 poly_div2((long long *) &accum.sigl);
236 poly_div2((long long *) &accum.sigl);
238 accum.sigh |= 0x80000000;
243 /* Convert to 64 bit signed-compatible */
244 accum.exp += bits + EXP_BIAS - 1;
246 reg_move(&accum, result);
253 /*--- poly_l2p1() -----------------------------------------------------------+
254 | Base 2 logarithm by a polynomial approximation. |
256 +---------------------------------------------------------------------------*/
258 poly_l2p1(FPU_REG * arg, FPU_REG * result)
262 FPU_REG arg_pl1, denom, accum, local_arg, poly_arg;
267 reg_add(arg, &CONST_1, &arg_pl1, FULL_PRECISION);
269 if ((arg_pl1.sign) | (arg_pl1.tag)) { /* We need a valid positive
273 reg_add(&CONST_1, &arg_pl1, &denom, FULL_PRECISION);
274 reg_div(arg, &denom, &local_arg, FULL_PRECISION);
275 local_arg.sign = 0; /* Make the sign positive */
277 /* Now we need to check that |local_arg| is less than 3-2*sqrt(2) =
278 * 0.17157.. = .0xafb0ccc0 * 2^-2 */
280 if (local_arg.exp >= EXP_BIAS - 3) {
281 if ((local_arg.exp > EXP_BIAS - 3) ||
282 (local_arg.sigh > (unsigned) 0xafb0ccc0)) {
283 /* The argument is large */
284 poly_l2(&arg_pl1, result);
288 /* Make a copy of local_arg */
289 reg_move(&local_arg, &poly_arg);
291 /* Get poly_arg bits aligned as required */
292 shrx((unsigned *) &(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));
294 mul64((long long *) &(poly_arg.sigl), (long long *) &(poly_arg.sigl), &Xsq);
297 /* Do the basic fixed point polynomial evaluation */
298 polynomial((u_int *) &accum.sigl, (unsigned *) &Xsq, lterms, HIPOWER - 1);
300 accum.tag = TW_Valid; /* set the tags to Valid */
301 accum.sign = SIGN_POS; /* and make accum positive */
303 /* make accum compatible and normalize */
304 accum.exp = EXP_BIAS - 1;
307 reg_u_mul(&local_arg, &accum, &accum, FULL_PRECISION);
309 reg_u_add(&local_arg, &accum, result, FULL_PRECISION);
311 /* Multiply the result by 2 */