float.c

   1 /* ----------------------------------------------------------------------- *
   2  *
   3  *   Copyright 1996-2009 The NASM Authors - All Rights Reserved
   4  *   See the file AUTHORS included with the NASM distribution for
   5  *   the specific copyright holders.
   6  *
   7  *   Redistribution and use in source and binary forms, with or without
   8  *   modification, are permitted provided that the following
   9  *   conditions are met:
  10  *
  11  *   * Redistributions of source code must retain the above copyright
  12  *     notice, this list of conditions and the following disclaimer.
  13  *   * Redistributions in binary form must reproduce the above
  14  *     copyright notice, this list of conditions and the following
  15  *     disclaimer in the documentation and/or other materials provided
  16  *     with the distribution.
  17  *
  18  *     THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
  19  *     CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
  20  *     INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
  21  *     MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  22  *     DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
  23  *     CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  24  *     SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  25  *     NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  26  *     LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  27  *     HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  28  *     CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
  29  *     OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
  30  *     EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  31  *
  32  * ----------------------------------------------------------------------- */
  33
  34 /*
  35  * float.c     floating-point constant support for the Netwide Assembler
  36  */
  37
  38 #include "compiler.h"
  39
  40 #include <ctype.h>
  41 #include <stdio.h>
  42 #include <stdlib.h>
  43 #include <string.h>
  44 #include <inttypes.h>
  45
  46 #include "nasm.h"
  47 #include "float.h"
  48
  49 /*
  50  * -----------------
  51  *  local variables
  52  * -----------------
  53  */
  54 static efunc error;
  55 static bool daz = false;        /* denormals as zero */
  56 static enum float_round rc = FLOAT_RC_NEAR;     /* rounding control */
  57
  58 /*
  59  * -----------
  60  *  constants
  61  * -----------
  62  */
  63
  64 /* "A limb is like a digit but bigger */
  65 typedef uint32_t fp_limb;
  66 typedef uint64_t fp_2limb;
  67
  68 #define LIMB_BITS       32
  69 #define LIMB_BYTES      (LIMB_BITS/8)
  70 #define LIMB_TOP_BIT    ((fp_limb)1 << (LIMB_BITS-1))
  71 #define LIMB_MASK       ((fp_limb)(~0))
  72 #define LIMB_ALL_BYTES  ((fp_limb)0x01010101)
  73 #define LIMB_BYTE(x)    ((x)*LIMB_ALL_BYTES)
  74
  75 /* 112 bits + 64 bits for accuracy + 16 bits for rounding */
  76 #define MANT_LIMBS 6
  77
  78 /* 52 digits fit in 176 bits because 10^53 > 2^176 > 10^52 */
  79 #define MANT_DIGITS 52
  80
  81 /* the format and the argument list depend on MANT_LIMBS */
  82 #define MANT_FMT "%08x_%08x_%08x_%08x_%08x_%08x"
  83 #define MANT_ARG SOME_ARG(mant, 0)
  84
  85 #define SOME_ARG(a,i) (a)[(i)+0], (a)[(i)+1], (a)[(i)+2], \
  86                       (a)[(i)+3], (a)[(i)+4], (a)[(i)+5]
  87
  88 /*
  89  * ---------------------------------------------------------------------------
  90  *  emit a printf()-like debug message... but only if DEBUG_FLOAT was defined
  91  * ---------------------------------------------------------------------------
  92  */
  93
  94 #ifdef DEBUG_FLOAT
  95 #define dprintf(x) printf x
  96 #else
  97 #define dprintf(x) do { } while (0)
  98 #endif
  99
 100 /*
 101  * ---------------------------------------------------------------------------
 102  *  multiply
 103  * ---------------------------------------------------------------------------
 104  */
 105 static int float_multiply(fp_limb *to, fp_limb *from)
 106 {
 107     fp_2limb temp[MANT_LIMBS * 2];
 108     int i, j;
 109
 110     /*
 111      * guaranteed that top bit of 'from' is set -- so we only have
 112      * to worry about _one_ bit shift to the left
 113      */
 114     dprintf(("%s=" MANT_FMT "\n", "mul1", SOME_ARG(to, 0)));
 115     dprintf(("%s=" MANT_FMT "\n", "mul2", SOME_ARG(from, 0)));
 116
 117     memset(temp, 0, sizeof temp);
 118
 119     for (i = 0; i < MANT_LIMBS; i++) {
 120         for (j = 0; j < MANT_LIMBS; j++) {
 121             fp_2limb n;
 122             n = (fp_2limb) to[i] * (fp_2limb) from[j];
 123             temp[i + j] += n >> LIMB_BITS;
 124             temp[i + j + 1] += (fp_limb)n;
 125         }
 126     }
 127
 128     for (i = MANT_LIMBS * 2; --i;) {
 129         temp[i - 1] += temp[i] >> LIMB_BITS;
 130         temp[i] &= LIMB_MASK;
 131     }
 132
 133     dprintf(("%s=" MANT_FMT "_" MANT_FMT "\n", "temp", SOME_ARG(temp, 0),
 134              SOME_ARG(temp, MANT_LIMBS)));
 135
 136     if (temp[0] & LIMB_TOP_BIT) {
 137         for (i = 0; i < MANT_LIMBS; i++) {
 138             to[i] = temp[i] & LIMB_MASK;
 139         }
 140         dprintf(("%s=" MANT_FMT " (%i)\n", "prod", SOME_ARG(to, 0), 0));
 141         return 0;
 142     } else {
 143         for (i = 0; i < MANT_LIMBS; i++) {
 144             to[i] = (temp[i] << 1) + !!(temp[i + 1] & LIMB_TOP_BIT);
 145         }
 146         dprintf(("%s=" MANT_FMT " (%i)\n", "prod", SOME_ARG(to, 0), -1));
 147         return -1;
 148     }
 149 }
 150
 151 /*
 152  * ---------------------------------------------------------------------------
 153  *  read an exponent; returns INT32_MAX on error
 154  * ---------------------------------------------------------------------------
 155  */
 156 static int32_t read_exponent(const char *string, int32_t max)
 157 {
 158     int32_t i = 0;
 159     bool neg = false;
 160
 161     if (*string == '+') {
 162         string++;
 163     } else if (*string == '-') {
 164         neg = true;
 165         string++;
 166     }
 167     while (*string) {
 168         if (*string >= '0' && *string <= '9') {
 169             i = (i * 10) + (*string - '0');
 170
 171             /*
 172              * To ensure that underflows and overflows are
 173              * handled properly we must avoid wraparounds of
 174              * the signed integer value that is used to hold
 175              * the exponent. Therefore we cap the exponent at
 176              * +/-5000, which is slightly more/less than
 177              * what's required for normal and denormal numbers
 178              * in single, double, and extended precision, but
 179              * sufficient to avoid signed integer wraparound.
 180              */
 181             if (i > max)
 182                 i = max;
 183         } else if (*string == '_') {
 184             /* do nothing */
 185         } else {
 186             error(ERR_NONFATAL|ERR_PASS1,
 187                   "invalid character in floating-point constant %s: '%c'",
 188                   "exponent", *string);
 189             return INT32_MAX;
 190         }
 191         string++;
 192     }
 193
 194     return neg ? -i : i;
 195 }
 196
 197 /*
 198  * ---------------------------------------------------------------------------
 199  *  convert
 200  * ---------------------------------------------------------------------------
 201  */
 202 static bool ieee_flconvert(const char *string, fp_limb *mant,
 203                            int32_t * exponent)
 204 {
 205     char digits[MANT_DIGITS];
 206     char *p, *q, *r;
 207     fp_limb mult[MANT_LIMBS], bit;
 208     fp_limb *m;
 209     int32_t tenpwr, twopwr;
 210     int32_t extratwos;
 211     bool started, seendot, warned;
 212
 213     warned = false;
 214     p = digits;
 215     tenpwr = 0;
 216     started = seendot = false;
 217
 218     while (*string && *string != 'E' && *string != 'e') {
 219         if (*string == '.') {
 220             if (!seendot) {
 221                 seendot = true;
 222             } else {
 223                 error(ERR_NONFATAL|ERR_PASS1,
 224                       "too many periods in floating-point constant");
 225                 return false;
 226             }
 227         } else if (*string >= '0' && *string <= '9') {
 228             if (*string == '0' && !started) {
 229                 if (seendot) {
 230                     tenpwr--;
 231                 }
 232             } else {
 233                 started = true;
 234                 if (p < digits + sizeof(digits)) {
 235                     *p++ = *string - '0';
 236                 } else {
 237                     if (!warned) {
 238                         error(ERR_WARNING|ERR_WARN_FL_TOOLONG|ERR_PASS1,
 239                               "floating-point constant significand contains "
 240                               "more than %i digits", MANT_DIGITS);
 241                         warned = true;
 242                     }
 243                 }
 244                 if (!seendot) {
 245                     tenpwr++;
 246                 }
 247             }
 248         } else if (*string == '_') {
 249             /* do nothing */
 250         } else {
 251             error(ERR_NONFATAL|ERR_PASS1,
 252                   "invalid character in floating-point constant %s: '%c'",
 253                   "significand", *string);
 254             return false;
 255         }
 256         string++;
 257     }
 258
 259     if (*string) {
 260         int32_t e;
 261
 262         string++;               /* eat the E */
 263         e = read_exponent(string, 5000);
 264         if (e == INT32_MAX)
 265             return false;
 266         tenpwr += e;
 267     }
 268
 269     /*
 270      * At this point, the memory interval [digits,p) contains a
 271      * series of decimal digits zzzzzzz, such that our number X
 272      * satisfies X = 0.zzzzzzz * 10^tenpwr.
 273      */
 274     q = digits;
 275     dprintf(("X = 0."));
 276     while (q < p) {
 277         dprintf(("%c", *q + '0'));
 278         q++;
 279     }
 280     dprintf((" * 10^%i\n", tenpwr));
 281
 282     /*
 283      * Now convert [digits,p) to our internal representation.
 284      */
 285     bit = LIMB_TOP_BIT;
 286     for (m = mant; m < mant + MANT_LIMBS; m++) {
 287         *m = 0;
 288     }
 289     m = mant;
 290     q = digits;
 291     started = false;
 292     twopwr = 0;
 293     while (m < mant + MANT_LIMBS) {
 294         fp_limb carry = 0;
 295         while (p > q && !p[-1]) {
 296             p--;
 297         }
 298         if (p <= q) {
 299             break;
 300         }
 301         for (r = p; r-- > q;) {
 302             int32_t i;
 303             i = 2 * *r + carry;
 304             if (i >= 10) {
 305                 carry = 1;
 306                 i -= 10;
 307             } else {
 308                 carry = 0;
 309             }
 310             *r = i;
 311         }
 312         if (carry) {
 313             *m |= bit;
 314             started = true;
 315         }
 316         if (started) {
 317             if (bit == 1) {
 318                 bit = LIMB_TOP_BIT;
 319                 m++;
 320             } else {
 321                 bit >>= 1;
 322             }
 323         } else {
 324             twopwr--;
 325         }
 326     }
 327     twopwr += tenpwr;
 328
 329     /*
 330      * At this point, the 'mant' array contains the first frac-
 331      * tional places of a base-2^16 real number which when mul-
 332      * tiplied by 2^twopwr and 5^tenpwr gives X.
 333      */
 334     dprintf(("X = " MANT_FMT " * 2^%i * 5^%i\n", MANT_ARG, twopwr,
 335              tenpwr));
 336
 337     /*
 338      * Now multiply 'mant' by 5^tenpwr.
 339      */
 340     if (tenpwr < 0) {           /* mult = 5^-1 = 0.2 */
 341         for (m = mult; m < mult + MANT_LIMBS - 1; m++) {
 342             *m = LIMB_BYTE(0xcc);
 343         }
 344         mult[MANT_LIMBS - 1] = LIMB_BYTE(0xcc)+1;
 345         extratwos = -2;
 346         tenpwr = -tenpwr;
 347
 348         /*
 349          * If tenpwr was 1000...000b, then it becomes 1000...000b. See
 350          * the "ANSI C" comment below for more details on that case.
 351          *
 352          * Because we already truncated tenpwr to +5000...-5000 inside
 353          * the exponent parsing code, this shouldn't happen though.
 354          */
 355     } else if (tenpwr > 0) {    /* mult = 5^+1 = 5.0 */
 356         mult[0] = (fp_limb)5 << (LIMB_BITS-3); /* 0xA000... */
 357         for (m = mult + 1; m < mult + MANT_LIMBS; m++) {
 358             *m = 0;
 359         }
 360         extratwos = 3;
 361     } else {
 362         extratwos = 0;
 363     }
 364     while (tenpwr) {
 365         dprintf(("loop=" MANT_FMT " * 2^%i * 5^%i (%i)\n", MANT_ARG,
 366                  twopwr, tenpwr, extratwos));
 367         if (tenpwr & 1) {
 368             dprintf(("mant*mult\n"));
 369             twopwr += extratwos + float_multiply(mant, mult);
 370         }
 371         dprintf(("mult*mult\n"));
 372         extratwos = extratwos * 2 + float_multiply(mult, mult);
 373         tenpwr >>= 1;
 374
 375         /*
 376          * In ANSI C, the result of right-shifting a signed integer is
 377          * considered implementation-specific. To ensure that the loop
 378          * terminates even if tenpwr was 1000...000b to begin with, we
 379          * manually clear the MSB, in case a 1 was shifted in.
 380          *
 381          * Because we already truncated tenpwr to +5000...-5000 inside
 382          * the exponent parsing code, this shouldn't matter; neverthe-
 383          * less it is the right thing to do here.
 384          */
 385         tenpwr &= (uint32_t) - 1 >> 1;
 386     }
 387
 388     /*
 389      * At this point, the 'mant' array contains the first frac-
 390      * tional places of a base-2^16 real number in [0.5,1) that
 391      * when multiplied by 2^twopwr gives X. Or it contains zero
 392      * of course. We are done.
 393      */
 394     *exponent = twopwr;
 395     return true;
 396 }
 397
 398 /*
 399  * ---------------------------------------------------------------------------
 400  *  operations of specific bits
 401  * ---------------------------------------------------------------------------
 402  */
 403
 404 /* Set a bit, using *bigendian* bit numbering (0 = MSB) */
 405 static void set_bit(fp_limb *mant, int bit)
 406 {
 407     mant[bit/LIMB_BITS] |= LIMB_TOP_BIT >> (bit & (LIMB_BITS-1));
 408 }
 409
 410 /* Test a single bit */
 411 static int test_bit(const fp_limb *mant, int bit)
 412 {
 413     return (mant[bit/LIMB_BITS] >> (~bit & (LIMB_BITS-1))) & 1;
 414 }
 415
 416 /* Report if the mantissa value is all zero */
 417 static bool is_zero(const fp_limb *mant)
 418 {
 419     int i;
 420
 421     for (i = 0; i < MANT_LIMBS; i++)
 422         if (mant[i])
 423             return false;
 424
 425     return true;
 426 }
 427
 428 /*
 429  * ---------------------------------------------------------------------------
 430  *  round a mantissa off after i words
 431  * ---------------------------------------------------------------------------
 432  */
 433
 434 #define ROUND_COLLECT_BITS                      \
 435     do {                                        \
 436         m = mant[i] & (2*bit-1);                \
 437         for (j = i+1; j < MANT_LIMBS; j++)      \
 438             m = m | mant[j];                    \
 439     } while (0)
 440
 441 #define ROUND_ABS_DOWN                          \
 442     do {                                        \
 443         mant[i] &= ~(bit-1);                    \
 444         for (j = i+1; j < MANT_LIMBS; j++)      \
 445             mant[j] = 0;                        \
 446         return false;                           \
 447     } while (0)
 448
 449 #define ROUND_ABS_UP                            \
 450     do {                                        \
 451         mant[i] = (mant[i] & ~(bit-1)) + bit;   \
 452         for (j = i+1; j < MANT_LIMBS; j++)      \
 453             mant[j] = 0;                        \
 454         while (i > 0 && !mant[i])               \
 455             ++mant[--i];                        \
 456         return !mant[0];                        \
 457     } while (0)
 458
 459 static bool ieee_round(bool minus, fp_limb *mant, int bits)
 460 {
 461     fp_limb m = 0;
 462     int32_t j;
 463     int i = bits / LIMB_BITS;
 464     int p = bits % LIMB_BITS;
 465     fp_limb bit = LIMB_TOP_BIT >> p;
 466
 467     if (rc == FLOAT_RC_NEAR) {
 468         if (mant[i] & bit) {
 469             mant[i] &= ~bit;
 470             ROUND_COLLECT_BITS;
 471             mant[i] |= bit;
 472             if (m) {
 473                 ROUND_ABS_UP;
 474             } else {
 475                 if (test_bit(mant, bits-1)) {
 476                     ROUND_ABS_UP;
 477                 } else {
 478                     ROUND_ABS_DOWN;
 479                 }
 480             }
 481         } else {
 482             ROUND_ABS_DOWN;
 483         }
 484     } else if (rc == FLOAT_RC_ZERO ||
 485                rc == (minus ? FLOAT_RC_UP : FLOAT_RC_DOWN)) {
 486         ROUND_ABS_DOWN;
 487     } else {
 488         /* rc == (minus ? FLOAT_RC_DOWN : FLOAT_RC_UP) */
 489         /* Round toward +/- infinity */
 490         ROUND_COLLECT_BITS;
 491         if (m) {
 492             ROUND_ABS_UP;
 493         } else {
 494             ROUND_ABS_DOWN;
 495         }
 496     }
 497     return false;
 498 }
 499
 500 /* Returns a value >= 16 if not a valid hex digit */
 501 static unsigned int hexval(char c)
 502 {
 503     unsigned int v = (unsigned char) c;
 504
 505     if (v >= '0' && v <= '9')
 506         return v - '0';
 507     else
 508         return (v|0x20) - 'a' + 10;
 509 }
 510
 511 /* Handle floating-point numbers with radix 2^bits and binary exponent */
 512 static bool ieee_flconvert_bin(const char *string, int bits,
 513                                fp_limb *mant, int32_t *exponent)
 514 {
 515     static const int log2tbl[16] =
 516         { -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3 };
 517     fp_limb mult[MANT_LIMBS + 1], *mp;
 518     int ms;
 519     int32_t twopwr;
 520     bool seendot, seendigit;
 521     unsigned char c;
 522     const int radix = 1 << bits;
 523     fp_limb v;
 524
 525     twopwr = 0;
 526     seendot = seendigit = false;
 527     ms = 0;
 528     mp = NULL;
 529
 530     memset(mult, 0, sizeof mult);
 531
 532     while ((c = *string++) != '\0') {
 533         if (c == '.') {
 534             if (!seendot)
 535                 seendot = true;
 536             else {
 537                 error(ERR_NONFATAL|ERR_PASS1,
 538                       "too many periods in floating-point constant");
 539                 return false;
 540             }
 541         } else if ((v = hexval(c)) < (unsigned int)radix) {
 542             if (!seendigit && v) {
 543                 int l = log2tbl[v];
 544
 545                 seendigit = true;
 546                 mp = mult;
 547                 ms = (LIMB_BITS-1)-l;
 548
 549                 twopwr = seendot ? twopwr-bits+l : l+1-bits;
 550             }
 551
 552             if (seendigit) {
 553                 if (ms <= 0) {
 554                     *mp |= v >> -ms;
 555                     mp++;
 556                     if (mp > &mult[MANT_LIMBS])
 557                         mp = &mult[MANT_LIMBS]; /* Guard slot */
 558                     ms += LIMB_BITS;
 559                 }
 560                 *mp |= v << ms;
 561                 ms -= bits;
 562
 563                 if (!seendot)
 564                     twopwr += bits;
 565             } else {
 566                 if (seendot)
 567                     twopwr -= bits;
 568             }
 569         } else if (c == 'p' || c == 'P') {
 570             int32_t e;
 571             e = read_exponent(string, 20000);
 572             if (e == INT32_MAX)
 573                 return false;
 574             twopwr += e;
 575             break;
 576         } else if (c == '_') {
 577             /* ignore */
 578         } else {
 579             error(ERR_NONFATAL|ERR_PASS1,
 580                   "floating-point constant: `%c' is invalid character", c);
 581             return false;
 582         }
 583     }
 584
 585     if (!seendigit) {
 586         memset(mant, 0, MANT_LIMBS*sizeof(fp_limb)); /* Zero */
 587         *exponent = 0;
 588     } else {
 589         memcpy(mant, mult, MANT_LIMBS*sizeof(fp_limb));
 590         *exponent = twopwr;
 591     }
 592
 593     return true;
 594 }
 595
 596 /*
 597  * Shift a mantissa to the right by i bits.
 598  */
 599 static void ieee_shr(fp_limb *mant, int i)
 600 {
 601     fp_limb n, m;
 602     int j = 0;
 603     int sr, sl, offs;
 604
 605     sr = i % LIMB_BITS; sl = LIMB_BITS-sr;
 606     offs = i/LIMB_BITS;
 607
 608     if (sr == 0) {
 609         if (offs)
 610             for (j = MANT_LIMBS-1; j >= offs; j--)
 611                 mant[j] = mant[j-offs];
 612     } else {
 613         n = mant[MANT_LIMBS-1-offs] >> sr;
 614         for (j = MANT_LIMBS-1; j > offs; j--) {
 615             m = mant[j-offs-1];
 616             mant[j] = (m << sl) | n;
 617             n = m >> sr;
 618         }
 619         mant[j--] = n;
 620     }
 621     while (j >= 0)
 622         mant[j--] = 0;
 623 }
 624
 625 /* Produce standard IEEE formats, with implicit or explicit integer
 626    bit; this makes the following assumptions:
 627
 628    - the sign bit is the MSB, followed by the exponent,
 629      followed by the integer bit if present.
 630    - the sign bit plus exponent fit in 16 bits.
 631    - the exponent bias is 2^(n-1)-1 for an n-bit exponent */
 632
 633 struct ieee_format {
 634     int bytes;
 635     int mantissa;               /* Fractional bits in the mantissa */
 636     int explicit;               /* Explicit integer */
 637     int exponent;               /* Bits in the exponent */
 638 };
 639
 640 /*
 641  * The 16- and 128-bit formats are expected to be in IEEE 754r.
 642  * AMD SSE5 uses the 16-bit format.
 643  *
 644  * The 32- and 64-bit formats are the original IEEE 754 formats.
 645  *
 646  * The 80-bit format is x87-specific, but widely used.
 647  *
 648  * The 8-bit format appears to be the consensus 8-bit floating-point
 649  * format.  It is apparently used in graphics applications.
 650  */
 651 static const struct ieee_format ieee_8   = {  1,   3, 0,  4 };
 652 static const struct ieee_format ieee_16  = {  2,  10, 0,  5 };
 653 static const struct ieee_format ieee_32  = {  4,  23, 0,  8 };
 654 static const struct ieee_format ieee_64  = {  8,  52, 0, 11 };
 655 static const struct ieee_format ieee_80  = { 10,  63, 1, 15 };
 656 static const struct ieee_format ieee_128 = { 16, 112, 0, 15 };
 657
 658 /* Types of values we can generate */
 659 enum floats {
 660     FL_ZERO,
 661     FL_DENORMAL,
 662     FL_NORMAL,
 663     FL_INFINITY,
 664     FL_QNAN,
 665     FL_SNAN
 666 };
 667
 668 static int to_packed_bcd(const char *str, const char *p,
 669                          int s, uint8_t *result,
 670                          const struct ieee_format *fmt)
 671 {
 672     int n = 0;
 673     char c;
 674     int tv = -1;
 675
 676     if (fmt != &ieee_80) {
 677         error(ERR_NONFATAL|ERR_PASS1,
 678               "packed BCD requires an 80-bit format");
 679         return 0;
 680     }
 681
 682     while (p >= str) {
 683         c = *p--;
 684         if (c >= '0' && c <= '9') {
 685             if (tv < 0) {
 686                 if (n == 9) {
 687                     error(ERR_WARNING|ERR_PASS1,
 688                           "packed BCD truncated to 18 digits");
 689                 }
 690                 tv = c-'0';
 691             } else {
 692                 if (n < 9)
 693                     *result++ = tv + ((c-'0') << 4);
 694                 n++;
 695                 tv = -1;
 696             }
 697         } else if (c == '_') {
 698             /* do nothing */
 699         } else {
 700             error(ERR_NONFATAL|ERR_PASS1,
 701                   "invalid character `%c' in packed BCD constant", c);
 702             return 0;
 703         }
 704     }
 705     if (tv >= 0) {
 706         if (n < 9)
 707             *result++ = tv;
 708         n++;
 709     }
 710     while (n < 9) {
 711         *result++ = 0;
 712         n++;
 713     }
 714     *result = (s < 0) ? 0x80 : 0;
 715
 716     return 1;                   /* success */
 717 }
 718
 719 static int to_float(const char *str, int s, uint8_t *result,
 720                     const struct ieee_format *fmt)
 721 {
 722     fp_limb mant[MANT_LIMBS];
 723     int32_t exponent = 0;
 724     const int32_t expmax = 1 << (fmt->exponent - 1);
 725     fp_limb one_mask = LIMB_TOP_BIT >>
 726         ((fmt->exponent+fmt->explicit) % LIMB_BITS);
 727     const int one_pos = (fmt->exponent+fmt->explicit)/LIMB_BITS;
 728     int i;
 729     int shift;
 730     enum floats type;
 731     bool ok;
 732     const bool minus = s < 0;
 733     const int bits = fmt->bytes * 8;
 734     const char *strend;
 735
 736     if (!str[0]) {
 737         error(ERR_PANIC,
 738               "internal errror: empty string passed to float_const");
 739         return 0;
 740     }
 741
 742     strend = strchr(str, '\0');
 743     if (strend[-1] == 'P' || strend[-1] == 'p')
 744         return to_packed_bcd(str, strend-2, s, result, fmt);
 745
 746     if (str[0] == '_') {
 747         /* Special tokens */
 748
 749         switch (str[2]) {
 750         case 'n':              /* __nan__ */
 751         case 'N':
 752         case 'q':              /* __qnan__ */
 753         case 'Q':
 754             type = FL_QNAN;
 755             break;
 756         case 's':              /* __snan__ */
 757         case 'S':
 758             type = FL_SNAN;
 759             break;
 760         case 'i':              /* __infinity__ */
 761         case 'I':
 762             type = FL_INFINITY;
 763             break;
 764         default:
 765             error(ERR_NONFATAL|ERR_PASS1,
 766                   "internal error: unknown FP constant token `%s'\n", str);
 767             type = FL_QNAN;
 768             break;
 769         }
 770     } else {
 771         if (str[0] == '0') {
 772             switch (str[1]) {
 773             case 'x': case 'X':
 774             case 'h': case 'H':
 775                 ok = ieee_flconvert_bin(str+2, 4, mant, &exponent);
 776                 break;
 777             case 'o': case 'O':
 778             case 'q': case 'Q':
 779                 ok = ieee_flconvert_bin(str+2, 3, mant, &exponent);
 780                 break;
 781             case 'b': case 'B':
 782             case 'y': case 'Y':
 783                 ok = ieee_flconvert_bin(str+2, 1, mant, &exponent);
 784                 break;
 785             case 'd': case 'D':
 786             case 't': case 'T':
 787                 ok = ieee_flconvert(str+2, mant, &exponent);
 788                 break;
 789             case 'p': case 'P':
 790                 return to_packed_bcd(str+2, strend-1, s, result, fmt);
 791             default:
 792                 /* Leading zero was just a zero? */
 793                 ok = ieee_flconvert(str, mant, &exponent);
 794                 break;
 795             }
 796         } else if (str[0] == '$') {
 797             ok = ieee_flconvert_bin(str+1, 4, mant, &exponent);
 798         } else {
 799             ok = ieee_flconvert(str, mant, &exponent);
 800         }
 801
 802         if (!ok) {
 803             type = FL_QNAN;
 804         } else if (mant[0] & LIMB_TOP_BIT) {
 805             /*
 806              * Non-zero.
 807              */
 808             exponent--;
 809             if (exponent >= 2 - expmax && exponent <= expmax) {
 810                 type = FL_NORMAL;
 811             } else if (exponent > 0) {
 812                 if (pass0 == 1)
 813                     error(ERR_WARNING|ERR_WARN_FL_OVERFLOW|ERR_PASS1,
 814                           "overflow in floating-point constant");
 815                 type = FL_INFINITY;
 816             } else {
 817                 /* underflow or denormal; the denormal code handles
 818                    actual underflow. */
 819                 type = FL_DENORMAL;
 820             }
 821         } else {
 822             /* Zero */
 823             type = FL_ZERO;
 824         }
 825     }
 826
 827     switch (type) {
 828     case FL_ZERO:
 829     zero:
 830         memset(mant, 0, sizeof mant);
 831         break;
 832
 833     case FL_DENORMAL:
 834     {
 835         shift = -(exponent + expmax - 2 - fmt->exponent)
 836             + fmt->explicit;
 837         ieee_shr(mant, shift);
 838         ieee_round(minus, mant, bits);
 839         if (mant[one_pos] & one_mask) {
 840             /* One's position is set, we rounded up into normal range */
 841             exponent = 1;
 842             if (!fmt->explicit)
 843                 mant[one_pos] &= ~one_mask;     /* remove explicit one */
 844             mant[0] |= exponent << (LIMB_BITS-1 - fmt->exponent);
 845         } else {
 846             if (daz || is_zero(mant)) {
 847                 /* Flush denormals to zero */
 848                 error(ERR_WARNING|ERR_WARN_FL_UNDERFLOW|ERR_PASS1,
 849                       "underflow in floating-point constant");
 850                 goto zero;
 851             } else {
 852                 error(ERR_WARNING|ERR_WARN_FL_DENORM|ERR_PASS1,
 853                       "denormal floating-point constant");
 854             }
 855         }
 856         break;
 857     }
 858
 859     case FL_NORMAL:
 860         exponent += expmax - 1;
 861         ieee_shr(mant, fmt->exponent+fmt->explicit);
 862         ieee_round(minus, mant, bits);
 863         /* did we scale up by one? */
 864         if (test_bit(mant, fmt->exponent+fmt->explicit-1)) {
 865             ieee_shr(mant, 1);
 866             exponent++;
 867             if (exponent >= (expmax << 1)-1) {
 868                     error(ERR_WARNING|ERR_WARN_FL_OVERFLOW|ERR_PASS1,
 869                           "overflow in floating-point constant");
 870                 type = FL_INFINITY;
 871                 goto overflow;
 872             }
 873         }
 874
 875         if (!fmt->explicit)
 876             mant[one_pos] &= ~one_mask; /* remove explicit one */
 877         mant[0] |= exponent << (LIMB_BITS-1 - fmt->exponent);
 878         break;
 879
 880     case FL_INFINITY:
 881     case FL_QNAN:
 882     case FL_SNAN:
 883     overflow:
 884         memset(mant, 0, sizeof mant);
 885         mant[0] = (((fp_limb)1 << fmt->exponent)-1)
 886             << (LIMB_BITS-1 - fmt->exponent);
 887         if (fmt->explicit)
 888             mant[one_pos] |= one_mask;
 889         if (type == FL_QNAN)
 890             set_bit(mant, fmt->exponent+fmt->explicit+1);
 891         else if (type == FL_SNAN)
 892             set_bit(mant, fmt->exponent+fmt->explicit+fmt->mantissa);
 893         break;
 894     }
 895
 896     mant[0] |= minus ? LIMB_TOP_BIT : 0;
 897
 898     for (i = fmt->bytes - 1; i >= 0; i--)
 899         *result++ = mant[i/LIMB_BYTES] >> (((LIMB_BYTES-1)-(i%LIMB_BYTES))*8);
 900
 901     return 1;                   /* success */
 902 }
 903
 904 int float_const(const char *number, int sign, uint8_t *result,
 905                 int bytes, efunc err)
 906 {
 907     error = err;
 908
 909     switch (bytes) {
 910     case 1:
 911         return to_float(number, sign, result, &ieee_8);
 912     case 2:
 913         return to_float(number, sign, result, &ieee_16);
 914     case 4:
 915         return to_float(number, sign, result, &ieee_32);
 916     case 8:
 917         return to_float(number, sign, result, &ieee_64);
 918     case 10:
 919         return to_float(number, sign, result, &ieee_80);
 920     case 16:
 921         return to_float(number, sign, result, &ieee_128);
 922     default:
 923         error(ERR_PANIC, "strange value %d passed to float_const", bytes);
 924         return 0;
 925     }
 926 }
 927
 928 /* Set floating-point options */
 929 int float_option(const char *option)
 930 {
 931     if (!nasm_stricmp(option, "daz")) {
 932         daz = true;
 933         return 0;
 934     } else if (!nasm_stricmp(option, "nodaz")) {
 935         daz = false;
 936         return 0;
 937     } else if (!nasm_stricmp(option, "near")) {
 938         rc = FLOAT_RC_NEAR;
 939         return 0;
 940     } else if (!nasm_stricmp(option, "down")) {
 941         rc = FLOAT_RC_DOWN;
 942         return 0;
 943     } else if (!nasm_stricmp(option, "up")) {
 944         rc = FLOAT_RC_UP;
 945         return 0;
 946     } else if (!nasm_stricmp(option, "zero")) {
 947         rc = FLOAT_RC_ZERO;
 948         return 0;
 949     } else if (!nasm_stricmp(option, "default")) {
 950         rc = FLOAT_RC_NEAR;
 951         daz = false;
 952         return 0;
 953     } else {
 954         return -1;              /* Unknown option */
 955     }
 956 }