Add comment explaining the meaning of encoding to data_out
[pspp-builds.git] / 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 #include <string.h>
25
26 #include <libpspp/assertion.h>
27 #include <libpspp/integer-format.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
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 }