float.c

   1 /* float.c     floating-point constant support for the Netwide Assembler
   2  *
   3  * The Netwide Assembler is copyright (C) 1996 Simon Tatham and
   4  * Julian Hall. All rights reserved. The software is
   5  * redistributable under the licence given in the file "Licence"
   6  * distributed in the NASM archive.
   7  *
   8  * initial version 13/ix/96 by Simon Tatham
   9  */
  10
  11 #include <ctype.h>
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #include <string.h>
  15 #include <inttypes.h>
  16
  17 #include "nasm.h"
  18
  19 #define TRUE 1
  20 #define FALSE 0
  21
  22 #define MANT_WORDS  10          /* 112 bits + 48 for accuracy == 160 */
  23 #define MANT_DIGITS 49          /* 50 digits don't fit in 160 bits */
  24
  25 /*
  26  * guaranteed top bit of from is set
  27  * => we only have to worry about _one_ bit shift to the left
  28  */
  29
  30 static int ieee_multiply(uint16_t *to, uint16_t *from)
  31 {
  32     uint32_t temp[MANT_WORDS * 2];
  33     int i, j;
  34
  35     for (i = 0; i < MANT_WORDS * 2; i++)
  36         temp[i] = 0;
  37
  38     for (i = 0; i < MANT_WORDS; i++)
  39         for (j = 0; j < MANT_WORDS; j++) {
  40             uint32_t n;
  41             n = (uint32_t)to[i] * (uint32_t)from[j];
  42             temp[i + j] += n >> 16;
  43             temp[i + j + 1] += n & 0xFFFF;
  44         }
  45
  46     for (i = MANT_WORDS * 2; --i;) {
  47         temp[i - 1] += temp[i] >> 16;
  48         temp[i] &= 0xFFFF;
  49     }
  50     if (temp[0] & 0x8000) {
  51         memcpy(to, temp, 2*MANT_WORDS);
  52         return 0;
  53     } else {
  54         for (i = 0; i < MANT_WORDS; i++)
  55             to[i] = (temp[i] << 1) + !!(temp[i + 1] & 0x8000);
  56         return -1;
  57     }
  58 }
  59
  60 static int hexval(char c)
  61 {
  62     if (c >= '0' && c <= '9')
  63         return c-'0';
  64     else if (c >= 'a' && c <= 'f')
  65         return c-'a'+10;
  66     else
  67         return c-'A'+10;
  68 }
  69
  70 static void ieee_flconvert_hex(char *string, uint16_t *mant,
  71                                int32_t *exponent, efunc error)
  72 {
  73     static const int log2tbl[16] =
  74         { -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3 };
  75     uint16_t mult[MANT_WORDS+1], *mp;
  76     int ms;
  77     int32_t twopwr;
  78     int seendot, seendigit;
  79     unsigned char c;
  80
  81     twopwr = 0;
  82     seendot = seendigit = 0;
  83
  84     memset(mult, 0, sizeof mult);
  85
  86     while ((c = *string++) != '\0') {
  87         if (c == '.') {
  88             if (!seendot)
  89                 seendot = TRUE;
  90             else {
  91                 error(ERR_NONFATAL,
  92                       "too many periods in floating-point constant");
  93                 return;
  94             }
  95         } else if (isxdigit(c)) {
  96             int v = hexval(c);
  97
  98             if (!seendigit && v) {
  99                 int l = log2tbl[v];
 100
 101                 seendigit = 1;
 102                 mp = mult;
 103                 ms = 15-l;
 104
 105                 twopwr = seendot ? twopwr-4+l : l-3;
 106             }
 107
 108             if (seendigit) {
 109                 if (ms <= 0) {
 110                     *mp |= v >> -ms;
 111                     mp++;
 112                     if (mp > &mult[MANT_WORDS])
 113                         mp = &mult[MANT_WORDS]; /* Guard slot */
 114                     ms += 16;
 115                 }
 116                 *mp |= v << ms;
 117                 ms -= 4;
 118
 119                 if (!seendot)
 120                     twopwr += 4;
 121             } else {
 122                 if (seendot)
 123                     twopwr -= 4;
 124             }
 125         } else if (c == 'p' || c == 'P') {
 126             twopwr += atoi(string);
 127             break;
 128         } else {
 129             error(ERR_NONFATAL,
 130                   "floating-point constant: `%c' is invalid character",
 131                   c);
 132             return;
 133         }
 134     }
 135
 136     if (!seendigit) {
 137         memset(mant, 0, 2*MANT_WORDS); /* Zero */
 138         *exponent = 0;
 139     } else {
 140         memcpy(mant, mult, 2*MANT_WORDS);
 141         *exponent = twopwr;
 142     }
 143 }
 144
 145 static void ieee_flconvert(char *string, uint16_t *mant,
 146                            int32_t *exponent, efunc error)
 147 {
 148     char digits[MANT_DIGITS];
 149     char *p, *q, *r;
 150     uint16_t mult[MANT_WORDS], bit;
 151     uint16_t *m;
 152     int32_t tenpwr, twopwr;
 153     int extratwos, started, seendot;
 154
 155     if (string[0] == '0' && (string[1] == 'x' || string[1] == 'X')) {
 156         ieee_flconvert_hex(string+2, mant, exponent, error);
 157         return;
 158     }
 159
 160     p = digits;
 161     tenpwr = 0;
 162     started = seendot = FALSE;
 163     while (*string && *string != 'E' && *string != 'e') {
 164         if (*string == '.') {
 165             if (!seendot)
 166                 seendot = TRUE;
 167             else {
 168                 error(ERR_NONFATAL,
 169                       "too many periods in floating-point constant");
 170                 return;
 171             }
 172         } else if (*string >= '0' && *string <= '9') {
 173             if (*string == '0' && !started) {
 174                 if (seendot)
 175                     tenpwr--;
 176             } else {
 177                 started = TRUE;
 178                 if (p < digits + sizeof(digits))
 179                     *p++ = *string - '0';
 180                 if (!seendot)
 181                     tenpwr++;
 182             }
 183         } else {
 184             error(ERR_NONFATAL,
 185                   "floating-point constant: `%c' is invalid character",
 186                   *string);
 187             return;
 188         }
 189         string++;
 190     }
 191     if (*string) {
 192         string++;               /* eat the E */
 193         tenpwr += atoi(string);
 194     }
 195
 196     /*
 197      * At this point, the memory interval [digits,p) contains a
 198      * series of decimal digits zzzzzzz such that our number X
 199      * satisfies
 200      *
 201      * X = 0.zzzzzzz * 10^tenpwr
 202      */
 203
 204     bit = 0x8000;
 205     for (m = mant; m < mant + MANT_WORDS; m++)
 206         *m = 0;
 207     m = mant;
 208     q = digits;
 209     started = FALSE;
 210     twopwr = 0;
 211     while (m < mant + MANT_WORDS) {
 212         uint16_t carry = 0;
 213         while (p > q && !p[-1])
 214             p--;
 215         if (p <= q)
 216             break;
 217         for (r = p; r-- > q;) {
 218             int i;
 219
 220             i = 2 * *r + carry;
 221             if (i >= 10)
 222                 carry = 1, i -= 10;
 223             else
 224                 carry = 0;
 225             *r = i;
 226         }
 227         if (carry)
 228             *m |= bit, started = TRUE;
 229         if (started) {
 230             if (bit == 1)
 231                 bit = 0x8000, m++;
 232             else
 233                 bit >>= 1;
 234         } else
 235             twopwr--;
 236     }
 237     twopwr += tenpwr;
 238
 239     /*
 240      * At this point the `mant' array contains the first six
 241      * fractional places of a base-2^16 real number, which when
 242      * multiplied by 2^twopwr and 5^tenpwr gives X. So now we
 243      * really do multiply by 5^tenpwr.
 244      */
 245
 246     if (tenpwr < 0) {
 247         for (m = mult; m < mult + MANT_WORDS; m++)
 248             *m = 0xCCCC;
 249         extratwos = -2;
 250         tenpwr = -tenpwr;
 251     } else if (tenpwr > 0) {
 252         mult[0] = 0xA000;
 253         for (m = mult + 1; m < mult + MANT_WORDS; m++)
 254             *m = 0;
 255         extratwos = 3;
 256     } else
 257         extratwos = 0;
 258     while (tenpwr) {
 259         if (tenpwr & 1)
 260             twopwr += extratwos + ieee_multiply(mant, mult);
 261         extratwos = extratwos * 2 + ieee_multiply(mult, mult);
 262         tenpwr >>= 1;
 263     }
 264
 265     /*
 266      * Conversion is done. The elements of `mant' contain the first
 267      * fractional places of a base-2^16 real number in [0.5,1)
 268      * which we can multiply by 2^twopwr to get X. Or, of course,
 269      * it contains zero.
 270      */
 271     *exponent = twopwr;
 272 }
 273
 274 /*
 275  * Shift a mantissa to the right by i (i < 16) bits.
 276  */
 277 static void ieee_shr(uint16_t *mant, int i)
 278 {
 279     uint16_t n = 0, m;
 280     int j;
 281
 282     for (j = 0; j < MANT_WORDS; j++) {
 283         m = (mant[j] << (16 - i)) & 0xFFFF;
 284         mant[j] = (mant[j] >> i) | n;
 285         n = m;
 286     }
 287 }
 288
 289 /*
 290  * Round a mantissa off after i words.
 291  */
 292 static int ieee_round(uint16_t *mant, int i)
 293 {
 294     if (mant[i] & 0x8000) {
 295         do {
 296             ++mant[--i];
 297             mant[i] &= 0xFFFF;
 298         } while (i > 0 && !mant[i]);
 299         return !i && !mant[i];
 300     }
 301     return 0;
 302 }
 303
 304 #define put(a,b) ( (*(a)=(b)), ((a)[1]=(b)>>8) )
 305
 306 /* Set a bit, using *bigendian* bit numbering (0 = MSB) */
 307 static void set_bit(uint16_t *mant, int bit)
 308 {
 309     mant[bit >> 4] |= 1 << (~bit & 15);
 310 }
 311
 312 /* Produce standard IEEE formats, with implicit "1" bit; this makes
 313    the following assumptions:
 314
 315    - the sign bit is the MSB, followed by the exponent.
 316    - the sign bit plus exponent fit in 16 bits.
 317    - the exponent bias is 2^(n-1)-1 for an n-bit exponent */
 318
 319 struct ieee_format {
 320     int words;
 321     int mantissa;               /* Bits in the mantissa */
 322     int exponent;               /* Bits in the exponent */
 323 };
 324
 325 static const struct ieee_format ieee_16  = { 1,  10,  5 };
 326 static const struct ieee_format ieee_32  = { 2,  23,  8 };
 327 static const struct ieee_format ieee_64  = { 4,  52, 11 };
 328 static const struct ieee_format ieee_128 = { 8, 112, 15 };
 329
 330 /* Produce all the standard IEEE formats: 16, 32, 64, and 128 bits */
 331 static int to_float(char *str, int32_t sign, uint8_t *result,
 332                     const struct ieee_format *fmt, efunc error)
 333 {
 334     uint16_t mant[MANT_WORDS], *mp;
 335     int32_t exponent;
 336     int32_t expmax = 1 << (fmt->exponent-1);
 337     uint16_t implicit_one = 0x8000 >> fmt->exponent;
 338     int i;
 339
 340     sign = (sign < 0 ? 0x8000L : 0L);
 341
 342     if (str[0] == '_') {
 343         /* NaN or Infinity */
 344         int32_t expmask = (1 << fmt->exponent)-1;
 345
 346         memset(mant, 0, sizeof mant);
 347         mant[0] = expmask << (15-fmt->exponent); /* Exponent: all bits one */
 348
 349         switch (str[2]) {
 350         case 'n':               /* __nan__ */
 351         case 'N':
 352         case 'q':               /* __qnan__ */
 353         case 'Q':
 354             set_bit(mant, fmt->exponent+1); /* Highest bit in mantissa */
 355             break;
 356         case 's':               /* __snan__ */
 357         case 'S':
 358             set_bit(mant, fmt->exponent+fmt->mantissa); /* Last bit */
 359             break;
 360         case 'i':               /* __infinity__ */
 361         case 'I':
 362             break;
 363         }
 364     } else {
 365         ieee_flconvert(str, mant, &exponent, error);
 366         if (mant[0] & 0x8000) {
 367             /*
 368              * Non-zero.
 369              */
 370             exponent--;
 371             if (exponent >= 2-expmax && exponent <= expmax) {
 372                 /*
 373                  * Normalised.
 374                  */
 375                 exponent += expmax;
 376                 ieee_shr(mant, fmt->exponent);
 377                 ieee_round(mant, fmt->words);
 378                 /* did we scale up by one? */
 379                 if (mant[0] & (implicit_one << 1)) {
 380                     ieee_shr(mant, 1);
 381                     exponent++;
 382                 }
 383
 384                 mant[0] &= (implicit_one-1);     /* remove leading one */
 385                 mant[0] |= exponent << (15 - fmt->exponent);
 386             } else if (exponent < 2-expmax &&
 387                        exponent >= 2-expmax-fmt->mantissa) {
 388                 /*
 389                  * Denormal.
 390                  */
 391                 int shift = -(exponent + expmax-2-fmt->exponent);
 392                 int sh = shift % 16, wds = shift / 16;
 393                 ieee_shr(mant, sh);
 394                 if (ieee_round(mant, fmt->words - wds)
 395                     || (sh > 0 && (mant[0] & (0x8000 >> (sh - 1))))) {
 396                     ieee_shr(mant, 1);
 397                     if (sh == 0)
 398                         mant[0] |= 0x8000;
 399                     exponent++;
 400                 }
 401
 402                 if (wds) {
 403                     for (i = fmt->words-1; i >= wds; i--)
 404                         mant[i] = mant[i-wds];
 405                     for (; i >= 0; i--)
 406                         mant[i] = 0;
 407                 }
 408             } else {
 409                 if (exponent > 0) {
 410                     error(ERR_NONFATAL, "overflow in floating-point constant");
 411                     return 0;
 412                 } else {
 413                     memset(mant, 0, 2*fmt->words);
 414                 }
 415             }
 416         } else {
 417             /* Zero */
 418             memset(mant, 0, 2*fmt->words);
 419         }
 420     }
 421
 422     mant[0] |= sign;
 423
 424     for (mp = &mant[fmt->words], i = 0; i < fmt->words; i++) {
 425         uint16_t m = *--mp;
 426         put(result, m);
 427         result += 2;
 428     }
 429
 430     return 1;                   /* success */
 431 }
 432
 433 /* 80-bit format with 64-bit mantissa *including an explicit integer 1*
 434    and 15-bit exponent. */
 435 static int to_ldoub(char *str, int32_t sign, uint8_t *result,
 436                     efunc error)
 437 {
 438     uint16_t mant[MANT_WORDS];
 439     int32_t exponent;
 440
 441     sign = (sign < 0 ? 0x8000L : 0L);
 442
 443     if (str[0] == '_') {
 444         uint16_t is_snan = 0, is_qnan = 0x8000;
 445         switch (str[2]) {
 446         case 'n':
 447         case 'N':
 448         case 'q':
 449         case 'Q':
 450             is_qnan = 0xc000;
 451             break;
 452         case 's':
 453         case 'S':
 454             is_snan = 1;
 455             break;
 456         case 'i':
 457         case 'I':
 458             break;
 459         }
 460         put(result + 0, is_snan);
 461         put(result + 2, 0);
 462         put(result + 4, 0);
 463         put(result + 6, is_qnan);
 464         put(result + 8, 0x7fff|sign);
 465         return 1;
 466     }
 467
 468     ieee_flconvert(str, mant, &exponent, error);
 469     if (mant[0] & 0x8000) {
 470         /*
 471          * Non-zero.
 472          */
 473         exponent--;
 474         if (exponent >= -16383 && exponent <= 16384) {
 475             /*
 476              * Normalised.
 477              */
 478             exponent += 16383;
 479             if (ieee_round(mant, 4))    /* did we scale up by one? */
 480                 ieee_shr(mant, 1), mant[0] |= 0x8000, exponent++;
 481             put(result + 0, mant[3]);
 482             put(result + 2, mant[2]);
 483             put(result + 4, mant[1]);
 484             put(result + 6, mant[0]);
 485             put(result + 8, exponent | sign);
 486         } else if (exponent < -16383 && exponent >= -16446) {
 487             /*
 488              * Denormal.
 489              */
 490             int shift = -(exponent + 16383);
 491             int sh = shift % 16, wds = shift / 16;
 492             ieee_shr(mant, sh);
 493             if (ieee_round(mant, 4 - wds)
 494                 || (sh > 0 && (mant[0] & (0x8000 >> (sh - 1))))) {
 495                 ieee_shr(mant, 1);
 496                 if (sh == 0)
 497                     mant[0] |= 0x8000;
 498                 exponent++;
 499             }
 500             put(result + 0, (wds <= 3 ? mant[3 - wds] : 0));
 501             put(result + 2, (wds <= 2 ? mant[2 - wds] : 0));
 502             put(result + 4, (wds <= 1 ? mant[1 - wds] : 0));
 503             put(result + 6, (wds == 0 ? mant[0] : 0));
 504             put(result + 8, sign);
 505         } else {
 506             if (exponent > 0) {
 507                 error(ERR_NONFATAL, "overflow in floating-point constant");
 508                 return 0;
 509             } else {
 510                 goto zero;
 511             }
 512         }
 513     } else {
 514         /*
 515          * Zero.
 516          */
 517     zero:
 518         put(result + 0, 0);
 519         put(result + 2, 0);
 520         put(result + 4, 0);
 521         put(result + 6, 0);
 522         put(result + 8, sign);
 523     }
 524     return 1;
 525 }
 526
 527 int float_const(char *number, int32_t sign, uint8_t *result, int bytes,
 528                 efunc error)
 529 {
 530     switch (bytes) {
 531     case 2:
 532         return to_float(number, sign, result, &ieee_16, error);
 533     case 4:
 534         return to_float(number, sign, result, &ieee_32, error);
 535     case 8:
 536         return to_float(number, sign, result, &ieee_64, error);
 537     case 10:
 538         return to_ldoub(number, sign, result, error);
 539     case 16:
 540         return to_float(number, sign, result, &ieee_128, error);
 541     default:
 542         error(ERR_PANIC, "strange value %d passed to float_const", bytes);
 543         return 0;
 544     }
 545 }