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