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