1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2006, 2011 Free Software Foundation, Inc.
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.
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.
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/>. */
19 #include "libpspp/float-format.h"
27 #include "libpspp/assertion.h"
28 #include "libpspp/integer-format.h"
31 /* Neutral intermediate representation for binary floating-point numbers. */
36 FINITE, /* Finite number (normalized or denormalized). */
37 INFINITE, /* Positive or negative infinity. */
38 NAN, /* Not a number. */
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. */
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.
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.
63 Unused for other classes. */
68 static void extract_number (enum float_format, const void *, struct fp *);
69 static void assemble_number (enum float_format, struct fp *, void *);
71 static inline uint32_t get_uint32 (const void *);
72 static inline uint64_t get_uint64 (const void *);
74 static inline void put_uint32 (uint32_t, void *);
75 static inline void put_uint64 (uint64_t, void *);
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. */
81 float_convert (enum float_format from, const void *src,
82 enum float_format to, void *dst)
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);
95 extract_number (from, src, &fp);
96 assemble_number (to, &fp, dst);
102 memmove (dst, src, float_get_size (from));
106 /* Converts SRC from format FROM to a native double and returns
109 float_get_double (enum float_format from, const void *src)
112 float_convert (from, src, FLOAT_NATIVE_DOUBLE, &dst);
116 /* Returns the number of bytes in a number in the given
119 float_get_size (enum float_format format)
123 case FLOAT_IEEE_SINGLE_LE:
124 case FLOAT_IEEE_SINGLE_BE:
129 case FLOAT_IEEE_DOUBLE_LE:
130 case FLOAT_IEEE_DOUBLE_BE:
137 return sizeof (struct fp);
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. */
153 float_identify (double expected_value, const void *number, size_t length,
154 enum float_format *best_guess)
156 /* Candidates for identification in order of decreasing
158 enum float_format candidates[] =
160 FLOAT_IEEE_SINGLE_LE,
161 FLOAT_IEEE_SINGLE_BE,
162 FLOAT_IEEE_DOUBLE_LE,
163 FLOAT_IEEE_DOUBLE_BE,
170 const size_t candidate_cnt = sizeof candidates / sizeof *candidates;
172 enum float_format *p;
176 for (p = candidates; p < candidates + candidate_cnt; p++)
177 if (float_get_size (*p) == length)
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)
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.) */
195 float_get_lowest (void)
202 assemble_number (FLOAT_NATIVE_DOUBLE, &fp, &x);
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)
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);
216 /* Returns NATIVE converted to a form that, when stored in
217 memory, will be in VAX-endian byte order. */
218 static inline uint32_t
219 native_to_vax32 (uint32_t native)
221 return native_to_be32 (((native & 0xff00ff00) >> 8) |
222 ((native & 0x00ff00ff) << 8));
225 /* Returns NATIVE converted to a form that, when stored in
226 memory, will be in VAX-endian byte order. */
227 static inline uint64_t
228 native_to_vax64 (uint64_t native)
230 return native_to_be64 (((native & UINT64_C(0xff00ff0000000000)) >> 40) |
231 ((native & UINT64_C(0x00ff00ff00000000)) >> 24) |
232 ((native & UINT64_C(0x00000000ff00ff00)) << 24) |
233 ((native & UINT64_C(0x0000000000ff00ff)) << 40));
236 /* Given VAX, obtained from memory in VAX-endian format, returns
238 static inline uint32_t
239 vax_to_native32 (uint32_t vax)
241 uint32_t be = be_to_native32 (vax);
242 return ((be & 0xff00ff00) >> 8) | ((be & 0x00ff00ff) << 8);
245 /* Given VAX, obtained from memory in VAX-endian format, returns
247 static inline uint64_t
248 vax_to_native64 (uint64_t vax)
250 uint64_t be = be_to_native64 (vax);
251 return (((be & UINT64_C(0xff00ff0000000000)) >> 40) |
252 ((be & UINT64_C(0x00ff00ff00000000)) >> 24) |
253 ((be & UINT64_C(0x00000000ff00ff00)) << 24) |
254 ((be & UINT64_C(0x0000000000ff00ff)) << 40));
257 static void extract_ieee (uint64_t, int exp_bits, int frac_bits, struct fp *);
258 static void extract_vax (uint64_t, int exp_bits, int frac_bits, struct fp *);
259 static void extract_z (uint64_t, int exp_bits, int frac_bits, struct fp *);
260 static void extract_hex (const char *, struct fp *);
262 /* Converts the number at BITS from format TYPE into neutral
265 extract_number (enum float_format type, const void *bits, struct fp *fp)
269 case FLOAT_IEEE_SINGLE_LE:
270 extract_ieee (le_to_native32 (get_uint32 (bits)), 8, 23, fp);
272 case FLOAT_IEEE_SINGLE_BE:
273 extract_ieee (be_to_native32 (get_uint32 (bits)), 8, 23, fp);
275 case FLOAT_IEEE_DOUBLE_LE:
276 extract_ieee (le_to_native64 (get_uint64 (bits)), 11, 52, fp);
278 case FLOAT_IEEE_DOUBLE_BE:
279 extract_ieee (be_to_native64 (get_uint64 (bits)), 11, 52, fp);
283 extract_vax (vax_to_native32 (get_uint32 (bits)), 8, 23, fp);
286 extract_vax (vax_to_native64 (get_uint64 (bits)), 8, 55, fp);
289 extract_vax (vax_to_native64 (get_uint64 (bits)), 11, 52, fp);
293 extract_z (be_to_native32 (get_uint32 (bits)), 7, 24, fp);
296 extract_z (be_to_native64 (get_uint64 (bits)), 7, 56, fp);
300 memcpy (fp, bits, sizeof *fp);
303 extract_hex (bits, fp);
307 assert (!(fp->class == FINITE && fp->fraction == 0));
310 /* Converts the IEEE format number in BITS, which has EXP_BITS of
311 exponent and FRAC_BITS of fraction, into neutral format at
314 extract_ieee (uint64_t bits, int exp_bits, int frac_bits, struct fp *fp)
316 const int bias = (1 << (exp_bits - 1)) - 1;
317 const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
318 const int max_raw_exp = (1 << exp_bits) - 1;
320 const uint64_t raw_frac = get_bits (bits, 0, frac_bits);
321 const int raw_exp = get_bits (bits, frac_bits, exp_bits);
322 const bool raw_sign = get_bits (bits, frac_bits + exp_bits, 1);
324 if (raw_sign && raw_exp == max_raw_exp - 1 && raw_frac == max_raw_frac - 1)
326 else if (raw_exp == max_raw_exp - 1 && raw_frac == max_raw_frac)
327 fp->class = raw_sign ? MISSING : HIGHEST;
328 else if (raw_exp == max_raw_exp)
331 fp->class = INFINITE;
335 fp->fraction = raw_frac << (64 - frac_bits);
338 else if (raw_exp == 0)
343 fp->exponent = 1 - bias;
344 fp->fraction = raw_frac << (64 - frac_bits);
352 fp->exponent = raw_exp - bias + 1;
353 fp->fraction = (raw_frac << (64 - frac_bits - 1)) | (UINT64_C(1) << 63);
356 fp->sign = raw_sign ? NEGATIVE : POSITIVE;
359 /* Converts the VAX format number in BITS, which has EXP_BITS of
360 exponent and FRAC_BITS of fraction, into neutral format at
363 extract_vax (uint64_t bits, int exp_bits, int frac_bits, struct fp *fp)
365 const int bias = 1 << (exp_bits - 1);
366 const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
367 const int max_raw_exp = (1 << exp_bits) - 1;
369 const uint64_t raw_frac = get_bits (bits, 0, frac_bits);
370 const int raw_exp = get_bits (bits, frac_bits, exp_bits);
371 const bool raw_sign = get_bits (bits, frac_bits + exp_bits, 1);
373 if (raw_sign && raw_exp == max_raw_exp && raw_frac == max_raw_frac - 1)
375 else if (raw_exp == max_raw_exp && raw_frac == max_raw_frac)
376 fp->class = raw_sign ? MISSING : HIGHEST;
377 else if (raw_exp == 0)
378 fp->class = raw_sign == 0 ? ZERO : RESERVED;
382 fp->fraction = (raw_frac << (64 - frac_bits - 1)) | (UINT64_C(1) << 63);
383 fp->exponent = raw_exp - bias;
386 fp->sign = raw_sign ? NEGATIVE : POSITIVE;
389 /* Converts the Z architecture format number in BITS, which has
390 EXP_BITS of exponent and FRAC_BITS of fraction, into neutral
393 extract_z (uint64_t bits, int exp_bits, int frac_bits, struct fp *fp)
395 const int bias = 1 << (exp_bits - 1);
396 const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
397 const int max_raw_exp = (1 << exp_bits) - 1;
399 uint64_t raw_frac = get_bits (bits, 0, frac_bits);
400 int raw_exp = get_bits (bits, frac_bits, exp_bits);
401 int raw_sign = get_bits (bits, frac_bits + exp_bits, 1);
403 fp->sign = raw_sign ? NEGATIVE : POSITIVE;
404 if (raw_exp == max_raw_exp && raw_frac == max_raw_frac)
405 fp->class = raw_sign ? MISSING : HIGHEST;
406 else if (raw_sign && raw_exp == max_raw_exp && raw_frac == max_raw_frac - 1)
408 else if (raw_frac != 0)
411 fp->fraction = raw_frac << (64 - frac_bits);
412 fp->exponent = (raw_exp - bias) * 4;
418 /* Returns the integer value of hex digit C. */
422 const char s[] = "0123456789abcdef";
423 const char *cp = strchr (s, tolower ((unsigned char) c));
429 /* Parses a hexadecimal floating point number string at S (useful
430 for testing) into neutral format at FP. */
432 extract_hex (const char *s, struct fp *fp)
442 if (!strcmp (s, "Infinity"))
443 fp->class = INFINITE;
444 else if (!strcmp (s, "Missing"))
446 else if (!strcmp (s, "Lowest"))
448 else if (!strcmp (s, "Highest"))
450 else if (!strcmp (s, "Reserved"))
451 fp->class = RESERVED;
456 if (!memcmp (s, "NaN:", 4))
470 for (; isxdigit ((unsigned char) *s); s++)
473 uint64_t digit = hexit_value (*s);
474 fp->fraction += digit << offset;
478 if (fp->class == FINITE)
480 if (fp->fraction == 0)
485 fp->exponent += strtol (s + 1, &tail, 10);
492 static uint64_t assemble_ieee (struct fp *, int exp_bits, int frac_bits);
493 static uint64_t assemble_vax (struct fp *, int exp_bits, int frac_bits);
494 static uint64_t assemble_z (struct fp *, int exp_bits, int frac_bits);
495 static void assemble_hex (struct fp *, void *);
497 /* Converts the neutral format floating point number in FP into
498 format TYPE at NUMBER. May modify FP as part of the
501 assemble_number (enum float_format type, struct fp *fp, void *number)
505 case FLOAT_IEEE_SINGLE_LE:
506 put_uint32 (native_to_le32 (assemble_ieee (fp, 8, 23)), number);
508 case FLOAT_IEEE_SINGLE_BE:
509 put_uint32 (native_to_be32 (assemble_ieee (fp, 8, 23)), number);
511 case FLOAT_IEEE_DOUBLE_LE:
512 put_uint64 (native_to_le64 (assemble_ieee (fp, 11, 52)), number);
514 case FLOAT_IEEE_DOUBLE_BE:
515 put_uint64 (native_to_be64 (assemble_ieee (fp, 11, 52)), number);
519 put_uint32 (native_to_vax32 (assemble_vax (fp, 8, 23)), number);
522 put_uint64 (native_to_vax64 (assemble_vax (fp, 8, 55)), number);
525 put_uint64 (native_to_vax64 (assemble_vax (fp, 11, 52)), number);
529 put_uint32 (native_to_be32 (assemble_z (fp, 7, 24)), number);
532 put_uint64 (native_to_be64 (assemble_z (fp, 7, 56)), number);
536 memcpy (number, fp, sizeof *fp);
539 assemble_hex (fp, number);
544 /* Rounds off FP's fraction to FRAC_BITS bits of precision.
545 Halfway values are rounded to even. */
547 normalize_and_round_fp (struct fp *fp, int frac_bits)
549 assert (fp->class == FINITE);
550 assert (fp->fraction != 0);
552 /* Make sure that the leading fraction bit is 1. */
553 while (!(fp->fraction & (UINT64_C(1) << 63)))
561 uint64_t last_frac_bit = UINT64_C(1) << (64 - frac_bits);
562 uint64_t decision_bit = last_frac_bit >> 1;
563 if (fp->fraction & decision_bit
564 && (fp->fraction & (decision_bit - 1)
565 || fp->fraction & last_frac_bit))
567 fp->fraction += last_frac_bit;
568 if ((fp->fraction >> 63) == 0)
570 fp->fraction = UINT64_C(1) << 63;
575 /* Mask off all but FRAC_BITS high-order bits.
576 If we rounded up, then these bits no longer have
577 meaningful values. */
578 fp->fraction &= ~(last_frac_bit - 1);
582 /* Converts the neutral format floating point number in FP into
583 IEEE format with EXP_BITS exponent bits and FRAC_BITS fraction
584 bits, and returns the value. */
586 assemble_ieee (struct fp *fp, int exp_bits, int frac_bits)
588 const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
590 const int bias = (1 << (exp_bits - 1)) - 1;
591 const int max_raw_exp = (1 << exp_bits) - 1;
592 const int min_norm_exp = 1 - bias;
593 const int min_denorm_exp = min_norm_exp - frac_bits;
594 const int max_norm_exp = max_raw_exp - 1 - bias;
600 raw_sign = fp->sign != POSITIVE;
605 normalize_and_round_fp (fp, frac_bits + 1);
606 if (fp->exponent - 1 > max_norm_exp)
608 /* Overflow to infinity. */
609 raw_exp = max_raw_exp;
612 else if (fp->exponent - 1 >= min_norm_exp)
615 raw_frac = (fp->fraction << 1) >> (64 - frac_bits);
616 raw_exp = (fp->exponent - 1) + bias;
618 else if (fp->exponent - 1 >= min_denorm_exp)
621 const int denorm_shift = min_norm_exp - fp->exponent;
622 raw_frac = (fp->fraction >> (64 - frac_bits)) >> denorm_shift;
627 /* Underflow to zero. */
635 raw_exp = max_raw_exp;
639 raw_frac = fp->fraction >> (64 - frac_bits);
642 raw_exp = max_raw_exp;
652 raw_exp = max_raw_exp - 1;
653 raw_frac = max_raw_frac;
658 raw_exp = max_raw_exp - 1;
659 raw_frac = max_raw_frac - 1;
664 raw_exp = max_raw_exp - 1;
665 raw_frac = max_raw_frac;
669 /* Convert to what processors commonly treat as signaling NAN. */
670 raw_frac = (UINT64_C(1) << frac_bits) - 1;
671 raw_exp = max_raw_exp;
678 return (((uint64_t) raw_sign << (frac_bits + exp_bits))
679 | ((uint64_t) raw_exp << frac_bits)
683 /* Converts the neutral format floating point number in FP into
684 VAX format with EXP_BITS exponent bits and FRAC_BITS fraction
685 bits, and returns the value. */
687 assemble_vax (struct fp *fp, int exp_bits, int frac_bits)
689 const int max_raw_exp = (1 << exp_bits) - 1;
690 const int bias = 1 << (exp_bits - 1);
691 const int min_finite_exp = 1 - bias;
692 const int max_finite_exp = max_raw_exp - bias;
693 const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
699 raw_sign = fp->sign != POSITIVE;
704 normalize_and_round_fp (fp, frac_bits + 1);
705 if (fp->exponent > max_finite_exp)
707 /* Overflow to reserved operand. */
712 else if (fp->exponent >= min_finite_exp)
715 raw_frac = (fp->fraction << 1) >> (64 - frac_bits);
716 raw_exp = fp->exponent + bias;
720 /* Underflow to zero. */
743 raw_exp = max_finite_exp + bias;
744 raw_frac = max_raw_frac;
749 raw_exp = max_finite_exp + bias;
750 raw_frac = max_raw_frac - 1;
755 raw_exp = max_finite_exp + bias;
756 raw_frac = max_raw_frac;
763 return (((uint64_t) raw_sign << (frac_bits + exp_bits))
764 | ((uint64_t) raw_exp << frac_bits)
768 /* Shift right until the exponent is a multiple of 4.
769 Rounding is not needed, because none of the formats we support
770 has more than 53 bits of precision. That is, we never shift a
771 1-bit off the right end of the fraction. */
773 normalize_hex_fp (struct fp *fp)
775 while (fp->exponent % 4)
782 /* Converts the neutral format floating point number in FP into Z
783 architecture format with EXP_BITS exponent bits and FRAC_BITS
784 fraction bits, and returns the value. */
786 assemble_z (struct fp *fp, int exp_bits, int frac_bits)
788 const int max_raw_exp = (1 << exp_bits) - 1;
789 const int bias = 1 << (exp_bits - 1);
790 const int max_norm_exp = (max_raw_exp - bias) * 4;
791 const int min_norm_exp = -bias * 4;
792 const int min_denorm_exp = min_norm_exp - (frac_bits - 1);
794 const uint64_t max_raw_frac = (UINT64_C(1) << frac_bits) - 1;
800 raw_sign = fp->sign != POSITIVE;
805 normalize_and_round_fp (fp, frac_bits);
806 normalize_hex_fp (fp);
807 if (fp->exponent > max_norm_exp)
809 /* Overflow to largest value. */
810 raw_exp = max_raw_exp;
811 raw_frac = max_raw_frac;
813 else if (fp->exponent >= min_norm_exp)
816 raw_frac = fp->fraction >> (64 - frac_bits);
817 raw_exp = (fp->exponent / 4) + bias;
819 else if (fp->exponent >= min_denorm_exp)
822 const int denorm_shift = min_norm_exp - fp->exponent;
823 raw_frac = (fp->fraction >> (64 - frac_bits)) >> denorm_shift;
828 /* Underflow to zero. */
835 /* Overflow to largest value. */
836 raw_exp = max_raw_exp;
837 raw_frac = max_raw_frac;
843 /* Treat all of these as zero, because Z architecture
844 doesn't have any reserved values. */
851 raw_exp = max_raw_exp;
852 raw_frac = max_raw_frac;
857 raw_exp = max_raw_exp;
858 raw_frac = max_raw_frac - 1;
863 raw_exp = max_raw_exp;
864 raw_frac = max_raw_frac;
871 return (((uint64_t) raw_sign << (frac_bits + exp_bits))
872 | ((uint64_t) raw_exp << frac_bits)
876 /* Converts the neutral format floating point number in FP into a
877 null-terminated human-readable hex string in OUTPUT. */
879 assemble_hex (struct fp *fp, void *output)
884 if (fp->sign == NEGATIVE)
890 normalize_and_round_fp (fp, 64);
891 normalize_hex_fp (fp);
892 assert (fp->fraction != 0);
895 while (fp->fraction != 0)
897 *s++ = (fp->fraction >> 60)["0123456789abcdef"];
900 if (fp->exponent != 0)
901 sprintf (s, "p%d", fp->exponent);
905 strcpy (s, "Infinity");
909 sprintf (s, "NaN:%016"PRIx64, fp->fraction);
917 strcpy (buffer, "Missing");
921 strcpy (buffer, "Lowest");
925 strcpy (buffer, "Highest");
929 strcpy (s, "Reserved");
933 strncpy (output, buffer, float_get_size (FLOAT_HEX));