sys/src/9/omap/fpi.c

   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         if(y->h || y->l)
 283                 i->l |= 0x01;
 284         renormalise(i);
 285 }
 286
 287 int
 288 fpicmp(Internal *x, Internal *y)
 289 {
 290         if(IsNaN(x) && IsNaN(y))
 291                 return 0;
 292         if(IsInfinity(x) && IsInfinity(y))
 293                 return y->s - x->s;
 294         if(x->e == y->e && x->h == y->h && x->l == y->l)
 295                 return y->s - x->s;
 296         if(x->e < y->e
 297            || (x->e == y->e && (x->h < y->h || (x->h == y->h && x->l < y->l))))
 298                 return y->s ? 1: -1;
 299         return x->s ? -1: 1;
 300 }