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