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)=(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     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             shr(mant, 11);
 235             round(mant, 4);
 236             if (mant[0] & 0x20)        /* did we scale up by one? */
 237                 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             shr(mant, sh);
 250             if (round(mant, 4-wds) || (sh>0 && (mant[0]&(0x8000>>(sh-1))))) {
 251                 shr(mant, 1);
 252                 if (sh==0)
 253                     mant[0] |= 0x8000;
 254                 exponent++;
 255             }
 256             put(result+6,(wds == 0 ? mant[0] : 0) | sign);
 257             put(result+4,(wds <= 1 ? mant[1-wds] : 0));
 258             put(result+2,(wds <= 2 ? mant[2-wds] : 0));
 259             put(result+0,(wds <= 3 ? mant[3-wds] : 0));
 260         } else {
 261             if (exponent > 0) {
 262                 error(ERR_NONFATAL, "overflow in floating-point constant");
 263                 return 0;
 264             } else
 265                 memset (result, 0, 8);
 266         }
 267     } else {
 268         /*
 269          * Zero.
 270          */
 271         memset (result, 0, 8);
 272     }
 273     return 1;                          /* success */
 274 }
 275
 276 static int to_float(char *str, long sign, unsigned char *result,
 277                     efunc error)
 278 {
 279     unsigned short mant[MANT_WORDS];
 280     long exponent;
 281
 282     sign = (sign < 0 ? 0x8000L : 0L);
 283
 284     flconvert (str, mant, &exponent, error);
 285     if (mant[0] & 0x8000) {
 286         /*
 287          * Non-zero.
 288          */
 289         exponent--;
 290         if (exponent >= -126 && exponent <= 128) {
 291             /*
 292              * Normalised.
 293              */
 294             exponent += 127;
 295             shr(mant, 8);
 296             round(mant, 2);
 297             if (mant[0] & 0x100)       /* did we scale up by one? */
 298                 shr(mant, 1), exponent++;
 299             mant[0] &= 0x7F;           /* remove leading one */
 300             put(result+2,(exponent << 7) | mant[0] | sign);
 301             put(result+0,mant[1]);
 302         } else if (exponent < -126 && exponent >= -149) {
 303             /*
 304              * Denormal.
 305              */
 306             int shift = -(exponent+118);
 307             int sh = shift % 16, wds = shift / 16;
 308             shr(mant, sh);
 309             if (round(mant, 2-wds) || (sh>0 && (mant[0]&(0x8000>>(sh-1))))) {
 310                 shr(mant, 1);
 311                 if (sh==0)
 312                     mant[0] |= 0x8000;
 313                 exponent++;
 314             }
 315             put(result+2,(wds == 0 ? mant[0] : 0) | sign);
 316             put(result+0,(wds <= 1 ? mant[1-wds] : 0));
 317         } else {
 318             if (exponent > 0) {
 319                 error(ERR_NONFATAL, "overflow in floating-point constant");
 320                 return 0;
 321             } else
 322                 memset (result, 0, 4);
 323         }
 324     } else {
 325         memset (result, 0, 4);
 326     }
 327     return 1;
 328 }
 329
 330 static int to_ldoub(char *str, long sign, unsigned char *result,
 331                     efunc error)
 332 {
 333     unsigned short mant[MANT_WORDS];
 334     long exponent;
 335
 336     sign = (sign < 0 ? 0x8000L : 0L);
 337
 338     flconvert (str, mant, &exponent, error);
 339     if (mant[0] & 0x8000) {
 340         /*
 341          * Non-zero.
 342          */
 343         exponent--;
 344         if (exponent >= -16383 && exponent <= 16384) {
 345             /*
 346              * Normalised.
 347              */
 348             exponent += 16383;
 349             if (round(mant, 4))        /* did we scale up by one? */
 350                 shr(mant, 1), mant[0] |= 0x8000, exponent++;
 351             put(result+8,exponent | sign);
 352             put(result+6,mant[0]);
 353             put(result+4,mant[1]);
 354             put(result+2,mant[2]);
 355             put(result+0,mant[3]);
 356         } else if (exponent < -16383 && exponent >= -16446) {
 357             /*
 358              * Denormal.
 359              */
 360             int shift = -(exponent+16383);
 361             int sh = shift % 16, wds = shift / 16;
 362             shr(mant, sh);
 363             if (round(mant, 4-wds) || (sh>0 && (mant[0]&(0x8000>>(sh-1))))) {
 364                 shr(mant, 1);
 365                 if (sh==0)
 366                     mant[0] |= 0x8000;
 367                 exponent++;
 368             }
 369             put(result+8,sign);
 370             put(result+6,(wds == 0 ? mant[0] : 0));
 371             put(result+4,(wds <= 1 ? mant[1-wds] : 0));
 372             put(result+2,(wds <= 2 ? mant[2-wds] : 0));
 373             put(result+0,(wds <= 3 ? mant[3-wds] : 0));
 374         } else {
 375             if (exponent > 0) {
 376                 error(ERR_NONFATAL, "overflow in floating-point constant");
 377                 return 0;
 378             } else
 379                 memset (result, 0, 10);
 380         }
 381     } else {
 382         /*
 383          * Zero.
 384          */
 385         memset (result, 0, 10);
 386     }
 387     return 1;
 388 }
 389
 390 int float_const (char *number, long sign, unsigned char *result, int bytes,
 391                  efunc error)
 392 {
 393     if (bytes == 4)
 394         return to_float (number, sign, result, error);
 395     else if (bytes == 8)
 396         return to_double (number, sign, result, error);
 397     else if (bytes == 10)
 398         return to_ldoub (number, sign, result, error);
 399     else {
 400         error(ERR_PANIC, "strange value %d passed to float_const", bytes);
 401         return 0;
 402     }
 403 }