src/libpspp/float-format.c

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