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 <stdio.h>
  12 #include <stdlib.h>
  13 #include <string.h>
  14 #include <inttypes.h>
  15
  16 #include "nasm.h"
  17
  18 #define TRUE 1
  19 #define FALSE 0
  20
  21 #define MANT_WORDS 6            /* 64 bits + 32 for accuracy == 96 */
  22 #define MANT_DIGITS 28          /* 29 digits don't fit in 96 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         for (i = 0; i < MANT_WORDS; i++)
  51             to[i] = temp[i] & 0xFFFF;
  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 void ieee_flconvert(char *string, uint16_t *mant,
  61                            int32_t *exponent, efunc error)
  62 {
  63     char digits[MANT_DIGITS];
  64     char *p, *q, *r;
  65     uint16_t mult[MANT_WORDS], bit;
  66     uint16_t *m;
  67     int32_t tenpwr, twopwr;
  68     int extratwos, started, seendot;
  69
  70     p = digits;
  71     tenpwr = 0;
  72     started = seendot = FALSE;
  73     while (*string && *string != 'E' && *string != 'e') {
  74         if (*string == '.') {
  75             if (!seendot)
  76                 seendot = TRUE;
  77             else {
  78                 error(ERR_NONFATAL,
  79                       "too many periods in floating-point constant");
  80                 return;
  81             }
  82         } else if (*string >= '0' && *string <= '9') {
  83             if (*string == '0' && !started) {
  84                 if (seendot)
  85                     tenpwr--;
  86             } else {
  87                 started = TRUE;
  88                 if (p < digits + sizeof(digits))
  89                     *p++ = *string - '0';
  90                 if (!seendot)
  91                     tenpwr++;
  92             }
  93         } else {
  94             error(ERR_NONFATAL,
  95                   "floating-point constant: `%c' is invalid character",
  96                   *string);
  97             return;
  98         }
  99         string++;
 100     }
 101     if (*string) {
 102         string++;               /* eat the E */
 103         tenpwr += atoi(string);
 104     }
 105
 106     /*
 107      * At this point, the memory interval [digits,p) contains a
 108      * series of decimal digits zzzzzzz such that our number X
 109      * satisfies
 110      *
 111      * X = 0.zzzzzzz * 10^tenpwr
 112      */
 113
 114     bit = 0x8000;
 115     for (m = mant; m < mant + MANT_WORDS; m++)
 116         *m = 0;
 117     m = mant;
 118     q = digits;
 119     started = FALSE;
 120     twopwr = 0;
 121     while (m < mant + MANT_WORDS) {
 122         uint16_t carry = 0;
 123         while (p > q && !p[-1])
 124             p--;
 125         if (p <= q)
 126             break;
 127         for (r = p; r-- > q;) {
 128             int i;
 129
 130             i = 2 * *r + carry;
 131             if (i >= 10)
 132                 carry = 1, i -= 10;
 133             else
 134                 carry = 0;
 135             *r = i;
 136         }
 137         if (carry)
 138             *m |= bit, started = TRUE;
 139         if (started) {
 140             if (bit == 1)
 141                 bit = 0x8000, m++;
 142             else
 143                 bit >>= 1;
 144         } else
 145             twopwr--;
 146     }
 147     twopwr += tenpwr;
 148
 149     /*
 150      * At this point the `mant' array contains the first six
 151      * fractional places of a base-2^16 real number, which when
 152      * multiplied by 2^twopwr and 5^tenpwr gives X. So now we
 153      * really do multiply by 5^tenpwr.
 154      */
 155
 156     if (tenpwr < 0) {
 157         for (m = mult; m < mult + MANT_WORDS; m++)
 158             *m = 0xCCCC;
 159         extratwos = -2;
 160         tenpwr = -tenpwr;
 161     } else if (tenpwr > 0) {
 162         mult[0] = 0xA000;
 163         for (m = mult + 1; m < mult + MANT_WORDS; m++)
 164             *m = 0;
 165         extratwos = 3;
 166     } else
 167         extratwos = 0;
 168     while (tenpwr) {
 169         if (tenpwr & 1)
 170             twopwr += extratwos + ieee_multiply(mant, mult);
 171         extratwos = extratwos * 2 + ieee_multiply(mult, mult);
 172         tenpwr >>= 1;
 173     }
 174
 175     /*
 176      * Conversion is done. The elements of `mant' contain the first
 177      * fractional places of a base-2^16 real number in [0.5,1)
 178      * which we can multiply by 2^twopwr to get X. Or, of course,
 179      * it contains zero.
 180      */
 181     *exponent = twopwr;
 182 }
 183
 184 /*
 185  * Shift a mantissa to the right by i (i < 16) bits.
 186  */
 187 static void ieee_shr(uint16_t *mant, int i)
 188 {
 189     uint16_t n = 0, m;
 190     int j;
 191
 192     for (j = 0; j < MANT_WORDS; j++) {
 193         m = (mant[j] << (16 - i)) & 0xFFFF;
 194         mant[j] = (mant[j] >> i) | n;
 195         n = m;
 196     }
 197 }
 198
 199 /*
 200  * Round a mantissa off after i words.
 201  */
 202 static int ieee_round(uint16_t *mant, int i)
 203 {
 204     if (mant[i] & 0x8000) {
 205         do {
 206             ++mant[--i];
 207             mant[i] &= 0xFFFF;
 208         } while (i > 0 && !mant[i]);
 209         return !i && !mant[i];
 210     }
 211     return 0;
 212 }
 213
 214 #define put(a,b) ( (*(a)=(b)), ((a)[1]=(b)>>8) )
 215
 216 static int to_double(char *str, int32_t sign, uint8_t *result,
 217                      efunc error)
 218 {
 219     uint16_t mant[MANT_WORDS];
 220     int32_t exponent;
 221
 222     sign = (sign < 0 ? 0x8000L : 0L);
 223
 224     ieee_flconvert(str, mant, &exponent, error);
 225     if (mant[0] & 0x8000) {
 226         /*
 227          * Non-zero.
 228          */
 229         exponent--;
 230         if (exponent >= -1022 && exponent <= 1024) {
 231             /*
 232              * Normalised.
 233              */
 234             exponent += 1023;
 235             ieee_shr(mant, 11);
 236             ieee_round(mant, 4);
 237             if (mant[0] & 0x20) /* did we scale up by one? */
 238                 ieee_shr(mant, 1), exponent++;
 239             mant[0] &= 0xF;     /* remove leading one */
 240             put(result + 6, (exponent << 4) | mant[0] | sign);
 241             put(result + 4, mant[1]);
 242             put(result + 2, mant[2]);
 243             put(result + 0, mant[3]);
 244         } else if (exponent < -1022 && exponent >= -1074) {
 245             /*
 246              * Denormal.
 247              */
 248             int shift = -(exponent + 1011);
 249             int sh = shift % 16, wds = shift / 16;
 250             ieee_shr(mant, sh);
 251             if (ieee_round(mant, 4 - wds)
 252                 || (sh > 0 && (mant[0] & (0x8000 >> (sh - 1))))) {
 253                 ieee_shr(mant, 1);
 254                 if (sh == 0)
 255                     mant[0] |= 0x8000;
 256                 exponent++;
 257             }
 258             put(result + 6, (wds == 0 ? mant[0] : 0) | sign);
 259             put(result + 4, (wds <= 1 ? mant[1 - wds] : 0));
 260             put(result + 2, (wds <= 2 ? mant[2 - wds] : 0));
 261             put(result + 0, (wds <= 3 ? mant[3 - wds] : 0));
 262         } else {
 263             if (exponent > 0) {
 264                 error(ERR_NONFATAL, "overflow in floating-point constant");
 265                 return 0;
 266             } else
 267                 memset(result, 0, 8);
 268         }
 269     } else {
 270         /*
 271          * Zero.
 272          */
 273         memset(result, 0, 8);
 274     }
 275     return 1;                   /* success */
 276 }
 277
 278 static int to_float(char *str, int32_t sign, uint8_t *result,
 279                     efunc error)
 280 {
 281     uint16_t mant[MANT_WORDS];
 282     int32_t exponent;
 283
 284     sign = (sign < 0 ? 0x8000L : 0L);
 285
 286     ieee_flconvert(str, mant, &exponent, error);
 287     if (mant[0] & 0x8000) {
 288         /*
 289          * Non-zero.
 290          */
 291         exponent--;
 292         if (exponent >= -126 && exponent <= 128) {
 293             /*
 294              * Normalised.
 295              */
 296             exponent += 127;
 297             ieee_shr(mant, 8);
 298             ieee_round(mant, 2);
 299             if (mant[0] & 0x100)        /* did we scale up by one? */
 300                 ieee_shr(mant, 1), exponent++;
 301             mant[0] &= 0x7F;    /* remove leading one */
 302             put(result + 2, (exponent << 7) | mant[0] | sign);
 303             put(result + 0, mant[1]);
 304         } else if (exponent < -126 && exponent >= -149) {
 305             /*
 306              * Denormal.
 307              */
 308             int shift = -(exponent + 118);
 309             int sh = shift % 16, wds = shift / 16;
 310             ieee_shr(mant, sh);
 311             if (ieee_round(mant, 2 - wds)
 312                 || (sh > 0 && (mant[0] & (0x8000 >> (sh - 1))))) {
 313                 ieee_shr(mant, 1);
 314                 if (sh == 0)
 315                     mant[0] |= 0x8000;
 316                 exponent++;
 317             }
 318             put(result + 2, (wds == 0 ? mant[0] : 0) | sign);
 319             put(result + 0, (wds <= 1 ? mant[1 - wds] : 0));
 320         } else {
 321             if (exponent > 0) {
 322                 error(ERR_NONFATAL, "overflow in floating-point constant");
 323                 return 0;
 324             } else
 325                 memset(result, 0, 4);
 326         }
 327     } else {
 328         memset(result, 0, 4);
 329     }
 330     return 1;
 331 }
 332
 333 static int to_ldoub(char *str, int32_t sign, uint8_t *result,
 334                     efunc error)
 335 {
 336     uint16_t mant[MANT_WORDS];
 337     int32_t exponent;
 338
 339     sign = (sign < 0 ? 0x8000L : 0L);
 340
 341     ieee_flconvert(str, mant, &exponent, error);
 342     if (mant[0] & 0x8000) {
 343         /*
 344          * Non-zero.
 345          */
 346         exponent--;
 347         if (exponent >= -16383 && exponent <= 16384) {
 348             /*
 349              * Normalised.
 350              */
 351             exponent += 16383;
 352             if (ieee_round(mant, 4))    /* did we scale up by one? */
 353                 ieee_shr(mant, 1), mant[0] |= 0x8000, exponent++;
 354             put(result + 8, exponent | sign);
 355             put(result + 6, mant[0]);
 356             put(result + 4, mant[1]);
 357             put(result + 2, mant[2]);
 358             put(result + 0, mant[3]);
 359         } else if (exponent < -16383 && exponent >= -16446) {
 360             /*
 361              * Denormal.
 362              */
 363             int shift = -(exponent + 16383);
 364             int sh = shift % 16, wds = shift / 16;
 365             ieee_shr(mant, sh);
 366             if (ieee_round(mant, 4 - wds)
 367                 || (sh > 0 && (mant[0] & (0x8000 >> (sh - 1))))) {
 368                 ieee_shr(mant, 1);
 369                 if (sh == 0)
 370                     mant[0] |= 0x8000;
 371                 exponent++;
 372             }
 373             put(result + 8, sign);
 374             put(result + 6, (wds == 0 ? mant[0] : 0));
 375             put(result + 4, (wds <= 1 ? mant[1 - wds] : 0));
 376             put(result + 2, (wds <= 2 ? mant[2 - wds] : 0));
 377             put(result + 0, (wds <= 3 ? mant[3 - wds] : 0));
 378         } else {
 379             if (exponent > 0) {
 380                 error(ERR_NONFATAL, "overflow in floating-point constant");
 381                 return 0;
 382             } else
 383                 memset(result, 0, 10);
 384         }
 385     } else {
 386         /*
 387          * Zero.
 388          */
 389         memset(result, 0, 10);
 390     }
 391     return 1;
 392 }
 393
 394 int float_const(char *number, int32_t sign, uint8_t *result, int bytes,
 395                 efunc error)
 396 {
 397     if (bytes == 4)
 398         return to_float(number, sign, result, error);
 399     else if (bytes == 8)
 400         return to_double(number, sign, result, error);
 401     else if (bytes == 10)
 402         return to_ldoub(number, sign, result, error);
 403     else {
 404         error(ERR_PANIC, "strange value %d passed to float_const", bytes);
 405         return 0;
 406     }
 407 }