i386/nasm/float.c

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