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