pintos-os.org Git - pspp/blob - src/libpspp/float-format.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2006, 2011 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include "libpspp/float-format.h"
  20
  21 #include <byteswap.h>
  22 #include <ctype.h>
  23 #include <inttypes.h>
  24 #include <stdlib.h>
  25 #include <string.h>
  26
  27 #include "libpspp/assertion.h"
  28 #include "libpspp/integer-format.h"
  29
  30 \f
  31 /* Neutral intermediate representation for binary floating-point numbers. */
  32 struct fp
  33   {
  34     enum
  35       {
  36         FINITE,         /* Finite number (normalized or denormalized). */
  37         INFINITE,       /* Positive or negative infinity. */
  38         NAN,            /* Not a number. */
  39
  40         ZERO,           /* Positive or negative zero. */
  41         MISSING,        /* System missing. */
  42         LOWEST,         /* LOWEST on e.g. missing values. */
  43         HIGHEST,        /* HIGHEST on e.g. missing values. */
  44         RESERVED        /* Special Vax representation. */
  45       }
  46     class;
  47
  48     enum
  49       {
  50         POSITIVE,
  51         NEGATIVE
  52       }
  53     sign;
  54
  55     /* class == FINITE: The number has value "fraction *
  56        2**exponent", considering bit 63 in fraction to be just
  57        right of the decimal point.
  58
  59        class == NAN: The fraction holds the significand, with its
  60        leftmost bit in bit 63, so that signaling and quiet NaN
  61        values can be preserved.
  62
  63        Unused for other classes. */
  64     uint64_t fraction;
  65     int exponent;
  66   };
  67
  68 static void extract_number (enum float_format, const void *, struct fp *);
  69 static void assemble_number (enum float_format, struct fp *, void *);
  70
  71 static inline uint16_t get_uint16 (const void *);
  72 static inline uint32_t get_uint32 (const void *);
  73 static inline uint64_t get_uint64 (const void *);
  74
  75 static inline void put_uint16 (uint16_t, void *);
  76 static inline void put_uint32 (uint32_t, void *);
  77 static inline void put_uint64 (uint64_t, void *);
  78
  79 /* Converts SRC from format FROM to format TO, storing the
  80    converted value into DST.
  81    SRC and DST are permitted to arbitrarily overlap. */
  82 void
  83 float_convert (enum float_format from, const void *src,
  84                enum float_format to, void *dst)
  85 {
  86   if (from != to)
  87     {
  88       if ((from == FLOAT_IEEE_SINGLE_LE || from == FLOAT_IEEE_SINGLE_BE)
  89           && (to == FLOAT_IEEE_SINGLE_LE || to == FLOAT_IEEE_SINGLE_BE))
  90         put_uint32 (bswap_32 (get_uint32 (src)), dst);
  91       else if ((from == FLOAT_IEEE_DOUBLE_LE || from == FLOAT_IEEE_DOUBLE_BE)
  92                && (to == FLOAT_IEEE_DOUBLE_LE || to == FLOAT_IEEE_DOUBLE_BE))
  93         put_uint64 (bswap_64 (get_uint64 (src)), dst);
  94       else
  95         {
  96           struct fp fp;
  97           extract_number (from, src, &fp);
  98           assemble_number (to, &fp, dst);
  99         }
 100     }
 101   else
 102     {
 103       if (src != dst)
 104         memmove (dst, src, float_get_size (from));
 105     }
 106 }
 107
 108 /* Converts SRC from format FROM to a native double and returns
 109    the double. */
 110 double
 111 float_get_double (enum float_format from, const void *src)
 112 {
 113   double dst;
 114   float_convert (from, src, FLOAT_NATIVE_DOUBLE, &dst);
 115   return dst;
 116 }
 117
 118 /* Returns the number of bytes in a number in the given
 119    FORMAT. */
 120 size_t
 121 float_get_size (enum float_format format)
 122 {
 123   switch (format)
 124     {
 125     case FLOAT_IEEE_SINGLE_LE:
 126     case FLOAT_IEEE_SINGLE_BE:
 127     case FLOAT_VAX_F:
 128     case FLOAT_Z_SHORT:
 129       return 4;
 130
 131     case FLOAT_IEEE_DOUBLE_LE:
 132     case FLOAT_IEEE_DOUBLE_BE:
 133     case FLOAT_VAX_D:
 134     case FLOAT_VAX_G:
 135     case FLOAT_Z_LONG:
 136       return 8;
 137
 138     case FLOAT_FP:
 139       return sizeof (struct fp);
 140
 141     case FLOAT_HEX:
 142       return 32;
 143     }
 144
 145   NOT_REACHED ();
 146 }
 147
 148 /* Attempts to identify the floating point format(s) in which the
 149    LENGTH bytes in NUMBER represent the given EXPECTED_VALUE.
 150    Returns the number of matches, which may be zero, one, or
 151    greater than one.  If a positive value is returned, then the
 152    most likely candidate format (based on how common the formats
 153    are in practice) is stored in *BEST_GUESS. */
 154 int
 155 float_identify (double expected_value, const void *number, size_t length,
 156                 enum float_format *best_guess)
 157 {
 158   /* Candidates for identification in order of decreasing
 159      preference. */
 160   enum float_format candidates[] =
 161     {
 162       FLOAT_IEEE_SINGLE_LE,
 163       FLOAT_IEEE_SINGLE_BE,
 164       FLOAT_IEEE_DOUBLE_LE,
 165       FLOAT_IEEE_DOUBLE_BE,
 166       FLOAT_VAX_F,
 167       FLOAT_VAX_D,
 168       FLOAT_VAX_G,
 169       FLOAT_Z_SHORT,
 170       FLOAT_Z_LONG,
 171     };
 172   const size_t candidate_cnt = sizeof candidates / sizeof *candidates;
 173
 174   enum float_format *p;
 175   int match_cnt;
 176
 177   match_cnt = 0;
 178   for (p = candidates; p < candidates + candidate_cnt; p++)
 179     if (float_get_size (*p) == length)
 180       {
 181         char tmp[8];
 182         assert (sizeof tmp >= float_get_size (*p));
 183         float_convert (FLOAT_NATIVE_DOUBLE, &expected_value, *p, tmp);
 184         if (!memcmp (tmp, number, length) && match_cnt++ == 0)
 185           *best_guess = *p;
 186       }
 187   return match_cnt;
 188 }
 189
 190 /* Returns the double value that is just greater than -DBL_MAX,
 191    which in PSPP syntax files is called LOWEST and used as the
 192    low end of numeric ranges that are supposed to be unbounded on
 193    the low end, as in the missing value set created by,
 194    e.g. MISSING VALUES X(LOWEST THRU 5).  (-DBL_MAX is used for
 195    SYSMIS so it is not available for LOWEST.) */
 196 double
 197 float_get_lowest (void)
 198 {
 199   struct fp fp;
 200   double x;
 201
 202   fp.class = LOWEST;
 203   fp.sign = POSITIVE;
 204   assemble_number (FLOAT_NATIVE_DOUBLE, &fp, &x);
 205   return x;
 206 }
 207 \f
 208 /* Returns CNT bits in X starting from the given bit OFS. */
 209 static inline uint64_t
 210 get_bits (uint64_t x, int ofs, int cnt)
 211 {
 212   assert (ofs >= 0 && ofs < 64);
 213   assert (cnt > 0 && cnt < 64);
 214   assert (ofs + cnt <= 64);
 215   return (x >> ofs) & ((UINT64_C(1) << cnt) - 1);
 216 }
 217
 218 /* Returns the 16-bit unsigned integer at P,
 219    which need not be aligned. */
 220 static inline uint16_t
 221 get_uint16 (const void *p)
 222 {
 223   uint16_t x;
 224   memcpy (&x, p, sizeof x);
 225   return x;
 226 }
 227
 228 /* Returns the 32-bit unsigned integer at P,
 229    which need not be aligned. */
 230 static inline uint32_t
 231 get_uint32 (const void *p)
 232 {
 233   uint32_t x;
 234   memcpy (&x, p, sizeof x);
 235   return x;
 236 }
 237
 238 /* Returns the 64-bit unsigned integer at P,
 239    which need not be aligned. */
 240 static inline uint64_t
 241 get_uint64 (const void *p)
 242 {
 243   uint64_t x;
 244   memcpy (&x, p, sizeof x);
 245   return x;
 246 }
 247
 248 /* Stores 16-bit unsigned integer X at P,
 249    which need not be aligned. */
 250 static inline void
 251 put_uint16 (uint16_t x, void *p)
 252 {
 253   memcpy (p, &x, sizeof x);
 254 }
 255
 256 /* Stores 32-bit unsigned integer X at P,
 257    which need not be aligned. */
 258 static inline void
 259 put_uint32 (uint32_t x, void *p)
 260 {
 261   memcpy (p, &x, sizeof x);
 262 }
 263
 264 /* Stores 64-bit unsigned integer X at P,
 265    which need not be aligned. */
 266 static inline void
 267 put_uint64 (uint64_t x, void *p)
 268 {
 269   memcpy (p, &x, sizeof x);
 270 }
 271
 272 /* Returns NATIVE converted to a form that, when stored in
 273    memory, will be in little-endian byte order. */
 274 static inline uint16_t
 275 native_to_le16 (uint16_t native)
 276 {
 277   return INTEGER_NATIVE == INTEGER_LSB_FIRST ? native : bswap_16 (native);
 278 }
 279
 280 /* Returns NATIVE converted to a form that, when stored in
 281    memory, will be in big-endian byte order. */
 282 static inline uint16_t
 283 native_to_be16 (uint16_t native)
 284 {
 285   return INTEGER_NATIVE == INTEGER_MSB_FIRST ? native : bswap_16 (native);
 286 }
 287
 288 /* Returns NATIVE converted to a form that, when stored in
 289    memory, will be in VAX-endian byte order. */
 290 static inline uint16_t
 291 native_to_vax16 (uint16_t native)
 292 {
 293   return native_to_le16 (native);
 294 }
 295
 296 /* Returns NATIVE converted to a form that, when stored in
 297    memory, will be in little-endian byte order. */
 298 static inline uint32_t
 299 native_to_le32 (uint32_t native)
 300 {
 301   return INTEGER_NATIVE == INTEGER_LSB_FIRST ? native : bswap_32 (native);
 302 }
 303
 304 /* Returns NATIVE converted to a form that, when stored in
 305    memory, will be in big-endian byte order. */
 306 static inline uint32_t
 307 native_to_be32 (uint32_t native)
 308 {
 309   return INTEGER_NATIVE == INTEGER_MSB_FIRST ? native : bswap_32 (native);
 310 }
 311
 312 /* Returns NATIVE converted to a form that, when stored in
 313    memory, will be in VAX-endian byte order. */
 314 static inline uint32_t
 315 native_to_vax32 (uint32_t native)
 316 {
 317   return native_to_be32 (((native & 0xff00ff00) >> 8) |
 318                          ((native & 0x00ff00ff) << 8));
 319 }
 320
 321 /* Returns NATIVE converted to a form that, when stored in
 322    memory, will be in little-endian byte order. */
 323 static inline uint64_t
 324 native_to_le64 (uint64_t native)
 325 {
 326   return INTEGER_NATIVE == INTEGER_LSB_FIRST ? native : bswap_64 (native);
 327 }
 328
 329 /* Returns NATIVE converted to a form that, when stored in
 330    memory, will be in big-endian byte order. */
 331 static inline uint64_t
 332 native_to_be64 (uint64_t native)
 333 {
 334   return INTEGER_NATIVE == INTEGER_MSB_FIRST ? native : bswap_64 (native);
 335 }
 336
 337 /* Returns NATIVE converted to a form that, when stored in
 338    memory, will be in VAX-endian byte order. */
 339 static inline uint64_t
 340 native_to_vax64 (uint64_t native)
 341 {
 342   return native_to_be64 (((native & UINT64_C(0xff00ff0000000000)) >> 40) |
 343                          ((native & UINT64_C(0x00ff00ff00000000)) >> 24) |
 344                          ((native & UINT64_C(0x00000000ff00ff00)) << 24) |
 345                          ((native & UINT64_C(0x0000000000ff00ff)) << 40));
 346 }
 347
 348 /* Given LE, obtained from memory in little-endian format,
 349    returns its value. */
 350 static inline uint16_t
 351 le_to_native16 (uint16_t le)
 352 {
 353   return INTEGER_NATIVE == INTEGER_LSB_FIRST ? le : bswap_16 (le);
 354 }
 355
 356 /* Given BE, obtained from memory in big-endian format, returns
 357    its value. */
 358 static inline uint16_t
 359 be_to_native16 (uint16_t be)
 360 {
 361   return INTEGER_NATIVE == INTEGER_MSB_FIRST ? be : bswap_16 (be);
 362 }
 363
 364 /* Given VAX, obtained from memory in VAX-endian format, returns
 365    its value. */
 366 static inline uint16_t
 367 vax_to_native16 (uint16_t vax)
 368 {
 369   return le_to_native16 (vax);
 370 }
 371
 372 /* Given LE, obtained from memory in little-endian format,
 373    returns its value. */
 374 static inline uint32_t
 375 le_to_native32 (uint32_t le)
 376 {
 377   return INTEGER_NATIVE == INTEGER_LSB_FIRST ? le : bswap_32 (le);
 378 }
 379
 380 /* Given BE, obtained from memory in big-endian format, returns
 381    its value. */
 382 static inline uint32_t
 383 be_to_native32 (uint32_t be)
 384 {
 385   return INTEGER_NATIVE == INTEGER_MSB_FIRST ? be : bswap_32 (be);
 386 }
 387
 388 /* Given VAX, obtained from memory in VAX-endian format, returns
 389    its value. */
 390 static inline uint32_t
 391 vax_to_native32 (uint32_t vax)
 392 {
 393   uint32_t be = be_to_native32 (vax);
 394   return ((be & 0xff00ff00) >> 8) | ((be & 0x00ff00ff) << 8);
 395 }
 396
 397 /* Given LE, obtained from memory in little-endian format,
 398    returns its value. */
 399 static inline uint64_t
 400 le_to_native64 (uint64_t le)
 401 {
 402   return INTEGER_NATIVE == INTEGER_LSB_FIRST ? le : bswap_64 (le);
 403 }
 404
 405 /* Given BE, obtained from memory in big-endian format, returns
 406    its value. */
 407 static inline uint64_t
 408 be_to_native64 (uint64_t be)
 409 {
 410   return INTEGER_NATIVE == INTEGER_MSB_FIRST ? be : bswap_64 (be);
 411 }
 412
 413 /* Given VAX, obtained from memory in VAX-endian format, returns
 414    its value. */
 415 static inline uint64_t
 416 vax_to_native64 (uint64_t vax)
 417 {
 418   uint64_t be = be_to_native64 (vax);
 419   return (((be & UINT64_C(0xff00ff0000000000)) >> 40) |
 420           ((be & UINT64_C(0x00ff00ff00000000)) >> 24) |
 421           ((be & UINT64_C(0x00000000ff00ff00)) << 24) |
 422           ((be & UINT64_C(0x0000000000ff00ff)) << 40));
 423 }
 424 \f
 425 static void extract_ieee (uint64_t, int exp_bits, int frac_bits, struct fp *);
 426 static void extract_vax (uint64_t, int exp_bits, int frac_bits, struct fp *);
 427 static void extract_z (uint64_t, int exp_bits, int frac_bits, struct fp *);
 428 static void extract_hex (const char *, struct fp *);
 429
 430 /* Converts the number at BITS from format TYPE into neutral
 431    format at FP. */
 432 static void
 433 extract_number (enum float_format type, const void *bits, struct fp *fp)
 434 {
 435   switch (type)
 436     {
 437     case FLOAT_IEEE_SINGLE_LE:
 438       extract_ieee (le_to_native32 (get_uint32 (bits)), 8, 23, fp);
 439       break;
 440     case FLOAT_IEEE_SINGLE_BE:
 441       extract_ieee (be_to_native32 (get_uint32 (bits)), 8, 23, fp);
 442       break;
 443     case FLOAT_IEEE_DOUBLE_LE:
 444       extract_ieee (le_to_native64 (get_uint64 (bits)), 11, 52, fp);
 445       break;
 446     case FLOAT_IEEE_DOUBLE_BE:
 447       extract_ieee (be_to_native64 (get_uint64 (bits)), 11, 52, fp);
 448       break;
 449
 450     case FLOAT_VAX_F:
 451       extract_vax (vax_to_native32 (get_uint32 (bits)), 8, 23, fp);
 452       break;
 453     case FLOAT_VAX_D:
 454       extract_vax (vax_to_native64 (get_uint64 (bits)), 8, 55, fp);
 455       break;
 456     case FLOAT_VAX_G:
 457       extract_vax (vax_to_native64 (get_uint64 (bits)), 11, 52, fp);
 458       break;
 459
 460     case FLOAT_Z_SHORT:
 461       extract_z (be_to_native32 (get_uint32 (bits)), 7, 24, fp);
 462       break;
 463     case FLOAT_Z_LONG:
 464       extract_z (be_to_native64 (get_uint64 (bits)), 7, 56, fp);
 465       break;
 466
 467     case FLOAT_FP:
 468       memcpy (fp, bits, sizeof *fp);
 469       break;
 470     case FLOAT_HEX:
 471       extract_hex (bits, fp);
 472       break;
 473     }
 474
 475   assert (!(fp->class == FINITE && fp->fraction == 0));
 476 }
 477
 478 /* Converts the IEEE format number in BITS, which has EXP_BITS of
 479    exponent and FRAC_BITS of fraction, into neutral format at
 480    FP. */
 481 static void
 482 extract_ieee (uint64_t bits, int exp_bits, int frac_bits, struct fp *fp)
 483 {
 484   const int bias = (1 << (exp_bits - 1)) - 1;
 485   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 486   const int max_raw_exp = (1 << exp_bits) - 1;
 487
 488   const uint64_t raw_frac = get_bits (bits, 0, frac_bits);
 489   const int raw_exp = get_bits (bits, frac_bits, exp_bits);
 490   const bool raw_sign = get_bits (bits, frac_bits + exp_bits, 1);
 491
 492   if (raw_sign && raw_exp == max_raw_exp - 1 && raw_frac == max_raw_frac - 1)
 493     fp->class = LOWEST;
 494   else if (raw_exp == max_raw_exp - 1 && raw_frac == max_raw_frac)
 495     fp->class = raw_sign ? MISSING : HIGHEST;
 496   else if (raw_exp == max_raw_exp)
 497     {
 498       if (raw_frac == 0)
 499         fp->class = INFINITE;
 500       else
 501         {
 502           fp->class = NAN;
 503           fp->fraction = raw_frac << (64 - frac_bits);
 504         }
 505     }
 506   else if (raw_exp == 0)
 507     {
 508       if (raw_frac != 0)
 509         {
 510           fp->class = FINITE;
 511           fp->exponent = 1 - bias;
 512           fp->fraction = raw_frac << (64 - frac_bits);
 513         }
 514       else
 515         fp->class = ZERO;
 516     }
 517   else
 518     {
 519       fp->class = FINITE;
 520       fp->exponent = raw_exp - bias + 1;
 521       fp->fraction = (raw_frac << (64 - frac_bits - 1)) | (UINT64_C(1) << 63);
 522     }
 523
 524   fp->sign = raw_sign ? NEGATIVE : POSITIVE;
 525 }
 526
 527 /* Converts the VAX format number in BITS, which has EXP_BITS of
 528    exponent and FRAC_BITS of fraction, into neutral format at
 529    FP. */
 530 static void
 531 extract_vax (uint64_t bits, int exp_bits, int frac_bits, struct fp *fp)
 532 {
 533   const int bias = 1 << (exp_bits - 1);
 534   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 535   const int max_raw_exp = (1 << exp_bits) - 1;
 536
 537   const uint64_t raw_frac = get_bits (bits, 0, frac_bits);
 538   const int raw_exp = get_bits (bits, frac_bits, exp_bits);
 539   const bool raw_sign = get_bits (bits, frac_bits + exp_bits, 1);
 540
 541   if (raw_sign && raw_exp == max_raw_exp && raw_frac == max_raw_frac - 1)
 542     fp->class = LOWEST;
 543   else if (raw_exp == max_raw_exp && raw_frac == max_raw_frac)
 544     fp->class = raw_sign ? MISSING : HIGHEST;
 545   else if (raw_exp == 0)
 546     fp->class = raw_sign == 0 ? ZERO : RESERVED;
 547   else
 548     {
 549       fp->class = FINITE;
 550       fp->fraction = (raw_frac << (64 - frac_bits - 1)) | (UINT64_C(1) << 63);
 551       fp->exponent = raw_exp - bias;
 552     }
 553
 554   fp->sign = raw_sign ? NEGATIVE : POSITIVE;
 555 }
 556
 557 /* Converts the Z architecture format number in BITS, which has
 558    EXP_BITS of exponent and FRAC_BITS of fraction, into neutral
 559    format at FP. */
 560 static void
 561 extract_z (uint64_t bits, int exp_bits, int frac_bits, struct fp *fp)
 562 {
 563   const int bias = 1 << (exp_bits - 1);
 564   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 565   const int max_raw_exp = (1 << exp_bits) - 1;
 566
 567   uint64_t raw_frac = get_bits (bits, 0, frac_bits);
 568   int raw_exp = get_bits (bits, frac_bits, exp_bits);
 569   int raw_sign = get_bits (bits, frac_bits + exp_bits, 1);
 570
 571   fp->sign = raw_sign ? NEGATIVE : POSITIVE;
 572   if (raw_exp == max_raw_exp && raw_frac == max_raw_frac)
 573     fp->class = raw_sign ? MISSING : HIGHEST;
 574   else if (raw_sign && raw_exp == max_raw_exp && raw_frac == max_raw_frac - 1)
 575     fp->class = LOWEST;
 576   else if (raw_frac != 0)
 577     {
 578       fp->class = FINITE;
 579       fp->fraction = raw_frac << (64 - frac_bits);
 580       fp->exponent = (raw_exp - bias) * 4;
 581     }
 582   else
 583     fp->class = ZERO;
 584 }
 585
 586 /* Returns the integer value of hex digit C. */
 587 static int
 588 hexit_value (int c)
 589 {
 590   const char s[] = "0123456789abcdef";
 591   const char *cp = strchr (s, tolower ((unsigned char) c));
 592
 593   assert (cp != NULL);
 594   return cp - s;
 595 }
 596
 597 /* Parses a hexadecimal floating point number string at S (useful
 598    for testing) into neutral format at FP. */
 599 static void
 600 extract_hex (const char *s, struct fp *fp)
 601 {
 602   if (*s == '-')
 603     {
 604       fp->sign = NEGATIVE;
 605       s++;
 606     }
 607   else
 608     fp->sign = POSITIVE;
 609
 610   if (!strcmp (s, "Infinity"))
 611     fp->class = INFINITE;
 612   else if (!strcmp (s, "Missing"))
 613     fp->class = MISSING;
 614   else if (!strcmp (s, "Lowest"))
 615     fp->class = LOWEST;
 616   else if (!strcmp (s, "Highest"))
 617     fp->class = HIGHEST;
 618   else if (!strcmp (s, "Reserved"))
 619     fp->class = RESERVED;
 620   else
 621     {
 622       int offset;
 623
 624       if (!memcmp (s, "NaN:", 4))
 625         {
 626           fp->class = NAN;
 627           s += 4;
 628         }
 629       else
 630         fp->class = FINITE;
 631
 632       if (*s == '.')
 633         s++;
 634
 635       fp->exponent = 0;
 636       fp->fraction = 0;
 637       offset = 60;
 638       for (; isxdigit ((unsigned char) *s); s++)
 639         if (offset >= 0)
 640           {
 641             uint64_t digit = hexit_value (*s);
 642             fp->fraction += digit << offset;
 643             offset -= 4;
 644           }
 645
 646       if (fp->class == FINITE)
 647         {
 648           if (fp->fraction == 0)
 649             fp->class = ZERO;
 650           else if (*s == 'p')
 651             {
 652               char *tail;
 653               fp->exponent += strtol (s + 1, &tail, 10);
 654               s = tail;
 655             }
 656         }
 657     }
 658 }
 659 \f
 660 static uint64_t assemble_ieee (struct fp *, int exp_bits, int frac_bits);
 661 static uint64_t assemble_vax (struct fp *, int exp_bits, int frac_bits);
 662 static uint64_t assemble_z (struct fp *, int exp_bits, int frac_bits);
 663 static void assemble_hex (struct fp *, void *);
 664
 665 /* Converts the neutral format floating point number in FP into
 666    format TYPE at NUMBER.  May modify FP as part of the
 667    process. */
 668 static void
 669 assemble_number (enum float_format type, struct fp *fp, void *number)
 670 {
 671   switch (type)
 672     {
 673     case FLOAT_IEEE_SINGLE_LE:
 674       put_uint32 (native_to_le32 (assemble_ieee (fp, 8, 23)), number);
 675       break;
 676     case FLOAT_IEEE_SINGLE_BE:
 677       put_uint32 (native_to_be32 (assemble_ieee (fp, 8, 23)), number);
 678       break;
 679     case FLOAT_IEEE_DOUBLE_LE:
 680       put_uint64 (native_to_le64 (assemble_ieee (fp, 11, 52)), number);
 681       break;
 682     case FLOAT_IEEE_DOUBLE_BE:
 683       put_uint64 (native_to_be64 (assemble_ieee (fp, 11, 52)), number);
 684       break;
 685
 686     case FLOAT_VAX_F:
 687       put_uint32 (native_to_vax32 (assemble_vax (fp, 8, 23)), number);
 688       break;
 689     case FLOAT_VAX_D:
 690       put_uint64 (native_to_vax64 (assemble_vax (fp, 8, 55)), number);
 691       break;
 692     case FLOAT_VAX_G:
 693       put_uint64 (native_to_vax64 (assemble_vax (fp, 11, 52)), number);
 694       break;
 695
 696     case FLOAT_Z_SHORT:
 697       put_uint32 (native_to_be32 (assemble_z (fp, 7, 24)), number);
 698       break;
 699     case FLOAT_Z_LONG:
 700       put_uint64 (native_to_be64 (assemble_z (fp, 7, 56)), number);
 701       break;
 702
 703     case FLOAT_FP:
 704       memcpy (number, fp, sizeof *fp);
 705       break;
 706     case FLOAT_HEX:
 707       assemble_hex (fp, number);
 708       break;
 709     }
 710 }
 711
 712 /* Rounds off FP's fraction to FRAC_BITS bits of precision.
 713    Halfway values are rounded to even. */
 714 static void
 715 normalize_and_round_fp (struct fp *fp, int frac_bits)
 716 {
 717   assert (fp->class == FINITE);
 718   assert (fp->fraction != 0);
 719
 720   /* Make sure that the leading fraction bit is 1. */
 721   while (!(fp->fraction & (UINT64_C(1) << 63)))
 722     {
 723       fp->fraction <<= 1;
 724       fp->exponent--;
 725     }
 726
 727   if (frac_bits < 64)
 728     {
 729       uint64_t last_frac_bit = UINT64_C(1) << (64 - frac_bits);
 730       uint64_t decision_bit = last_frac_bit >> 1;
 731       if (fp->fraction & decision_bit
 732           && (fp->fraction & (decision_bit - 1)
 733               || fp->fraction & last_frac_bit))
 734         {
 735           fp->fraction += last_frac_bit;
 736           if ((fp->fraction >> 63) == 0)
 737             {
 738               fp->fraction = UINT64_C(1) << 63;
 739               fp->exponent++;
 740             }
 741         }
 742
 743       /* Mask off all but FRAC_BITS high-order bits.
 744          If we rounded up, then these bits no longer have
 745          meaningful values. */
 746       fp->fraction &= ~(last_frac_bit - 1);
 747     }
 748 }
 749
 750 /* Converts the neutral format floating point number in FP into
 751    IEEE format with EXP_BITS exponent bits and FRAC_BITS fraction
 752    bits, and returns the value. */
 753 static uint64_t
 754 assemble_ieee (struct fp *fp, int exp_bits, int frac_bits)
 755 {
 756   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 757
 758   const int bias = (1 << (exp_bits - 1)) - 1;
 759   const int max_raw_exp = (1 << exp_bits) - 1;
 760   const int min_norm_exp = 1 - bias;
 761   const int min_denorm_exp = min_norm_exp - frac_bits;
 762   const int max_norm_exp = max_raw_exp - 1 - bias;
 763
 764   uint64_t raw_frac;
 765   int raw_exp;
 766   bool raw_sign;
 767
 768   raw_sign = fp->sign != POSITIVE;
 769
 770   switch (fp->class)
 771     {
 772     case FINITE:
 773       normalize_and_round_fp (fp, frac_bits + 1);
 774       if (fp->exponent - 1 > max_norm_exp)
 775         {
 776           /* Overflow to infinity. */
 777           raw_exp = max_raw_exp;
 778           raw_frac = 0;
 779         }
 780       else if (fp->exponent - 1 >= min_norm_exp)
 781         {
 782           /* Normal. */
 783           raw_frac = (fp->fraction << 1) >> (64 - frac_bits);
 784           raw_exp = (fp->exponent - 1) + bias;
 785         }
 786       else if (fp->exponent - 1 >= min_denorm_exp)
 787         {
 788           /* Denormal. */
 789           const int denorm_shift = min_norm_exp - fp->exponent;
 790           raw_frac = (fp->fraction >> (64 - frac_bits)) >> denorm_shift;
 791           raw_exp = 0;
 792         }
 793       else
 794         {
 795           /* Underflow to zero. */
 796           raw_frac = 0;
 797           raw_exp = 0;
 798         }
 799       break;
 800
 801     case INFINITE:
 802       raw_frac = 0;
 803       raw_exp = max_raw_exp;
 804       break;
 805
 806     case NAN:
 807       raw_frac = fp->fraction >> (64 - frac_bits);
 808       if (raw_frac == 0)
 809         raw_frac = 1;
 810       raw_exp = max_raw_exp;
 811       break;
 812
 813     case ZERO:
 814       raw_frac = 0;
 815       raw_exp = 0;
 816       break;
 817
 818     case MISSING:
 819       raw_sign = 1;
 820       raw_exp = max_raw_exp - 1;
 821       raw_frac = max_raw_frac;
 822       break;
 823
 824     case LOWEST:
 825       raw_sign = 1;
 826       raw_exp = max_raw_exp - 1;
 827       raw_frac = max_raw_frac - 1;
 828       break;
 829
 830     case HIGHEST:
 831       raw_sign = 0;
 832       raw_exp = max_raw_exp - 1;
 833       raw_frac = max_raw_frac;
 834       break;
 835
 836     case RESERVED:
 837       /* Convert to what processors commonly treat as signaling NAN. */
 838       raw_frac = (UINT64_C(1) << frac_bits) - 1;
 839       raw_exp = max_raw_exp;
 840       break;
 841
 842     default:
 843       NOT_REACHED ();
 844     }
 845
 846   return (((uint64_t) raw_sign << (frac_bits + exp_bits))
 847           | ((uint64_t) raw_exp << frac_bits)
 848           | raw_frac);
 849 }
 850
 851 /* Converts the neutral format floating point number in FP into
 852    VAX format with EXP_BITS exponent bits and FRAC_BITS fraction
 853    bits, and returns the value. */
 854 static uint64_t
 855 assemble_vax (struct fp *fp, int exp_bits, int frac_bits)
 856 {
 857   const int max_raw_exp = (1 << exp_bits) - 1;
 858   const int bias = 1 << (exp_bits - 1);
 859   const int min_finite_exp = 1 - bias;
 860   const int max_finite_exp = max_raw_exp - bias;
 861   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 862
 863   uint64_t raw_frac;
 864   int raw_exp;
 865   bool raw_sign;
 866
 867   raw_sign = fp->sign != POSITIVE;
 868
 869   switch (fp->class)
 870     {
 871     case FINITE:
 872       normalize_and_round_fp (fp, frac_bits + 1);
 873       if (fp->exponent > max_finite_exp)
 874         {
 875           /* Overflow to reserved operand. */
 876           raw_sign = 1;
 877           raw_exp = 0;
 878           raw_frac = 0;
 879         }
 880       else if (fp->exponent >= min_finite_exp)
 881         {
 882           /* Finite. */
 883           raw_frac = (fp->fraction << 1) >> (64 - frac_bits);
 884           raw_exp = fp->exponent + bias;
 885         }
 886       else
 887         {
 888           /* Underflow to zero. */
 889           raw_sign = 0;
 890           raw_frac = 0;
 891           raw_exp = 0;
 892         }
 893       break;
 894
 895     case INFINITE:
 896     case NAN:
 897     case RESERVED:
 898       raw_sign = 1;
 899       raw_exp = 0;
 900       raw_frac = 0;
 901       break;
 902
 903     case ZERO:
 904       raw_sign = 0;
 905       raw_frac = 0;
 906       raw_exp = 0;
 907       break;
 908
 909     case MISSING:
 910       raw_sign = 1;
 911       raw_exp = max_finite_exp + bias;
 912       raw_frac = max_raw_frac;
 913       break;
 914
 915     case LOWEST:
 916       raw_sign = 1;
 917       raw_exp = max_finite_exp + bias;
 918       raw_frac = max_raw_frac - 1;
 919       break;
 920
 921     case HIGHEST:
 922       raw_sign = 0;
 923       raw_exp = max_finite_exp + bias;
 924       raw_frac = max_raw_frac;
 925       break;
 926
 927     default:
 928       NOT_REACHED ();
 929     }
 930
 931   return (((uint64_t) raw_sign << (frac_bits + exp_bits))
 932           | ((uint64_t) raw_exp << frac_bits)
 933           | raw_frac);
 934 }
 935
 936 /* Shift right until the exponent is a multiple of 4.
 937    Rounding is not needed, because none of the formats we support
 938    has more than 53 bits of precision.  That is, we never shift a
 939    1-bit off the right end of the fraction.  */
 940 static void
 941 normalize_hex_fp (struct fp *fp)
 942 {
 943   while (fp->exponent % 4)
 944     {
 945       fp->fraction >>= 1;
 946       fp->exponent++;
 947     }
 948 }
 949
 950 /* Converts the neutral format floating point number in FP into Z
 951    architecture format with EXP_BITS exponent bits and FRAC_BITS
 952    fraction bits, and returns the value. */
 953 static uint64_t
 954 assemble_z (struct fp *fp, int exp_bits, int frac_bits)
 955 {
 956   const int max_raw_exp = (1 << exp_bits) - 1;
 957   const int bias = 1 << (exp_bits - 1);
 958   const int max_norm_exp = (max_raw_exp - bias) * 4;
 959   const int min_norm_exp = -bias * 4;
 960   const int min_denorm_exp = min_norm_exp - (frac_bits - 1);
 961
 962   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 963
 964   uint64_t raw_frac;
 965   int raw_exp;
 966   int raw_sign;
 967
 968   raw_sign = fp->sign != POSITIVE;
 969
 970   switch (fp->class)
 971     {
 972     case FINITE:
 973       normalize_and_round_fp (fp, frac_bits);
 974       normalize_hex_fp (fp);
 975       if (fp->exponent > max_norm_exp)
 976         {
 977           /* Overflow to largest value. */
 978           raw_exp = max_raw_exp;
 979           raw_frac = max_raw_frac;
 980         }
 981       else if (fp->exponent >= min_norm_exp)
 982         {
 983           /* Normal. */
 984           raw_frac = fp->fraction >> (64 - frac_bits);
 985           raw_exp = (fp->exponent / 4) + bias;
 986         }
 987       else if (fp->exponent >= min_denorm_exp)
 988         {
 989           /* Denormal. */
 990           const int denorm_shift = min_norm_exp - fp->exponent;
 991           raw_frac = (fp->fraction >> (64 - frac_bits)) >> denorm_shift;
 992           raw_exp = 0;
 993         }
 994       else
 995         {
 996           /* Underflow to zero. */
 997           raw_frac = 0;
 998           raw_exp = 0;
 999         }
1000       break;
1001
1002     case INFINITE:
1003       /* Overflow to largest value. */
1004       raw_exp = max_raw_exp;
1005       raw_frac = max_raw_frac;
1006       break;
1007
1008     case NAN:
1009     case RESERVED:
1010     case ZERO:
1011       /* Treat all of these as zero, because Z architecture
1012          doesn't have any reserved values. */
1013       raw_exp = 0;
1014       raw_frac = 0;
1015       break;
1016
1017     case MISSING:
1018       raw_sign = 1;
1019       raw_exp = max_raw_exp;
1020       raw_frac = max_raw_frac;
1021       break;
1022
1023     case LOWEST:
1024       raw_sign = 1;
1025       raw_exp = max_raw_exp;
1026       raw_frac = max_raw_frac - 1;
1027       break;
1028
1029     case HIGHEST:
1030       raw_sign = 0;
1031       raw_exp = max_raw_exp;
1032       raw_frac = max_raw_frac;
1033       break;
1034
1035     default:
1036       NOT_REACHED ();
1037     }
1038
1039   return (((uint64_t) raw_sign << (frac_bits + exp_bits))
1040           | ((uint64_t) raw_exp << frac_bits)
1041           | raw_frac);
1042 }
1043
1044 /* Converts the neutral format floating point number in FP into a
1045    null-terminated human-readable hex string in OUTPUT. */
1046 static void
1047 assemble_hex (struct fp *fp, void *output)
1048 {
1049   char buffer[32];
1050   char *s = buffer;
1051
1052   if (fp->sign == NEGATIVE)
1053     *s++ = '-';
1054
1055   switch (fp->class)
1056     {
1057     case FINITE:
1058       normalize_and_round_fp (fp, 64);
1059       normalize_hex_fp (fp);
1060       assert (fp->fraction != 0);
1061
1062       *s++ = '.';
1063       while (fp->fraction != 0)
1064         {
1065           *s++ = (fp->fraction >> 60)["0123456789abcdef"];
1066           fp->fraction <<= 4;
1067         }
1068       if (fp->exponent != 0)
1069         sprintf (s, "p%d", fp->exponent);
1070       break;
1071
1072     case INFINITE:
1073       strcpy (s, "Infinity");
1074       break;
1075
1076     case NAN:
1077       sprintf (s, "NaN:%016"PRIx64, fp->fraction);
1078       break;
1079
1080     case ZERO:
1081       strcpy (s, "0");
1082       break;
1083
1084     case MISSING:
1085       strcpy (buffer, "Missing");
1086       break;
1087
1088     case LOWEST:
1089       strcpy (buffer, "Lowest");
1090       break;
1091
1092     case HIGHEST:
1093       strcpy (buffer, "Highest");
1094       break;
1095
1096     case RESERVED:
1097       strcpy (s, "Reserved");
1098       break;
1099     }
1100
1101   strncpy (output, buffer, float_get_size (FLOAT_HEX));
1102 }