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