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
fpiround(Internal * i)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
matchexponents(Internal * x,Internal * y)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
shift(Internal * i)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
normalise(Internal * i)74 normalise(Internal *i)
75 {
76 while((i->h & HiddenBit) == 0)
77 shift(i);
78 }
79
80 static void
renormalise(Internal * i)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
fpinormalise(Internal * x)97 fpinormalise(Internal *x)
98 {
99 if(!IsWeird(x) && !IsZero(x))
100 normalise(x);
101 }
102
103 void
fpiadd(Internal * x,Internal * y,Internal * i)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
fpisub(Internal * x,Internal * y,Internal * i)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
fpimul(Internal * x,Internal * y,Internal * i)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
fpidiv(Internal * x,Internal * y,Internal * i)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
fpicmp(Internal * x,Internal * y)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 }
306