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

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