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 candidate_cnt = sizeof candidates / sizeof *candidates;
 171
 172   enum float_format *p;
 173   int match_cnt;
 174
 175   match_cnt = 0;
 176   for (p = candidates; p < candidates + candidate_cnt; 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) && match_cnt++ == 0)
 183           *best_guess = *p;
 184       }
 185   return match_cnt;
 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 CNT bits in X starting from the given bit OFS. */
 207 static inline uint64_t
 208 get_bits (uint64_t x, int ofs, int cnt)
 209 {
 210   assert (ofs >= 0 && ofs < 64);
 211   assert (cnt > 0 && cnt < 64);
 212   assert (ofs + cnt <= 64);
 213   return (x >> ofs) & ((UINT64_C(1) << cnt) - 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             {
 484               char *tail;
 485               fp->exponent += strtol (s + 1, &tail, 10);
 486               s = tail;
 487             }
 488         }
 489     }
 490 }
 491 \f
 492 static uint64_t assemble_ieee (struct fp *, int exp_bits, int frac_bits);
 493 static uint64_t assemble_vax (struct fp *, int exp_bits, int frac_bits);
 494 static uint64_t assemble_z (struct fp *, int exp_bits, int frac_bits);
 495 static void assemble_hex (struct fp *, void *);
 496
 497 /* Converts the neutral format floating point number in FP into
 498    format TYPE at NUMBER.  May modify FP as part of the
 499    process. */
 500 static void
 501 assemble_number (enum float_format type, struct fp *fp, void *number)
 502 {
 503   switch (type)
 504     {
 505     case FLOAT_IEEE_SINGLE_LE:
 506       put_uint32 (native_to_le32 (assemble_ieee (fp, 8, 23)), number);
 507       break;
 508     case FLOAT_IEEE_SINGLE_BE:
 509       put_uint32 (native_to_be32 (assemble_ieee (fp, 8, 23)), number);
 510       break;
 511     case FLOAT_IEEE_DOUBLE_LE:
 512       put_uint64 (native_to_le64 (assemble_ieee (fp, 11, 52)), number);
 513       break;
 514     case FLOAT_IEEE_DOUBLE_BE:
 515       put_uint64 (native_to_be64 (assemble_ieee (fp, 11, 52)), number);
 516       break;
 517
 518     case FLOAT_VAX_F:
 519       put_uint32 (native_to_vax32 (assemble_vax (fp, 8, 23)), number);
 520       break;
 521     case FLOAT_VAX_D:
 522       put_uint64 (native_to_vax64 (assemble_vax (fp, 8, 55)), number);
 523       break;
 524     case FLOAT_VAX_G:
 525       put_uint64 (native_to_vax64 (assemble_vax (fp, 11, 52)), number);
 526       break;
 527
 528     case FLOAT_Z_SHORT:
 529       put_uint32 (native_to_be32 (assemble_z (fp, 7, 24)), number);
 530       break;
 531     case FLOAT_Z_LONG:
 532       put_uint64 (native_to_be64 (assemble_z (fp, 7, 56)), number);
 533       break;
 534
 535     case FLOAT_FP:
 536       memcpy (number, fp, sizeof *fp);
 537       break;
 538     case FLOAT_HEX:
 539       assemble_hex (fp, number);
 540       break;
 541     }
 542 }
 543
 544 /* Rounds off FP's fraction to FRAC_BITS bits of precision.
 545    Halfway values are rounded to even. */
 546 static void
 547 normalize_and_round_fp (struct fp *fp, int frac_bits)
 548 {
 549   assert (fp->class == FINITE);
 550   assert (fp->fraction != 0);
 551
 552   /* Make sure that the leading fraction bit is 1. */
 553   while (!(fp->fraction & (UINT64_C(1) << 63)))
 554     {
 555       fp->fraction <<= 1;
 556       fp->exponent--;
 557     }
 558
 559   if (frac_bits < 64)
 560     {
 561       uint64_t last_frac_bit = UINT64_C(1) << (64 - frac_bits);
 562       uint64_t decision_bit = last_frac_bit >> 1;
 563       if (fp->fraction & decision_bit
 564           && (fp->fraction & (decision_bit - 1)
 565               || fp->fraction & last_frac_bit))
 566         {
 567           fp->fraction += last_frac_bit;
 568           if ((fp->fraction >> 63) == 0)
 569             {
 570               fp->fraction = UINT64_C(1) << 63;
 571               fp->exponent++;
 572             }
 573         }
 574
 575       /* Mask off all but FRAC_BITS high-order bits.
 576          If we rounded up, then these bits no longer have
 577          meaningful values. */
 578       fp->fraction &= ~(last_frac_bit - 1);
 579     }
 580 }
 581
 582 /* Converts the neutral format floating point number in FP into
 583    IEEE format with EXP_BITS exponent bits and FRAC_BITS fraction
 584    bits, and returns the value. */
 585 static uint64_t
 586 assemble_ieee (struct fp *fp, int exp_bits, int frac_bits)
 587 {
 588   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 589
 590   const int bias = (1 << (exp_bits - 1)) - 1;
 591   const int max_raw_exp = (1 << exp_bits) - 1;
 592   const int min_norm_exp = 1 - bias;
 593   const int min_denorm_exp = min_norm_exp - frac_bits;
 594   const int max_norm_exp = max_raw_exp - 1 - bias;
 595
 596   uint64_t raw_frac;
 597   int raw_exp;
 598   bool raw_sign;
 599
 600   raw_sign = fp->sign != POSITIVE;
 601
 602   switch (fp->class)
 603     {
 604     case FINITE:
 605       normalize_and_round_fp (fp, frac_bits + 1);
 606       if (fp->exponent - 1 > max_norm_exp)
 607         {
 608           /* Overflow to infinity. */
 609           raw_exp = max_raw_exp;
 610           raw_frac = 0;
 611         }
 612       else if (fp->exponent - 1 >= min_norm_exp)
 613         {
 614           /* Normal. */
 615           raw_frac = (fp->fraction << 1) >> (64 - frac_bits);
 616           raw_exp = (fp->exponent - 1) + bias;
 617         }
 618       else if (fp->exponent - 1 >= min_denorm_exp)
 619         {
 620           /* Denormal. */
 621           const int denorm_shift = min_norm_exp - fp->exponent;
 622           raw_frac = (fp->fraction >> (64 - frac_bits)) >> denorm_shift;
 623           raw_exp = 0;
 624         }
 625       else
 626         {
 627           /* Underflow to zero. */
 628           raw_frac = 0;
 629           raw_exp = 0;
 630         }
 631       break;
 632
 633     case INFINITE:
 634       raw_frac = 0;
 635       raw_exp = max_raw_exp;
 636       break;
 637
 638     case NAN:
 639       raw_frac = fp->fraction >> (64 - frac_bits);
 640       if (raw_frac == 0)
 641         raw_frac = 1;
 642       raw_exp = max_raw_exp;
 643       break;
 644
 645     case ZERO:
 646       raw_frac = 0;
 647       raw_exp = 0;
 648       break;
 649
 650     case MISSING:
 651       raw_sign = 1;
 652       raw_exp = max_raw_exp - 1;
 653       raw_frac = max_raw_frac;
 654       break;
 655
 656     case LOWEST:
 657       raw_sign = 1;
 658       raw_exp = max_raw_exp - 1;
 659       raw_frac = max_raw_frac - 1;
 660       break;
 661
 662     case HIGHEST:
 663       raw_sign = 0;
 664       raw_exp = max_raw_exp - 1;
 665       raw_frac = max_raw_frac;
 666       break;
 667
 668     case RESERVED:
 669       /* Convert to what processors commonly treat as signaling NAN. */
 670       raw_frac = (UINT64_C(1) << frac_bits) - 1;
 671       raw_exp = max_raw_exp;
 672       break;
 673
 674     default:
 675       NOT_REACHED ();
 676     }
 677
 678   return (((uint64_t) raw_sign << (frac_bits + exp_bits))
 679           | ((uint64_t) raw_exp << frac_bits)
 680           | raw_frac);
 681 }
 682
 683 /* Converts the neutral format floating point number in FP into
 684    VAX format with EXP_BITS exponent bits and FRAC_BITS fraction
 685    bits, and returns the value. */
 686 static uint64_t
 687 assemble_vax (struct fp *fp, int exp_bits, int frac_bits)
 688 {
 689   const int max_raw_exp = (1 << exp_bits) - 1;
 690   const int bias = 1 << (exp_bits - 1);
 691   const int min_finite_exp = 1 - bias;
 692   const int max_finite_exp = max_raw_exp - bias;
 693   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 694
 695   uint64_t raw_frac;
 696   int raw_exp;
 697   bool raw_sign;
 698
 699   raw_sign = fp->sign != POSITIVE;
 700
 701   switch (fp->class)
 702     {
 703     case FINITE:
 704       normalize_and_round_fp (fp, frac_bits + 1);
 705       if (fp->exponent > max_finite_exp)
 706         {
 707           /* Overflow to reserved operand. */
 708           raw_sign = 1;
 709           raw_exp = 0;
 710           raw_frac = 0;
 711         }
 712       else if (fp->exponent >= min_finite_exp)
 713         {
 714           /* Finite. */
 715           raw_frac = (fp->fraction << 1) >> (64 - frac_bits);
 716           raw_exp = fp->exponent + bias;
 717         }
 718       else
 719         {
 720           /* Underflow to zero. */
 721           raw_sign = 0;
 722           raw_frac = 0;
 723           raw_exp = 0;
 724         }
 725       break;
 726
 727     case INFINITE:
 728     case NAN:
 729     case RESERVED:
 730       raw_sign = 1;
 731       raw_exp = 0;
 732       raw_frac = 0;
 733       break;
 734
 735     case ZERO:
 736       raw_sign = 0;
 737       raw_frac = 0;
 738       raw_exp = 0;
 739       break;
 740
 741     case MISSING:
 742       raw_sign = 1;
 743       raw_exp = max_finite_exp + bias;
 744       raw_frac = max_raw_frac;
 745       break;
 746
 747     case LOWEST:
 748       raw_sign = 1;
 749       raw_exp = max_finite_exp + bias;
 750       raw_frac = max_raw_frac - 1;
 751       break;
 752
 753     case HIGHEST:
 754       raw_sign = 0;
 755       raw_exp = max_finite_exp + bias;
 756       raw_frac = max_raw_frac;
 757       break;
 758
 759     default:
 760       NOT_REACHED ();
 761     }
 762
 763   return (((uint64_t) raw_sign << (frac_bits + exp_bits))
 764           | ((uint64_t) raw_exp << frac_bits)
 765           | raw_frac);
 766 }
 767
 768 /* Shift right until the exponent is a multiple of 4.
 769    Rounding is not needed, because none of the formats we support
 770    has more than 53 bits of precision.  That is, we never shift a
 771    1-bit off the right end of the fraction.  */
 772 static void
 773 normalize_hex_fp (struct fp *fp)
 774 {
 775   while (fp->exponent % 4)
 776     {
 777       fp->fraction >>= 1;
 778       fp->exponent++;
 779     }
 780 }
 781
 782 /* Converts the neutral format floating point number in FP into Z
 783    architecture format with EXP_BITS exponent bits and FRAC_BITS
 784    fraction bits, and returns the value. */
 785 static uint64_t
 786 assemble_z (struct fp *fp, int exp_bits, int frac_bits)
 787 {
 788   const int max_raw_exp = (1 << exp_bits) - 1;
 789   const int bias = 1 << (exp_bits - 1);
 790   const int max_norm_exp = (max_raw_exp - bias) * 4;
 791   const int min_norm_exp = -bias * 4;
 792   const int min_denorm_exp = min_norm_exp - (frac_bits - 1);
 793
 794   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 795
 796   uint64_t raw_frac;
 797   int raw_exp;
 798   int raw_sign;
 799
 800   raw_sign = fp->sign != POSITIVE;
 801
 802   switch (fp->class)
 803     {
 804     case FINITE:
 805       normalize_and_round_fp (fp, frac_bits);
 806       normalize_hex_fp (fp);
 807       if (fp->exponent > max_norm_exp)
 808         {
 809           /* Overflow to largest value. */
 810           raw_exp = max_raw_exp;
 811           raw_frac = max_raw_frac;
 812         }
 813       else if (fp->exponent >= min_norm_exp)
 814         {
 815           /* Normal. */
 816           raw_frac = fp->fraction >> (64 - frac_bits);
 817           raw_exp = (fp->exponent / 4) + bias;
 818         }
 819       else if (fp->exponent >= min_denorm_exp)
 820         {
 821           /* Denormal. */
 822           const int denorm_shift = min_norm_exp - fp->exponent;
 823           raw_frac = (fp->fraction >> (64 - frac_bits)) >> denorm_shift;
 824           raw_exp = 0;
 825         }
 826       else
 827         {
 828           /* Underflow to zero. */
 829           raw_frac = 0;
 830           raw_exp = 0;
 831         }
 832       break;
 833
 834     case INFINITE:
 835       /* Overflow to largest value. */
 836       raw_exp = max_raw_exp;
 837       raw_frac = max_raw_frac;
 838       break;
 839
 840     case NAN:
 841     case RESERVED:
 842     case ZERO:
 843       /* Treat all of these as zero, because Z architecture
 844          doesn't have any reserved values. */
 845       raw_exp = 0;
 846       raw_frac = 0;
 847       break;
 848
 849     case MISSING:
 850       raw_sign = 1;
 851       raw_exp = max_raw_exp;
 852       raw_frac = max_raw_frac;
 853       break;
 854
 855     case LOWEST:
 856       raw_sign = 1;
 857       raw_exp = max_raw_exp;
 858       raw_frac = max_raw_frac - 1;
 859       break;
 860
 861     case HIGHEST:
 862       raw_sign = 0;
 863       raw_exp = max_raw_exp;
 864       raw_frac = max_raw_frac;
 865       break;
 866
 867     default:
 868       NOT_REACHED ();
 869     }
 870
 871   return (((uint64_t) raw_sign << (frac_bits + exp_bits))
 872           | ((uint64_t) raw_exp << frac_bits)
 873           | raw_frac);
 874 }
 875
 876 /* Converts the neutral format floating point number in FP into a
 877    null-terminated human-readable hex string in OUTPUT. */
 878 static void
 879 assemble_hex (struct fp *fp, void *output)
 880 {
 881   char buffer[32];
 882   char *s = buffer;
 883
 884   if (fp->sign == NEGATIVE)
 885     *s++ = '-';
 886
 887   switch (fp->class)
 888     {
 889     case FINITE:
 890       normalize_and_round_fp (fp, 64);
 891       normalize_hex_fp (fp);
 892       assert (fp->fraction != 0);
 893
 894       *s++ = '.';
 895       while (fp->fraction != 0)
 896         {
 897           *s++ = (fp->fraction >> 60)["0123456789abcdef"];
 898           fp->fraction <<= 4;
 899         }
 900       if (fp->exponent != 0)
 901         sprintf (s, "p%d", fp->exponent);
 902       break;
 903
 904     case INFINITE:
 905       strcpy (s, "Infinity");
 906       break;
 907
 908     case NAN:
 909       sprintf (s, "NaN:%016"PRIx64, fp->fraction);
 910       break;
 911
 912     case ZERO:
 913       strcpy (s, "0");
 914       break;
 915
 916     case MISSING:
 917       strcpy (buffer, "Missing");
 918       break;
 919
 920     case LOWEST:
 921       strcpy (buffer, "Lowest");
 922       break;
 923
 924     case HIGHEST:
 925       strcpy (buffer, "Highest");
 926       break;
 927
 928     case RESERVED:
 929       strcpy (s, "Reserved");
 930       break;
 931     }
 932
 933   strncpy (output, buffer, float_get_size (FLOAT_HEX));
 934 }