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