Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
/*
2 
 * Floating Point Interpreter.
3 
 * shamelessly stolen from an original by ark.
4 
 *
5 
 * NB: the Internal arguments to fpisub and fpidiv are reversed from
6 
 * what you might naively expect: they compute y-x and y/x, respectively.
7 
 */
8 
#include "fpi.h"
9 
 
10 
void
11 
fpiround(Internal *i)
12 
{
13 
	unsigned long guard;
14 
 
15 
	guard = i->l & GuardMask;
16 
	i->l &= ~GuardMask;
17 
	if(guard > (LsBit>>1) || (guard == (LsBit>>1) && (i->l & LsBit))){
18 
		i->l += LsBit;
19 
		if(i->l & CarryBit){
20 
			i->l &= ~CarryBit;
21 
			i->h++;
22 
			if(i->h & CarryBit){
23 
				if (i->h & 0x01)
24 
					i->l |= CarryBit;
25 
				i->l >>= 1;
26 
				i->h >>= 1;
27 
				i->e++;
28 
			}
29 
		}
30 
	}
31 
}
32 
 
33 
static void
34 
matchexponents(Internal *x, Internal *y)
35 
{
36 
	int count;
37 
 
38 
	count = y->e - x->e;
39 
	x->e = y->e;
40 
	if(count >= 2*FractBits){
41 
		x->l = x->l || x->h;
42 
		x->h = 0;
43 
		return;
44 
	}
45 
	if(count >= FractBits){
46 
		count -= FractBits;
47 
		x->l = x->h|(x->l != 0);
48 
		x->h = 0;
49 
	}
50 
	while(count > 0){
51 
		count--;
52 
		if(x->h & 0x01)
53 
			x->l |= CarryBit;
54 
		if(x->l & 0x01)
55 
			x->l |= 2;
56 
		x->l >>= 1;
57 
		x->h >>= 1;
58 
	}
59 
}
60 
 
61 
static void
62 
shift(Internal *i)
63 
{
64 
	i->e--;
65 
	i->h <<= 1;
66 
	i->l <<= 1;
67 
	if(i->l & CarryBit){
68 
		i->l &= ~CarryBit;
69 
		i->h |= 0x01;
70 
	}
71 
}
72 
 
73 
static void
74 
normalise(Internal *i)
75 
{
76 
	while((i->h & HiddenBit) == 0)
77 
		shift(i);
78 
}
79 
 
80 
static void
81 
renormalise(Internal *i)
82 
{
83 
	if(i->e < -2 * FractBits)
84 
		i->e = -2 * FractBits;
85 
	while(i->e < 1){
86 
		i->e++;
87 
		if(i->h & 0x01)
88 
			i->l |= CarryBit;
89 
		i->h >>= 1;
90 
		i->l = (i->l>>1)|(i->l & 0x01);
91 
	}
92 
	if(i->e >= ExpInfinity)
93 
		SetInfinity(i);
94 
}
95 
 
96 
void
97 
fpinormalise(Internal *x)
98 
{
99 
	if(!IsWeird(x) && !IsZero(x))
100 
		normalise(x);
101 
}
102 
 
103 
void
104 
fpiadd(Internal *x, Internal *y, Internal *i)
105 
{
106 
	Internal *t;
107 
 
108 
	i->s = x->s;
109 
	if(IsWeird(x) || IsWeird(y)){
110 
		if(IsNaN(x) || IsNaN(y))
111 
			SetQNaN(i);
112 
		else
113 
			SetInfinity(i);
114 
		return;
115 
	}
116 
	if(x->e > y->e){
117 
		t = x;
118 
		x = y;
119 
		y = t;
120 
	}
121 
	matchexponents(x, y);
122 
	i->e = x->e;
123 
	i->h = x->h + y->h;
124 
	i->l = x->l + y->l;
125 
	if(i->l & CarryBit){
126 
		i->h++;
127 
		i->l &= ~CarryBit;
128 
	}
129 
	if(i->h & (HiddenBit<<1)){
130 
		if(i->h & 0x01)
131 
			i->l |= CarryBit;
132 
		i->l = (i->l>>1)|(i->l & 0x01);
133 
		i->h >>= 1;
134 
		i->e++;
135 
	}
136 
	if(IsWeird(i))
137 
		SetInfinity(i);
138 
}
139 
 
140 
void
141 
fpisub(Internal *x, Internal *y, Internal *i)
142 
{
143 
	Internal *t;
144 
 
145 
	if(y->e < x->e
146 
	   || (y->e == x->e && (y->h < x->h || (y->h == x->h && y->l < x->l)))){
147 
		t = x;
148 
		x = y;
149 
		y = t;
150 
	}
151 
	i->s = y->s;
152 
	if(IsNaN(y)){
153 
		SetQNaN(i);
154 
		return;
155 
	}
156 
	if(IsInfinity(y)){
157 
		if(IsInfinity(x))
158 
			SetQNaN(i);
159 
		else
160 
			SetInfinity(i);
161 
		return;
162 
	}
163 
	matchexponents(x, y);
164 
	i->e = y->e;
165 
	i->h = y->h - x->h;
166 
	i->l = y->l - x->l;
167 
	if(i->l < 0){
168 
		i->l += CarryBit;
169 
		i->h--;
170 
	}
171 
	if(i->h == 0 && i->l == 0)
172 
		SetZero(i);
173 
	else while(i->e > 1 && (i->h & HiddenBit) == 0)
174 
		shift(i);
175 
}
176 
 
177 
#define	CHUNK		(FractBits/2)
178 
#define	CMASK		((1<<CHUNK)-1)
179 
#define	HI(x)		((short)((x)>>CHUNK) & CMASK)
180 
#define	LO(x)		((short)(x) & CMASK)
181 
#define	SPILL(x)	((x)>>CHUNK)
182 
#define	M(x, y)		((long)a[x]*(long)b[y])
183 
#define	C(h, l)		(((long)((h) & CMASK)<<CHUNK)|((l) & CMASK))
184 
 
185 
void
186 
fpimul(Internal *x, Internal *y, Internal *i)
187 
{
188 
	long a[4], b[4], c[7], f[4];
189 
 
190 
	i->s = x->s^y->s;
191 
	if(IsWeird(x) || IsWeird(y)){
192 
		if(IsNaN(x) || IsNaN(y) || IsZero(x) || IsZero(y))
193 
			SetQNaN(i);
194 
		else
195 
			SetInfinity(i);
196 
		return;
197 
	}
198 
	else if(IsZero(x) || IsZero(y)){
199 
		SetZero(i);
200 
		return;
201 
	}
202 
	normalise(x);
203 
	normalise(y);
204 
	i->e = x->e + y->e - (ExpBias - 1);
205 
 
206 
	a[0] = HI(x->h); b[0] = HI(y->h);
207 
	a[1] = LO(x->h); b[1] = LO(y->h);
208 
	a[2] = HI(x->l); b[2] = HI(y->l);
209 
	a[3] = LO(x->l); b[3] = LO(y->l);
210 
 
211 
	c[6] =                               M(3, 3);
212 
	c[5] =                     M(2, 3) + M(3, 2) + SPILL(c[6]);
213 
	c[4] =           M(1, 3) + M(2, 2) + M(3, 1) + SPILL(c[5]);
214 
	c[3] = M(0, 3) + M(1, 2) + M(2, 1) + M(3, 0) + SPILL(c[4]);
215 
	c[2] = M(0, 2) + M(1, 1) + M(2, 0)           + SPILL(c[3]);
216 
	c[1] = M(0, 1) + M(1, 0)                     + SPILL(c[2]);
217 
	c[0] = M(0, 0)                               + SPILL(c[1]);
218 
 
219 
	f[0] = c[0];
220 
	f[1] = C(c[1], c[2]);
221 
	f[2] = C(c[3], c[4]);
222 
	f[3] = C(c[5], c[6]);
223 
 
224 
	if((f[0] & HiddenBit) == 0){
225 
		f[0] <<= 1;
226 
		f[1] <<= 1;
227 
		f[2] <<= 1;
228 
		f[3] <<= 1;
229 
		if(f[1] & CarryBit){
230 
			f[0] |= 1;
231 
			f[1] &= ~CarryBit;
232 
		}
233 
		if(f[2] & CarryBit){
234 
			f[1] |= 1;
235 
			f[2] &= ~CarryBit;
236 
		}
237 
		if(f[3] & CarryBit){
238 
			f[2] |= 1;
239 
			f[3] &= ~CarryBit;
240 
		}
241 
		i->e--;
242 
	}
243 
	i->h = f[0];
244 
	i->l = f[1];
245 
	if(f[2] || f[3])
246 
		i->l |= 1;
247 
	renormalise(i);
248 
}
249 
 
250 
void
251 
fpidiv(Internal *x, Internal *y, Internal *i)
252 
{
253 
	i->s = x->s^y->s;
254 
	if(IsNaN(x) || IsNaN(y)
255 
	   || (IsInfinity(x) && IsInfinity(y)) || (IsZero(x) && IsZero(y))){
256 
		SetQNaN(i);
257 
		return;
258 
	}
259 
	else if(IsZero(x) || IsInfinity(y)){
260 
		SetInfinity(i);
261 
		return;
262 
	}
263 
	else if(IsInfinity(x) || IsZero(y)){
264 
		SetZero(i);
265 
		return;
266 
	}
267 
	normalise(x);
268 
	normalise(y);
269 
	i->h = 0;
270 
	i->l = 0;
271 
	i->e = y->e - x->e + (ExpBias + 2*FractBits - 1);
272 
	do{
273 
		if(y->h > x->h || (y->h == x->h && y->l >= x->l)){
274 
			i->l |= 0x01;
275 
			y->h -= x->h;
276 
			y->l -= x->l;
277 
			if(y->l < 0){
278 
				y->l += CarryBit;
279 
				y->h--;
280 
			}
281 
		}
282 
		shift(y);
283 
		shift(i);
284 
	}while ((i->h & HiddenBit) == 0);
285 
	if(y->h || y->l)
286 
		i->l |= 0x01;
287 
	renormalise(i);
288 
}
289 
 
290 
int
291 
fpicmp(Internal *x, Internal *y)
292 
{
293 
	if(IsNaN(x) && IsNaN(y))
294 
		return 0;
295 
	if(IsInfinity(x) && IsInfinity(y))
296 
		return y->s - x->s;
297 
	if(IsZero(x) && IsZero(y))
298 
		return 0;
299 
	if(x->e == y->e && x->h == y->h && x->l == y->l)
300 
		return y->s - x->s;
301 
	if(x->e < y->e
302 
	   || (x->e == y->e && (x->h < y->h || (x->h == y->h && x->l < y->l))))
303 
		return y->s ? 1: -1;
304 
	return x->s ? -1: 1;
305 
}