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