1 /*
2 * poly_l2.c
3 *
4 * Compute the base 2 log of a FPU_REG, using a polynomial approximation.
5 *
6 *
7 * Copyright (C) 1992, 1993 W. Metzenthen, 22 Parker St, Ormond,
8 * Vic 3163, Australia.
9 * E-mail apm233m@vaxc.cc.monash.edu.au
10 * All rights reserved.
11 *
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 operating system. Any other use is not permitted
15 * under this copyright.
16 *
17 * Redistribution and use in source and binary forms, with or without
18 * modification, are permitted provided that the following conditions
19 * are met:
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
24 * either:
25 * a) an offer to provide the source code for a nominal distribution
26 * fee, or
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
30 * downloaded.
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
35 * permission.
36 *
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.
47 *
48 */
49
50
51 #include "exception.h"
52 #include "reg_constant.h"
53 #include "fpu_emu.h"
54 #include "control_w.h"
55
56
57
58 #define HIPOWER 9
59 static unsigned short lterms[HIPOWER][4] =
60 {
61 /* Ideal computation with these coeffs gives about
62 64.6 bit rel accuracy. */
63 { 0xe177, 0xb82f, 0x7652, 0x7154 },
64 { 0xee0f, 0xe80f, 0x2770, 0x7b1c },
65 { 0x0fc0, 0xbe87, 0xb143, 0x49dd },
66 { 0x78b9, 0xdadd, 0xec54, 0x34c2 },
67 { 0x003a, 0x5de9, 0x628b, 0x2909 },
68 { 0x5588, 0xed16, 0x4abf, 0x2193 },
69 { 0xb461, 0x85f7, 0x347a, 0x1c6a },
70 { 0x0975, 0x87b3, 0xd5bf, 0x1876 },
71 { 0xe85c, 0xcec9, 0x84e7, 0x187d }
72 };
73
74
75
76
77 /*--- poly_l2() -------------------------------------------------------------+
78 | Base 2 logarithm by a polynomial approximation. |
79 +---------------------------------------------------------------------------*/
poly_l2(FPU_REG * arg,FPU_REG * result)80 void poly_l2(FPU_REG *arg, FPU_REG *result)
81 {
82 short exponent;
83 char zero; /* flag for an Xx == 0 */
84 unsigned short bits, shift;
85 long long Xsq;
86 FPU_REG accum, denom, num, Xx;
87
88
89 exponent = arg->exp - EXP_BIAS;
90
91 accum.tag = TW_Valid; /* set the tags to Valid */
92
93 if ( arg->sigh > (unsigned)0xb504f334 )
94 {
95 /* This is good enough for the computation of the polynomial
96 sum, but actually results in a loss of precision for
97 the computation of Xx. This will matter only if exponent
98 becomes zero. */
99 exponent++;
100 accum.sign = 1; /* sign to negative */
101 num.exp = EXP_BIAS; /* needed to prevent errors in div routine */
102 reg_u_div(&CONST_1, arg, &num, FULL_PRECISION);
103 }
104 else
105 {
106 accum.sign = 0; /* set the sign to positive */
107 num.sigl = arg->sigl; /* copy the mantissa */
108 num.sigh = arg->sigh;
109 }
110
111
112 /* shift num left, lose the ms bit */
113 num.sigh <<= 1;
114 if ( num.sigl & 0x80000000 ) num.sigh |= 1;
115 num.sigl <<= 1;
116
117 denom.sigl = num.sigl;
118 denom.sigh = num.sigh;
119 poly_div4((long long *)&(denom.sigl));
120 denom.sigh += 0x80000000; /* set the msb */
121 Xx.exp = EXP_BIAS; /* needed to prevent errors in div routine */
122 reg_u_div(&num, &denom, &Xx, FULL_PRECISION);
123
124 zero = !(Xx.sigh | Xx.sigl);
125
126 mul64((long long *)&Xx.sigl, (long long *)&Xx.sigl, &Xsq);
127 poly_div16(&Xsq);
128
129 accum.exp = -1; /* exponent of accum */
130
131 /* Do the basic fixed point polynomial evaluation */
132 polynomial((unsigned *)&accum.sigl, (unsigned *)&Xsq, lterms, HIPOWER-1);
133
134 if ( !exponent )
135 {
136 /* If the exponent is zero, then we would lose precision by
137 sticking to fixed point computation here */
138 /* We need to re-compute Xx because of loss of precision. */
139 FPU_REG lXx;
140 char sign;
141
142 sign = accum.sign;
143 accum.sign = 0;
144
145 /* make accum compatible and normalize */
146 accum.exp = EXP_BIAS + accum.exp;
147 normalize(&accum);
148
149 if ( zero )
150 {
151 reg_move(&CONST_Z, result);
152 }
153 else
154 {
155 /* we need to re-compute lXx to better accuracy */
156 num.tag = TW_Valid; /* set the tags to Vaild */
157 num.sign = 0; /* set the sign to positive */
158 num.exp = EXP_BIAS - 1;
159 if ( sign )
160 {
161 /* The argument is of the form 1-x */
162 /* Use 1-1/(1-x) = x/(1-x) */
163 *((long long *)&num.sigl) = - *((long long *)&(arg->sigl));
164 normalize(&num);
165 reg_div(&num, arg, &num, FULL_PRECISION);
166 }
167 else
168 {
169 normalize(&num);
170 }
171
172 denom.tag = TW_Valid; /* set the tags to Valid */
173 denom.sign = SIGN_POS; /* set the sign to positive */
174 denom.exp = EXP_BIAS;
175
176 reg_div(&num, &denom, &lXx, FULL_PRECISION);
177
178 reg_u_mul(&lXx, &accum, &accum, FULL_PRECISION);
179
180 reg_u_add(&lXx, &accum, result, FULL_PRECISION);
181
182 normalize(result);
183 }
184
185 result->sign = sign;
186 return;
187 }
188
189 mul64((long long *)&accum.sigl,
190 (long long *)&Xx.sigl, (long long *)&accum.sigl);
191
192 *((long long *)(&accum.sigl)) += *((long long *)(&Xx.sigl));
193
194 if ( Xx.sigh > accum.sigh )
195 {
196 /* There was an overflow */
197
198 poly_div2((long long *)&accum.sigl);
199 accum.sigh |= 0x80000000;
200 accum.exp++;
201 }
202
203 /* When we add the exponent to the accum result later, we will
204 require that their signs are the same. Here we ensure that
205 this is so. */
206 if ( exponent && ((exponent < 0) ^ (accum.sign)) )
207 {
208 /* signs are different */
209
210 accum.sign = !accum.sign;
211
212 /* An exceptional case is when accum is zero */
213 if ( accum.sigl | accum.sigh )
214 {
215 /* find 1-accum */
216 /* Shift to get exponent == 0 */
217 if ( accum.exp < 0 )
218 {
219 poly_div2((long long *)&accum.sigl);
220 accum.exp++;
221 }
222 /* Just negate, but throw away the sign */
223 *((long long *)&(accum.sigl)) = - *((long long *)&(accum.sigl));
224 if ( exponent < 0 )
225 exponent++;
226 else
227 exponent--;
228 }
229 }
230
231 shift = exponent >= 0 ? exponent : -exponent ;
232 bits = 0;
233 if ( shift )
234 {
235 if ( accum.exp )
236 {
237 accum.exp++;
238 poly_div2((long long *)&accum.sigl);
239 }
240 while ( shift )
241 {
242 poly_div2((long long *)&accum.sigl);
243 if ( shift & 1)
244 accum.sigh |= 0x80000000;
245 shift >>= 1;
246 bits++;
247 }
248 }
249
250 /* Convert to 64 bit signed-compatible */
251 accum.exp += bits + EXP_BIAS - 1;
252
253 reg_move(&accum, result);
254 normalize(result);
255
256 return;
257 }
258
259
260 /*--- poly_l2p1() -----------------------------------------------------------+
261 | Base 2 logarithm by a polynomial approximation. |
262 | log2(x+1) |
263 +---------------------------------------------------------------------------*/
poly_l2p1(FPU_REG * arg,FPU_REG * result)264 int poly_l2p1(FPU_REG *arg, FPU_REG *result)
265 {
266 char sign = 0;
267 long long Xsq;
268 FPU_REG arg_pl1, denom, accum, local_arg, poly_arg;
269
270
271 sign = arg->sign;
272
273 reg_add(arg, &CONST_1, &arg_pl1, FULL_PRECISION);
274
275 if ( (arg_pl1.sign) | (arg_pl1.tag) )
276 { /* We need a valid positive number! */
277 return 1;
278 }
279
280 reg_add(&CONST_1, &arg_pl1, &denom, FULL_PRECISION);
281 reg_div(arg, &denom, &local_arg, FULL_PRECISION);
282 local_arg.sign = 0; /* Make the sign positive */
283
284 /* Now we need to check that |local_arg| is less than
285 3-2*sqrt(2) = 0.17157.. = .0xafb0ccc0 * 2^-2 */
286
287 if ( local_arg.exp >= EXP_BIAS - 3 )
288 {
289 if ( (local_arg.exp > EXP_BIAS - 3) ||
290 (local_arg.sigh > (unsigned)0xafb0ccc0) )
291 {
292 /* The argument is large */
293 poly_l2(&arg_pl1, result); return 0;
294 }
295 }
296
297 /* Make a copy of local_arg */
298 reg_move(&local_arg, &poly_arg);
299
300 /* Get poly_arg bits aligned as required */
301 shrx((unsigned *)&(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));
302
303 mul64((long long *)&(poly_arg.sigl), (long long *)&(poly_arg.sigl), &Xsq);
304 poly_div16(&Xsq);
305
306 /* Do the basic fixed point polynomial evaluation */
307 polynomial(&(accum.sigl), (unsigned *)&Xsq, lterms, HIPOWER-1);
308
309 accum.tag = TW_Valid; /* set the tags to Valid */
310 accum.sign = SIGN_POS; /* and make accum positive */
311
312 /* make accum compatible and normalize */
313 accum.exp = EXP_BIAS - 1;
314 normalize(&accum);
315
316 reg_u_mul(&local_arg, &accum, &accum, FULL_PRECISION);
317
318 reg_u_add(&local_arg, &accum, result, FULL_PRECISION);
319
320 /* Multiply the result by 2 */
321 result->exp++;
322
323 result->sign = sign;
324
325 return 0;
326 }
327