![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to poly_l2.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [MW Coherent from dump] / coherent / b / kernel / emulator|
Return to poly_l2.c CVS log | Up to [MW Coherent from dump] / coherent / b / kernel / emulator |
1.1 root 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 [email protected] |
8: | |
9: | |
10: +---------------------------------------------------------------------------*/
11:
12: #include "fpu_system.h"
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)
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)
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: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.