src/libpspp/float-format.c

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