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 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 flconvert(char *string, unsigned short *mant, long *exponent,
  60                       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 + multiply (mant, mult);
 170         extratwos = extratwos * 2 + 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 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 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)[0]=(unsigned char)((b) & 0xFF)), \
 214                    ((a)[1]=(unsigned char)(((b) >> 8) & 0xFF)) )
 215
 216 static int to_double(char *str, long sign, unsigned char *result,
 217                      efunc error)
 218 {
 219     unsigned short mant[MANT_WORDS];
 220     long exponent;
 221
 222     sign = (sign < 0 ? 0x8000L : 0L);
 223
 224     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             shr(mant, 11);
 236             round(mant, 4);
 237             if (mant[0] & 0x20)        /* did we scale up by one? */
 238                 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             shr(mant, sh);
 251             if (round(mant, 4-wds) || (sh>0 && (mant[0]&(0x8000>>(sh-1))))) {
 252                 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     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             shr(mant, 8);
 297             round(mant, 2);
 298             if (mant[0] & 0x100)       /* did we scale up by one? */
 299                 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             shr(mant, sh);
 310             if (round(mant, 2-wds) || (sh>0 && (mant[0]&(0x8000>>(sh-1))))) {
 311                 shr(mant, 1);
 312                 if (sh==0)
 313                     mant[0] |= 0x8000;
 314                 exponent++;
 315             }
 316             put(result+2,(wds == 0 ? mant[0] : 0) | sign);
 317             put(result+0,(wds <= 1 ? mant[1-wds] : 0));
 318         } else {
 319             if (exponent > 0) {
 320                 error(ERR_NONFATAL, "overflow in floating-point constant");
 321                 return 0;
 322             } else
 323                 memset (result, 0, 4);
 324         }
 325     } else {
 326         memset (result, 0, 4);
 327     }
 328     return 1;
 329 }
 330
 331 static int to_ldoub(char *str, long sign, unsigned char *result,
 332                     efunc error)
 333 {
 334     unsigned short mant[MANT_WORDS];
 335     long exponent;
 336
 337     sign = (sign < 0 ? 0x8000L : 0L);
 338
 339     flconvert (str, mant, &exponent, error);
 340     if (mant[0] & 0x8000) {
 341         /*
 342          * Non-zero.
 343          */
 344         exponent--;
 345         if (exponent >= -16383 && exponent <= 16384) {
 346             /*
 347              * Normalised.
 348              */
 349             exponent += 16383;
 350             if (round(mant, 4))        /* did we scale up by one? */
 351                 shr(mant, 1), mant[0] |= 0x8000, exponent++;
 352             put(result+8,exponent | sign);
 353             put(result+6,mant[0]);
 354             put(result+4,mant[1]);
 355             put(result+2,mant[2]);
 356             put(result+0,mant[3]);
 357         } else if (exponent < -16383 && exponent >= -16446) {
 358             /*
 359              * Denormal.
 360              */
 361             int shift = -(exponent+16383);
 362             int sh = shift % 16, wds = shift / 16;
 363             shr(mant, sh);
 364             if (round(mant, 4-wds) || (sh>0 && (mant[0]&(0x8000>>(sh-1))))) {
 365                 shr(mant, 1);
 366                 if (sh==0)
 367                     mant[0] |= 0x8000;
 368                 exponent++;
 369             }
 370             put(result+8,sign);
 371             put(result+6,(wds == 0 ? mant[0] : 0));
 372             put(result+4,(wds <= 1 ? mant[1-wds] : 0));
 373             put(result+2,(wds <= 2 ? mant[2-wds] : 0));
 374             put(result+0,(wds <= 3 ? mant[3-wds] : 0));
 375         } else {
 376             if (exponent > 0) {
 377                 error(ERR_NONFATAL, "overflow in floating-point constant");
 378                 return 0;
 379             } else
 380                 memset (result, 0, 10);
 381         }
 382     } else {
 383         /*
 384          * Zero.
 385          */
 386         memset (result, 0, 10);
 387     }
 388     return 1;
 389 }
 390
 391 int float_const (char *number, long sign, unsigned char *result, int bytes,
 392                  efunc error)
 393 {
 394     if (bytes == 4)
 395         return to_float (number, sign, result, error);
 396     else if (bytes == 8)
 397         return to_double (number, sign, result, error);
 398     else if (bytes == 10)
 399         return to_ldoub (number, sign, result, error);
 400     else {
 401         error(ERR_PANIC, "strange value %d passed to float_const", bytes);
 402         return 0;
 403     }
 404 }