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