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 the 32-bit unsigned integer at P,
 217    which need not be aligned. */
 218 static inline uint32_t
 219 get_uint32 (const void *p)
 220 {
 221   uint32_t x;
 222   memcpy (&x, p, sizeof x);
 223   return x;
 224 }
 225
 226 /* Returns the 64-bit unsigned integer at P,
 227    which need not be aligned. */
 228 static inline uint64_t
 229 get_uint64 (const void *p)
 230 {
 231   uint64_t x;
 232   memcpy (&x, p, sizeof x);
 233   return x;
 234 }
 235
 236 /* Stores 32-bit unsigned integer X at P,
 237    which need not be aligned. */
 238 static inline void
 239 put_uint32 (uint32_t x, void *p)
 240 {
 241   memcpy (p, &x, sizeof x);
 242 }
 243
 244 /* Stores 64-bit unsigned integer X at P,
 245    which need not be aligned. */
 246 static inline void
 247 put_uint64 (uint64_t x, void *p)
 248 {
 249   memcpy (p, &x, sizeof x);
 250 }
 251
 252 /* Returns NATIVE converted to a form that, when stored in
 253    memory, will be in little-endian byte order. */
 254 static inline uint32_t
 255 native_to_le32 (uint32_t native)
 256 {
 257   return INTEGER_NATIVE == INTEGER_LSB_FIRST ? native : bswap_32 (native);
 258 }
 259
 260 /* Returns NATIVE converted to a form that, when stored in
 261    memory, will be in big-endian byte order. */
 262 static inline uint32_t
 263 native_to_be32 (uint32_t native)
 264 {
 265   return INTEGER_NATIVE == INTEGER_MSB_FIRST ? native : bswap_32 (native);
 266 }
 267
 268 /* Returns NATIVE converted to a form that, when stored in
 269    memory, will be in VAX-endian byte order. */
 270 static inline uint32_t
 271 native_to_vax32 (uint32_t native)
 272 {
 273   return native_to_be32 (((native & 0xff00ff00) >> 8) |
 274                          ((native & 0x00ff00ff) << 8));
 275 }
 276
 277 /* Returns NATIVE converted to a form that, when stored in
 278    memory, will be in little-endian byte order. */
 279 static inline uint64_t
 280 native_to_le64 (uint64_t native)
 281 {
 282   return INTEGER_NATIVE == INTEGER_LSB_FIRST ? native : bswap_64 (native);
 283 }
 284
 285 /* Returns NATIVE converted to a form that, when stored in
 286    memory, will be in big-endian byte order. */
 287 static inline uint64_t
 288 native_to_be64 (uint64_t native)
 289 {
 290   return INTEGER_NATIVE == INTEGER_MSB_FIRST ? native : bswap_64 (native);
 291 }
 292
 293 /* Returns NATIVE converted to a form that, when stored in
 294    memory, will be in VAX-endian byte order. */
 295 static inline uint64_t
 296 native_to_vax64 (uint64_t native)
 297 {
 298   return native_to_be64 (((native & UINT64_C(0xff00ff0000000000)) >> 40) |
 299                          ((native & UINT64_C(0x00ff00ff00000000)) >> 24) |
 300                          ((native & UINT64_C(0x00000000ff00ff00)) << 24) |
 301                          ((native & UINT64_C(0x0000000000ff00ff)) << 40));
 302 }
 303
 304 /* Given LE, obtained from memory in little-endian format,
 305    returns its value. */
 306 static inline uint32_t
 307 le_to_native32 (uint32_t le)
 308 {
 309   return INTEGER_NATIVE == INTEGER_LSB_FIRST ? le : bswap_32 (le);
 310 }
 311
 312 /* Given BE, obtained from memory in big-endian format, returns
 313    its value. */
 314 static inline uint32_t
 315 be_to_native32 (uint32_t be)
 316 {
 317   return INTEGER_NATIVE == INTEGER_MSB_FIRST ? be : bswap_32 (be);
 318 }
 319
 320 /* Given VAX, obtained from memory in VAX-endian format, returns
 321    its value. */
 322 static inline uint32_t
 323 vax_to_native32 (uint32_t vax)
 324 {
 325   uint32_t be = be_to_native32 (vax);
 326   return ((be & 0xff00ff00) >> 8) | ((be & 0x00ff00ff) << 8);
 327 }
 328
 329 /* Given LE, obtained from memory in little-endian format,
 330    returns its value. */
 331 static inline uint64_t
 332 le_to_native64 (uint64_t le)
 333 {
 334   return INTEGER_NATIVE == INTEGER_LSB_FIRST ? le : bswap_64 (le);
 335 }
 336
 337 /* Given BE, obtained from memory in big-endian format, returns
 338    its value. */
 339 static inline uint64_t
 340 be_to_native64 (uint64_t be)
 341 {
 342   return INTEGER_NATIVE == INTEGER_MSB_FIRST ? be : bswap_64 (be);
 343 }
 344
 345 /* Given VAX, obtained from memory in VAX-endian format, returns
 346    its value. */
 347 static inline uint64_t
 348 vax_to_native64 (uint64_t vax)
 349 {
 350   uint64_t be = be_to_native64 (vax);
 351   return (((be & UINT64_C(0xff00ff0000000000)) >> 40) |
 352           ((be & UINT64_C(0x00ff00ff00000000)) >> 24) |
 353           ((be & UINT64_C(0x00000000ff00ff00)) << 24) |
 354           ((be & UINT64_C(0x0000000000ff00ff)) << 40));
 355 }
 356 \f
 357 static void extract_ieee (uint64_t, int exp_bits, int frac_bits, struct fp *);
 358 static void extract_vax (uint64_t, int exp_bits, int frac_bits, struct fp *);
 359 static void extract_z (uint64_t, int exp_bits, int frac_bits, struct fp *);
 360 static void extract_hex (const char *, struct fp *);
 361
 362 /* Converts the number at BITS from format TYPE into neutral
 363    format at FP. */
 364 static void
 365 extract_number (enum float_format type, const void *bits, struct fp *fp)
 366 {
 367   switch (type)
 368     {
 369     case FLOAT_IEEE_SINGLE_LE:
 370       extract_ieee (le_to_native32 (get_uint32 (bits)), 8, 23, fp);
 371       break;
 372     case FLOAT_IEEE_SINGLE_BE:
 373       extract_ieee (be_to_native32 (get_uint32 (bits)), 8, 23, fp);
 374       break;
 375     case FLOAT_IEEE_DOUBLE_LE:
 376       extract_ieee (le_to_native64 (get_uint64 (bits)), 11, 52, fp);
 377       break;
 378     case FLOAT_IEEE_DOUBLE_BE:
 379       extract_ieee (be_to_native64 (get_uint64 (bits)), 11, 52, fp);
 380       break;
 381
 382     case FLOAT_VAX_F:
 383       extract_vax (vax_to_native32 (get_uint32 (bits)), 8, 23, fp);
 384       break;
 385     case FLOAT_VAX_D:
 386       extract_vax (vax_to_native64 (get_uint64 (bits)), 8, 55, fp);
 387       break;
 388     case FLOAT_VAX_G:
 389       extract_vax (vax_to_native64 (get_uint64 (bits)), 11, 52, fp);
 390       break;
 391
 392     case FLOAT_Z_SHORT:
 393       extract_z (be_to_native32 (get_uint32 (bits)), 7, 24, fp);
 394       break;
 395     case FLOAT_Z_LONG:
 396       extract_z (be_to_native64 (get_uint64 (bits)), 7, 56, fp);
 397       break;
 398
 399     case FLOAT_FP:
 400       memcpy (fp, bits, sizeof *fp);
 401       break;
 402     case FLOAT_HEX:
 403       extract_hex (bits, fp);
 404       break;
 405     }
 406
 407   assert (!(fp->class == FINITE && fp->fraction == 0));
 408 }
 409
 410 /* Converts the IEEE format number in BITS, which has EXP_BITS of
 411    exponent and FRAC_BITS of fraction, into neutral format at
 412    FP. */
 413 static void
 414 extract_ieee (uint64_t bits, int exp_bits, int frac_bits, struct fp *fp)
 415 {
 416   const int bias = (1 << (exp_bits - 1)) - 1;
 417   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 418   const int max_raw_exp = (1 << exp_bits) - 1;
 419
 420   const uint64_t raw_frac = get_bits (bits, 0, frac_bits);
 421   const int raw_exp = get_bits (bits, frac_bits, exp_bits);
 422   const bool raw_sign = get_bits (bits, frac_bits + exp_bits, 1);
 423
 424   if (raw_sign && raw_exp == max_raw_exp - 1 && raw_frac == max_raw_frac - 1)
 425     fp->class = LOWEST;
 426   else if (raw_exp == max_raw_exp - 1 && raw_frac == max_raw_frac)
 427     fp->class = raw_sign ? MISSING : HIGHEST;
 428   else if (raw_exp == max_raw_exp)
 429     {
 430       if (raw_frac == 0)
 431         fp->class = INFINITE;
 432       else
 433         {
 434           fp->class = NAN;
 435           fp->fraction = raw_frac << (64 - frac_bits);
 436         }
 437     }
 438   else if (raw_exp == 0)
 439     {
 440       if (raw_frac != 0)
 441         {
 442           fp->class = FINITE;
 443           fp->exponent = 1 - bias;
 444           fp->fraction = raw_frac << (64 - frac_bits);
 445         }
 446       else
 447         fp->class = ZERO;
 448     }
 449   else
 450     {
 451       fp->class = FINITE;
 452       fp->exponent = raw_exp - bias + 1;
 453       fp->fraction = (raw_frac << (64 - frac_bits - 1)) | (UINT64_C(1) << 63);
 454     }
 455
 456   fp->sign = raw_sign ? NEGATIVE : POSITIVE;
 457 }
 458
 459 /* Converts the VAX format number in BITS, which has EXP_BITS of
 460    exponent and FRAC_BITS of fraction, into neutral format at
 461    FP. */
 462 static void
 463 extract_vax (uint64_t bits, int exp_bits, int frac_bits, struct fp *fp)
 464 {
 465   const int bias = 1 << (exp_bits - 1);
 466   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 467   const int max_raw_exp = (1 << exp_bits) - 1;
 468
 469   const uint64_t raw_frac = get_bits (bits, 0, frac_bits);
 470   const int raw_exp = get_bits (bits, frac_bits, exp_bits);
 471   const bool raw_sign = get_bits (bits, frac_bits + exp_bits, 1);
 472
 473   if (raw_sign && raw_exp == max_raw_exp && raw_frac == max_raw_frac - 1)
 474     fp->class = LOWEST;
 475   else if (raw_exp == max_raw_exp && raw_frac == max_raw_frac)
 476     fp->class = raw_sign ? MISSING : HIGHEST;
 477   else if (raw_exp == 0)
 478     fp->class = raw_sign == 0 ? ZERO : RESERVED;
 479   else
 480     {
 481       fp->class = FINITE;
 482       fp->fraction = (raw_frac << (64 - frac_bits - 1)) | (UINT64_C(1) << 63);
 483       fp->exponent = raw_exp - bias;
 484     }
 485
 486   fp->sign = raw_sign ? NEGATIVE : POSITIVE;
 487 }
 488
 489 /* Converts the Z architecture format number in BITS, which has
 490    EXP_BITS of exponent and FRAC_BITS of fraction, into neutral
 491    format at FP. */
 492 static void
 493 extract_z (uint64_t bits, int exp_bits, int frac_bits, struct fp *fp)
 494 {
 495   const int bias = 1 << (exp_bits - 1);
 496   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 497   const int max_raw_exp = (1 << exp_bits) - 1;
 498
 499   uint64_t raw_frac = get_bits (bits, 0, frac_bits);
 500   int raw_exp = get_bits (bits, frac_bits, exp_bits);
 501   int raw_sign = get_bits (bits, frac_bits + exp_bits, 1);
 502
 503   fp->sign = raw_sign ? NEGATIVE : POSITIVE;
 504   if (raw_exp == max_raw_exp && raw_frac == max_raw_frac)
 505     fp->class = raw_sign ? MISSING : HIGHEST;
 506   else if (raw_sign && raw_exp == max_raw_exp && raw_frac == max_raw_frac - 1)
 507     fp->class = LOWEST;
 508   else if (raw_frac != 0)
 509     {
 510       fp->class = FINITE;
 511       fp->fraction = raw_frac << (64 - frac_bits);
 512       fp->exponent = (raw_exp - bias) * 4;
 513     }
 514   else
 515     fp->class = ZERO;
 516 }
 517
 518 /* Returns the integer value of hex digit C. */
 519 static int
 520 hexit_value (int c)
 521 {
 522   const char s[] = "0123456789abcdef";
 523   const char *cp = strchr (s, tolower ((unsigned char) c));
 524
 525   assert (cp != NULL);
 526   return cp - s;
 527 }
 528
 529 /* Parses a hexadecimal floating point number string at S (useful
 530    for testing) into neutral format at FP. */
 531 static void
 532 extract_hex (const char *s, struct fp *fp)
 533 {
 534   if (*s == '-')
 535     {
 536       fp->sign = NEGATIVE;
 537       s++;
 538     }
 539   else
 540     fp->sign = POSITIVE;
 541
 542   if (!strcmp (s, "Infinity"))
 543     fp->class = INFINITE;
 544   else if (!strcmp (s, "Missing"))
 545     fp->class = MISSING;
 546   else if (!strcmp (s, "Lowest"))
 547     fp->class = LOWEST;
 548   else if (!strcmp (s, "Highest"))
 549     fp->class = HIGHEST;
 550   else if (!strcmp (s, "Reserved"))
 551     fp->class = RESERVED;
 552   else
 553     {
 554       int offset;
 555
 556       if (!memcmp (s, "NaN:", 4))
 557         {
 558           fp->class = NAN;
 559           s += 4;
 560         }
 561       else
 562         fp->class = FINITE;
 563
 564       if (*s == '.')
 565         s++;
 566
 567       fp->exponent = 0;
 568       fp->fraction = 0;
 569       offset = 60;
 570       for (; isxdigit ((unsigned char) *s); s++)
 571         if (offset >= 0)
 572           {
 573             uint64_t digit = hexit_value (*s);
 574             fp->fraction += digit << offset;
 575             offset -= 4;
 576           }
 577
 578       if (fp->class == FINITE)
 579         {
 580           if (fp->fraction == 0)
 581             fp->class = ZERO;
 582           else if (*s == 'p')
 583             {
 584               char *tail;
 585               fp->exponent += strtol (s + 1, &tail, 10);
 586               s = tail;
 587             }
 588         }
 589     }
 590 }
 591 \f
 592 static uint64_t assemble_ieee (struct fp *, int exp_bits, int frac_bits);
 593 static uint64_t assemble_vax (struct fp *, int exp_bits, int frac_bits);
 594 static uint64_t assemble_z (struct fp *, int exp_bits, int frac_bits);
 595 static void assemble_hex (struct fp *, void *);
 596
 597 /* Converts the neutral format floating point number in FP into
 598    format TYPE at NUMBER.  May modify FP as part of the
 599    process. */
 600 static void
 601 assemble_number (enum float_format type, struct fp *fp, void *number)
 602 {
 603   switch (type)
 604     {
 605     case FLOAT_IEEE_SINGLE_LE:
 606       put_uint32 (native_to_le32 (assemble_ieee (fp, 8, 23)), number);
 607       break;
 608     case FLOAT_IEEE_SINGLE_BE:
 609       put_uint32 (native_to_be32 (assemble_ieee (fp, 8, 23)), number);
 610       break;
 611     case FLOAT_IEEE_DOUBLE_LE:
 612       put_uint64 (native_to_le64 (assemble_ieee (fp, 11, 52)), number);
 613       break;
 614     case FLOAT_IEEE_DOUBLE_BE:
 615       put_uint64 (native_to_be64 (assemble_ieee (fp, 11, 52)), number);
 616       break;
 617
 618     case FLOAT_VAX_F:
 619       put_uint32 (native_to_vax32 (assemble_vax (fp, 8, 23)), number);
 620       break;
 621     case FLOAT_VAX_D:
 622       put_uint64 (native_to_vax64 (assemble_vax (fp, 8, 55)), number);
 623       break;
 624     case FLOAT_VAX_G:
 625       put_uint64 (native_to_vax64 (assemble_vax (fp, 11, 52)), number);
 626       break;
 627
 628     case FLOAT_Z_SHORT:
 629       put_uint32 (native_to_be32 (assemble_z (fp, 7, 24)), number);
 630       break;
 631     case FLOAT_Z_LONG:
 632       put_uint64 (native_to_be64 (assemble_z (fp, 7, 56)), number);
 633       break;
 634
 635     case FLOAT_FP:
 636       memcpy (number, fp, sizeof *fp);
 637       break;
 638     case FLOAT_HEX:
 639       assemble_hex (fp, number);
 640       break;
 641     }
 642 }
 643
 644 /* Rounds off FP's fraction to FRAC_BITS bits of precision.
 645    Halfway values are rounded to even. */
 646 static void
 647 normalize_and_round_fp (struct fp *fp, int frac_bits)
 648 {
 649   assert (fp->class == FINITE);
 650   assert (fp->fraction != 0);
 651
 652   /* Make sure that the leading fraction bit is 1. */
 653   while (!(fp->fraction & (UINT64_C(1) << 63)))
 654     {
 655       fp->fraction <<= 1;
 656       fp->exponent--;
 657     }
 658
 659   if (frac_bits < 64)
 660     {
 661       uint64_t last_frac_bit = UINT64_C(1) << (64 - frac_bits);
 662       uint64_t decision_bit = last_frac_bit >> 1;
 663       if (fp->fraction & decision_bit
 664           && (fp->fraction & (decision_bit - 1)
 665               || fp->fraction & last_frac_bit))
 666         {
 667           fp->fraction += last_frac_bit;
 668           if ((fp->fraction >> 63) == 0)
 669             {
 670               fp->fraction = UINT64_C(1) << 63;
 671               fp->exponent++;
 672             }
 673         }
 674
 675       /* Mask off all but FRAC_BITS high-order bits.
 676          If we rounded up, then these bits no longer have
 677          meaningful values. */
 678       fp->fraction &= ~(last_frac_bit - 1);
 679     }
 680 }
 681
 682 /* Converts the neutral format floating point number in FP into
 683    IEEE format with EXP_BITS exponent bits and FRAC_BITS fraction
 684    bits, and returns the value. */
 685 static uint64_t
 686 assemble_ieee (struct fp *fp, int exp_bits, int frac_bits)
 687 {
 688   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 689
 690   const int bias = (1 << (exp_bits - 1)) - 1;
 691   const int max_raw_exp = (1 << exp_bits) - 1;
 692   const int min_norm_exp = 1 - bias;
 693   const int min_denorm_exp = min_norm_exp - frac_bits;
 694   const int max_norm_exp = max_raw_exp - 1 - bias;
 695
 696   uint64_t raw_frac;
 697   int raw_exp;
 698   bool raw_sign;
 699
 700   raw_sign = fp->sign != POSITIVE;
 701
 702   switch (fp->class)
 703     {
 704     case FINITE:
 705       normalize_and_round_fp (fp, frac_bits + 1);
 706       if (fp->exponent - 1 > max_norm_exp)
 707         {
 708           /* Overflow to infinity. */
 709           raw_exp = max_raw_exp;
 710           raw_frac = 0;
 711         }
 712       else if (fp->exponent - 1 >= min_norm_exp)
 713         {
 714           /* Normal. */
 715           raw_frac = (fp->fraction << 1) >> (64 - frac_bits);
 716           raw_exp = (fp->exponent - 1) + bias;
 717         }
 718       else if (fp->exponent - 1 >= min_denorm_exp)
 719         {
 720           /* Denormal. */
 721           const int denorm_shift = min_norm_exp - fp->exponent;
 722           raw_frac = (fp->fraction >> (64 - frac_bits)) >> denorm_shift;
 723           raw_exp = 0;
 724         }
 725       else
 726         {
 727           /* Underflow to zero. */
 728           raw_frac = 0;
 729           raw_exp = 0;
 730         }
 731       break;
 732
 733     case INFINITE:
 734       raw_frac = 0;
 735       raw_exp = max_raw_exp;
 736       break;
 737
 738     case NAN:
 739       raw_frac = fp->fraction >> (64 - frac_bits);
 740       if (raw_frac == 0)
 741         raw_frac = 1;
 742       raw_exp = max_raw_exp;
 743       break;
 744
 745     case ZERO:
 746       raw_frac = 0;
 747       raw_exp = 0;
 748       break;
 749
 750     case MISSING:
 751       raw_sign = 1;
 752       raw_exp = max_raw_exp - 1;
 753       raw_frac = max_raw_frac;
 754       break;
 755
 756     case LOWEST:
 757       raw_sign = 1;
 758       raw_exp = max_raw_exp - 1;
 759       raw_frac = max_raw_frac - 1;
 760       break;
 761
 762     case HIGHEST:
 763       raw_sign = 0;
 764       raw_exp = max_raw_exp - 1;
 765       raw_frac = max_raw_frac;
 766       break;
 767
 768     case RESERVED:
 769       /* Convert to what processors commonly treat as signaling NAN. */
 770       raw_frac = (UINT64_C(1) << frac_bits) - 1;
 771       raw_exp = max_raw_exp;
 772       break;
 773
 774     default:
 775       NOT_REACHED ();
 776     }
 777
 778   return (((uint64_t) raw_sign << (frac_bits + exp_bits))
 779           | ((uint64_t) raw_exp << frac_bits)
 780           | raw_frac);
 781 }
 782
 783 /* Converts the neutral format floating point number in FP into
 784    VAX format with EXP_BITS exponent bits and FRAC_BITS fraction
 785    bits, and returns the value. */
 786 static uint64_t
 787 assemble_vax (struct fp *fp, int exp_bits, int frac_bits)
 788 {
 789   const int max_raw_exp = (1 << exp_bits) - 1;
 790   const int bias = 1 << (exp_bits - 1);
 791   const int min_finite_exp = 1 - bias;
 792   const int max_finite_exp = max_raw_exp - bias;
 793   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 794
 795   uint64_t raw_frac;
 796   int raw_exp;
 797   bool raw_sign;
 798
 799   raw_sign = fp->sign != POSITIVE;
 800
 801   switch (fp->class)
 802     {
 803     case FINITE:
 804       normalize_and_round_fp (fp, frac_bits + 1);
 805       if (fp->exponent > max_finite_exp)
 806         {
 807           /* Overflow to reserved operand. */
 808           raw_sign = 1;
 809           raw_exp = 0;
 810           raw_frac = 0;
 811         }
 812       else if (fp->exponent >= min_finite_exp)
 813         {
 814           /* Finite. */
 815           raw_frac = (fp->fraction << 1) >> (64 - frac_bits);
 816           raw_exp = fp->exponent + bias;
 817         }
 818       else
 819         {
 820           /* Underflow to zero. */
 821           raw_sign = 0;
 822           raw_frac = 0;
 823           raw_exp = 0;
 824         }
 825       break;
 826
 827     case INFINITE:
 828     case NAN:
 829     case RESERVED:
 830       raw_sign = 1;
 831       raw_exp = 0;
 832       raw_frac = 0;
 833       break;
 834
 835     case ZERO:
 836       raw_sign = 0;
 837       raw_frac = 0;
 838       raw_exp = 0;
 839       break;
 840
 841     case MISSING:
 842       raw_sign = 1;
 843       raw_exp = max_finite_exp + bias;
 844       raw_frac = max_raw_frac;
 845       break;
 846
 847     case LOWEST:
 848       raw_sign = 1;
 849       raw_exp = max_finite_exp + bias;
 850       raw_frac = max_raw_frac - 1;
 851       break;
 852
 853     case HIGHEST:
 854       raw_sign = 0;
 855       raw_exp = max_finite_exp + bias;
 856       raw_frac = max_raw_frac;
 857       break;
 858
 859     default:
 860       NOT_REACHED ();
 861     }
 862
 863   return (((uint64_t) raw_sign << (frac_bits + exp_bits))
 864           | ((uint64_t) raw_exp << frac_bits)
 865           | raw_frac);
 866 }
 867
 868 /* Shift right until the exponent is a multiple of 4.
 869    Rounding is not needed, because none of the formats we support
 870    has more than 53 bits of precision.  That is, we never shift a
 871    1-bit off the right end of the fraction.  */
 872 static void
 873 normalize_hex_fp (struct fp *fp)
 874 {
 875   while (fp->exponent % 4)
 876     {
 877       fp->fraction >>= 1;
 878       fp->exponent++;
 879     }
 880 }
 881
 882 /* Converts the neutral format floating point number in FP into Z
 883    architecture format with EXP_BITS exponent bits and FRAC_BITS
 884    fraction bits, and returns the value. */
 885 static uint64_t
 886 assemble_z (struct fp *fp, int exp_bits, int frac_bits)
 887 {
 888   const int max_raw_exp = (1 << exp_bits) - 1;
 889   const int bias = 1 << (exp_bits - 1);
 890   const int max_norm_exp = (max_raw_exp - bias) * 4;
 891   const int min_norm_exp = -bias * 4;
 892   const int min_denorm_exp = min_norm_exp - (frac_bits - 1);
 893
 894   const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
 895
 896   uint64_t raw_frac;
 897   int raw_exp;
 898   int raw_sign;
 899
 900   raw_sign = fp->sign != POSITIVE;
 901
 902   switch (fp->class)
 903     {
 904     case FINITE:
 905       normalize_and_round_fp (fp, frac_bits);
 906       normalize_hex_fp (fp);
 907       if (fp->exponent > max_norm_exp)
 908         {
 909           /* Overflow to largest value. */
 910           raw_exp = max_raw_exp;
 911           raw_frac = max_raw_frac;
 912         }
 913       else if (fp->exponent >= min_norm_exp)
 914         {
 915           /* Normal. */
 916           raw_frac = fp->fraction >> (64 - frac_bits);
 917           raw_exp = (fp->exponent / 4) + bias;
 918         }
 919       else if (fp->exponent >= min_denorm_exp)
 920         {
 921           /* Denormal. */
 922           const int denorm_shift = min_norm_exp - fp->exponent;
 923           raw_frac = (fp->fraction >> (64 - frac_bits)) >> denorm_shift;
 924           raw_exp = 0;
 925         }
 926       else
 927         {
 928           /* Underflow to zero. */
 929           raw_frac = 0;
 930           raw_exp = 0;
 931         }
 932       break;
 933
 934     case INFINITE:
 935       /* Overflow to largest value. */
 936       raw_exp = max_raw_exp;
 937       raw_frac = max_raw_frac;
 938       break;
 939
 940     case NAN:
 941     case RESERVED:
 942     case ZERO:
 943       /* Treat all of these as zero, because Z architecture
 944          doesn't have any reserved values. */
 945       raw_exp = 0;
 946       raw_frac = 0;
 947       break;
 948
 949     case MISSING:
 950       raw_sign = 1;
 951       raw_exp = max_raw_exp;
 952       raw_frac = max_raw_frac;
 953       break;
 954
 955     case LOWEST:
 956       raw_sign = 1;
 957       raw_exp = max_raw_exp;
 958       raw_frac = max_raw_frac - 1;
 959       break;
 960
 961     case HIGHEST:
 962       raw_sign = 0;
 963       raw_exp = max_raw_exp;
 964       raw_frac = max_raw_frac;
 965       break;
 966
 967     default:
 968       NOT_REACHED ();
 969     }
 970
 971   return (((uint64_t) raw_sign << (frac_bits + exp_bits))
 972           | ((uint64_t) raw_exp << frac_bits)
 973           | raw_frac);
 974 }
 975
 976 /* Converts the neutral format floating point number in FP into a
 977    null-terminated human-readable hex string in OUTPUT. */
 978 static void
 979 assemble_hex (struct fp *fp, void *output)
 980 {
 981   char buffer[32];
 982   char *s = buffer;
 983
 984   if (fp->sign == NEGATIVE)
 985     *s++ = '-';
 986
 987   switch (fp->class)
 988     {
 989     case FINITE:
 990       normalize_and_round_fp (fp, 64);
 991       normalize_hex_fp (fp);
 992       assert (fp->fraction != 0);
 993
 994       *s++ = '.';
 995       while (fp->fraction != 0)
 996         {
 997           *s++ = (fp->fraction >> 60)["0123456789abcdef"];
 998           fp->fraction <<= 4;
 999         }
1000       if (fp->exponent != 0)
1001         sprintf (s, "p%d", fp->exponent);
1002       break;
1003
1004     case INFINITE:
1005       strcpy (s, "Infinity");
1006       break;
1007
1008     case NAN:
1009       sprintf (s, "NaN:%016"PRIx64, fp->fraction);
1010       break;
1011
1012     case ZERO:
1013       strcpy (s, "0");
1014       break;
1015
1016     case MISSING:
1017       strcpy (buffer, "Missing");
1018       break;
1019
1020     case LOWEST:
1021       strcpy (buffer, "Lowest");
1022       break;
1023
1024     case HIGHEST:
1025       strcpy (buffer, "Highest");
1026       break;
1027
1028     case RESERVED:
1029       strcpy (s, "Reserved");
1030       break;
1031     }
1032
1033   strncpy (output, buffer, float_get_size (FLOAT_HEX));
1034 }