src/libpspp/float-format.c

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