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 /* Produce standard IEEE formats, with implicit "1" bit; this makes
 307    the following assumptions:
 308
 309    - the sign bit is the MSB, followed by the exponent.
 310    - the sign bit plus exponent fit in 16 bits.
 311    - the exponent bias is 2^(n-1)-1 for an n-bit exponent */
 312
 313 struct ieee_format {
 314     int words;
 315     int mantissa;               /* Bits in the mantissa */
 316     int exponent;               /* Bits in the exponent */
 317 };
 318
 319 static const struct ieee_format ieee_16  = { 1,  10,  5 };
 320 static const struct ieee_format ieee_32  = { 2,  23,  8 };
 321 static const struct ieee_format ieee_64  = { 4,  52, 11 };
 322 static const struct ieee_format ieee_128 = { 8, 112, 15 };
 323
 324 /* Produce all the standard IEEE formats: 16, 32, 64, and 128 bits */
 325 static int to_float(char *str, int32_t sign, uint8_t *result,
 326                     const struct ieee_format *fmt, efunc error)
 327 {
 328     uint16_t mant[MANT_WORDS], *mp;
 329     int32_t exponent;
 330     int32_t expmax = 1 << (fmt->exponent-1);
 331     uint16_t implicit_one = 0x8000 >> fmt->exponent;
 332     int i;
 333
 334     sign = (sign < 0 ? 0x8000L : 0L);
 335
 336     ieee_flconvert(str, mant, &exponent, error);
 337     if (mant[0] & 0x8000) {
 338         /*
 339          * Non-zero.
 340          */
 341         exponent--;
 342         if (exponent >= 2-expmax && exponent <= expmax) {
 343             /*
 344              * Normalised.
 345              */
 346             exponent += expmax;
 347             ieee_shr(mant, fmt->exponent);
 348             ieee_round(mant, fmt->words);
 349             /* did we scale up by one? */
 350             if (mant[0] & (implicit_one << 1)) {
 351                 ieee_shr(mant, 1);
 352                 exponent++;
 353             }
 354
 355             mant[0] &= (implicit_one-1);     /* remove leading one */
 356             mant[0] |= exponent << (15 - fmt->exponent);
 357         } else if (exponent < 2-expmax && exponent >= 2-expmax-fmt->mantissa) {
 358             /*
 359              * Denormal.
 360              */
 361             int shift = -(exponent + expmax-2-fmt->exponent);
 362             int sh = shift % 16, wds = shift / 16;
 363             ieee_shr(mant, sh);
 364             if (ieee_round(mant, fmt->words - wds)
 365                 || (sh > 0 && (mant[0] & (0x8000 >> (sh - 1))))) {
 366                 ieee_shr(mant, 1);
 367                 if (sh == 0)
 368                     mant[0] |= 0x8000;
 369                 exponent++;
 370             }
 371
 372             if (wds) {
 373                 for (i = fmt->words-1; i >= wds; i--)
 374                     mant[i] = mant[i-wds];
 375                 for (; i >= 0; i--)
 376                     mant[i] = 0;
 377             }
 378         } else {
 379             if (exponent > 0) {
 380                 error(ERR_NONFATAL, "overflow in floating-point constant");
 381                 return 0;
 382             } else {
 383                 memset(mant, 0, 2*fmt->words);
 384             }
 385         }
 386     } else {
 387         /* Zero */
 388         memset(mant, 0, 2*fmt->words);
 389     }
 390
 391     mant[0] |= sign;
 392
 393     for (mp = &mant[fmt->words], i = 0; i < fmt->words; i++) {
 394         uint16_t m = *--mp;
 395         put(result, m);
 396         result += 2;
 397     }
 398
 399     return 1;                   /* success */
 400 }
 401
 402 /* 80-bit format with 64-bit mantissa *including an explicit integer 1*
 403    and 15-bit exponent. */
 404 static int to_ldoub(char *str, int32_t sign, uint8_t *result,
 405                     efunc error)
 406 {
 407     uint16_t mant[MANT_WORDS];
 408     int32_t exponent;
 409
 410     sign = (sign < 0 ? 0x8000L : 0L);
 411
 412     ieee_flconvert(str, mant, &exponent, error);
 413     if (mant[0] & 0x8000) {
 414         /*
 415          * Non-zero.
 416          */
 417         exponent--;
 418         if (exponent >= -16383 && exponent <= 16384) {
 419             /*
 420              * Normalised.
 421              */
 422             exponent += 16383;
 423             if (ieee_round(mant, 4))    /* did we scale up by one? */
 424                 ieee_shr(mant, 1), mant[0] |= 0x8000, exponent++;
 425             put(result + 8, exponent | sign);
 426             put(result + 6, mant[0]);
 427             put(result + 4, mant[1]);
 428             put(result + 2, mant[2]);
 429             put(result + 0, mant[3]);
 430         } else if (exponent < -16383 && exponent >= -16446) {
 431             /*
 432              * Denormal.
 433              */
 434             int shift = -(exponent + 16383);
 435             int sh = shift % 16, wds = shift / 16;
 436             ieee_shr(mant, sh);
 437             if (ieee_round(mant, 4 - wds)
 438                 || (sh > 0 && (mant[0] & (0x8000 >> (sh - 1))))) {
 439                 ieee_shr(mant, 1);
 440                 if (sh == 0)
 441                     mant[0] |= 0x8000;
 442                 exponent++;
 443             }
 444             put(result + 8, sign);
 445             put(result + 6, (wds == 0 ? mant[0] : 0));
 446             put(result + 4, (wds <= 1 ? mant[1 - wds] : 0));
 447             put(result + 2, (wds <= 2 ? mant[2 - wds] : 0));
 448             put(result + 0, (wds <= 3 ? mant[3 - wds] : 0));
 449         } else {
 450             if (exponent > 0) {
 451                 error(ERR_NONFATAL, "overflow in floating-point constant");
 452                 return 0;
 453             } else
 454                 memset(result, 0, 10);
 455         }
 456     } else {
 457         /*
 458          * Zero.
 459          */
 460         memset(result, 0, 10);
 461     }
 462     return 1;
 463 }
 464
 465 int float_const(char *number, int32_t sign, uint8_t *result, int bytes,
 466                 efunc error)
 467 {
 468     switch (bytes) {
 469     case 2:
 470         return to_float(number, sign, result, &ieee_16, error);
 471     case 4:
 472         return to_float(number, sign, result, &ieee_32, error);
 473     case 8:
 474         return to_float(number, sign, result, &ieee_64, error);
 475     case 10:
 476         return to_ldoub(number, sign, result, error);
 477     case 16:
 478         return to_float(number, sign, result, &ieee_128, error);
 479     default:
 480         error(ERR_PANIC, "strange value %d passed to float_const", bytes);
 481         return 0;
 482     }
 483 }