src/libpspp/float-format.c

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