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