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