xref: /openbmc/qemu/fpu/softfloat.c (revision 4921d0a7)
1 /*
2  * QEMU float support
3  *
4  * The code in this source file is derived from release 2a of the SoftFloat
5  * IEC/IEEE Floating-point Arithmetic Package. Those parts of the code (and
6  * some later contributions) are provided under that license, as detailed below.
7  * It has subsequently been modified by contributors to the QEMU Project,
8  * so some portions are provided under:
9  *  the SoftFloat-2a license
10  *  the BSD license
11  *  GPL-v2-or-later
12  *
13  * Any future contributions to this file after December 1st 2014 will be
14  * taken to be licensed under the Softfloat-2a license unless specifically
15  * indicated otherwise.
16  */
17 
18 /*
19 ===============================================================================
20 This C source file is part of the SoftFloat IEC/IEEE Floating-point
21 Arithmetic Package, Release 2a.
22 
23 Written by John R. Hauser.  This work was made possible in part by the
24 International Computer Science Institute, located at Suite 600, 1947 Center
25 Street, Berkeley, California 94704.  Funding was partially provided by the
26 National Science Foundation under grant MIP-9311980.  The original version
27 of this code was written as part of a project to build a fixed-point vector
28 processor in collaboration with the University of California at Berkeley,
29 overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
30 is available through the Web page `http://HTTP.CS.Berkeley.EDU/~jhauser/
31 arithmetic/SoftFloat.html'.
32 
33 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort
34 has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
35 TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO
36 PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
37 AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
38 
39 Derivative works are acceptable, even for commercial purposes, so long as
40 (1) they include prominent notice that the work is derivative, and (2) they
41 include prominent notice akin to these four paragraphs for those parts of
42 this code that are retained.
43 
44 ===============================================================================
45 */
46 
47 /* BSD licensing:
48  * Copyright (c) 2006, Fabrice Bellard
49  * All rights reserved.
50  *
51  * Redistribution and use in source and binary forms, with or without
52  * modification, are permitted provided that the following conditions are met:
53  *
54  * 1. Redistributions of source code must retain the above copyright notice,
55  * this list of conditions and the following disclaimer.
56  *
57  * 2. Redistributions in binary form must reproduce the above copyright notice,
58  * this list of conditions and the following disclaimer in the documentation
59  * and/or other materials provided with the distribution.
60  *
61  * 3. Neither the name of the copyright holder nor the names of its contributors
62  * may be used to endorse or promote products derived from this software without
63  * specific prior written permission.
64  *
65  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
66  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
67  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
68  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
69  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
70  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
71  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
72  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
73  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
74  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
75  * THE POSSIBILITY OF SUCH DAMAGE.
76  */
77 
78 /* Portions of this work are licensed under the terms of the GNU GPL,
79  * version 2 or later. See the COPYING file in the top-level directory.
80  */
81 
82 /* softfloat (and in particular the code in softfloat-specialize.h) is
83  * target-dependent and needs the TARGET_* macros.
84  */
85 #include "qemu/osdep.h"
86 #include <math.h>
87 #include "qemu/bitops.h"
88 #include "fpu/softfloat.h"
89 
90 /* We only need stdlib for abort() */
91 
92 /*----------------------------------------------------------------------------
93 | Primitive arithmetic functions, including multi-word arithmetic, and
94 | division and square root approximations.  (Can be specialized to target if
95 | desired.)
96 *----------------------------------------------------------------------------*/
97 #include "fpu/softfloat-macros.h"
98 
99 /*
100  * Hardfloat
101  *
102  * Fast emulation of guest FP instructions is challenging for two reasons.
103  * First, FP instruction semantics are similar but not identical, particularly
104  * when handling NaNs. Second, emulating at reasonable speed the guest FP
105  * exception flags is not trivial: reading the host's flags register with a
106  * feclearexcept & fetestexcept pair is slow [slightly slower than soft-fp],
107  * and trapping on every FP exception is not fast nor pleasant to work with.
108  *
109  * We address these challenges by leveraging the host FPU for a subset of the
110  * operations. To do this we expand on the idea presented in this paper:
111  *
112  * Guo, Yu-Chuan, et al. "Translating the ARM Neon and VFP instructions in a
113  * binary translator." Software: Practice and Experience 46.12 (2016):1591-1615.
114  *
115  * The idea is thus to leverage the host FPU to (1) compute FP operations
116  * and (2) identify whether FP exceptions occurred while avoiding
117  * expensive exception flag register accesses.
118  *
119  * An important optimization shown in the paper is that given that exception
120  * flags are rarely cleared by the guest, we can avoid recomputing some flags.
121  * This is particularly useful for the inexact flag, which is very frequently
122  * raised in floating-point workloads.
123  *
124  * We optimize the code further by deferring to soft-fp whenever FP exception
125  * detection might get hairy. Two examples: (1) when at least one operand is
126  * denormal/inf/NaN; (2) when operands are not guaranteed to lead to a 0 result
127  * and the result is < the minimum normal.
128  */
129 #define GEN_INPUT_FLUSH__NOCHECK(name, soft_t)                          \
130     static inline void name(soft_t *a, float_status *s)                 \
131     {                                                                   \
132         if (unlikely(soft_t ## _is_denormal(*a))) {                     \
133             *a = soft_t ## _set_sign(soft_t ## _zero,                   \
134                                      soft_t ## _is_neg(*a));            \
135             float_raise(float_flag_input_denormal, s);                  \
136         }                                                               \
137     }
138 
139 GEN_INPUT_FLUSH__NOCHECK(float32_input_flush__nocheck, float32)
140 GEN_INPUT_FLUSH__NOCHECK(float64_input_flush__nocheck, float64)
141 #undef GEN_INPUT_FLUSH__NOCHECK
142 
143 #define GEN_INPUT_FLUSH1(name, soft_t)                  \
144     static inline void name(soft_t *a, float_status *s) \
145     {                                                   \
146         if (likely(!s->flush_inputs_to_zero)) {         \
147             return;                                     \
148         }                                               \
149         soft_t ## _input_flush__nocheck(a, s);          \
150     }
151 
152 GEN_INPUT_FLUSH1(float32_input_flush1, float32)
153 GEN_INPUT_FLUSH1(float64_input_flush1, float64)
154 #undef GEN_INPUT_FLUSH1
155 
156 #define GEN_INPUT_FLUSH2(name, soft_t)                                  \
157     static inline void name(soft_t *a, soft_t *b, float_status *s)      \
158     {                                                                   \
159         if (likely(!s->flush_inputs_to_zero)) {                         \
160             return;                                                     \
161         }                                                               \
162         soft_t ## _input_flush__nocheck(a, s);                          \
163         soft_t ## _input_flush__nocheck(b, s);                          \
164     }
165 
166 GEN_INPUT_FLUSH2(float32_input_flush2, float32)
167 GEN_INPUT_FLUSH2(float64_input_flush2, float64)
168 #undef GEN_INPUT_FLUSH2
169 
170 #define GEN_INPUT_FLUSH3(name, soft_t)                                  \
171     static inline void name(soft_t *a, soft_t *b, soft_t *c, float_status *s) \
172     {                                                                   \
173         if (likely(!s->flush_inputs_to_zero)) {                         \
174             return;                                                     \
175         }                                                               \
176         soft_t ## _input_flush__nocheck(a, s);                          \
177         soft_t ## _input_flush__nocheck(b, s);                          \
178         soft_t ## _input_flush__nocheck(c, s);                          \
179     }
180 
181 GEN_INPUT_FLUSH3(float32_input_flush3, float32)
182 GEN_INPUT_FLUSH3(float64_input_flush3, float64)
183 #undef GEN_INPUT_FLUSH3
184 
185 /*
186  * Choose whether to use fpclassify or float32/64_* primitives in the generated
187  * hardfloat functions. Each combination of number of inputs and float size
188  * gets its own value.
189  */
190 #if defined(__x86_64__)
191 # define QEMU_HARDFLOAT_1F32_USE_FP 0
192 # define QEMU_HARDFLOAT_1F64_USE_FP 1
193 # define QEMU_HARDFLOAT_2F32_USE_FP 0
194 # define QEMU_HARDFLOAT_2F64_USE_FP 1
195 # define QEMU_HARDFLOAT_3F32_USE_FP 0
196 # define QEMU_HARDFLOAT_3F64_USE_FP 1
197 #else
198 # define QEMU_HARDFLOAT_1F32_USE_FP 0
199 # define QEMU_HARDFLOAT_1F64_USE_FP 0
200 # define QEMU_HARDFLOAT_2F32_USE_FP 0
201 # define QEMU_HARDFLOAT_2F64_USE_FP 0
202 # define QEMU_HARDFLOAT_3F32_USE_FP 0
203 # define QEMU_HARDFLOAT_3F64_USE_FP 0
204 #endif
205 
206 /*
207  * QEMU_HARDFLOAT_USE_ISINF chooses whether to use isinf() over
208  * float{32,64}_is_infinity when !USE_FP.
209  * On x86_64/aarch64, using the former over the latter can yield a ~6% speedup.
210  * On power64 however, using isinf() reduces fp-bench performance by up to 50%.
211  */
212 #if defined(__x86_64__) || defined(__aarch64__)
213 # define QEMU_HARDFLOAT_USE_ISINF   1
214 #else
215 # define QEMU_HARDFLOAT_USE_ISINF   0
216 #endif
217 
218 /*
219  * Some targets clear the FP flags before most FP operations. This prevents
220  * the use of hardfloat, since hardfloat relies on the inexact flag being
221  * already set.
222  */
223 #if defined(TARGET_PPC) || defined(__FAST_MATH__)
224 # if defined(__FAST_MATH__)
225 #  warning disabling hardfloat due to -ffast-math: hardfloat requires an exact \
226     IEEE implementation
227 # endif
228 # define QEMU_NO_HARDFLOAT 1
229 # define QEMU_SOFTFLOAT_ATTR QEMU_FLATTEN
230 #else
231 # define QEMU_NO_HARDFLOAT 0
232 # define QEMU_SOFTFLOAT_ATTR QEMU_FLATTEN __attribute__((noinline))
233 #endif
234 
235 static inline bool can_use_fpu(const float_status *s)
236 {
237     if (QEMU_NO_HARDFLOAT) {
238         return false;
239     }
240     return likely(s->float_exception_flags & float_flag_inexact &&
241                   s->float_rounding_mode == float_round_nearest_even);
242 }
243 
244 /*
245  * Hardfloat generation functions. Each operation can have two flavors:
246  * either using softfloat primitives (e.g. float32_is_zero_or_normal) for
247  * most condition checks, or native ones (e.g. fpclassify).
248  *
249  * The flavor is chosen by the callers. Instead of using macros, we rely on the
250  * compiler to propagate constants and inline everything into the callers.
251  *
252  * We only generate functions for operations with two inputs, since only
253  * these are common enough to justify consolidating them into common code.
254  */
255 
256 typedef union {
257     float32 s;
258     float h;
259 } union_float32;
260 
261 typedef union {
262     float64 s;
263     double h;
264 } union_float64;
265 
266 typedef bool (*f32_check_fn)(union_float32 a, union_float32 b);
267 typedef bool (*f64_check_fn)(union_float64 a, union_float64 b);
268 
269 typedef float32 (*soft_f32_op2_fn)(float32 a, float32 b, float_status *s);
270 typedef float64 (*soft_f64_op2_fn)(float64 a, float64 b, float_status *s);
271 typedef float   (*hard_f32_op2_fn)(float a, float b);
272 typedef double  (*hard_f64_op2_fn)(double a, double b);
273 
274 /* 2-input is-zero-or-normal */
275 static inline bool f32_is_zon2(union_float32 a, union_float32 b)
276 {
277     if (QEMU_HARDFLOAT_2F32_USE_FP) {
278         /*
279          * Not using a temp variable for consecutive fpclassify calls ends up
280          * generating faster code.
281          */
282         return (fpclassify(a.h) == FP_NORMAL || fpclassify(a.h) == FP_ZERO) &&
283                (fpclassify(b.h) == FP_NORMAL || fpclassify(b.h) == FP_ZERO);
284     }
285     return float32_is_zero_or_normal(a.s) &&
286            float32_is_zero_or_normal(b.s);
287 }
288 
289 static inline bool f64_is_zon2(union_float64 a, union_float64 b)
290 {
291     if (QEMU_HARDFLOAT_2F64_USE_FP) {
292         return (fpclassify(a.h) == FP_NORMAL || fpclassify(a.h) == FP_ZERO) &&
293                (fpclassify(b.h) == FP_NORMAL || fpclassify(b.h) == FP_ZERO);
294     }
295     return float64_is_zero_or_normal(a.s) &&
296            float64_is_zero_or_normal(b.s);
297 }
298 
299 /* 3-input is-zero-or-normal */
300 static inline
301 bool f32_is_zon3(union_float32 a, union_float32 b, union_float32 c)
302 {
303     if (QEMU_HARDFLOAT_3F32_USE_FP) {
304         return (fpclassify(a.h) == FP_NORMAL || fpclassify(a.h) == FP_ZERO) &&
305                (fpclassify(b.h) == FP_NORMAL || fpclassify(b.h) == FP_ZERO) &&
306                (fpclassify(c.h) == FP_NORMAL || fpclassify(c.h) == FP_ZERO);
307     }
308     return float32_is_zero_or_normal(a.s) &&
309            float32_is_zero_or_normal(b.s) &&
310            float32_is_zero_or_normal(c.s);
311 }
312 
313 static inline
314 bool f64_is_zon3(union_float64 a, union_float64 b, union_float64 c)
315 {
316     if (QEMU_HARDFLOAT_3F64_USE_FP) {
317         return (fpclassify(a.h) == FP_NORMAL || fpclassify(a.h) == FP_ZERO) &&
318                (fpclassify(b.h) == FP_NORMAL || fpclassify(b.h) == FP_ZERO) &&
319                (fpclassify(c.h) == FP_NORMAL || fpclassify(c.h) == FP_ZERO);
320     }
321     return float64_is_zero_or_normal(a.s) &&
322            float64_is_zero_or_normal(b.s) &&
323            float64_is_zero_or_normal(c.s);
324 }
325 
326 static inline bool f32_is_inf(union_float32 a)
327 {
328     if (QEMU_HARDFLOAT_USE_ISINF) {
329         return isinf(a.h);
330     }
331     return float32_is_infinity(a.s);
332 }
333 
334 static inline bool f64_is_inf(union_float64 a)
335 {
336     if (QEMU_HARDFLOAT_USE_ISINF) {
337         return isinf(a.h);
338     }
339     return float64_is_infinity(a.s);
340 }
341 
342 static inline float32
343 float32_gen2(float32 xa, float32 xb, float_status *s,
344              hard_f32_op2_fn hard, soft_f32_op2_fn soft,
345              f32_check_fn pre, f32_check_fn post)
346 {
347     union_float32 ua, ub, ur;
348 
349     ua.s = xa;
350     ub.s = xb;
351 
352     if (unlikely(!can_use_fpu(s))) {
353         goto soft;
354     }
355 
356     float32_input_flush2(&ua.s, &ub.s, s);
357     if (unlikely(!pre(ua, ub))) {
358         goto soft;
359     }
360 
361     ur.h = hard(ua.h, ub.h);
362     if (unlikely(f32_is_inf(ur))) {
363         float_raise(float_flag_overflow, s);
364     } else if (unlikely(fabsf(ur.h) <= FLT_MIN) && post(ua, ub)) {
365         goto soft;
366     }
367     return ur.s;
368 
369  soft:
370     return soft(ua.s, ub.s, s);
371 }
372 
373 static inline float64
374 float64_gen2(float64 xa, float64 xb, float_status *s,
375              hard_f64_op2_fn hard, soft_f64_op2_fn soft,
376              f64_check_fn pre, f64_check_fn post)
377 {
378     union_float64 ua, ub, ur;
379 
380     ua.s = xa;
381     ub.s = xb;
382 
383     if (unlikely(!can_use_fpu(s))) {
384         goto soft;
385     }
386 
387     float64_input_flush2(&ua.s, &ub.s, s);
388     if (unlikely(!pre(ua, ub))) {
389         goto soft;
390     }
391 
392     ur.h = hard(ua.h, ub.h);
393     if (unlikely(f64_is_inf(ur))) {
394         float_raise(float_flag_overflow, s);
395     } else if (unlikely(fabs(ur.h) <= DBL_MIN) && post(ua, ub)) {
396         goto soft;
397     }
398     return ur.s;
399 
400  soft:
401     return soft(ua.s, ub.s, s);
402 }
403 
404 /*
405  * Classify a floating point number. Everything above float_class_qnan
406  * is a NaN so cls >= float_class_qnan is any NaN.
407  */
408 
409 typedef enum __attribute__ ((__packed__)) {
410     float_class_unclassified,
411     float_class_zero,
412     float_class_normal,
413     float_class_inf,
414     float_class_qnan,  /* all NaNs from here */
415     float_class_snan,
416 } FloatClass;
417 
418 #define float_cmask(bit)  (1u << (bit))
419 
420 enum {
421     float_cmask_zero    = float_cmask(float_class_zero),
422     float_cmask_normal  = float_cmask(float_class_normal),
423     float_cmask_inf     = float_cmask(float_class_inf),
424     float_cmask_qnan    = float_cmask(float_class_qnan),
425     float_cmask_snan    = float_cmask(float_class_snan),
426 
427     float_cmask_infzero = float_cmask_zero | float_cmask_inf,
428     float_cmask_anynan  = float_cmask_qnan | float_cmask_snan,
429 };
430 
431 /* Flags for parts_minmax. */
432 enum {
433     /* Set for minimum; clear for maximum. */
434     minmax_ismin = 1,
435     /* Set for the IEEE 754-2008 minNum() and maxNum() operations. */
436     minmax_isnum = 2,
437     /* Set for the IEEE 754-2008 minNumMag() and minNumMag() operations. */
438     minmax_ismag = 4,
439     /*
440      * Set for the IEEE 754-2019 minimumNumber() and maximumNumber()
441      * operations.
442      */
443     minmax_isnumber = 8,
444 };
445 
446 /* Simple helpers for checking if, or what kind of, NaN we have */
447 static inline __attribute__((unused)) bool is_nan(FloatClass c)
448 {
449     return unlikely(c >= float_class_qnan);
450 }
451 
452 static inline __attribute__((unused)) bool is_snan(FloatClass c)
453 {
454     return c == float_class_snan;
455 }
456 
457 static inline __attribute__((unused)) bool is_qnan(FloatClass c)
458 {
459     return c == float_class_qnan;
460 }
461 
462 /*
463  * Structure holding all of the decomposed parts of a float.
464  * The exponent is unbiased and the fraction is normalized.
465  *
466  * The fraction words are stored in big-endian word ordering,
467  * so that truncation from a larger format to a smaller format
468  * can be done simply by ignoring subsequent elements.
469  */
470 
471 typedef struct {
472     FloatClass cls;
473     bool sign;
474     int32_t exp;
475     union {
476         /* Routines that know the structure may reference the singular name. */
477         uint64_t frac;
478         /*
479          * Routines expanded with multiple structures reference "hi" and "lo"
480          * depending on the operation.  In FloatParts64, "hi" and "lo" are
481          * both the same word and aliased here.
482          */
483         uint64_t frac_hi;
484         uint64_t frac_lo;
485     };
486 } FloatParts64;
487 
488 typedef struct {
489     FloatClass cls;
490     bool sign;
491     int32_t exp;
492     uint64_t frac_hi;
493     uint64_t frac_lo;
494 } FloatParts128;
495 
496 typedef struct {
497     FloatClass cls;
498     bool sign;
499     int32_t exp;
500     uint64_t frac_hi;
501     uint64_t frac_hm;  /* high-middle */
502     uint64_t frac_lm;  /* low-middle */
503     uint64_t frac_lo;
504 } FloatParts256;
505 
506 /* These apply to the most significant word of each FloatPartsN. */
507 #define DECOMPOSED_BINARY_POINT    63
508 #define DECOMPOSED_IMPLICIT_BIT    (1ull << DECOMPOSED_BINARY_POINT)
509 
510 /* Structure holding all of the relevant parameters for a format.
511  *   exp_size: the size of the exponent field
512  *   exp_bias: the offset applied to the exponent field
513  *   exp_max: the maximum normalised exponent
514  *   frac_size: the size of the fraction field
515  *   frac_shift: shift to normalise the fraction with DECOMPOSED_BINARY_POINT
516  * The following are computed based the size of fraction
517  *   round_mask: bits below lsb which must be rounded
518  * The following optional modifiers are available:
519  *   arm_althp: handle ARM Alternative Half Precision
520  *   m68k_denormal: explicit integer bit for extended precision may be 1
521  */
522 typedef struct {
523     int exp_size;
524     int exp_bias;
525     int exp_re_bias;
526     int exp_max;
527     int frac_size;
528     int frac_shift;
529     bool arm_althp;
530     bool m68k_denormal;
531     uint64_t round_mask;
532 } FloatFmt;
533 
534 /* Expand fields based on the size of exponent and fraction */
535 #define FLOAT_PARAMS_(E)                                \
536     .exp_size       = E,                                \
537     .exp_bias       = ((1 << E) - 1) >> 1,              \
538     .exp_re_bias    = (1 << (E - 1)) + (1 << (E - 2)),  \
539     .exp_max        = (1 << E) - 1
540 
541 #define FLOAT_PARAMS(E, F)                              \
542     FLOAT_PARAMS_(E),                                   \
543     .frac_size      = F,                                \
544     .frac_shift     = (-F - 1) & 63,                    \
545     .round_mask     = (1ull << ((-F - 1) & 63)) - 1
546 
547 static const FloatFmt float16_params = {
548     FLOAT_PARAMS(5, 10)
549 };
550 
551 static const FloatFmt float16_params_ahp = {
552     FLOAT_PARAMS(5, 10),
553     .arm_althp = true
554 };
555 
556 static const FloatFmt bfloat16_params = {
557     FLOAT_PARAMS(8, 7)
558 };
559 
560 static const FloatFmt float32_params = {
561     FLOAT_PARAMS(8, 23)
562 };
563 
564 static const FloatFmt float64_params = {
565     FLOAT_PARAMS(11, 52)
566 };
567 
568 static const FloatFmt float128_params = {
569     FLOAT_PARAMS(15, 112)
570 };
571 
572 #define FLOATX80_PARAMS(R)              \
573     FLOAT_PARAMS_(15),                  \
574     .frac_size = R == 64 ? 63 : R,      \
575     .frac_shift = 0,                    \
576     .round_mask = R == 64 ? -1 : (1ull << ((-R - 1) & 63)) - 1
577 
578 static const FloatFmt floatx80_params[3] = {
579     [floatx80_precision_s] = { FLOATX80_PARAMS(23) },
580     [floatx80_precision_d] = { FLOATX80_PARAMS(52) },
581     [floatx80_precision_x] = {
582         FLOATX80_PARAMS(64),
583 #ifdef TARGET_M68K
584         .m68k_denormal = true,
585 #endif
586     },
587 };
588 
589 /* Unpack a float to parts, but do not canonicalize.  */
590 static void unpack_raw64(FloatParts64 *r, const FloatFmt *fmt, uint64_t raw)
591 {
592     const int f_size = fmt->frac_size;
593     const int e_size = fmt->exp_size;
594 
595     *r = (FloatParts64) {
596         .cls = float_class_unclassified,
597         .sign = extract64(raw, f_size + e_size, 1),
598         .exp = extract64(raw, f_size, e_size),
599         .frac = extract64(raw, 0, f_size)
600     };
601 }
602 
603 static void QEMU_FLATTEN float16_unpack_raw(FloatParts64 *p, float16 f)
604 {
605     unpack_raw64(p, &float16_params, f);
606 }
607 
608 static void QEMU_FLATTEN bfloat16_unpack_raw(FloatParts64 *p, bfloat16 f)
609 {
610     unpack_raw64(p, &bfloat16_params, f);
611 }
612 
613 static void QEMU_FLATTEN float32_unpack_raw(FloatParts64 *p, float32 f)
614 {
615     unpack_raw64(p, &float32_params, f);
616 }
617 
618 static void QEMU_FLATTEN float64_unpack_raw(FloatParts64 *p, float64 f)
619 {
620     unpack_raw64(p, &float64_params, f);
621 }
622 
623 static void QEMU_FLATTEN floatx80_unpack_raw(FloatParts128 *p, floatx80 f)
624 {
625     *p = (FloatParts128) {
626         .cls = float_class_unclassified,
627         .sign = extract32(f.high, 15, 1),
628         .exp = extract32(f.high, 0, 15),
629         .frac_hi = f.low
630     };
631 }
632 
633 static void QEMU_FLATTEN float128_unpack_raw(FloatParts128 *p, float128 f)
634 {
635     const int f_size = float128_params.frac_size - 64;
636     const int e_size = float128_params.exp_size;
637 
638     *p = (FloatParts128) {
639         .cls = float_class_unclassified,
640         .sign = extract64(f.high, f_size + e_size, 1),
641         .exp = extract64(f.high, f_size, e_size),
642         .frac_hi = extract64(f.high, 0, f_size),
643         .frac_lo = f.low,
644     };
645 }
646 
647 /* Pack a float from parts, but do not canonicalize.  */
648 static uint64_t pack_raw64(const FloatParts64 *p, const FloatFmt *fmt)
649 {
650     const int f_size = fmt->frac_size;
651     const int e_size = fmt->exp_size;
652     uint64_t ret;
653 
654     ret = (uint64_t)p->sign << (f_size + e_size);
655     ret = deposit64(ret, f_size, e_size, p->exp);
656     ret = deposit64(ret, 0, f_size, p->frac);
657     return ret;
658 }
659 
660 static float16 QEMU_FLATTEN float16_pack_raw(const FloatParts64 *p)
661 {
662     return make_float16(pack_raw64(p, &float16_params));
663 }
664 
665 static bfloat16 QEMU_FLATTEN bfloat16_pack_raw(const FloatParts64 *p)
666 {
667     return pack_raw64(p, &bfloat16_params);
668 }
669 
670 static float32 QEMU_FLATTEN float32_pack_raw(const FloatParts64 *p)
671 {
672     return make_float32(pack_raw64(p, &float32_params));
673 }
674 
675 static float64 QEMU_FLATTEN float64_pack_raw(const FloatParts64 *p)
676 {
677     return make_float64(pack_raw64(p, &float64_params));
678 }
679 
680 static float128 QEMU_FLATTEN float128_pack_raw(const FloatParts128 *p)
681 {
682     const int f_size = float128_params.frac_size - 64;
683     const int e_size = float128_params.exp_size;
684     uint64_t hi;
685 
686     hi = (uint64_t)p->sign << (f_size + e_size);
687     hi = deposit64(hi, f_size, e_size, p->exp);
688     hi = deposit64(hi, 0, f_size, p->frac_hi);
689     return make_float128(hi, p->frac_lo);
690 }
691 
692 /*----------------------------------------------------------------------------
693 | Functions and definitions to determine:  (1) whether tininess for underflow
694 | is detected before or after rounding by default, (2) what (if anything)
695 | happens when exceptions are raised, (3) how signaling NaNs are distinguished
696 | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
697 | are propagated from function inputs to output.  These details are target-
698 | specific.
699 *----------------------------------------------------------------------------*/
700 #include "softfloat-specialize.c.inc"
701 
702 #define PARTS_GENERIC_64_128(NAME, P) \
703     _Generic((P), FloatParts64 *: parts64_##NAME, \
704                   FloatParts128 *: parts128_##NAME)
705 
706 #define PARTS_GENERIC_64_128_256(NAME, P) \
707     _Generic((P), FloatParts64 *: parts64_##NAME, \
708                   FloatParts128 *: parts128_##NAME, \
709                   FloatParts256 *: parts256_##NAME)
710 
711 #define parts_default_nan(P, S)    PARTS_GENERIC_64_128(default_nan, P)(P, S)
712 #define parts_silence_nan(P, S)    PARTS_GENERIC_64_128(silence_nan, P)(P, S)
713 
714 static void parts64_return_nan(FloatParts64 *a, float_status *s);
715 static void parts128_return_nan(FloatParts128 *a, float_status *s);
716 
717 #define parts_return_nan(P, S)     PARTS_GENERIC_64_128(return_nan, P)(P, S)
718 
719 static FloatParts64 *parts64_pick_nan(FloatParts64 *a, FloatParts64 *b,
720                                       float_status *s);
721 static FloatParts128 *parts128_pick_nan(FloatParts128 *a, FloatParts128 *b,
722                                         float_status *s);
723 
724 #define parts_pick_nan(A, B, S)    PARTS_GENERIC_64_128(pick_nan, A)(A, B, S)
725 
726 static FloatParts64 *parts64_pick_nan_muladd(FloatParts64 *a, FloatParts64 *b,
727                                              FloatParts64 *c, float_status *s,
728                                              int ab_mask, int abc_mask);
729 static FloatParts128 *parts128_pick_nan_muladd(FloatParts128 *a,
730                                                FloatParts128 *b,
731                                                FloatParts128 *c,
732                                                float_status *s,
733                                                int ab_mask, int abc_mask);
734 
735 #define parts_pick_nan_muladd(A, B, C, S, ABM, ABCM) \
736     PARTS_GENERIC_64_128(pick_nan_muladd, A)(A, B, C, S, ABM, ABCM)
737 
738 static void parts64_canonicalize(FloatParts64 *p, float_status *status,
739                                  const FloatFmt *fmt);
740 static void parts128_canonicalize(FloatParts128 *p, float_status *status,
741                                   const FloatFmt *fmt);
742 
743 #define parts_canonicalize(A, S, F) \
744     PARTS_GENERIC_64_128(canonicalize, A)(A, S, F)
745 
746 static void parts64_uncanon_normal(FloatParts64 *p, float_status *status,
747                                    const FloatFmt *fmt);
748 static void parts128_uncanon_normal(FloatParts128 *p, float_status *status,
749                                     const FloatFmt *fmt);
750 
751 #define parts_uncanon_normal(A, S, F) \
752     PARTS_GENERIC_64_128(uncanon_normal, A)(A, S, F)
753 
754 static void parts64_uncanon(FloatParts64 *p, float_status *status,
755                             const FloatFmt *fmt);
756 static void parts128_uncanon(FloatParts128 *p, float_status *status,
757                              const FloatFmt *fmt);
758 
759 #define parts_uncanon(A, S, F) \
760     PARTS_GENERIC_64_128(uncanon, A)(A, S, F)
761 
762 static void parts64_add_normal(FloatParts64 *a, FloatParts64 *b);
763 static void parts128_add_normal(FloatParts128 *a, FloatParts128 *b);
764 static void parts256_add_normal(FloatParts256 *a, FloatParts256 *b);
765 
766 #define parts_add_normal(A, B) \
767     PARTS_GENERIC_64_128_256(add_normal, A)(A, B)
768 
769 static bool parts64_sub_normal(FloatParts64 *a, FloatParts64 *b);
770 static bool parts128_sub_normal(FloatParts128 *a, FloatParts128 *b);
771 static bool parts256_sub_normal(FloatParts256 *a, FloatParts256 *b);
772 
773 #define parts_sub_normal(A, B) \
774     PARTS_GENERIC_64_128_256(sub_normal, A)(A, B)
775 
776 static FloatParts64 *parts64_addsub(FloatParts64 *a, FloatParts64 *b,
777                                     float_status *s, bool subtract);
778 static FloatParts128 *parts128_addsub(FloatParts128 *a, FloatParts128 *b,
779                                       float_status *s, bool subtract);
780 
781 #define parts_addsub(A, B, S, Z) \
782     PARTS_GENERIC_64_128(addsub, A)(A, B, S, Z)
783 
784 static FloatParts64 *parts64_mul(FloatParts64 *a, FloatParts64 *b,
785                                  float_status *s);
786 static FloatParts128 *parts128_mul(FloatParts128 *a, FloatParts128 *b,
787                                    float_status *s);
788 
789 #define parts_mul(A, B, S) \
790     PARTS_GENERIC_64_128(mul, A)(A, B, S)
791 
792 static FloatParts64 *parts64_muladd(FloatParts64 *a, FloatParts64 *b,
793                                     FloatParts64 *c, int flags,
794                                     float_status *s);
795 static FloatParts128 *parts128_muladd(FloatParts128 *a, FloatParts128 *b,
796                                       FloatParts128 *c, int flags,
797                                       float_status *s);
798 
799 #define parts_muladd(A, B, C, Z, S) \
800     PARTS_GENERIC_64_128(muladd, A)(A, B, C, Z, S)
801 
802 static FloatParts64 *parts64_div(FloatParts64 *a, FloatParts64 *b,
803                                  float_status *s);
804 static FloatParts128 *parts128_div(FloatParts128 *a, FloatParts128 *b,
805                                    float_status *s);
806 
807 #define parts_div(A, B, S) \
808     PARTS_GENERIC_64_128(div, A)(A, B, S)
809 
810 static FloatParts64 *parts64_modrem(FloatParts64 *a, FloatParts64 *b,
811                                     uint64_t *mod_quot, float_status *s);
812 static FloatParts128 *parts128_modrem(FloatParts128 *a, FloatParts128 *b,
813                                       uint64_t *mod_quot, float_status *s);
814 
815 #define parts_modrem(A, B, Q, S) \
816     PARTS_GENERIC_64_128(modrem, A)(A, B, Q, S)
817 
818 static void parts64_sqrt(FloatParts64 *a, float_status *s, const FloatFmt *f);
819 static void parts128_sqrt(FloatParts128 *a, float_status *s, const FloatFmt *f);
820 
821 #define parts_sqrt(A, S, F) \
822     PARTS_GENERIC_64_128(sqrt, A)(A, S, F)
823 
824 static bool parts64_round_to_int_normal(FloatParts64 *a, FloatRoundMode rm,
825                                         int scale, int frac_size);
826 static bool parts128_round_to_int_normal(FloatParts128 *a, FloatRoundMode r,
827                                          int scale, int frac_size);
828 
829 #define parts_round_to_int_normal(A, R, C, F) \
830     PARTS_GENERIC_64_128(round_to_int_normal, A)(A, R, C, F)
831 
832 static void parts64_round_to_int(FloatParts64 *a, FloatRoundMode rm,
833                                  int scale, float_status *s,
834                                  const FloatFmt *fmt);
835 static void parts128_round_to_int(FloatParts128 *a, FloatRoundMode r,
836                                   int scale, float_status *s,
837                                   const FloatFmt *fmt);
838 
839 #define parts_round_to_int(A, R, C, S, F) \
840     PARTS_GENERIC_64_128(round_to_int, A)(A, R, C, S, F)
841 
842 static int64_t parts64_float_to_sint(FloatParts64 *p, FloatRoundMode rmode,
843                                      int scale, int64_t min, int64_t max,
844                                      float_status *s);
845 static int64_t parts128_float_to_sint(FloatParts128 *p, FloatRoundMode rmode,
846                                      int scale, int64_t min, int64_t max,
847                                      float_status *s);
848 
849 #define parts_float_to_sint(P, R, Z, MN, MX, S) \
850     PARTS_GENERIC_64_128(float_to_sint, P)(P, R, Z, MN, MX, S)
851 
852 static uint64_t parts64_float_to_uint(FloatParts64 *p, FloatRoundMode rmode,
853                                       int scale, uint64_t max,
854                                       float_status *s);
855 static uint64_t parts128_float_to_uint(FloatParts128 *p, FloatRoundMode rmode,
856                                        int scale, uint64_t max,
857                                        float_status *s);
858 
859 #define parts_float_to_uint(P, R, Z, M, S) \
860     PARTS_GENERIC_64_128(float_to_uint, P)(P, R, Z, M, S)
861 
862 static int64_t parts64_float_to_sint_modulo(FloatParts64 *p,
863                                             FloatRoundMode rmode,
864                                             int bitsm1, float_status *s);
865 static int64_t parts128_float_to_sint_modulo(FloatParts128 *p,
866                                              FloatRoundMode rmode,
867                                              int bitsm1, float_status *s);
868 
869 #define parts_float_to_sint_modulo(P, R, M, S) \
870     PARTS_GENERIC_64_128(float_to_sint_modulo, P)(P, R, M, S)
871 
872 static void parts64_sint_to_float(FloatParts64 *p, int64_t a,
873                                   int scale, float_status *s);
874 static void parts128_sint_to_float(FloatParts128 *p, int64_t a,
875                                    int scale, float_status *s);
876 
877 #define parts_float_to_sint(P, R, Z, MN, MX, S) \
878     PARTS_GENERIC_64_128(float_to_sint, P)(P, R, Z, MN, MX, S)
879 
880 #define parts_sint_to_float(P, I, Z, S) \
881     PARTS_GENERIC_64_128(sint_to_float, P)(P, I, Z, S)
882 
883 static void parts64_uint_to_float(FloatParts64 *p, uint64_t a,
884                                   int scale, float_status *s);
885 static void parts128_uint_to_float(FloatParts128 *p, uint64_t a,
886                                    int scale, float_status *s);
887 
888 #define parts_uint_to_float(P, I, Z, S) \
889     PARTS_GENERIC_64_128(uint_to_float, P)(P, I, Z, S)
890 
891 static FloatParts64 *parts64_minmax(FloatParts64 *a, FloatParts64 *b,
892                                     float_status *s, int flags);
893 static FloatParts128 *parts128_minmax(FloatParts128 *a, FloatParts128 *b,
894                                       float_status *s, int flags);
895 
896 #define parts_minmax(A, B, S, F) \
897     PARTS_GENERIC_64_128(minmax, A)(A, B, S, F)
898 
899 static FloatRelation parts64_compare(FloatParts64 *a, FloatParts64 *b,
900                                      float_status *s, bool q);
901 static FloatRelation parts128_compare(FloatParts128 *a, FloatParts128 *b,
902                                       float_status *s, bool q);
903 
904 #define parts_compare(A, B, S, Q) \
905     PARTS_GENERIC_64_128(compare, A)(A, B, S, Q)
906 
907 static void parts64_scalbn(FloatParts64 *a, int n, float_status *s);
908 static void parts128_scalbn(FloatParts128 *a, int n, float_status *s);
909 
910 #define parts_scalbn(A, N, S) \
911     PARTS_GENERIC_64_128(scalbn, A)(A, N, S)
912 
913 static void parts64_log2(FloatParts64 *a, float_status *s, const FloatFmt *f);
914 static void parts128_log2(FloatParts128 *a, float_status *s, const FloatFmt *f);
915 
916 #define parts_log2(A, S, F) \
917     PARTS_GENERIC_64_128(log2, A)(A, S, F)
918 
919 /*
920  * Helper functions for softfloat-parts.c.inc, per-size operations.
921  */
922 
923 #define FRAC_GENERIC_64_128(NAME, P) \
924     _Generic((P), FloatParts64 *: frac64_##NAME, \
925                   FloatParts128 *: frac128_##NAME)
926 
927 #define FRAC_GENERIC_64_128_256(NAME, P) \
928     _Generic((P), FloatParts64 *: frac64_##NAME, \
929                   FloatParts128 *: frac128_##NAME, \
930                   FloatParts256 *: frac256_##NAME)
931 
932 static bool frac64_add(FloatParts64 *r, FloatParts64 *a, FloatParts64 *b)
933 {
934     return uadd64_overflow(a->frac, b->frac, &r->frac);
935 }
936 
937 static bool frac128_add(FloatParts128 *r, FloatParts128 *a, FloatParts128 *b)
938 {
939     bool c = 0;
940     r->frac_lo = uadd64_carry(a->frac_lo, b->frac_lo, &c);
941     r->frac_hi = uadd64_carry(a->frac_hi, b->frac_hi, &c);
942     return c;
943 }
944 
945 static bool frac256_add(FloatParts256 *r, FloatParts256 *a, FloatParts256 *b)
946 {
947     bool c = 0;
948     r->frac_lo = uadd64_carry(a->frac_lo, b->frac_lo, &c);
949     r->frac_lm = uadd64_carry(a->frac_lm, b->frac_lm, &c);
950     r->frac_hm = uadd64_carry(a->frac_hm, b->frac_hm, &c);
951     r->frac_hi = uadd64_carry(a->frac_hi, b->frac_hi, &c);
952     return c;
953 }
954 
955 #define frac_add(R, A, B)  FRAC_GENERIC_64_128_256(add, R)(R, A, B)
956 
957 static bool frac64_addi(FloatParts64 *r, FloatParts64 *a, uint64_t c)
958 {
959     return uadd64_overflow(a->frac, c, &r->frac);
960 }
961 
962 static bool frac128_addi(FloatParts128 *r, FloatParts128 *a, uint64_t c)
963 {
964     c = uadd64_overflow(a->frac_lo, c, &r->frac_lo);
965     return uadd64_overflow(a->frac_hi, c, &r->frac_hi);
966 }
967 
968 #define frac_addi(R, A, C)  FRAC_GENERIC_64_128(addi, R)(R, A, C)
969 
970 static void frac64_allones(FloatParts64 *a)
971 {
972     a->frac = -1;
973 }
974 
975 static void frac128_allones(FloatParts128 *a)
976 {
977     a->frac_hi = a->frac_lo = -1;
978 }
979 
980 #define frac_allones(A)  FRAC_GENERIC_64_128(allones, A)(A)
981 
982 static FloatRelation frac64_cmp(FloatParts64 *a, FloatParts64 *b)
983 {
984     return (a->frac == b->frac ? float_relation_equal
985             : a->frac < b->frac ? float_relation_less
986             : float_relation_greater);
987 }
988 
989 static FloatRelation frac128_cmp(FloatParts128 *a, FloatParts128 *b)
990 {
991     uint64_t ta = a->frac_hi, tb = b->frac_hi;
992     if (ta == tb) {
993         ta = a->frac_lo, tb = b->frac_lo;
994         if (ta == tb) {
995             return float_relation_equal;
996         }
997     }
998     return ta < tb ? float_relation_less : float_relation_greater;
999 }
1000 
1001 #define frac_cmp(A, B)  FRAC_GENERIC_64_128(cmp, A)(A, B)
1002 
1003 static void frac64_clear(FloatParts64 *a)
1004 {
1005     a->frac = 0;
1006 }
1007 
1008 static void frac128_clear(FloatParts128 *a)
1009 {
1010     a->frac_hi = a->frac_lo = 0;
1011 }
1012 
1013 #define frac_clear(A)  FRAC_GENERIC_64_128(clear, A)(A)
1014 
1015 static bool frac64_div(FloatParts64 *a, FloatParts64 *b)
1016 {
1017     uint64_t n1, n0, r, q;
1018     bool ret;
1019 
1020     /*
1021      * We want a 2*N / N-bit division to produce exactly an N-bit
1022      * result, so that we do not lose any precision and so that we
1023      * do not have to renormalize afterward.  If A.frac < B.frac,
1024      * then division would produce an (N-1)-bit result; shift A left
1025      * by one to produce the an N-bit result, and return true to
1026      * decrement the exponent to match.
1027      *
1028      * The udiv_qrnnd algorithm that we're using requires normalization,
1029      * i.e. the msb of the denominator must be set, which is already true.
1030      */
1031     ret = a->frac < b->frac;
1032     if (ret) {
1033         n0 = a->frac;
1034         n1 = 0;
1035     } else {
1036         n0 = a->frac >> 1;
1037         n1 = a->frac << 63;
1038     }
1039     q = udiv_qrnnd(&r, n0, n1, b->frac);
1040 
1041     /* Set lsb if there is a remainder, to set inexact. */
1042     a->frac = q | (r != 0);
1043 
1044     return ret;
1045 }
1046 
1047 static bool frac128_div(FloatParts128 *a, FloatParts128 *b)
1048 {
1049     uint64_t q0, q1, a0, a1, b0, b1;
1050     uint64_t r0, r1, r2, r3, t0, t1, t2, t3;
1051     bool ret = false;
1052 
1053     a0 = a->frac_hi, a1 = a->frac_lo;
1054     b0 = b->frac_hi, b1 = b->frac_lo;
1055 
1056     ret = lt128(a0, a1, b0, b1);
1057     if (!ret) {
1058         a1 = shr_double(a0, a1, 1);
1059         a0 = a0 >> 1;
1060     }
1061 
1062     /* Use 128/64 -> 64 division as estimate for 192/128 -> 128 division. */
1063     q0 = estimateDiv128To64(a0, a1, b0);
1064 
1065     /*
1066      * Estimate is high because B1 was not included (unless B1 == 0).
1067      * Reduce quotient and increase remainder until remainder is non-negative.
1068      * This loop will execute 0 to 2 times.
1069      */
1070     mul128By64To192(b0, b1, q0, &t0, &t1, &t2);
1071     sub192(a0, a1, 0, t0, t1, t2, &r0, &r1, &r2);
1072     while (r0 != 0) {
1073         q0--;
1074         add192(r0, r1, r2, 0, b0, b1, &r0, &r1, &r2);
1075     }
1076 
1077     /* Repeat using the remainder, producing a second word of quotient. */
1078     q1 = estimateDiv128To64(r1, r2, b0);
1079     mul128By64To192(b0, b1, q1, &t1, &t2, &t3);
1080     sub192(r1, r2, 0, t1, t2, t3, &r1, &r2, &r3);
1081     while (r1 != 0) {
1082         q1--;
1083         add192(r1, r2, r3, 0, b0, b1, &r1, &r2, &r3);
1084     }
1085 
1086     /* Any remainder indicates inexact; set sticky bit. */
1087     q1 |= (r2 | r3) != 0;
1088 
1089     a->frac_hi = q0;
1090     a->frac_lo = q1;
1091     return ret;
1092 }
1093 
1094 #define frac_div(A, B)  FRAC_GENERIC_64_128(div, A)(A, B)
1095 
1096 static bool frac64_eqz(FloatParts64 *a)
1097 {
1098     return a->frac == 0;
1099 }
1100 
1101 static bool frac128_eqz(FloatParts128 *a)
1102 {
1103     return (a->frac_hi | a->frac_lo) == 0;
1104 }
1105 
1106 #define frac_eqz(A)  FRAC_GENERIC_64_128(eqz, A)(A)
1107 
1108 static void frac64_mulw(FloatParts128 *r, FloatParts64 *a, FloatParts64 *b)
1109 {
1110     mulu64(&r->frac_lo, &r->frac_hi, a->frac, b->frac);
1111 }
1112 
1113 static void frac128_mulw(FloatParts256 *r, FloatParts128 *a, FloatParts128 *b)
1114 {
1115     mul128To256(a->frac_hi, a->frac_lo, b->frac_hi, b->frac_lo,
1116                 &r->frac_hi, &r->frac_hm, &r->frac_lm, &r->frac_lo);
1117 }
1118 
1119 #define frac_mulw(R, A, B)  FRAC_GENERIC_64_128(mulw, A)(R, A, B)
1120 
1121 static void frac64_neg(FloatParts64 *a)
1122 {
1123     a->frac = -a->frac;
1124 }
1125 
1126 static void frac128_neg(FloatParts128 *a)
1127 {
1128     bool c = 0;
1129     a->frac_lo = usub64_borrow(0, a->frac_lo, &c);
1130     a->frac_hi = usub64_borrow(0, a->frac_hi, &c);
1131 }
1132 
1133 static void frac256_neg(FloatParts256 *a)
1134 {
1135     bool c = 0;
1136     a->frac_lo = usub64_borrow(0, a->frac_lo, &c);
1137     a->frac_lm = usub64_borrow(0, a->frac_lm, &c);
1138     a->frac_hm = usub64_borrow(0, a->frac_hm, &c);
1139     a->frac_hi = usub64_borrow(0, a->frac_hi, &c);
1140 }
1141 
1142 #define frac_neg(A)  FRAC_GENERIC_64_128_256(neg, A)(A)
1143 
1144 static int frac64_normalize(FloatParts64 *a)
1145 {
1146     if (a->frac) {
1147         int shift = clz64(a->frac);
1148         a->frac <<= shift;
1149         return shift;
1150     }
1151     return 64;
1152 }
1153 
1154 static int frac128_normalize(FloatParts128 *a)
1155 {
1156     if (a->frac_hi) {
1157         int shl = clz64(a->frac_hi);
1158         a->frac_hi = shl_double(a->frac_hi, a->frac_lo, shl);
1159         a->frac_lo <<= shl;
1160         return shl;
1161     } else if (a->frac_lo) {
1162         int shl = clz64(a->frac_lo);
1163         a->frac_hi = a->frac_lo << shl;
1164         a->frac_lo = 0;
1165         return shl + 64;
1166     }
1167     return 128;
1168 }
1169 
1170 static int frac256_normalize(FloatParts256 *a)
1171 {
1172     uint64_t a0 = a->frac_hi, a1 = a->frac_hm;
1173     uint64_t a2 = a->frac_lm, a3 = a->frac_lo;
1174     int ret, shl;
1175 
1176     if (likely(a0)) {
1177         shl = clz64(a0);
1178         if (shl == 0) {
1179             return 0;
1180         }
1181         ret = shl;
1182     } else {
1183         if (a1) {
1184             ret = 64;
1185             a0 = a1, a1 = a2, a2 = a3, a3 = 0;
1186         } else if (a2) {
1187             ret = 128;
1188             a0 = a2, a1 = a3, a2 = 0, a3 = 0;
1189         } else if (a3) {
1190             ret = 192;
1191             a0 = a3, a1 = 0, a2 = 0, a3 = 0;
1192         } else {
1193             ret = 256;
1194             a0 = 0, a1 = 0, a2 = 0, a3 = 0;
1195             goto done;
1196         }
1197         shl = clz64(a0);
1198         if (shl == 0) {
1199             goto done;
1200         }
1201         ret += shl;
1202     }
1203 
1204     a0 = shl_double(a0, a1, shl);
1205     a1 = shl_double(a1, a2, shl);
1206     a2 = shl_double(a2, a3, shl);
1207     a3 <<= shl;
1208 
1209  done:
1210     a->frac_hi = a0;
1211     a->frac_hm = a1;
1212     a->frac_lm = a2;
1213     a->frac_lo = a3;
1214     return ret;
1215 }
1216 
1217 #define frac_normalize(A)  FRAC_GENERIC_64_128_256(normalize, A)(A)
1218 
1219 static void frac64_modrem(FloatParts64 *a, FloatParts64 *b, uint64_t *mod_quot)
1220 {
1221     uint64_t a0, a1, b0, t0, t1, q, quot;
1222     int exp_diff = a->exp - b->exp;
1223     int shift;
1224 
1225     a0 = a->frac;
1226     a1 = 0;
1227 
1228     if (exp_diff < -1) {
1229         if (mod_quot) {
1230             *mod_quot = 0;
1231         }
1232         return;
1233     }
1234     if (exp_diff == -1) {
1235         a0 >>= 1;
1236         exp_diff = 0;
1237     }
1238 
1239     b0 = b->frac;
1240     quot = q = b0 <= a0;
1241     if (q) {
1242         a0 -= b0;
1243     }
1244 
1245     exp_diff -= 64;
1246     while (exp_diff > 0) {
1247         q = estimateDiv128To64(a0, a1, b0);
1248         q = q > 2 ? q - 2 : 0;
1249         mul64To128(b0, q, &t0, &t1);
1250         sub128(a0, a1, t0, t1, &a0, &a1);
1251         shortShift128Left(a0, a1, 62, &a0, &a1);
1252         exp_diff -= 62;
1253         quot = (quot << 62) + q;
1254     }
1255 
1256     exp_diff += 64;
1257     if (exp_diff > 0) {
1258         q = estimateDiv128To64(a0, a1, b0);
1259         q = q > 2 ? (q - 2) >> (64 - exp_diff) : 0;
1260         mul64To128(b0, q << (64 - exp_diff), &t0, &t1);
1261         sub128(a0, a1, t0, t1, &a0, &a1);
1262         shortShift128Left(0, b0, 64 - exp_diff, &t0, &t1);
1263         while (le128(t0, t1, a0, a1)) {
1264             ++q;
1265             sub128(a0, a1, t0, t1, &a0, &a1);
1266         }
1267         quot = (exp_diff < 64 ? quot << exp_diff : 0) + q;
1268     } else {
1269         t0 = b0;
1270         t1 = 0;
1271     }
1272 
1273     if (mod_quot) {
1274         *mod_quot = quot;
1275     } else {
1276         sub128(t0, t1, a0, a1, &t0, &t1);
1277         if (lt128(t0, t1, a0, a1) ||
1278             (eq128(t0, t1, a0, a1) && (q & 1))) {
1279             a0 = t0;
1280             a1 = t1;
1281             a->sign = !a->sign;
1282         }
1283     }
1284 
1285     if (likely(a0)) {
1286         shift = clz64(a0);
1287         shortShift128Left(a0, a1, shift, &a0, &a1);
1288     } else if (likely(a1)) {
1289         shift = clz64(a1);
1290         a0 = a1 << shift;
1291         a1 = 0;
1292         shift += 64;
1293     } else {
1294         a->cls = float_class_zero;
1295         return;
1296     }
1297 
1298     a->exp = b->exp + exp_diff - shift;
1299     a->frac = a0 | (a1 != 0);
1300 }
1301 
1302 static void frac128_modrem(FloatParts128 *a, FloatParts128 *b,
1303                            uint64_t *mod_quot)
1304 {
1305     uint64_t a0, a1, a2, b0, b1, t0, t1, t2, q, quot;
1306     int exp_diff = a->exp - b->exp;
1307     int shift;
1308 
1309     a0 = a->frac_hi;
1310     a1 = a->frac_lo;
1311     a2 = 0;
1312 
1313     if (exp_diff < -1) {
1314         if (mod_quot) {
1315             *mod_quot = 0;
1316         }
1317         return;
1318     }
1319     if (exp_diff == -1) {
1320         shift128Right(a0, a1, 1, &a0, &a1);
1321         exp_diff = 0;
1322     }
1323 
1324     b0 = b->frac_hi;
1325     b1 = b->frac_lo;
1326 
1327     quot = q = le128(b0, b1, a0, a1);
1328     if (q) {
1329         sub128(a0, a1, b0, b1, &a0, &a1);
1330     }
1331 
1332     exp_diff -= 64;
1333     while (exp_diff > 0) {
1334         q = estimateDiv128To64(a0, a1, b0);
1335         q = q > 4 ? q - 4 : 0;
1336         mul128By64To192(b0, b1, q, &t0, &t1, &t2);
1337         sub192(a0, a1, a2, t0, t1, t2, &a0, &a1, &a2);
1338         shortShift192Left(a0, a1, a2, 61, &a0, &a1, &a2);
1339         exp_diff -= 61;
1340         quot = (quot << 61) + q;
1341     }
1342 
1343     exp_diff += 64;
1344     if (exp_diff > 0) {
1345         q = estimateDiv128To64(a0, a1, b0);
1346         q = q > 4 ? (q - 4) >> (64 - exp_diff) : 0;
1347         mul128By64To192(b0, b1, q << (64 - exp_diff), &t0, &t1, &t2);
1348         sub192(a0, a1, a2, t0, t1, t2, &a0, &a1, &a2);
1349         shortShift192Left(0, b0, b1, 64 - exp_diff, &t0, &t1, &t2);
1350         while (le192(t0, t1, t2, a0, a1, a2)) {
1351             ++q;
1352             sub192(a0, a1, a2, t0, t1, t2, &a0, &a1, &a2);
1353         }
1354         quot = (exp_diff < 64 ? quot << exp_diff : 0) + q;
1355     } else {
1356         t0 = b0;
1357         t1 = b1;
1358         t2 = 0;
1359     }
1360 
1361     if (mod_quot) {
1362         *mod_quot = quot;
1363     } else {
1364         sub192(t0, t1, t2, a0, a1, a2, &t0, &t1, &t2);
1365         if (lt192(t0, t1, t2, a0, a1, a2) ||
1366             (eq192(t0, t1, t2, a0, a1, a2) && (q & 1))) {
1367             a0 = t0;
1368             a1 = t1;
1369             a2 = t2;
1370             a->sign = !a->sign;
1371         }
1372     }
1373 
1374     if (likely(a0)) {
1375         shift = clz64(a0);
1376         shortShift192Left(a0, a1, a2, shift, &a0, &a1, &a2);
1377     } else if (likely(a1)) {
1378         shift = clz64(a1);
1379         shortShift128Left(a1, a2, shift, &a0, &a1);
1380         a2 = 0;
1381         shift += 64;
1382     } else if (likely(a2)) {
1383         shift = clz64(a2);
1384         a0 = a2 << shift;
1385         a1 = a2 = 0;
1386         shift += 128;
1387     } else {
1388         a->cls = float_class_zero;
1389         return;
1390     }
1391 
1392     a->exp = b->exp + exp_diff - shift;
1393     a->frac_hi = a0;
1394     a->frac_lo = a1 | (a2 != 0);
1395 }
1396 
1397 #define frac_modrem(A, B, Q)  FRAC_GENERIC_64_128(modrem, A)(A, B, Q)
1398 
1399 static void frac64_shl(FloatParts64 *a, int c)
1400 {
1401     a->frac <<= c;
1402 }
1403 
1404 static void frac128_shl(FloatParts128 *a, int c)
1405 {
1406     uint64_t a0 = a->frac_hi, a1 = a->frac_lo;
1407 
1408     if (c & 64) {
1409         a0 = a1, a1 = 0;
1410     }
1411 
1412     c &= 63;
1413     if (c) {
1414         a0 = shl_double(a0, a1, c);
1415         a1 = a1 << c;
1416     }
1417 
1418     a->frac_hi = a0;
1419     a->frac_lo = a1;
1420 }
1421 
1422 #define frac_shl(A, C)  FRAC_GENERIC_64_128(shl, A)(A, C)
1423 
1424 static void frac64_shr(FloatParts64 *a, int c)
1425 {
1426     a->frac >>= c;
1427 }
1428 
1429 static void frac128_shr(FloatParts128 *a, int c)
1430 {
1431     uint64_t a0 = a->frac_hi, a1 = a->frac_lo;
1432 
1433     if (c & 64) {
1434         a1 = a0, a0 = 0;
1435     }
1436 
1437     c &= 63;
1438     if (c) {
1439         a1 = shr_double(a0, a1, c);
1440         a0 = a0 >> c;
1441     }
1442 
1443     a->frac_hi = a0;
1444     a->frac_lo = a1;
1445 }
1446 
1447 #define frac_shr(A, C)  FRAC_GENERIC_64_128(shr, A)(A, C)
1448 
1449 static void frac64_shrjam(FloatParts64 *a, int c)
1450 {
1451     uint64_t a0 = a->frac;
1452 
1453     if (likely(c != 0)) {
1454         if (likely(c < 64)) {
1455             a0 = (a0 >> c) | (shr_double(a0, 0, c) != 0);
1456         } else {
1457             a0 = a0 != 0;
1458         }
1459         a->frac = a0;
1460     }
1461 }
1462 
1463 static void frac128_shrjam(FloatParts128 *a, int c)
1464 {
1465     uint64_t a0 = a->frac_hi, a1 = a->frac_lo;
1466     uint64_t sticky = 0;
1467 
1468     if (unlikely(c == 0)) {
1469         return;
1470     } else if (likely(c < 64)) {
1471         /* nothing */
1472     } else if (likely(c < 128)) {
1473         sticky = a1;
1474         a1 = a0;
1475         a0 = 0;
1476         c &= 63;
1477         if (c == 0) {
1478             goto done;
1479         }
1480     } else {
1481         sticky = a0 | a1;
1482         a0 = a1 = 0;
1483         goto done;
1484     }
1485 
1486     sticky |= shr_double(a1, 0, c);
1487     a1 = shr_double(a0, a1, c);
1488     a0 = a0 >> c;
1489 
1490  done:
1491     a->frac_lo = a1 | (sticky != 0);
1492     a->frac_hi = a0;
1493 }
1494 
1495 static void frac256_shrjam(FloatParts256 *a, int c)
1496 {
1497     uint64_t a0 = a->frac_hi, a1 = a->frac_hm;
1498     uint64_t a2 = a->frac_lm, a3 = a->frac_lo;
1499     uint64_t sticky = 0;
1500 
1501     if (unlikely(c == 0)) {
1502         return;
1503     } else if (likely(c < 64)) {
1504         /* nothing */
1505     } else if (likely(c < 256)) {
1506         if (unlikely(c & 128)) {
1507             sticky |= a2 | a3;
1508             a3 = a1, a2 = a0, a1 = 0, a0 = 0;
1509         }
1510         if (unlikely(c & 64)) {
1511             sticky |= a3;
1512             a3 = a2, a2 = a1, a1 = a0, a0 = 0;
1513         }
1514         c &= 63;
1515         if (c == 0) {
1516             goto done;
1517         }
1518     } else {
1519         sticky = a0 | a1 | a2 | a3;
1520         a0 = a1 = a2 = a3 = 0;
1521         goto done;
1522     }
1523 
1524     sticky |= shr_double(a3, 0, c);
1525     a3 = shr_double(a2, a3, c);
1526     a2 = shr_double(a1, a2, c);
1527     a1 = shr_double(a0, a1, c);
1528     a0 = a0 >> c;
1529 
1530  done:
1531     a->frac_lo = a3 | (sticky != 0);
1532     a->frac_lm = a2;
1533     a->frac_hm = a1;
1534     a->frac_hi = a0;
1535 }
1536 
1537 #define frac_shrjam(A, C)  FRAC_GENERIC_64_128_256(shrjam, A)(A, C)
1538 
1539 static bool frac64_sub(FloatParts64 *r, FloatParts64 *a, FloatParts64 *b)
1540 {
1541     return usub64_overflow(a->frac, b->frac, &r->frac);
1542 }
1543 
1544 static bool frac128_sub(FloatParts128 *r, FloatParts128 *a, FloatParts128 *b)
1545 {
1546     bool c = 0;
1547     r->frac_lo = usub64_borrow(a->frac_lo, b->frac_lo, &c);
1548     r->frac_hi = usub64_borrow(a->frac_hi, b->frac_hi, &c);
1549     return c;
1550 }
1551 
1552 static bool frac256_sub(FloatParts256 *r, FloatParts256 *a, FloatParts256 *b)
1553 {
1554     bool c = 0;
1555     r->frac_lo = usub64_borrow(a->frac_lo, b->frac_lo, &c);
1556     r->frac_lm = usub64_borrow(a->frac_lm, b->frac_lm, &c);
1557     r->frac_hm = usub64_borrow(a->frac_hm, b->frac_hm, &c);
1558     r->frac_hi = usub64_borrow(a->frac_hi, b->frac_hi, &c);
1559     return c;
1560 }
1561 
1562 #define frac_sub(R, A, B)  FRAC_GENERIC_64_128_256(sub, R)(R, A, B)
1563 
1564 static void frac64_truncjam(FloatParts64 *r, FloatParts128 *a)
1565 {
1566     r->frac = a->frac_hi | (a->frac_lo != 0);
1567 }
1568 
1569 static void frac128_truncjam(FloatParts128 *r, FloatParts256 *a)
1570 {
1571     r->frac_hi = a->frac_hi;
1572     r->frac_lo = a->frac_hm | ((a->frac_lm | a->frac_lo) != 0);
1573 }
1574 
1575 #define frac_truncjam(R, A)  FRAC_GENERIC_64_128(truncjam, R)(R, A)
1576 
1577 static void frac64_widen(FloatParts128 *r, FloatParts64 *a)
1578 {
1579     r->frac_hi = a->frac;
1580     r->frac_lo = 0;
1581 }
1582 
1583 static void frac128_widen(FloatParts256 *r, FloatParts128 *a)
1584 {
1585     r->frac_hi = a->frac_hi;
1586     r->frac_hm = a->frac_lo;
1587     r->frac_lm = 0;
1588     r->frac_lo = 0;
1589 }
1590 
1591 #define frac_widen(A, B)  FRAC_GENERIC_64_128(widen, B)(A, B)
1592 
1593 /*
1594  * Reciprocal sqrt table.  1 bit of exponent, 6-bits of mantessa.
1595  * From https://git.musl-libc.org/cgit/musl/tree/src/math/sqrt_data.c
1596  * and thus MIT licenced.
1597  */
1598 static const uint16_t rsqrt_tab[128] = {
1599     0xb451, 0xb2f0, 0xb196, 0xb044, 0xaef9, 0xadb6, 0xac79, 0xab43,
1600     0xaa14, 0xa8eb, 0xa7c8, 0xa6aa, 0xa592, 0xa480, 0xa373, 0xa26b,
1601     0xa168, 0xa06a, 0x9f70, 0x9e7b, 0x9d8a, 0x9c9d, 0x9bb5, 0x9ad1,
1602     0x99f0, 0x9913, 0x983a, 0x9765, 0x9693, 0x95c4, 0x94f8, 0x9430,
1603     0x936b, 0x92a9, 0x91ea, 0x912e, 0x9075, 0x8fbe, 0x8f0a, 0x8e59,
1604     0x8daa, 0x8cfe, 0x8c54, 0x8bac, 0x8b07, 0x8a64, 0x89c4, 0x8925,
1605     0x8889, 0x87ee, 0x8756, 0x86c0, 0x862b, 0x8599, 0x8508, 0x8479,
1606     0x83ec, 0x8361, 0x82d8, 0x8250, 0x81c9, 0x8145, 0x80c2, 0x8040,
1607     0xff02, 0xfd0e, 0xfb25, 0xf947, 0xf773, 0xf5aa, 0xf3ea, 0xf234,
1608     0xf087, 0xeee3, 0xed47, 0xebb3, 0xea27, 0xe8a3, 0xe727, 0xe5b2,
1609     0xe443, 0xe2dc, 0xe17a, 0xe020, 0xdecb, 0xdd7d, 0xdc34, 0xdaf1,
1610     0xd9b3, 0xd87b, 0xd748, 0xd61a, 0xd4f1, 0xd3cd, 0xd2ad, 0xd192,
1611     0xd07b, 0xcf69, 0xce5b, 0xcd51, 0xcc4a, 0xcb48, 0xca4a, 0xc94f,
1612     0xc858, 0xc764, 0xc674, 0xc587, 0xc49d, 0xc3b7, 0xc2d4, 0xc1f4,
1613     0xc116, 0xc03c, 0xbf65, 0xbe90, 0xbdbe, 0xbcef, 0xbc23, 0xbb59,
1614     0xba91, 0xb9cc, 0xb90a, 0xb84a, 0xb78c, 0xb6d0, 0xb617, 0xb560,
1615 };
1616 
1617 #define partsN(NAME)   glue(glue(glue(parts,N),_),NAME)
1618 #define FloatPartsN    glue(FloatParts,N)
1619 #define FloatPartsW    glue(FloatParts,W)
1620 
1621 #define N 64
1622 #define W 128
1623 
1624 #include "softfloat-parts-addsub.c.inc"
1625 #include "softfloat-parts.c.inc"
1626 
1627 #undef  N
1628 #undef  W
1629 #define N 128
1630 #define W 256
1631 
1632 #include "softfloat-parts-addsub.c.inc"
1633 #include "softfloat-parts.c.inc"
1634 
1635 #undef  N
1636 #undef  W
1637 #define N            256
1638 
1639 #include "softfloat-parts-addsub.c.inc"
1640 
1641 #undef  N
1642 #undef  W
1643 #undef  partsN
1644 #undef  FloatPartsN
1645 #undef  FloatPartsW
1646 
1647 /*
1648  * Pack/unpack routines with a specific FloatFmt.
1649  */
1650 
1651 static void float16a_unpack_canonical(FloatParts64 *p, float16 f,
1652                                       float_status *s, const FloatFmt *params)
1653 {
1654     float16_unpack_raw(p, f);
1655     parts_canonicalize(p, s, params);
1656 }
1657 
1658 static void float16_unpack_canonical(FloatParts64 *p, float16 f,
1659                                      float_status *s)
1660 {
1661     float16a_unpack_canonical(p, f, s, &float16_params);
1662 }
1663 
1664 static void bfloat16_unpack_canonical(FloatParts64 *p, bfloat16 f,
1665                                       float_status *s)
1666 {
1667     bfloat16_unpack_raw(p, f);
1668     parts_canonicalize(p, s, &bfloat16_params);
1669 }
1670 
1671 static float16 float16a_round_pack_canonical(FloatParts64 *p,
1672                                              float_status *s,
1673                                              const FloatFmt *params)
1674 {
1675     parts_uncanon(p, s, params);
1676     return float16_pack_raw(p);
1677 }
1678 
1679 static float16 float16_round_pack_canonical(FloatParts64 *p,
1680                                             float_status *s)
1681 {
1682     return float16a_round_pack_canonical(p, s, &float16_params);
1683 }
1684 
1685 static bfloat16 bfloat16_round_pack_canonical(FloatParts64 *p,
1686                                               float_status *s)
1687 {
1688     parts_uncanon(p, s, &bfloat16_params);
1689     return bfloat16_pack_raw(p);
1690 }
1691 
1692 static void float32_unpack_canonical(FloatParts64 *p, float32 f,
1693                                      float_status *s)
1694 {
1695     float32_unpack_raw(p, f);
1696     parts_canonicalize(p, s, &float32_params);
1697 }
1698 
1699 static float32 float32_round_pack_canonical(FloatParts64 *p,
1700                                             float_status *s)
1701 {
1702     parts_uncanon(p, s, &float32_params);
1703     return float32_pack_raw(p);
1704 }
1705 
1706 static void float64_unpack_canonical(FloatParts64 *p, float64 f,
1707                                      float_status *s)
1708 {
1709     float64_unpack_raw(p, f);
1710     parts_canonicalize(p, s, &float64_params);
1711 }
1712 
1713 static float64 float64_round_pack_canonical(FloatParts64 *p,
1714                                             float_status *s)
1715 {
1716     parts_uncanon(p, s, &float64_params);
1717     return float64_pack_raw(p);
1718 }
1719 
1720 static float64 float64r32_round_pack_canonical(FloatParts64 *p,
1721                                                float_status *s)
1722 {
1723     parts_uncanon(p, s, &float32_params);
1724 
1725     /*
1726      * In parts_uncanon, we placed the fraction for float32 at the lsb.
1727      * We need to adjust the fraction higher so that the least N bits are
1728      * zero, and the fraction is adjacent to the float64 implicit bit.
1729      */
1730     switch (p->cls) {
1731     case float_class_normal:
1732         if (unlikely(p->exp == 0)) {
1733             /*
1734              * The result is denormal for float32, but can be represented
1735              * in normalized form for float64.  Adjust, per canonicalize.
1736              */
1737             int shift = frac_normalize(p);
1738             p->exp = (float32_params.frac_shift -
1739                       float32_params.exp_bias - shift + 1 +
1740                       float64_params.exp_bias);
1741             frac_shr(p, float64_params.frac_shift);
1742         } else {
1743             frac_shl(p, float32_params.frac_shift - float64_params.frac_shift);
1744             p->exp += float64_params.exp_bias - float32_params.exp_bias;
1745         }
1746         break;
1747     case float_class_snan:
1748     case float_class_qnan:
1749         frac_shl(p, float32_params.frac_shift - float64_params.frac_shift);
1750         p->exp = float64_params.exp_max;
1751         break;
1752     case float_class_inf:
1753         p->exp = float64_params.exp_max;
1754         break;
1755     case float_class_zero:
1756         break;
1757     default:
1758         g_assert_not_reached();
1759     }
1760 
1761     return float64_pack_raw(p);
1762 }
1763 
1764 static void float128_unpack_canonical(FloatParts128 *p, float128 f,
1765                                       float_status *s)
1766 {
1767     float128_unpack_raw(p, f);
1768     parts_canonicalize(p, s, &float128_params);
1769 }
1770 
1771 static float128 float128_round_pack_canonical(FloatParts128 *p,
1772                                               float_status *s)
1773 {
1774     parts_uncanon(p, s, &float128_params);
1775     return float128_pack_raw(p);
1776 }
1777 
1778 /* Returns false if the encoding is invalid. */
1779 static bool floatx80_unpack_canonical(FloatParts128 *p, floatx80 f,
1780                                       float_status *s)
1781 {
1782     /* Ensure rounding precision is set before beginning. */
1783     switch (s->floatx80_rounding_precision) {
1784     case floatx80_precision_x:
1785     case floatx80_precision_d:
1786     case floatx80_precision_s:
1787         break;
1788     default:
1789         g_assert_not_reached();
1790     }
1791 
1792     if (unlikely(floatx80_invalid_encoding(f))) {
1793         float_raise(float_flag_invalid, s);
1794         return false;
1795     }
1796 
1797     floatx80_unpack_raw(p, f);
1798 
1799     if (likely(p->exp != floatx80_params[floatx80_precision_x].exp_max)) {
1800         parts_canonicalize(p, s, &floatx80_params[floatx80_precision_x]);
1801     } else {
1802         /* The explicit integer bit is ignored, after invalid checks. */
1803         p->frac_hi &= MAKE_64BIT_MASK(0, 63);
1804         p->cls = (p->frac_hi == 0 ? float_class_inf
1805                   : parts_is_snan_frac(p->frac_hi, s)
1806                   ? float_class_snan : float_class_qnan);
1807     }
1808     return true;
1809 }
1810 
1811 static floatx80 floatx80_round_pack_canonical(FloatParts128 *p,
1812                                               float_status *s)
1813 {
1814     const FloatFmt *fmt = &floatx80_params[s->floatx80_rounding_precision];
1815     uint64_t frac;
1816     int exp;
1817 
1818     switch (p->cls) {
1819     case float_class_normal:
1820         if (s->floatx80_rounding_precision == floatx80_precision_x) {
1821             parts_uncanon_normal(p, s, fmt);
1822             frac = p->frac_hi;
1823             exp = p->exp;
1824         } else {
1825             FloatParts64 p64;
1826 
1827             p64.sign = p->sign;
1828             p64.exp = p->exp;
1829             frac_truncjam(&p64, p);
1830             parts_uncanon_normal(&p64, s, fmt);
1831             frac = p64.frac;
1832             exp = p64.exp;
1833         }
1834         if (exp != fmt->exp_max) {
1835             break;
1836         }
1837         /* rounded to inf -- fall through to set frac correctly */
1838 
1839     case float_class_inf:
1840         /* x86 and m68k differ in the setting of the integer bit. */
1841         frac = floatx80_infinity_low;
1842         exp = fmt->exp_max;
1843         break;
1844 
1845     case float_class_zero:
1846         frac = 0;
1847         exp = 0;
1848         break;
1849 
1850     case float_class_snan:
1851     case float_class_qnan:
1852         /* NaNs have the integer bit set. */
1853         frac = p->frac_hi | (1ull << 63);
1854         exp = fmt->exp_max;
1855         break;
1856 
1857     default:
1858         g_assert_not_reached();
1859     }
1860 
1861     return packFloatx80(p->sign, exp, frac);
1862 }
1863 
1864 /*
1865  * Addition and subtraction
1866  */
1867 
1868 static float16 QEMU_FLATTEN
1869 float16_addsub(float16 a, float16 b, float_status *status, bool subtract)
1870 {
1871     FloatParts64 pa, pb, *pr;
1872 
1873     float16_unpack_canonical(&pa, a, status);
1874     float16_unpack_canonical(&pb, b, status);
1875     pr = parts_addsub(&pa, &pb, status, subtract);
1876 
1877     return float16_round_pack_canonical(pr, status);
1878 }
1879 
1880 float16 float16_add(float16 a, float16 b, float_status *status)
1881 {
1882     return float16_addsub(a, b, status, false);
1883 }
1884 
1885 float16 float16_sub(float16 a, float16 b, float_status *status)
1886 {
1887     return float16_addsub(a, b, status, true);
1888 }
1889 
1890 static float32 QEMU_SOFTFLOAT_ATTR
1891 soft_f32_addsub(float32 a, float32 b, float_status *status, bool subtract)
1892 {
1893     FloatParts64 pa, pb, *pr;
1894 
1895     float32_unpack_canonical(&pa, a, status);
1896     float32_unpack_canonical(&pb, b, status);
1897     pr = parts_addsub(&pa, &pb, status, subtract);
1898 
1899     return float32_round_pack_canonical(pr, status);
1900 }
1901 
1902 static float32 soft_f32_add(float32 a, float32 b, float_status *status)
1903 {
1904     return soft_f32_addsub(a, b, status, false);
1905 }
1906 
1907 static float32 soft_f32_sub(float32 a, float32 b, float_status *status)
1908 {
1909     return soft_f32_addsub(a, b, status, true);
1910 }
1911 
1912 static float64 QEMU_SOFTFLOAT_ATTR
1913 soft_f64_addsub(float64 a, float64 b, float_status *status, bool subtract)
1914 {
1915     FloatParts64 pa, pb, *pr;
1916 
1917     float64_unpack_canonical(&pa, a, status);
1918     float64_unpack_canonical(&pb, b, status);
1919     pr = parts_addsub(&pa, &pb, status, subtract);
1920 
1921     return float64_round_pack_canonical(pr, status);
1922 }
1923 
1924 static float64 soft_f64_add(float64 a, float64 b, float_status *status)
1925 {
1926     return soft_f64_addsub(a, b, status, false);
1927 }
1928 
1929 static float64 soft_f64_sub(float64 a, float64 b, float_status *status)
1930 {
1931     return soft_f64_addsub(a, b, status, true);
1932 }
1933 
1934 static float hard_f32_add(float a, float b)
1935 {
1936     return a + b;
1937 }
1938 
1939 static float hard_f32_sub(float a, float b)
1940 {
1941     return a - b;
1942 }
1943 
1944 static double hard_f64_add(double a, double b)
1945 {
1946     return a + b;
1947 }
1948 
1949 static double hard_f64_sub(double a, double b)
1950 {
1951     return a - b;
1952 }
1953 
1954 static bool f32_addsubmul_post(union_float32 a, union_float32 b)
1955 {
1956     if (QEMU_HARDFLOAT_2F32_USE_FP) {
1957         return !(fpclassify(a.h) == FP_ZERO && fpclassify(b.h) == FP_ZERO);
1958     }
1959     return !(float32_is_zero(a.s) && float32_is_zero(b.s));
1960 }
1961 
1962 static bool f64_addsubmul_post(union_float64 a, union_float64 b)
1963 {
1964     if (QEMU_HARDFLOAT_2F64_USE_FP) {
1965         return !(fpclassify(a.h) == FP_ZERO && fpclassify(b.h) == FP_ZERO);
1966     } else {
1967         return !(float64_is_zero(a.s) && float64_is_zero(b.s));
1968     }
1969 }
1970 
1971 static float32 float32_addsub(float32 a, float32 b, float_status *s,
1972                               hard_f32_op2_fn hard, soft_f32_op2_fn soft)
1973 {
1974     return float32_gen2(a, b, s, hard, soft,
1975                         f32_is_zon2, f32_addsubmul_post);
1976 }
1977 
1978 static float64 float64_addsub(float64 a, float64 b, float_status *s,
1979                               hard_f64_op2_fn hard, soft_f64_op2_fn soft)
1980 {
1981     return float64_gen2(a, b, s, hard, soft,
1982                         f64_is_zon2, f64_addsubmul_post);
1983 }
1984 
1985 float32 QEMU_FLATTEN
1986 float32_add(float32 a, float32 b, float_status *s)
1987 {
1988     return float32_addsub(a, b, s, hard_f32_add, soft_f32_add);
1989 }
1990 
1991 float32 QEMU_FLATTEN
1992 float32_sub(float32 a, float32 b, float_status *s)
1993 {
1994     return float32_addsub(a, b, s, hard_f32_sub, soft_f32_sub);
1995 }
1996 
1997 float64 QEMU_FLATTEN
1998 float64_add(float64 a, float64 b, float_status *s)
1999 {
2000     return float64_addsub(a, b, s, hard_f64_add, soft_f64_add);
2001 }
2002 
2003 float64 QEMU_FLATTEN
2004 float64_sub(float64 a, float64 b, float_status *s)
2005 {
2006     return float64_addsub(a, b, s, hard_f64_sub, soft_f64_sub);
2007 }
2008 
2009 static float64 float64r32_addsub(float64 a, float64 b, float_status *status,
2010                                  bool subtract)
2011 {
2012     FloatParts64 pa, pb, *pr;
2013 
2014     float64_unpack_canonical(&pa, a, status);
2015     float64_unpack_canonical(&pb, b, status);
2016     pr = parts_addsub(&pa, &pb, status, subtract);
2017 
2018     return float64r32_round_pack_canonical(pr, status);
2019 }
2020 
2021 float64 float64r32_add(float64 a, float64 b, float_status *status)
2022 {
2023     return float64r32_addsub(a, b, status, false);
2024 }
2025 
2026 float64 float64r32_sub(float64 a, float64 b, float_status *status)
2027 {
2028     return float64r32_addsub(a, b, status, true);
2029 }
2030 
2031 static bfloat16 QEMU_FLATTEN
2032 bfloat16_addsub(bfloat16 a, bfloat16 b, float_status *status, bool subtract)
2033 {
2034     FloatParts64 pa, pb, *pr;
2035 
2036     bfloat16_unpack_canonical(&pa, a, status);
2037     bfloat16_unpack_canonical(&pb, b, status);
2038     pr = parts_addsub(&pa, &pb, status, subtract);
2039 
2040     return bfloat16_round_pack_canonical(pr, status);
2041 }
2042 
2043 bfloat16 bfloat16_add(bfloat16 a, bfloat16 b, float_status *status)
2044 {
2045     return bfloat16_addsub(a, b, status, false);
2046 }
2047 
2048 bfloat16 bfloat16_sub(bfloat16 a, bfloat16 b, float_status *status)
2049 {
2050     return bfloat16_addsub(a, b, status, true);
2051 }
2052 
2053 static float128 QEMU_FLATTEN
2054 float128_addsub(float128 a, float128 b, float_status *status, bool subtract)
2055 {
2056     FloatParts128 pa, pb, *pr;
2057 
2058     float128_unpack_canonical(&pa, a, status);
2059     float128_unpack_canonical(&pb, b, status);
2060     pr = parts_addsub(&pa, &pb, status, subtract);
2061 
2062     return float128_round_pack_canonical(pr, status);
2063 }
2064 
2065 float128 float128_add(float128 a, float128 b, float_status *status)
2066 {
2067     return float128_addsub(a, b, status, false);
2068 }
2069 
2070 float128 float128_sub(float128 a, float128 b, float_status *status)
2071 {
2072     return float128_addsub(a, b, status, true);
2073 }
2074 
2075 static floatx80 QEMU_FLATTEN
2076 floatx80_addsub(floatx80 a, floatx80 b, float_status *status, bool subtract)
2077 {
2078     FloatParts128 pa, pb, *pr;
2079 
2080     if (!floatx80_unpack_canonical(&pa, a, status) ||
2081         !floatx80_unpack_canonical(&pb, b, status)) {
2082         return floatx80_default_nan(status);
2083     }
2084 
2085     pr = parts_addsub(&pa, &pb, status, subtract);
2086     return floatx80_round_pack_canonical(pr, status);
2087 }
2088 
2089 floatx80 floatx80_add(floatx80 a, floatx80 b, float_status *status)
2090 {
2091     return floatx80_addsub(a, b, status, false);
2092 }
2093 
2094 floatx80 floatx80_sub(floatx80 a, floatx80 b, float_status *status)
2095 {
2096     return floatx80_addsub(a, b, status, true);
2097 }
2098 
2099 /*
2100  * Multiplication
2101  */
2102 
2103 float16 QEMU_FLATTEN float16_mul(float16 a, float16 b, float_status *status)
2104 {
2105     FloatParts64 pa, pb, *pr;
2106 
2107     float16_unpack_canonical(&pa, a, status);
2108     float16_unpack_canonical(&pb, b, status);
2109     pr = parts_mul(&pa, &pb, status);
2110 
2111     return float16_round_pack_canonical(pr, status);
2112 }
2113 
2114 static float32 QEMU_SOFTFLOAT_ATTR
2115 soft_f32_mul(float32 a, float32 b, float_status *status)
2116 {
2117     FloatParts64 pa, pb, *pr;
2118 
2119     float32_unpack_canonical(&pa, a, status);
2120     float32_unpack_canonical(&pb, b, status);
2121     pr = parts_mul(&pa, &pb, status);
2122 
2123     return float32_round_pack_canonical(pr, status);
2124 }
2125 
2126 static float64 QEMU_SOFTFLOAT_ATTR
2127 soft_f64_mul(float64 a, float64 b, float_status *status)
2128 {
2129     FloatParts64 pa, pb, *pr;
2130 
2131     float64_unpack_canonical(&pa, a, status);
2132     float64_unpack_canonical(&pb, b, status);
2133     pr = parts_mul(&pa, &pb, status);
2134 
2135     return float64_round_pack_canonical(pr, status);
2136 }
2137 
2138 static float hard_f32_mul(float a, float b)
2139 {
2140     return a * b;
2141 }
2142 
2143 static double hard_f64_mul(double a, double b)
2144 {
2145     return a * b;
2146 }
2147 
2148 float32 QEMU_FLATTEN
2149 float32_mul(float32 a, float32 b, float_status *s)
2150 {
2151     return float32_gen2(a, b, s, hard_f32_mul, soft_f32_mul,
2152                         f32_is_zon2, f32_addsubmul_post);
2153 }
2154 
2155 float64 QEMU_FLATTEN
2156 float64_mul(float64 a, float64 b, float_status *s)
2157 {
2158     return float64_gen2(a, b, s, hard_f64_mul, soft_f64_mul,
2159                         f64_is_zon2, f64_addsubmul_post);
2160 }
2161 
2162 float64 float64r32_mul(float64 a, float64 b, float_status *status)
2163 {
2164     FloatParts64 pa, pb, *pr;
2165 
2166     float64_unpack_canonical(&pa, a, status);
2167     float64_unpack_canonical(&pb, b, status);
2168     pr = parts_mul(&pa, &pb, status);
2169 
2170     return float64r32_round_pack_canonical(pr, status);
2171 }
2172 
2173 bfloat16 QEMU_FLATTEN
2174 bfloat16_mul(bfloat16 a, bfloat16 b, float_status *status)
2175 {
2176     FloatParts64 pa, pb, *pr;
2177 
2178     bfloat16_unpack_canonical(&pa, a, status);
2179     bfloat16_unpack_canonical(&pb, b, status);
2180     pr = parts_mul(&pa, &pb, status);
2181 
2182     return bfloat16_round_pack_canonical(pr, status);
2183 }
2184 
2185 float128 QEMU_FLATTEN
2186 float128_mul(float128 a, float128 b, float_status *status)
2187 {
2188     FloatParts128 pa, pb, *pr;
2189 
2190     float128_unpack_canonical(&pa, a, status);
2191     float128_unpack_canonical(&pb, b, status);
2192     pr = parts_mul(&pa, &pb, status);
2193 
2194     return float128_round_pack_canonical(pr, status);
2195 }
2196 
2197 floatx80 QEMU_FLATTEN
2198 floatx80_mul(floatx80 a, floatx80 b, float_status *status)
2199 {
2200     FloatParts128 pa, pb, *pr;
2201 
2202     if (!floatx80_unpack_canonical(&pa, a, status) ||
2203         !floatx80_unpack_canonical(&pb, b, status)) {
2204         return floatx80_default_nan(status);
2205     }
2206 
2207     pr = parts_mul(&pa, &pb, status);
2208     return floatx80_round_pack_canonical(pr, status);
2209 }
2210 
2211 /*
2212  * Fused multiply-add
2213  */
2214 
2215 float16 QEMU_FLATTEN float16_muladd(float16 a, float16 b, float16 c,
2216                                     int flags, float_status *status)
2217 {
2218     FloatParts64 pa, pb, pc, *pr;
2219 
2220     float16_unpack_canonical(&pa, a, status);
2221     float16_unpack_canonical(&pb, b, status);
2222     float16_unpack_canonical(&pc, c, status);
2223     pr = parts_muladd(&pa, &pb, &pc, flags, status);
2224 
2225     return float16_round_pack_canonical(pr, status);
2226 }
2227 
2228 static float32 QEMU_SOFTFLOAT_ATTR
2229 soft_f32_muladd(float32 a, float32 b, float32 c, int flags,
2230                 float_status *status)
2231 {
2232     FloatParts64 pa, pb, pc, *pr;
2233 
2234     float32_unpack_canonical(&pa, a, status);
2235     float32_unpack_canonical(&pb, b, status);
2236     float32_unpack_canonical(&pc, c, status);
2237     pr = parts_muladd(&pa, &pb, &pc, flags, status);
2238 
2239     return float32_round_pack_canonical(pr, status);
2240 }
2241 
2242 static float64 QEMU_SOFTFLOAT_ATTR
2243 soft_f64_muladd(float64 a, float64 b, float64 c, int flags,
2244                 float_status *status)
2245 {
2246     FloatParts64 pa, pb, pc, *pr;
2247 
2248     float64_unpack_canonical(&pa, a, status);
2249     float64_unpack_canonical(&pb, b, status);
2250     float64_unpack_canonical(&pc, c, status);
2251     pr = parts_muladd(&pa, &pb, &pc, flags, status);
2252 
2253     return float64_round_pack_canonical(pr, status);
2254 }
2255 
2256 static bool force_soft_fma;
2257 
2258 float32 QEMU_FLATTEN
2259 float32_muladd(float32 xa, float32 xb, float32 xc, int flags, float_status *s)
2260 {
2261     union_float32 ua, ub, uc, ur;
2262 
2263     ua.s = xa;
2264     ub.s = xb;
2265     uc.s = xc;
2266 
2267     if (unlikely(!can_use_fpu(s))) {
2268         goto soft;
2269     }
2270     if (unlikely(flags & float_muladd_halve_result)) {
2271         goto soft;
2272     }
2273 
2274     float32_input_flush3(&ua.s, &ub.s, &uc.s, s);
2275     if (unlikely(!f32_is_zon3(ua, ub, uc))) {
2276         goto soft;
2277     }
2278 
2279     if (unlikely(force_soft_fma)) {
2280         goto soft;
2281     }
2282 
2283     /*
2284      * When (a || b) == 0, there's no need to check for under/over flow,
2285      * since we know the addend is (normal || 0) and the product is 0.
2286      */
2287     if (float32_is_zero(ua.s) || float32_is_zero(ub.s)) {
2288         union_float32 up;
2289         bool prod_sign;
2290 
2291         prod_sign = float32_is_neg(ua.s) ^ float32_is_neg(ub.s);
2292         prod_sign ^= !!(flags & float_muladd_negate_product);
2293         up.s = float32_set_sign(float32_zero, prod_sign);
2294 
2295         if (flags & float_muladd_negate_c) {
2296             uc.h = -uc.h;
2297         }
2298         ur.h = up.h + uc.h;
2299     } else {
2300         union_float32 ua_orig = ua;
2301         union_float32 uc_orig = uc;
2302 
2303         if (flags & float_muladd_negate_product) {
2304             ua.h = -ua.h;
2305         }
2306         if (flags & float_muladd_negate_c) {
2307             uc.h = -uc.h;
2308         }
2309 
2310         ur.h = fmaf(ua.h, ub.h, uc.h);
2311 
2312         if (unlikely(f32_is_inf(ur))) {
2313             float_raise(float_flag_overflow, s);
2314         } else if (unlikely(fabsf(ur.h) <= FLT_MIN)) {
2315             ua = ua_orig;
2316             uc = uc_orig;
2317             goto soft;
2318         }
2319     }
2320     if (flags & float_muladd_negate_result) {
2321         return float32_chs(ur.s);
2322     }
2323     return ur.s;
2324 
2325  soft:
2326     return soft_f32_muladd(ua.s, ub.s, uc.s, flags, s);
2327 }
2328 
2329 float64 QEMU_FLATTEN
2330 float64_muladd(float64 xa, float64 xb, float64 xc, int flags, float_status *s)
2331 {
2332     union_float64 ua, ub, uc, ur;
2333 
2334     ua.s = xa;
2335     ub.s = xb;
2336     uc.s = xc;
2337 
2338     if (unlikely(!can_use_fpu(s))) {
2339         goto soft;
2340     }
2341     if (unlikely(flags & float_muladd_halve_result)) {
2342         goto soft;
2343     }
2344 
2345     float64_input_flush3(&ua.s, &ub.s, &uc.s, s);
2346     if (unlikely(!f64_is_zon3(ua, ub, uc))) {
2347         goto soft;
2348     }
2349 
2350     if (unlikely(force_soft_fma)) {
2351         goto soft;
2352     }
2353 
2354     /*
2355      * When (a || b) == 0, there's no need to check for under/over flow,
2356      * since we know the addend is (normal || 0) and the product is 0.
2357      */
2358     if (float64_is_zero(ua.s) || float64_is_zero(ub.s)) {
2359         union_float64 up;
2360         bool prod_sign;
2361 
2362         prod_sign = float64_is_neg(ua.s) ^ float64_is_neg(ub.s);
2363         prod_sign ^= !!(flags & float_muladd_negate_product);
2364         up.s = float64_set_sign(float64_zero, prod_sign);
2365 
2366         if (flags & float_muladd_negate_c) {
2367             uc.h = -uc.h;
2368         }
2369         ur.h = up.h + uc.h;
2370     } else {
2371         union_float64 ua_orig = ua;
2372         union_float64 uc_orig = uc;
2373 
2374         if (flags & float_muladd_negate_product) {
2375             ua.h = -ua.h;
2376         }
2377         if (flags & float_muladd_negate_c) {
2378             uc.h = -uc.h;
2379         }
2380 
2381         ur.h = fma(ua.h, ub.h, uc.h);
2382 
2383         if (unlikely(f64_is_inf(ur))) {
2384             float_raise(float_flag_overflow, s);
2385         } else if (unlikely(fabs(ur.h) <= FLT_MIN)) {
2386             ua = ua_orig;
2387             uc = uc_orig;
2388             goto soft;
2389         }
2390     }
2391     if (flags & float_muladd_negate_result) {
2392         return float64_chs(ur.s);
2393     }
2394     return ur.s;
2395 
2396  soft:
2397     return soft_f64_muladd(ua.s, ub.s, uc.s, flags, s);
2398 }
2399 
2400 float64 float64r32_muladd(float64 a, float64 b, float64 c,
2401                           int flags, float_status *status)
2402 {
2403     FloatParts64 pa, pb, pc, *pr;
2404 
2405     float64_unpack_canonical(&pa, a, status);
2406     float64_unpack_canonical(&pb, b, status);
2407     float64_unpack_canonical(&pc, c, status);
2408     pr = parts_muladd(&pa, &pb, &pc, flags, status);
2409 
2410     return float64r32_round_pack_canonical(pr, status);
2411 }
2412 
2413 bfloat16 QEMU_FLATTEN bfloat16_muladd(bfloat16 a, bfloat16 b, bfloat16 c,
2414                                       int flags, float_status *status)
2415 {
2416     FloatParts64 pa, pb, pc, *pr;
2417 
2418     bfloat16_unpack_canonical(&pa, a, status);
2419     bfloat16_unpack_canonical(&pb, b, status);
2420     bfloat16_unpack_canonical(&pc, c, status);
2421     pr = parts_muladd(&pa, &pb, &pc, flags, status);
2422 
2423     return bfloat16_round_pack_canonical(pr, status);
2424 }
2425 
2426 float128 QEMU_FLATTEN float128_muladd(float128 a, float128 b, float128 c,
2427                                       int flags, float_status *status)
2428 {
2429     FloatParts128 pa, pb, pc, *pr;
2430 
2431     float128_unpack_canonical(&pa, a, status);
2432     float128_unpack_canonical(&pb, b, status);
2433     float128_unpack_canonical(&pc, c, status);
2434     pr = parts_muladd(&pa, &pb, &pc, flags, status);
2435 
2436     return float128_round_pack_canonical(pr, status);
2437 }
2438 
2439 /*
2440  * Division
2441  */
2442 
2443 float16 float16_div(float16 a, float16 b, float_status *status)
2444 {
2445     FloatParts64 pa, pb, *pr;
2446 
2447     float16_unpack_canonical(&pa, a, status);
2448     float16_unpack_canonical(&pb, b, status);
2449     pr = parts_div(&pa, &pb, status);
2450 
2451     return float16_round_pack_canonical(pr, status);
2452 }
2453 
2454 static float32 QEMU_SOFTFLOAT_ATTR
2455 soft_f32_div(float32 a, float32 b, float_status *status)
2456 {
2457     FloatParts64 pa, pb, *pr;
2458 
2459     float32_unpack_canonical(&pa, a, status);
2460     float32_unpack_canonical(&pb, b, status);
2461     pr = parts_div(&pa, &pb, status);
2462 
2463     return float32_round_pack_canonical(pr, status);
2464 }
2465 
2466 static float64 QEMU_SOFTFLOAT_ATTR
2467 soft_f64_div(float64 a, float64 b, float_status *status)
2468 {
2469     FloatParts64 pa, pb, *pr;
2470 
2471     float64_unpack_canonical(&pa, a, status);
2472     float64_unpack_canonical(&pb, b, status);
2473     pr = parts_div(&pa, &pb, status);
2474 
2475     return float64_round_pack_canonical(pr, status);
2476 }
2477 
2478 static float hard_f32_div(float a, float b)
2479 {
2480     return a / b;
2481 }
2482 
2483 static double hard_f64_div(double a, double b)
2484 {
2485     return a / b;
2486 }
2487 
2488 static bool f32_div_pre(union_float32 a, union_float32 b)
2489 {
2490     if (QEMU_HARDFLOAT_2F32_USE_FP) {
2491         return (fpclassify(a.h) == FP_NORMAL || fpclassify(a.h) == FP_ZERO) &&
2492                fpclassify(b.h) == FP_NORMAL;
2493     }
2494     return float32_is_zero_or_normal(a.s) && float32_is_normal(b.s);
2495 }
2496 
2497 static bool f64_div_pre(union_float64 a, union_float64 b)
2498 {
2499     if (QEMU_HARDFLOAT_2F64_USE_FP) {
2500         return (fpclassify(a.h) == FP_NORMAL || fpclassify(a.h) == FP_ZERO) &&
2501                fpclassify(b.h) == FP_NORMAL;
2502     }
2503     return float64_is_zero_or_normal(a.s) && float64_is_normal(b.s);
2504 }
2505 
2506 static bool f32_div_post(union_float32 a, union_float32 b)
2507 {
2508     if (QEMU_HARDFLOAT_2F32_USE_FP) {
2509         return fpclassify(a.h) != FP_ZERO;
2510     }
2511     return !float32_is_zero(a.s);
2512 }
2513 
2514 static bool f64_div_post(union_float64 a, union_float64 b)
2515 {
2516     if (QEMU_HARDFLOAT_2F64_USE_FP) {
2517         return fpclassify(a.h) != FP_ZERO;
2518     }
2519     return !float64_is_zero(a.s);
2520 }
2521 
2522 float32 QEMU_FLATTEN
2523 float32_div(float32 a, float32 b, float_status *s)
2524 {
2525     return float32_gen2(a, b, s, hard_f32_div, soft_f32_div,
2526                         f32_div_pre, f32_div_post);
2527 }
2528 
2529 float64 QEMU_FLATTEN
2530 float64_div(float64 a, float64 b, float_status *s)
2531 {
2532     return float64_gen2(a, b, s, hard_f64_div, soft_f64_div,
2533                         f64_div_pre, f64_div_post);
2534 }
2535 
2536 float64 float64r32_div(float64 a, float64 b, float_status *status)
2537 {
2538     FloatParts64 pa, pb, *pr;
2539 
2540     float64_unpack_canonical(&pa, a, status);
2541     float64_unpack_canonical(&pb, b, status);
2542     pr = parts_div(&pa, &pb, status);
2543 
2544     return float64r32_round_pack_canonical(pr, status);
2545 }
2546 
2547 bfloat16 QEMU_FLATTEN
2548 bfloat16_div(bfloat16 a, bfloat16 b, float_status *status)
2549 {
2550     FloatParts64 pa, pb, *pr;
2551 
2552     bfloat16_unpack_canonical(&pa, a, status);
2553     bfloat16_unpack_canonical(&pb, b, status);
2554     pr = parts_div(&pa, &pb, status);
2555 
2556     return bfloat16_round_pack_canonical(pr, status);
2557 }
2558 
2559 float128 QEMU_FLATTEN
2560 float128_div(float128 a, float128 b, float_status *status)
2561 {
2562     FloatParts128 pa, pb, *pr;
2563 
2564     float128_unpack_canonical(&pa, a, status);
2565     float128_unpack_canonical(&pb, b, status);
2566     pr = parts_div(&pa, &pb, status);
2567 
2568     return float128_round_pack_canonical(pr, status);
2569 }
2570 
2571 floatx80 floatx80_div(floatx80 a, floatx80 b, float_status *status)
2572 {
2573     FloatParts128 pa, pb, *pr;
2574 
2575     if (!floatx80_unpack_canonical(&pa, a, status) ||
2576         !floatx80_unpack_canonical(&pb, b, status)) {
2577         return floatx80_default_nan(status);
2578     }
2579 
2580     pr = parts_div(&pa, &pb, status);
2581     return floatx80_round_pack_canonical(pr, status);
2582 }
2583 
2584 /*
2585  * Remainder
2586  */
2587 
2588 float32 float32_rem(float32 a, float32 b, float_status *status)
2589 {
2590     FloatParts64 pa, pb, *pr;
2591 
2592     float32_unpack_canonical(&pa, a, status);
2593     float32_unpack_canonical(&pb, b, status);
2594     pr = parts_modrem(&pa, &pb, NULL, status);
2595 
2596     return float32_round_pack_canonical(pr, status);
2597 }
2598 
2599 float64 float64_rem(float64 a, float64 b, float_status *status)
2600 {
2601     FloatParts64 pa, pb, *pr;
2602 
2603     float64_unpack_canonical(&pa, a, status);
2604     float64_unpack_canonical(&pb, b, status);
2605     pr = parts_modrem(&pa, &pb, NULL, status);
2606 
2607     return float64_round_pack_canonical(pr, status);
2608 }
2609 
2610 float128 float128_rem(float128 a, float128 b, float_status *status)
2611 {
2612     FloatParts128 pa, pb, *pr;
2613 
2614     float128_unpack_canonical(&pa, a, status);
2615     float128_unpack_canonical(&pb, b, status);
2616     pr = parts_modrem(&pa, &pb, NULL, status);
2617 
2618     return float128_round_pack_canonical(pr, status);
2619 }
2620 
2621 /*
2622  * Returns the remainder of the extended double-precision floating-point value
2623  * `a' with respect to the corresponding value `b'.
2624  * If 'mod' is false, the operation is performed according to the IEC/IEEE
2625  * Standard for Binary Floating-Point Arithmetic.  If 'mod' is true, return
2626  * the remainder based on truncating the quotient toward zero instead and
2627  * *quotient is set to the low 64 bits of the absolute value of the integer
2628  * quotient.
2629  */
2630 floatx80 floatx80_modrem(floatx80 a, floatx80 b, bool mod,
2631                          uint64_t *quotient, float_status *status)
2632 {
2633     FloatParts128 pa, pb, *pr;
2634 
2635     *quotient = 0;
2636     if (!floatx80_unpack_canonical(&pa, a, status) ||
2637         !floatx80_unpack_canonical(&pb, b, status)) {
2638         return floatx80_default_nan(status);
2639     }
2640     pr = parts_modrem(&pa, &pb, mod ? quotient : NULL, status);
2641 
2642     return floatx80_round_pack_canonical(pr, status);
2643 }
2644 
2645 floatx80 floatx80_rem(floatx80 a, floatx80 b, float_status *status)
2646 {
2647     uint64_t quotient;
2648     return floatx80_modrem(a, b, false, &quotient, status);
2649 }
2650 
2651 floatx80 floatx80_mod(floatx80 a, floatx80 b, float_status *status)
2652 {
2653     uint64_t quotient;
2654     return floatx80_modrem(a, b, true, &quotient, status);
2655 }
2656 
2657 /*
2658  * Float to Float conversions
2659  *
2660  * Returns the result of converting one float format to another. The
2661  * conversion is performed according to the IEC/IEEE Standard for
2662  * Binary Floating-Point Arithmetic.
2663  *
2664  * Usually this only needs to take care of raising invalid exceptions
2665  * and handling the conversion on NaNs.
2666  */
2667 
2668 static void parts_float_to_ahp(FloatParts64 *a, float_status *s)
2669 {
2670     switch (a->cls) {
2671     case float_class_snan:
2672         float_raise(float_flag_invalid_snan, s);
2673         /* fall through */
2674     case float_class_qnan:
2675         /*
2676          * There is no NaN in the destination format.  Raise Invalid
2677          * and return a zero with the sign of the input NaN.
2678          */
2679         float_raise(float_flag_invalid, s);
2680         a->cls = float_class_zero;
2681         break;
2682 
2683     case float_class_inf:
2684         /*
2685          * There is no Inf in the destination format.  Raise Invalid
2686          * and return the maximum normal with the correct sign.
2687          */
2688         float_raise(float_flag_invalid, s);
2689         a->cls = float_class_normal;
2690         a->exp = float16_params_ahp.exp_max;
2691         a->frac = MAKE_64BIT_MASK(float16_params_ahp.frac_shift,
2692                                   float16_params_ahp.frac_size + 1);
2693         break;
2694 
2695     case float_class_normal:
2696     case float_class_zero:
2697         break;
2698 
2699     default:
2700         g_assert_not_reached();
2701     }
2702 }
2703 
2704 static void parts64_float_to_float(FloatParts64 *a, float_status *s)
2705 {
2706     if (is_nan(a->cls)) {
2707         parts_return_nan(a, s);
2708     }
2709 }
2710 
2711 static void parts128_float_to_float(FloatParts128 *a, float_status *s)
2712 {
2713     if (is_nan(a->cls)) {
2714         parts_return_nan(a, s);
2715     }
2716 }
2717 
2718 #define parts_float_to_float(P, S) \
2719     PARTS_GENERIC_64_128(float_to_float, P)(P, S)
2720 
2721 static void parts_float_to_float_narrow(FloatParts64 *a, FloatParts128 *b,
2722                                         float_status *s)
2723 {
2724     a->cls = b->cls;
2725     a->sign = b->sign;
2726     a->exp = b->exp;
2727 
2728     if (a->cls == float_class_normal) {
2729         frac_truncjam(a, b);
2730     } else if (is_nan(a->cls)) {
2731         /* Discard the low bits of the NaN. */
2732         a->frac = b->frac_hi;
2733         parts_return_nan(a, s);
2734     }
2735 }
2736 
2737 static void parts_float_to_float_widen(FloatParts128 *a, FloatParts64 *b,
2738                                        float_status *s)
2739 {
2740     a->cls = b->cls;
2741     a->sign = b->sign;
2742     a->exp = b->exp;
2743     frac_widen(a, b);
2744 
2745     if (is_nan(a->cls)) {
2746         parts_return_nan(a, s);
2747     }
2748 }
2749 
2750 float32 float16_to_float32(float16 a, bool ieee, float_status *s)
2751 {
2752     const FloatFmt *fmt16 = ieee ? &float16_params : &float16_params_ahp;
2753     FloatParts64 p;
2754 
2755     float16a_unpack_canonical(&p, a, s, fmt16);
2756     parts_float_to_float(&p, s);
2757     return float32_round_pack_canonical(&p, s);
2758 }
2759 
2760 float64 float16_to_float64(float16 a, bool ieee, float_status *s)
2761 {
2762     const FloatFmt *fmt16 = ieee ? &float16_params : &float16_params_ahp;
2763     FloatParts64 p;
2764 
2765     float16a_unpack_canonical(&p, a, s, fmt16);
2766     parts_float_to_float(&p, s);
2767     return float64_round_pack_canonical(&p, s);
2768 }
2769 
2770 float16 float32_to_float16(float32 a, bool ieee, float_status *s)
2771 {
2772     FloatParts64 p;
2773     const FloatFmt *fmt;
2774 
2775     float32_unpack_canonical(&p, a, s);
2776     if (ieee) {
2777         parts_float_to_float(&p, s);
2778         fmt = &float16_params;
2779     } else {
2780         parts_float_to_ahp(&p, s);
2781         fmt = &float16_params_ahp;
2782     }
2783     return float16a_round_pack_canonical(&p, s, fmt);
2784 }
2785 
2786 static float64 QEMU_SOFTFLOAT_ATTR
2787 soft_float32_to_float64(float32 a, float_status *s)
2788 {
2789     FloatParts64 p;
2790 
2791     float32_unpack_canonical(&p, a, s);
2792     parts_float_to_float(&p, s);
2793     return float64_round_pack_canonical(&p, s);
2794 }
2795 
2796 float64 float32_to_float64(float32 a, float_status *s)
2797 {
2798     if (likely(float32_is_normal(a))) {
2799         /* Widening conversion can never produce inexact results.  */
2800         union_float32 uf;
2801         union_float64 ud;
2802         uf.s = a;
2803         ud.h = uf.h;
2804         return ud.s;
2805     } else if (float32_is_zero(a)) {
2806         return float64_set_sign(float64_zero, float32_is_neg(a));
2807     } else {
2808         return soft_float32_to_float64(a, s);
2809     }
2810 }
2811 
2812 float16 float64_to_float16(float64 a, bool ieee, float_status *s)
2813 {
2814     FloatParts64 p;
2815     const FloatFmt *fmt;
2816 
2817     float64_unpack_canonical(&p, a, s);
2818     if (ieee) {
2819         parts_float_to_float(&p, s);
2820         fmt = &float16_params;
2821     } else {
2822         parts_float_to_ahp(&p, s);
2823         fmt = &float16_params_ahp;
2824     }
2825     return float16a_round_pack_canonical(&p, s, fmt);
2826 }
2827 
2828 float32 float64_to_float32(float64 a, float_status *s)
2829 {
2830     FloatParts64 p;
2831 
2832     float64_unpack_canonical(&p, a, s);
2833     parts_float_to_float(&p, s);
2834     return float32_round_pack_canonical(&p, s);
2835 }
2836 
2837 float32 bfloat16_to_float32(bfloat16 a, float_status *s)
2838 {
2839     FloatParts64 p;
2840 
2841     bfloat16_unpack_canonical(&p, a, s);
2842     parts_float_to_float(&p, s);
2843     return float32_round_pack_canonical(&p, s);
2844 }
2845 
2846 float64 bfloat16_to_float64(bfloat16 a, float_status *s)
2847 {
2848     FloatParts64 p;
2849 
2850     bfloat16_unpack_canonical(&p, a, s);
2851     parts_float_to_float(&p, s);
2852     return float64_round_pack_canonical(&p, s);
2853 }
2854 
2855 bfloat16 float32_to_bfloat16(float32 a, float_status *s)
2856 {
2857     FloatParts64 p;
2858 
2859     float32_unpack_canonical(&p, a, s);
2860     parts_float_to_float(&p, s);
2861     return bfloat16_round_pack_canonical(&p, s);
2862 }
2863 
2864 bfloat16 float64_to_bfloat16(float64 a, float_status *s)
2865 {
2866     FloatParts64 p;
2867 
2868     float64_unpack_canonical(&p, a, s);
2869     parts_float_to_float(&p, s);
2870     return bfloat16_round_pack_canonical(&p, s);
2871 }
2872 
2873 float32 float128_to_float32(float128 a, float_status *s)
2874 {
2875     FloatParts64 p64;
2876     FloatParts128 p128;
2877 
2878     float128_unpack_canonical(&p128, a, s);
2879     parts_float_to_float_narrow(&p64, &p128, s);
2880     return float32_round_pack_canonical(&p64, s);
2881 }
2882 
2883 float64 float128_to_float64(float128 a, float_status *s)
2884 {
2885     FloatParts64 p64;
2886     FloatParts128 p128;
2887 
2888     float128_unpack_canonical(&p128, a, s);
2889     parts_float_to_float_narrow(&p64, &p128, s);
2890     return float64_round_pack_canonical(&p64, s);
2891 }
2892 
2893 float128 float32_to_float128(float32 a, float_status *s)
2894 {
2895     FloatParts64 p64;
2896     FloatParts128 p128;
2897 
2898     float32_unpack_canonical(&p64, a, s);
2899     parts_float_to_float_widen(&p128, &p64, s);
2900     return float128_round_pack_canonical(&p128, s);
2901 }
2902 
2903 float128 float64_to_float128(float64 a, float_status *s)
2904 {
2905     FloatParts64 p64;
2906     FloatParts128 p128;
2907 
2908     float64_unpack_canonical(&p64, a, s);
2909     parts_float_to_float_widen(&p128, &p64, s);
2910     return float128_round_pack_canonical(&p128, s);
2911 }
2912 
2913 float32 floatx80_to_float32(floatx80 a, float_status *s)
2914 {
2915     FloatParts64 p64;
2916     FloatParts128 p128;
2917 
2918     if (floatx80_unpack_canonical(&p128, a, s)) {
2919         parts_float_to_float_narrow(&p64, &p128, s);
2920     } else {
2921         parts_default_nan(&p64, s);
2922     }
2923     return float32_round_pack_canonical(&p64, s);
2924 }
2925 
2926 float64 floatx80_to_float64(floatx80 a, float_status *s)
2927 {
2928     FloatParts64 p64;
2929     FloatParts128 p128;
2930 
2931     if (floatx80_unpack_canonical(&p128, a, s)) {
2932         parts_float_to_float_narrow(&p64, &p128, s);
2933     } else {
2934         parts_default_nan(&p64, s);
2935     }
2936     return float64_round_pack_canonical(&p64, s);
2937 }
2938 
2939 float128 floatx80_to_float128(floatx80 a, float_status *s)
2940 {
2941     FloatParts128 p;
2942 
2943     if (floatx80_unpack_canonical(&p, a, s)) {
2944         parts_float_to_float(&p, s);
2945     } else {
2946         parts_default_nan(&p, s);
2947     }
2948     return float128_round_pack_canonical(&p, s);
2949 }
2950 
2951 floatx80 float32_to_floatx80(float32 a, float_status *s)
2952 {
2953     FloatParts64 p64;
2954     FloatParts128 p128;
2955 
2956     float32_unpack_canonical(&p64, a, s);
2957     parts_float_to_float_widen(&p128, &p64, s);
2958     return floatx80_round_pack_canonical(&p128, s);
2959 }
2960 
2961 floatx80 float64_to_floatx80(float64 a, float_status *s)
2962 {
2963     FloatParts64 p64;
2964     FloatParts128 p128;
2965 
2966     float64_unpack_canonical(&p64, a, s);
2967     parts_float_to_float_widen(&p128, &p64, s);
2968     return floatx80_round_pack_canonical(&p128, s);
2969 }
2970 
2971 floatx80 float128_to_floatx80(float128 a, float_status *s)
2972 {
2973     FloatParts128 p;
2974 
2975     float128_unpack_canonical(&p, a, s);
2976     parts_float_to_float(&p, s);
2977     return floatx80_round_pack_canonical(&p, s);
2978 }
2979 
2980 /*
2981  * Round to integral value
2982  */
2983 
2984 float16 float16_round_to_int(float16 a, float_status *s)
2985 {
2986     FloatParts64 p;
2987 
2988     float16_unpack_canonical(&p, a, s);
2989     parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float16_params);
2990     return float16_round_pack_canonical(&p, s);
2991 }
2992 
2993 float32 float32_round_to_int(float32 a, float_status *s)
2994 {
2995     FloatParts64 p;
2996 
2997     float32_unpack_canonical(&p, a, s);
2998     parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float32_params);
2999     return float32_round_pack_canonical(&p, s);
3000 }
3001 
3002 float64 float64_round_to_int(float64 a, float_status *s)
3003 {
3004     FloatParts64 p;
3005 
3006     float64_unpack_canonical(&p, a, s);
3007     parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float64_params);
3008     return float64_round_pack_canonical(&p, s);
3009 }
3010 
3011 bfloat16 bfloat16_round_to_int(bfloat16 a, float_status *s)
3012 {
3013     FloatParts64 p;
3014 
3015     bfloat16_unpack_canonical(&p, a, s);
3016     parts_round_to_int(&p, s->float_rounding_mode, 0, s, &bfloat16_params);
3017     return bfloat16_round_pack_canonical(&p, s);
3018 }
3019 
3020 float128 float128_round_to_int(float128 a, float_status *s)
3021 {
3022     FloatParts128 p;
3023 
3024     float128_unpack_canonical(&p, a, s);
3025     parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float128_params);
3026     return float128_round_pack_canonical(&p, s);
3027 }
3028 
3029 floatx80 floatx80_round_to_int(floatx80 a, float_status *status)
3030 {
3031     FloatParts128 p;
3032 
3033     if (!floatx80_unpack_canonical(&p, a, status)) {
3034         return floatx80_default_nan(status);
3035     }
3036 
3037     parts_round_to_int(&p, status->float_rounding_mode, 0, status,
3038                        &floatx80_params[status->floatx80_rounding_precision]);
3039     return floatx80_round_pack_canonical(&p, status);
3040 }
3041 
3042 /*
3043  * Floating-point to signed integer conversions
3044  */
3045 
3046 int8_t float16_to_int8_scalbn(float16 a, FloatRoundMode rmode, int scale,
3047                               float_status *s)
3048 {
3049     FloatParts64 p;
3050 
3051     float16_unpack_canonical(&p, a, s);
3052     return parts_float_to_sint(&p, rmode, scale, INT8_MIN, INT8_MAX, s);
3053 }
3054 
3055 int16_t float16_to_int16_scalbn(float16 a, FloatRoundMode rmode, int scale,
3056                                 float_status *s)
3057 {
3058     FloatParts64 p;
3059 
3060     float16_unpack_canonical(&p, a, s);
3061     return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
3062 }
3063 
3064 int32_t float16_to_int32_scalbn(float16 a, FloatRoundMode rmode, int scale,
3065                                 float_status *s)
3066 {
3067     FloatParts64 p;
3068 
3069     float16_unpack_canonical(&p, a, s);
3070     return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
3071 }
3072 
3073 int64_t float16_to_int64_scalbn(float16 a, FloatRoundMode rmode, int scale,
3074                                 float_status *s)
3075 {
3076     FloatParts64 p;
3077 
3078     float16_unpack_canonical(&p, a, s);
3079     return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
3080 }
3081 
3082 int16_t float32_to_int16_scalbn(float32 a, FloatRoundMode rmode, int scale,
3083                                 float_status *s)
3084 {
3085     FloatParts64 p;
3086 
3087     float32_unpack_canonical(&p, a, s);
3088     return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
3089 }
3090 
3091 int32_t float32_to_int32_scalbn(float32 a, FloatRoundMode rmode, int scale,
3092                                 float_status *s)
3093 {
3094     FloatParts64 p;
3095 
3096     float32_unpack_canonical(&p, a, s);
3097     return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
3098 }
3099 
3100 int64_t float32_to_int64_scalbn(float32 a, FloatRoundMode rmode, int scale,
3101                                 float_status *s)
3102 {
3103     FloatParts64 p;
3104 
3105     float32_unpack_canonical(&p, a, s);
3106     return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
3107 }
3108 
3109 int16_t float64_to_int16_scalbn(float64 a, FloatRoundMode rmode, int scale,
3110                                 float_status *s)
3111 {
3112     FloatParts64 p;
3113 
3114     float64_unpack_canonical(&p, a, s);
3115     return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
3116 }
3117 
3118 int32_t float64_to_int32_scalbn(float64 a, FloatRoundMode rmode, int scale,
3119                                 float_status *s)
3120 {
3121     FloatParts64 p;
3122 
3123     float64_unpack_canonical(&p, a, s);
3124     return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
3125 }
3126 
3127 int64_t float64_to_int64_scalbn(float64 a, FloatRoundMode rmode, int scale,
3128                                 float_status *s)
3129 {
3130     FloatParts64 p;
3131 
3132     float64_unpack_canonical(&p, a, s);
3133     return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
3134 }
3135 
3136 int8_t bfloat16_to_int8_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
3137                                float_status *s)
3138 {
3139     FloatParts64 p;
3140 
3141     bfloat16_unpack_canonical(&p, a, s);
3142     return parts_float_to_sint(&p, rmode, scale, INT8_MIN, INT8_MAX, s);
3143 }
3144 
3145 int16_t bfloat16_to_int16_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
3146                                  float_status *s)
3147 {
3148     FloatParts64 p;
3149 
3150     bfloat16_unpack_canonical(&p, a, s);
3151     return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
3152 }
3153 
3154 int32_t bfloat16_to_int32_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
3155                                  float_status *s)
3156 {
3157     FloatParts64 p;
3158 
3159     bfloat16_unpack_canonical(&p, a, s);
3160     return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
3161 }
3162 
3163 int64_t bfloat16_to_int64_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
3164                                  float_status *s)
3165 {
3166     FloatParts64 p;
3167 
3168     bfloat16_unpack_canonical(&p, a, s);
3169     return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
3170 }
3171 
3172 static int32_t float128_to_int32_scalbn(float128 a, FloatRoundMode rmode,
3173                                         int scale, float_status *s)
3174 {
3175     FloatParts128 p;
3176 
3177     float128_unpack_canonical(&p, a, s);
3178     return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
3179 }
3180 
3181 static int64_t float128_to_int64_scalbn(float128 a, FloatRoundMode rmode,
3182                                         int scale, float_status *s)
3183 {
3184     FloatParts128 p;
3185 
3186     float128_unpack_canonical(&p, a, s);
3187     return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
3188 }
3189 
3190 static Int128 float128_to_int128_scalbn(float128 a, FloatRoundMode rmode,
3191                                         int scale, float_status *s)
3192 {
3193     int flags = 0;
3194     Int128 r;
3195     FloatParts128 p;
3196 
3197     float128_unpack_canonical(&p, a, s);
3198 
3199     switch (p.cls) {
3200     case float_class_snan:
3201         flags |= float_flag_invalid_snan;
3202         /* fall through */
3203     case float_class_qnan:
3204         flags |= float_flag_invalid;
3205         r = UINT128_MAX;
3206         break;
3207 
3208     case float_class_inf:
3209         flags = float_flag_invalid | float_flag_invalid_cvti;
3210         r = p.sign ? INT128_MIN : INT128_MAX;
3211         break;
3212 
3213     case float_class_zero:
3214         return int128_zero();
3215 
3216     case float_class_normal:
3217         if (parts_round_to_int_normal(&p, rmode, scale, 128 - 2)) {
3218             flags = float_flag_inexact;
3219         }
3220 
3221         if (p.exp < 127) {
3222             int shift = 127 - p.exp;
3223             r = int128_urshift(int128_make128(p.frac_lo, p.frac_hi), shift);
3224             if (p.sign) {
3225                 r = int128_neg(r);
3226             }
3227         } else if (p.exp == 127 && p.sign && p.frac_lo == 0 &&
3228                    p.frac_hi == DECOMPOSED_IMPLICIT_BIT) {
3229             r = INT128_MIN;
3230         } else {
3231             flags = float_flag_invalid | float_flag_invalid_cvti;
3232             r = p.sign ? INT128_MIN : INT128_MAX;
3233         }
3234         break;
3235 
3236     default:
3237         g_assert_not_reached();
3238     }
3239 
3240     float_raise(flags, s);
3241     return r;
3242 }
3243 
3244 static int32_t floatx80_to_int32_scalbn(floatx80 a, FloatRoundMode rmode,
3245                                         int scale, float_status *s)
3246 {
3247     FloatParts128 p;
3248 
3249     if (!floatx80_unpack_canonical(&p, a, s)) {
3250         parts_default_nan(&p, s);
3251     }
3252     return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
3253 }
3254 
3255 static int64_t floatx80_to_int64_scalbn(floatx80 a, FloatRoundMode rmode,
3256                                         int scale, float_status *s)
3257 {
3258     FloatParts128 p;
3259 
3260     if (!floatx80_unpack_canonical(&p, a, s)) {
3261         parts_default_nan(&p, s);
3262     }
3263     return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
3264 }
3265 
3266 int8_t float16_to_int8(float16 a, float_status *s)
3267 {
3268     return float16_to_int8_scalbn(a, s->float_rounding_mode, 0, s);
3269 }
3270 
3271 int16_t float16_to_int16(float16 a, float_status *s)
3272 {
3273     return float16_to_int16_scalbn(a, s->float_rounding_mode, 0, s);
3274 }
3275 
3276 int32_t float16_to_int32(float16 a, float_status *s)
3277 {
3278     return float16_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3279 }
3280 
3281 int64_t float16_to_int64(float16 a, float_status *s)
3282 {
3283     return float16_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3284 }
3285 
3286 int16_t float32_to_int16(float32 a, float_status *s)
3287 {
3288     return float32_to_int16_scalbn(a, s->float_rounding_mode, 0, s);
3289 }
3290 
3291 int32_t float32_to_int32(float32 a, float_status *s)
3292 {
3293     return float32_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3294 }
3295 
3296 int64_t float32_to_int64(float32 a, float_status *s)
3297 {
3298     return float32_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3299 }
3300 
3301 int16_t float64_to_int16(float64 a, float_status *s)
3302 {
3303     return float64_to_int16_scalbn(a, s->float_rounding_mode, 0, s);
3304 }
3305 
3306 int32_t float64_to_int32(float64 a, float_status *s)
3307 {
3308     return float64_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3309 }
3310 
3311 int64_t float64_to_int64(float64 a, float_status *s)
3312 {
3313     return float64_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3314 }
3315 
3316 int32_t float128_to_int32(float128 a, float_status *s)
3317 {
3318     return float128_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3319 }
3320 
3321 int64_t float128_to_int64(float128 a, float_status *s)
3322 {
3323     return float128_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3324 }
3325 
3326 Int128 float128_to_int128(float128 a, float_status *s)
3327 {
3328     return float128_to_int128_scalbn(a, s->float_rounding_mode, 0, s);
3329 }
3330 
3331 int32_t floatx80_to_int32(floatx80 a, float_status *s)
3332 {
3333     return floatx80_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3334 }
3335 
3336 int64_t floatx80_to_int64(floatx80 a, float_status *s)
3337 {
3338     return floatx80_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3339 }
3340 
3341 int16_t float16_to_int16_round_to_zero(float16 a, float_status *s)
3342 {
3343     return float16_to_int16_scalbn(a, float_round_to_zero, 0, s);
3344 }
3345 
3346 int32_t float16_to_int32_round_to_zero(float16 a, float_status *s)
3347 {
3348     return float16_to_int32_scalbn(a, float_round_to_zero, 0, s);
3349 }
3350 
3351 int64_t float16_to_int64_round_to_zero(float16 a, float_status *s)
3352 {
3353     return float16_to_int64_scalbn(a, float_round_to_zero, 0, s);
3354 }
3355 
3356 int16_t float32_to_int16_round_to_zero(float32 a, float_status *s)
3357 {
3358     return float32_to_int16_scalbn(a, float_round_to_zero, 0, s);
3359 }
3360 
3361 int32_t float32_to_int32_round_to_zero(float32 a, float_status *s)
3362 {
3363     return float32_to_int32_scalbn(a, float_round_to_zero, 0, s);
3364 }
3365 
3366 int64_t float32_to_int64_round_to_zero(float32 a, float_status *s)
3367 {
3368     return float32_to_int64_scalbn(a, float_round_to_zero, 0, s);
3369 }
3370 
3371 int16_t float64_to_int16_round_to_zero(float64 a, float_status *s)
3372 {
3373     return float64_to_int16_scalbn(a, float_round_to_zero, 0, s);
3374 }
3375 
3376 int32_t float64_to_int32_round_to_zero(float64 a, float_status *s)
3377 {
3378     return float64_to_int32_scalbn(a, float_round_to_zero, 0, s);
3379 }
3380 
3381 int64_t float64_to_int64_round_to_zero(float64 a, float_status *s)
3382 {
3383     return float64_to_int64_scalbn(a, float_round_to_zero, 0, s);
3384 }
3385 
3386 int32_t float128_to_int32_round_to_zero(float128 a, float_status *s)
3387 {
3388     return float128_to_int32_scalbn(a, float_round_to_zero, 0, s);
3389 }
3390 
3391 int64_t float128_to_int64_round_to_zero(float128 a, float_status *s)
3392 {
3393     return float128_to_int64_scalbn(a, float_round_to_zero, 0, s);
3394 }
3395 
3396 Int128 float128_to_int128_round_to_zero(float128 a, float_status *s)
3397 {
3398     return float128_to_int128_scalbn(a, float_round_to_zero, 0, s);
3399 }
3400 
3401 int32_t floatx80_to_int32_round_to_zero(floatx80 a, float_status *s)
3402 {
3403     return floatx80_to_int32_scalbn(a, float_round_to_zero, 0, s);
3404 }
3405 
3406 int64_t floatx80_to_int64_round_to_zero(floatx80 a, float_status *s)
3407 {
3408     return floatx80_to_int64_scalbn(a, float_round_to_zero, 0, s);
3409 }
3410 
3411 int8_t bfloat16_to_int8(bfloat16 a, float_status *s)
3412 {
3413     return bfloat16_to_int8_scalbn(a, s->float_rounding_mode, 0, s);
3414 }
3415 
3416 int16_t bfloat16_to_int16(bfloat16 a, float_status *s)
3417 {
3418     return bfloat16_to_int16_scalbn(a, s->float_rounding_mode, 0, s);
3419 }
3420 
3421 int32_t bfloat16_to_int32(bfloat16 a, float_status *s)
3422 {
3423     return bfloat16_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3424 }
3425 
3426 int64_t bfloat16_to_int64(bfloat16 a, float_status *s)
3427 {
3428     return bfloat16_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3429 }
3430 
3431 int8_t bfloat16_to_int8_round_to_zero(bfloat16 a, float_status *s)
3432 {
3433     return bfloat16_to_int8_scalbn(a, float_round_to_zero, 0, s);
3434 }
3435 
3436 int16_t bfloat16_to_int16_round_to_zero(bfloat16 a, float_status *s)
3437 {
3438     return bfloat16_to_int16_scalbn(a, float_round_to_zero, 0, s);
3439 }
3440 
3441 int32_t bfloat16_to_int32_round_to_zero(bfloat16 a, float_status *s)
3442 {
3443     return bfloat16_to_int32_scalbn(a, float_round_to_zero, 0, s);
3444 }
3445 
3446 int64_t bfloat16_to_int64_round_to_zero(bfloat16 a, float_status *s)
3447 {
3448     return bfloat16_to_int64_scalbn(a, float_round_to_zero, 0, s);
3449 }
3450 
3451 int32_t float64_to_int32_modulo(float64 a, FloatRoundMode rmode,
3452                                 float_status *s)
3453 {
3454     FloatParts64 p;
3455 
3456     float64_unpack_canonical(&p, a, s);
3457     return parts_float_to_sint_modulo(&p, rmode, 31, s);
3458 }
3459 
3460 int64_t float64_to_int64_modulo(float64 a, FloatRoundMode rmode,
3461                                 float_status *s)
3462 {
3463     FloatParts64 p;
3464 
3465     float64_unpack_canonical(&p, a, s);
3466     return parts_float_to_sint_modulo(&p, rmode, 63, s);
3467 }
3468 
3469 /*
3470  * Floating-point to unsigned integer conversions
3471  */
3472 
3473 uint8_t float16_to_uint8_scalbn(float16 a, FloatRoundMode rmode, int scale,
3474                                 float_status *s)
3475 {
3476     FloatParts64 p;
3477 
3478     float16_unpack_canonical(&p, a, s);
3479     return parts_float_to_uint(&p, rmode, scale, UINT8_MAX, s);
3480 }
3481 
3482 uint16_t float16_to_uint16_scalbn(float16 a, FloatRoundMode rmode, int scale,
3483                                   float_status *s)
3484 {
3485     FloatParts64 p;
3486 
3487     float16_unpack_canonical(&p, a, s);
3488     return parts_float_to_uint(&p, rmode, scale, UINT16_MAX, s);
3489 }
3490 
3491 uint32_t float16_to_uint32_scalbn(float16 a, FloatRoundMode rmode, int scale,
3492                                   float_status *s)
3493 {
3494     FloatParts64 p;
3495 
3496     float16_unpack_canonical(&p, a, s);
3497     return parts_float_to_uint(&p, rmode, scale, UINT32_MAX, s);
3498 }
3499 
3500 uint64_t float16_to_uint64_scalbn(float16 a, FloatRoundMode rmode, int scale,
3501                                   float_status *s)
3502 {
3503     FloatParts64 p;
3504 
3505     float16_unpack_canonical(&p, a, s);
3506     return parts_float_to_uint(&p, rmode, scale, UINT64_MAX, s);
3507 }
3508 
3509 uint16_t float32_to_uint16_scalbn(float32 a, FloatRoundMode rmode, int scale,
3510                                   float_status *s)
3511 {
3512     FloatParts64 p;
3513 
3514     float32_unpack_canonical(&p, a, s);
3515     return parts_float_to_uint(&p, rmode, scale, UINT16_MAX, s);
3516 }
3517 
3518 uint32_t float32_to_uint32_scalbn(float32 a, FloatRoundMode rmode, int scale,
3519                                   float_status *s)
3520 {
3521     FloatParts64 p;
3522 
3523     float32_unpack_canonical(&p, a, s);
3524     return parts_float_to_uint(&p, rmode, scale, UINT32_MAX, s);
3525 }
3526 
3527 uint64_t float32_to_uint64_scalbn(float32 a, FloatRoundMode rmode, int scale,
3528                                   float_status *s)
3529 {
3530     FloatParts64 p;
3531 
3532     float32_unpack_canonical(&p, a, s);
3533     return parts_float_to_uint(&p, rmode, scale, UINT64_MAX, s);
3534 }
3535 
3536 uint16_t float64_to_uint16_scalbn(float64 a, FloatRoundMode rmode, int scale,
3537                                   float_status *s)
3538 {
3539     FloatParts64 p;
3540 
3541     float64_unpack_canonical(&p, a, s);
3542     return parts_float_to_uint(&p, rmode, scale, UINT16_MAX, s);
3543 }
3544 
3545 uint32_t float64_to_uint32_scalbn(float64 a, FloatRoundMode rmode, int scale,
3546                                   float_status *s)
3547 {
3548     FloatParts64 p;
3549 
3550     float64_unpack_canonical(&p, a, s);
3551     return parts_float_to_uint(&p, rmode, scale, UINT32_MAX, s);
3552 }
3553 
3554 uint64_t float64_to_uint64_scalbn(float64 a, FloatRoundMode rmode, int scale,
3555                                   float_status *s)
3556 {
3557     FloatParts64 p;
3558 
3559     float64_unpack_canonical(&p, a, s);
3560     return parts_float_to_uint(&p, rmode, scale, UINT64_MAX, s);
3561 }
3562 
3563 uint8_t bfloat16_to_uint8_scalbn(bfloat16 a, FloatRoundMode rmode,
3564                                  int scale, float_status *s)
3565 {
3566     FloatParts64 p;
3567 
3568     bfloat16_unpack_canonical(&p, a, s);
3569     return parts_float_to_uint(&p, rmode, scale, UINT8_MAX, s);
3570 }
3571 
3572 uint16_t bfloat16_to_uint16_scalbn(bfloat16 a, FloatRoundMode rmode,
3573                                    int scale, float_status *s)
3574 {
3575     FloatParts64 p;
3576 
3577     bfloat16_unpack_canonical(&p, a, s);
3578     return parts_float_to_uint(&p, rmode, scale, UINT16_MAX, s);
3579 }
3580 
3581 uint32_t bfloat16_to_uint32_scalbn(bfloat16 a, FloatRoundMode rmode,
3582                                    int scale, float_status *s)
3583 {
3584     FloatParts64 p;
3585 
3586     bfloat16_unpack_canonical(&p, a, s);
3587     return parts_float_to_uint(&p, rmode, scale, UINT32_MAX, s);
3588 }
3589 
3590 uint64_t bfloat16_to_uint64_scalbn(bfloat16 a, FloatRoundMode rmode,
3591                                    int scale, float_status *s)
3592 {
3593     FloatParts64 p;
3594 
3595     bfloat16_unpack_canonical(&p, a, s);
3596     return parts_float_to_uint(&p, rmode, scale, UINT64_MAX, s);
3597 }
3598 
3599 static uint32_t float128_to_uint32_scalbn(float128 a, FloatRoundMode rmode,
3600                                           int scale, float_status *s)
3601 {
3602     FloatParts128 p;
3603 
3604     float128_unpack_canonical(&p, a, s);
3605     return parts_float_to_uint(&p, rmode, scale, UINT32_MAX, s);
3606 }
3607 
3608 static uint64_t float128_to_uint64_scalbn(float128 a, FloatRoundMode rmode,
3609                                           int scale, float_status *s)
3610 {
3611     FloatParts128 p;
3612 
3613     float128_unpack_canonical(&p, a, s);
3614     return parts_float_to_uint(&p, rmode, scale, UINT64_MAX, s);
3615 }
3616 
3617 static Int128 float128_to_uint128_scalbn(float128 a, FloatRoundMode rmode,
3618                                          int scale, float_status *s)
3619 {
3620     int flags = 0;
3621     Int128 r;
3622     FloatParts128 p;
3623 
3624     float128_unpack_canonical(&p, a, s);
3625 
3626     switch (p.cls) {
3627     case float_class_snan:
3628         flags |= float_flag_invalid_snan;
3629         /* fall through */
3630     case float_class_qnan:
3631         flags |= float_flag_invalid;
3632         r = UINT128_MAX;
3633         break;
3634 
3635     case float_class_inf:
3636         flags = float_flag_invalid | float_flag_invalid_cvti;
3637         r = p.sign ? int128_zero() : UINT128_MAX;
3638         break;
3639 
3640     case float_class_zero:
3641         return int128_zero();
3642 
3643     case float_class_normal:
3644         if (parts_round_to_int_normal(&p, rmode, scale, 128 - 2)) {
3645             flags = float_flag_inexact;
3646             if (p.cls == float_class_zero) {
3647                 r = int128_zero();
3648                 break;
3649             }
3650         }
3651 
3652         if (p.sign) {
3653             flags = float_flag_invalid | float_flag_invalid_cvti;
3654             r = int128_zero();
3655         } else if (p.exp <= 127) {
3656             int shift = 127 - p.exp;
3657             r = int128_urshift(int128_make128(p.frac_lo, p.frac_hi), shift);
3658         } else {
3659             flags = float_flag_invalid | float_flag_invalid_cvti;
3660             r = UINT128_MAX;
3661         }
3662         break;
3663 
3664     default:
3665         g_assert_not_reached();
3666     }
3667 
3668     float_raise(flags, s);
3669     return r;
3670 }
3671 
3672 uint8_t float16_to_uint8(float16 a, float_status *s)
3673 {
3674     return float16_to_uint8_scalbn(a, s->float_rounding_mode, 0, s);
3675 }
3676 
3677 uint16_t float16_to_uint16(float16 a, float_status *s)
3678 {
3679     return float16_to_uint16_scalbn(a, s->float_rounding_mode, 0, s);
3680 }
3681 
3682 uint32_t float16_to_uint32(float16 a, float_status *s)
3683 {
3684     return float16_to_uint32_scalbn(a, s->float_rounding_mode, 0, s);
3685 }
3686 
3687 uint64_t float16_to_uint64(float16 a, float_status *s)
3688 {
3689     return float16_to_uint64_scalbn(a, s->float_rounding_mode, 0, s);
3690 }
3691 
3692 uint16_t float32_to_uint16(float32 a, float_status *s)
3693 {
3694     return float32_to_uint16_scalbn(a, s->float_rounding_mode, 0, s);
3695 }
3696 
3697 uint32_t float32_to_uint32(float32 a, float_status *s)
3698 {
3699     return float32_to_uint32_scalbn(a, s->float_rounding_mode, 0, s);
3700 }
3701 
3702 uint64_t float32_to_uint64(float32 a, float_status *s)
3703 {
3704     return float32_to_uint64_scalbn(a, s->float_rounding_mode, 0, s);
3705 }
3706 
3707 uint16_t float64_to_uint16(float64 a, float_status *s)
3708 {
3709     return float64_to_uint16_scalbn(a, s->float_rounding_mode, 0, s);
3710 }
3711 
3712 uint32_t float64_to_uint32(float64 a, float_status *s)
3713 {
3714     return float64_to_uint32_scalbn(a, s->float_rounding_mode, 0, s);
3715 }
3716 
3717 uint64_t float64_to_uint64(float64 a, float_status *s)
3718 {
3719     return float64_to_uint64_scalbn(a, s->float_rounding_mode, 0, s);
3720 }
3721 
3722 uint32_t float128_to_uint32(float128 a, float_status *s)
3723 {
3724     return float128_to_uint32_scalbn(a, s->float_rounding_mode, 0, s);
3725 }
3726 
3727 uint64_t float128_to_uint64(float128 a, float_status *s)
3728 {
3729     return float128_to_uint64_scalbn(a, s->float_rounding_mode, 0, s);
3730 }
3731 
3732 Int128 float128_to_uint128(float128 a, float_status *s)
3733 {
3734     return float128_to_uint128_scalbn(a, s->float_rounding_mode, 0, s);
3735 }
3736 
3737 uint16_t float16_to_uint16_round_to_zero(float16 a, float_status *s)
3738 {
3739     return float16_to_uint16_scalbn(a, float_round_to_zero, 0, s);
3740 }
3741 
3742 uint32_t float16_to_uint32_round_to_zero(float16 a, float_status *s)
3743 {
3744     return float16_to_uint32_scalbn(a, float_round_to_zero, 0, s);
3745 }
3746 
3747 uint64_t float16_to_uint64_round_to_zero(float16 a, float_status *s)
3748 {
3749     return float16_to_uint64_scalbn(a, float_round_to_zero, 0, s);
3750 }
3751 
3752 uint16_t float32_to_uint16_round_to_zero(float32 a, float_status *s)
3753 {
3754     return float32_to_uint16_scalbn(a, float_round_to_zero, 0, s);
3755 }
3756 
3757 uint32_t float32_to_uint32_round_to_zero(float32 a, float_status *s)
3758 {
3759     return float32_to_uint32_scalbn(a, float_round_to_zero, 0, s);
3760 }
3761 
3762 uint64_t float32_to_uint64_round_to_zero(float32 a, float_status *s)
3763 {
3764     return float32_to_uint64_scalbn(a, float_round_to_zero, 0, s);
3765 }
3766 
3767 uint16_t float64_to_uint16_round_to_zero(float64 a, float_status *s)
3768 {
3769     return float64_to_uint16_scalbn(a, float_round_to_zero, 0, s);
3770 }
3771 
3772 uint32_t float64_to_uint32_round_to_zero(float64 a, float_status *s)
3773 {
3774     return float64_to_uint32_scalbn(a, float_round_to_zero, 0, s);
3775 }
3776 
3777 uint64_t float64_to_uint64_round_to_zero(float64 a, float_status *s)
3778 {
3779     return float64_to_uint64_scalbn(a, float_round_to_zero, 0, s);
3780 }
3781 
3782 uint32_t float128_to_uint32_round_to_zero(float128 a, float_status *s)
3783 {
3784     return float128_to_uint32_scalbn(a, float_round_to_zero, 0, s);
3785 }
3786 
3787 uint64_t float128_to_uint64_round_to_zero(float128 a, float_status *s)
3788 {
3789     return float128_to_uint64_scalbn(a, float_round_to_zero, 0, s);
3790 }
3791 
3792 Int128 float128_to_uint128_round_to_zero(float128 a, float_status *s)
3793 {
3794     return float128_to_uint128_scalbn(a, float_round_to_zero, 0, s);
3795 }
3796 
3797 uint8_t bfloat16_to_uint8(bfloat16 a, float_status *s)
3798 {
3799     return bfloat16_to_uint8_scalbn(a, s->float_rounding_mode, 0, s);
3800 }
3801 
3802 uint16_t bfloat16_to_uint16(bfloat16 a, float_status *s)
3803 {
3804     return bfloat16_to_uint16_scalbn(a, s->float_rounding_mode, 0, s);
3805 }
3806 
3807 uint32_t bfloat16_to_uint32(bfloat16 a, float_status *s)
3808 {
3809     return bfloat16_to_uint32_scalbn(a, s->float_rounding_mode, 0, s);
3810 }
3811 
3812 uint64_t bfloat16_to_uint64(bfloat16 a, float_status *s)
3813 {
3814     return bfloat16_to_uint64_scalbn(a, s->float_rounding_mode, 0, s);
3815 }
3816 
3817 uint8_t bfloat16_to_uint8_round_to_zero(bfloat16 a, float_status *s)
3818 {
3819     return bfloat16_to_uint8_scalbn(a, float_round_to_zero, 0, s);
3820 }
3821 
3822 uint16_t bfloat16_to_uint16_round_to_zero(bfloat16 a, float_status *s)
3823 {
3824     return bfloat16_to_uint16_scalbn(a, float_round_to_zero, 0, s);
3825 }
3826 
3827 uint32_t bfloat16_to_uint32_round_to_zero(bfloat16 a, float_status *s)
3828 {
3829     return bfloat16_to_uint32_scalbn(a, float_round_to_zero, 0, s);
3830 }
3831 
3832 uint64_t bfloat16_to_uint64_round_to_zero(bfloat16 a, float_status *s)
3833 {
3834     return bfloat16_to_uint64_scalbn(a, float_round_to_zero, 0, s);
3835 }
3836 
3837 /*
3838  * Signed integer to floating-point conversions
3839  */
3840 
3841 float16 int64_to_float16_scalbn(int64_t a, int scale, float_status *status)
3842 {
3843     FloatParts64 p;
3844 
3845     parts_sint_to_float(&p, a, scale, status);
3846     return float16_round_pack_canonical(&p, status);
3847 }
3848 
3849 float16 int32_to_float16_scalbn(int32_t a, int scale, float_status *status)
3850 {
3851     return int64_to_float16_scalbn(a, scale, status);
3852 }
3853 
3854 float16 int16_to_float16_scalbn(int16_t a, int scale, float_status *status)
3855 {
3856     return int64_to_float16_scalbn(a, scale, status);
3857 }
3858 
3859 float16 int64_to_float16(int64_t a, float_status *status)
3860 {
3861     return int64_to_float16_scalbn(a, 0, status);
3862 }
3863 
3864 float16 int32_to_float16(int32_t a, float_status *status)
3865 {
3866     return int64_to_float16_scalbn(a, 0, status);
3867 }
3868 
3869 float16 int16_to_float16(int16_t a, float_status *status)
3870 {
3871     return int64_to_float16_scalbn(a, 0, status);
3872 }
3873 
3874 float16 int8_to_float16(int8_t a, float_status *status)
3875 {
3876     return int64_to_float16_scalbn(a, 0, status);
3877 }
3878 
3879 float32 int64_to_float32_scalbn(int64_t a, int scale, float_status *status)
3880 {
3881     FloatParts64 p;
3882 
3883     /* Without scaling, there are no overflow concerns. */
3884     if (likely(scale == 0) && can_use_fpu(status)) {
3885         union_float32 ur;
3886         ur.h = a;
3887         return ur.s;
3888     }
3889 
3890     parts64_sint_to_float(&p, a, scale, status);
3891     return float32_round_pack_canonical(&p, status);
3892 }
3893 
3894 float32 int32_to_float32_scalbn(int32_t a, int scale, float_status *status)
3895 {
3896     return int64_to_float32_scalbn(a, scale, status);
3897 }
3898 
3899 float32 int16_to_float32_scalbn(int16_t a, int scale, float_status *status)
3900 {
3901     return int64_to_float32_scalbn(a, scale, status);
3902 }
3903 
3904 float32 int64_to_float32(int64_t a, float_status *status)
3905 {
3906     return int64_to_float32_scalbn(a, 0, status);
3907 }
3908 
3909 float32 int32_to_float32(int32_t a, float_status *status)
3910 {
3911     return int64_to_float32_scalbn(a, 0, status);
3912 }
3913 
3914 float32 int16_to_float32(int16_t a, float_status *status)
3915 {
3916     return int64_to_float32_scalbn(a, 0, status);
3917 }
3918 
3919 float64 int64_to_float64_scalbn(int64_t a, int scale, float_status *status)
3920 {
3921     FloatParts64 p;
3922 
3923     /* Without scaling, there are no overflow concerns. */
3924     if (likely(scale == 0) && can_use_fpu(status)) {
3925         union_float64 ur;
3926         ur.h = a;
3927         return ur.s;
3928     }
3929 
3930     parts_sint_to_float(&p, a, scale, status);
3931     return float64_round_pack_canonical(&p, status);
3932 }
3933 
3934 float64 int32_to_float64_scalbn(int32_t a, int scale, float_status *status)
3935 {
3936     return int64_to_float64_scalbn(a, scale, status);
3937 }
3938 
3939 float64 int16_to_float64_scalbn(int16_t a, int scale, float_status *status)
3940 {
3941     return int64_to_float64_scalbn(a, scale, status);
3942 }
3943 
3944 float64 int64_to_float64(int64_t a, float_status *status)
3945 {
3946     return int64_to_float64_scalbn(a, 0, status);
3947 }
3948 
3949 float64 int32_to_float64(int32_t a, float_status *status)
3950 {
3951     return int64_to_float64_scalbn(a, 0, status);
3952 }
3953 
3954 float64 int16_to_float64(int16_t a, float_status *status)
3955 {
3956     return int64_to_float64_scalbn(a, 0, status);
3957 }
3958 
3959 bfloat16 int64_to_bfloat16_scalbn(int64_t a, int scale, float_status *status)
3960 {
3961     FloatParts64 p;
3962 
3963     parts_sint_to_float(&p, a, scale, status);
3964     return bfloat16_round_pack_canonical(&p, status);
3965 }
3966 
3967 bfloat16 int32_to_bfloat16_scalbn(int32_t a, int scale, float_status *status)
3968 {
3969     return int64_to_bfloat16_scalbn(a, scale, status);
3970 }
3971 
3972 bfloat16 int16_to_bfloat16_scalbn(int16_t a, int scale, float_status *status)
3973 {
3974     return int64_to_bfloat16_scalbn(a, scale, status);
3975 }
3976 
3977 bfloat16 int8_to_bfloat16_scalbn(int8_t a, int scale, float_status *status)
3978 {
3979     return int64_to_bfloat16_scalbn(a, scale, status);
3980 }
3981 
3982 bfloat16 int64_to_bfloat16(int64_t a, float_status *status)
3983 {
3984     return int64_to_bfloat16_scalbn(a, 0, status);
3985 }
3986 
3987 bfloat16 int32_to_bfloat16(int32_t a, float_status *status)
3988 {
3989     return int64_to_bfloat16_scalbn(a, 0, status);
3990 }
3991 
3992 bfloat16 int16_to_bfloat16(int16_t a, float_status *status)
3993 {
3994     return int64_to_bfloat16_scalbn(a, 0, status);
3995 }
3996 
3997 bfloat16 int8_to_bfloat16(int8_t a, float_status *status)
3998 {
3999     return int64_to_bfloat16_scalbn(a, 0, status);
4000 }
4001 
4002 float128 int128_to_float128(Int128 a, float_status *status)
4003 {
4004     FloatParts128 p = { };
4005     int shift;
4006 
4007     if (int128_nz(a)) {
4008         p.cls = float_class_normal;
4009         if (!int128_nonneg(a)) {
4010             p.sign = true;
4011             a = int128_neg(a);
4012         }
4013 
4014         shift = clz64(int128_gethi(a));
4015         if (shift == 64) {
4016             shift += clz64(int128_getlo(a));
4017         }
4018 
4019         p.exp = 127 - shift;
4020         a = int128_lshift(a, shift);
4021 
4022         p.frac_hi = int128_gethi(a);
4023         p.frac_lo = int128_getlo(a);
4024     } else {
4025         p.cls = float_class_zero;
4026     }
4027 
4028     return float128_round_pack_canonical(&p, status);
4029 }
4030 
4031 float128 int64_to_float128(int64_t a, float_status *status)
4032 {
4033     FloatParts128 p;
4034 
4035     parts_sint_to_float(&p, a, 0, status);
4036     return float128_round_pack_canonical(&p, status);
4037 }
4038 
4039 float128 int32_to_float128(int32_t a, float_status *status)
4040 {
4041     return int64_to_float128(a, status);
4042 }
4043 
4044 floatx80 int64_to_floatx80(int64_t a, float_status *status)
4045 {
4046     FloatParts128 p;
4047 
4048     parts_sint_to_float(&p, a, 0, status);
4049     return floatx80_round_pack_canonical(&p, status);
4050 }
4051 
4052 floatx80 int32_to_floatx80(int32_t a, float_status *status)
4053 {
4054     return int64_to_floatx80(a, status);
4055 }
4056 
4057 /*
4058  * Unsigned Integer to floating-point conversions
4059  */
4060 
4061 float16 uint64_to_float16_scalbn(uint64_t a, int scale, float_status *status)
4062 {
4063     FloatParts64 p;
4064 
4065     parts_uint_to_float(&p, a, scale, status);
4066     return float16_round_pack_canonical(&p, status);
4067 }
4068 
4069 float16 uint32_to_float16_scalbn(uint32_t a, int scale, float_status *status)
4070 {
4071     return uint64_to_float16_scalbn(a, scale, status);
4072 }
4073 
4074 float16 uint16_to_float16_scalbn(uint16_t a, int scale, float_status *status)
4075 {
4076     return uint64_to_float16_scalbn(a, scale, status);
4077 }
4078 
4079 float16 uint64_to_float16(uint64_t a, float_status *status)
4080 {
4081     return uint64_to_float16_scalbn(a, 0, status);
4082 }
4083 
4084 float16 uint32_to_float16(uint32_t a, float_status *status)
4085 {
4086     return uint64_to_float16_scalbn(a, 0, status);
4087 }
4088 
4089 float16 uint16_to_float16(uint16_t a, float_status *status)
4090 {
4091     return uint64_to_float16_scalbn(a, 0, status);
4092 }
4093 
4094 float16 uint8_to_float16(uint8_t a, float_status *status)
4095 {
4096     return uint64_to_float16_scalbn(a, 0, status);
4097 }
4098 
4099 float32 uint64_to_float32_scalbn(uint64_t a, int scale, float_status *status)
4100 {
4101     FloatParts64 p;
4102 
4103     /* Without scaling, there are no overflow concerns. */
4104     if (likely(scale == 0) && can_use_fpu(status)) {
4105         union_float32 ur;
4106         ur.h = a;
4107         return ur.s;
4108     }
4109 
4110     parts_uint_to_float(&p, a, scale, status);
4111     return float32_round_pack_canonical(&p, status);
4112 }
4113 
4114 float32 uint32_to_float32_scalbn(uint32_t a, int scale, float_status *status)
4115 {
4116     return uint64_to_float32_scalbn(a, scale, status);
4117 }
4118 
4119 float32 uint16_to_float32_scalbn(uint16_t a, int scale, float_status *status)
4120 {
4121     return uint64_to_float32_scalbn(a, scale, status);
4122 }
4123 
4124 float32 uint64_to_float32(uint64_t a, float_status *status)
4125 {
4126     return uint64_to_float32_scalbn(a, 0, status);
4127 }
4128 
4129 float32 uint32_to_float32(uint32_t a, float_status *status)
4130 {
4131     return uint64_to_float32_scalbn(a, 0, status);
4132 }
4133 
4134 float32 uint16_to_float32(uint16_t a, float_status *status)
4135 {
4136     return uint64_to_float32_scalbn(a, 0, status);
4137 }
4138 
4139 float64 uint64_to_float64_scalbn(uint64_t a, int scale, float_status *status)
4140 {
4141     FloatParts64 p;
4142 
4143     /* Without scaling, there are no overflow concerns. */
4144     if (likely(scale == 0) && can_use_fpu(status)) {
4145         union_float64 ur;
4146         ur.h = a;
4147         return ur.s;
4148     }
4149 
4150     parts_uint_to_float(&p, a, scale, status);
4151     return float64_round_pack_canonical(&p, status);
4152 }
4153 
4154 float64 uint32_to_float64_scalbn(uint32_t a, int scale, float_status *status)
4155 {
4156     return uint64_to_float64_scalbn(a, scale, status);
4157 }
4158 
4159 float64 uint16_to_float64_scalbn(uint16_t a, int scale, float_status *status)
4160 {
4161     return uint64_to_float64_scalbn(a, scale, status);
4162 }
4163 
4164 float64 uint64_to_float64(uint64_t a, float_status *status)
4165 {
4166     return uint64_to_float64_scalbn(a, 0, status);
4167 }
4168 
4169 float64 uint32_to_float64(uint32_t a, float_status *status)
4170 {
4171     return uint64_to_float64_scalbn(a, 0, status);
4172 }
4173 
4174 float64 uint16_to_float64(uint16_t a, float_status *status)
4175 {
4176     return uint64_to_float64_scalbn(a, 0, status);
4177 }
4178 
4179 bfloat16 uint64_to_bfloat16_scalbn(uint64_t a, int scale, float_status *status)
4180 {
4181     FloatParts64 p;
4182 
4183     parts_uint_to_float(&p, a, scale, status);
4184     return bfloat16_round_pack_canonical(&p, status);
4185 }
4186 
4187 bfloat16 uint32_to_bfloat16_scalbn(uint32_t a, int scale, float_status *status)
4188 {
4189     return uint64_to_bfloat16_scalbn(a, scale, status);
4190 }
4191 
4192 bfloat16 uint16_to_bfloat16_scalbn(uint16_t a, int scale, float_status *status)
4193 {
4194     return uint64_to_bfloat16_scalbn(a, scale, status);
4195 }
4196 
4197 bfloat16 uint8_to_bfloat16_scalbn(uint8_t a, int scale, float_status *status)
4198 {
4199     return uint64_to_bfloat16_scalbn(a, scale, status);
4200 }
4201 
4202 bfloat16 uint64_to_bfloat16(uint64_t a, float_status *status)
4203 {
4204     return uint64_to_bfloat16_scalbn(a, 0, status);
4205 }
4206 
4207 bfloat16 uint32_to_bfloat16(uint32_t a, float_status *status)
4208 {
4209     return uint64_to_bfloat16_scalbn(a, 0, status);
4210 }
4211 
4212 bfloat16 uint16_to_bfloat16(uint16_t a, float_status *status)
4213 {
4214     return uint64_to_bfloat16_scalbn(a, 0, status);
4215 }
4216 
4217 bfloat16 uint8_to_bfloat16(uint8_t a, float_status *status)
4218 {
4219     return uint64_to_bfloat16_scalbn(a, 0, status);
4220 }
4221 
4222 float128 uint64_to_float128(uint64_t a, float_status *status)
4223 {
4224     FloatParts128 p;
4225 
4226     parts_uint_to_float(&p, a, 0, status);
4227     return float128_round_pack_canonical(&p, status);
4228 }
4229 
4230 float128 uint128_to_float128(Int128 a, float_status *status)
4231 {
4232     FloatParts128 p = { };
4233     int shift;
4234 
4235     if (int128_nz(a)) {
4236         p.cls = float_class_normal;
4237 
4238         shift = clz64(int128_gethi(a));
4239         if (shift == 64) {
4240             shift += clz64(int128_getlo(a));
4241         }
4242 
4243         p.exp = 127 - shift;
4244         a = int128_lshift(a, shift);
4245 
4246         p.frac_hi = int128_gethi(a);
4247         p.frac_lo = int128_getlo(a);
4248     } else {
4249         p.cls = float_class_zero;
4250     }
4251 
4252     return float128_round_pack_canonical(&p, status);
4253 }
4254 
4255 /*
4256  * Minimum and maximum
4257  */
4258 
4259 static float16 float16_minmax(float16 a, float16 b, float_status *s, int flags)
4260 {
4261     FloatParts64 pa, pb, *pr;
4262 
4263     float16_unpack_canonical(&pa, a, s);
4264     float16_unpack_canonical(&pb, b, s);
4265     pr = parts_minmax(&pa, &pb, s, flags);
4266 
4267     return float16_round_pack_canonical(pr, s);
4268 }
4269 
4270 static bfloat16 bfloat16_minmax(bfloat16 a, bfloat16 b,
4271                                 float_status *s, int flags)
4272 {
4273     FloatParts64 pa, pb, *pr;
4274 
4275     bfloat16_unpack_canonical(&pa, a, s);
4276     bfloat16_unpack_canonical(&pb, b, s);
4277     pr = parts_minmax(&pa, &pb, s, flags);
4278 
4279     return bfloat16_round_pack_canonical(pr, s);
4280 }
4281 
4282 static float32 float32_minmax(float32 a, float32 b, float_status *s, int flags)
4283 {
4284     FloatParts64 pa, pb, *pr;
4285 
4286     float32_unpack_canonical(&pa, a, s);
4287     float32_unpack_canonical(&pb, b, s);
4288     pr = parts_minmax(&pa, &pb, s, flags);
4289 
4290     return float32_round_pack_canonical(pr, s);
4291 }
4292 
4293 static float64 float64_minmax(float64 a, float64 b, float_status *s, int flags)
4294 {
4295     FloatParts64 pa, pb, *pr;
4296 
4297     float64_unpack_canonical(&pa, a, s);
4298     float64_unpack_canonical(&pb, b, s);
4299     pr = parts_minmax(&pa, &pb, s, flags);
4300 
4301     return float64_round_pack_canonical(pr, s);
4302 }
4303 
4304 static float128 float128_minmax(float128 a, float128 b,
4305                                 float_status *s, int flags)
4306 {
4307     FloatParts128 pa, pb, *pr;
4308 
4309     float128_unpack_canonical(&pa, a, s);
4310     float128_unpack_canonical(&pb, b, s);
4311     pr = parts_minmax(&pa, &pb, s, flags);
4312 
4313     return float128_round_pack_canonical(pr, s);
4314 }
4315 
4316 #define MINMAX_1(type, name, flags) \
4317     type type##_##name(type a, type b, float_status *s) \
4318     { return type##_minmax(a, b, s, flags); }
4319 
4320 #define MINMAX_2(type) \
4321     MINMAX_1(type, max, 0)                                                \
4322     MINMAX_1(type, maxnum, minmax_isnum)                                  \
4323     MINMAX_1(type, maxnummag, minmax_isnum | minmax_ismag)                \
4324     MINMAX_1(type, maximum_number, minmax_isnumber)                       \
4325     MINMAX_1(type, min, minmax_ismin)                                     \
4326     MINMAX_1(type, minnum, minmax_ismin | minmax_isnum)                   \
4327     MINMAX_1(type, minnummag, minmax_ismin | minmax_isnum | minmax_ismag) \
4328     MINMAX_1(type, minimum_number, minmax_ismin | minmax_isnumber)        \
4329 
4330 MINMAX_2(float16)
4331 MINMAX_2(bfloat16)
4332 MINMAX_2(float32)
4333 MINMAX_2(float64)
4334 MINMAX_2(float128)
4335 
4336 #undef MINMAX_1
4337 #undef MINMAX_2
4338 
4339 /*
4340  * Floating point compare
4341  */
4342 
4343 static FloatRelation QEMU_FLATTEN
4344 float16_do_compare(float16 a, float16 b, float_status *s, bool is_quiet)
4345 {
4346     FloatParts64 pa, pb;
4347 
4348     float16_unpack_canonical(&pa, a, s);
4349     float16_unpack_canonical(&pb, b, s);
4350     return parts_compare(&pa, &pb, s, is_quiet);
4351 }
4352 
4353 FloatRelation float16_compare(float16 a, float16 b, float_status *s)
4354 {
4355     return float16_do_compare(a, b, s, false);
4356 }
4357 
4358 FloatRelation float16_compare_quiet(float16 a, float16 b, float_status *s)
4359 {
4360     return float16_do_compare(a, b, s, true);
4361 }
4362 
4363 static FloatRelation QEMU_SOFTFLOAT_ATTR
4364 float32_do_compare(float32 a, float32 b, float_status *s, bool is_quiet)
4365 {
4366     FloatParts64 pa, pb;
4367 
4368     float32_unpack_canonical(&pa, a, s);
4369     float32_unpack_canonical(&pb, b, s);
4370     return parts_compare(&pa, &pb, s, is_quiet);
4371 }
4372 
4373 static FloatRelation QEMU_FLATTEN
4374 float32_hs_compare(float32 xa, float32 xb, float_status *s, bool is_quiet)
4375 {
4376     union_float32 ua, ub;
4377 
4378     ua.s = xa;
4379     ub.s = xb;
4380 
4381     if (QEMU_NO_HARDFLOAT) {
4382         goto soft;
4383     }
4384 
4385     float32_input_flush2(&ua.s, &ub.s, s);
4386     if (isgreaterequal(ua.h, ub.h)) {
4387         if (isgreater(ua.h, ub.h)) {
4388             return float_relation_greater;
4389         }
4390         return float_relation_equal;
4391     }
4392     if (likely(isless(ua.h, ub.h))) {
4393         return float_relation_less;
4394     }
4395     /*
4396      * The only condition remaining is unordered.
4397      * Fall through to set flags.
4398      */
4399  soft:
4400     return float32_do_compare(ua.s, ub.s, s, is_quiet);
4401 }
4402 
4403 FloatRelation float32_compare(float32 a, float32 b, float_status *s)
4404 {
4405     return float32_hs_compare(a, b, s, false);
4406 }
4407 
4408 FloatRelation float32_compare_quiet(float32 a, float32 b, float_status *s)
4409 {
4410     return float32_hs_compare(a, b, s, true);
4411 }
4412 
4413 static FloatRelation QEMU_SOFTFLOAT_ATTR
4414 float64_do_compare(float64 a, float64 b, float_status *s, bool is_quiet)
4415 {
4416     FloatParts64 pa, pb;
4417 
4418     float64_unpack_canonical(&pa, a, s);
4419     float64_unpack_canonical(&pb, b, s);
4420     return parts_compare(&pa, &pb, s, is_quiet);
4421 }
4422 
4423 static FloatRelation QEMU_FLATTEN
4424 float64_hs_compare(float64 xa, float64 xb, float_status *s, bool is_quiet)
4425 {
4426     union_float64 ua, ub;
4427 
4428     ua.s = xa;
4429     ub.s = xb;
4430 
4431     if (QEMU_NO_HARDFLOAT) {
4432         goto soft;
4433     }
4434 
4435     float64_input_flush2(&ua.s, &ub.s, s);
4436     if (isgreaterequal(ua.h, ub.h)) {
4437         if (isgreater(ua.h, ub.h)) {
4438             return float_relation_greater;
4439         }
4440         return float_relation_equal;
4441     }
4442     if (likely(isless(ua.h, ub.h))) {
4443         return float_relation_less;
4444     }
4445     /*
4446      * The only condition remaining is unordered.
4447      * Fall through to set flags.
4448      */
4449  soft:
4450     return float64_do_compare(ua.s, ub.s, s, is_quiet);
4451 }
4452 
4453 FloatRelation float64_compare(float64 a, float64 b, float_status *s)
4454 {
4455     return float64_hs_compare(a, b, s, false);
4456 }
4457 
4458 FloatRelation float64_compare_quiet(float64 a, float64 b, float_status *s)
4459 {
4460     return float64_hs_compare(a, b, s, true);
4461 }
4462 
4463 static FloatRelation QEMU_FLATTEN
4464 bfloat16_do_compare(bfloat16 a, bfloat16 b, float_status *s, bool is_quiet)
4465 {
4466     FloatParts64 pa, pb;
4467 
4468     bfloat16_unpack_canonical(&pa, a, s);
4469     bfloat16_unpack_canonical(&pb, b, s);
4470     return parts_compare(&pa, &pb, s, is_quiet);
4471 }
4472 
4473 FloatRelation bfloat16_compare(bfloat16 a, bfloat16 b, float_status *s)
4474 {
4475     return bfloat16_do_compare(a, b, s, false);
4476 }
4477 
4478 FloatRelation bfloat16_compare_quiet(bfloat16 a, bfloat16 b, float_status *s)
4479 {
4480     return bfloat16_do_compare(a, b, s, true);
4481 }
4482 
4483 static FloatRelation QEMU_FLATTEN
4484 float128_do_compare(float128 a, float128 b, float_status *s, bool is_quiet)
4485 {
4486     FloatParts128 pa, pb;
4487 
4488     float128_unpack_canonical(&pa, a, s);
4489     float128_unpack_canonical(&pb, b, s);
4490     return parts_compare(&pa, &pb, s, is_quiet);
4491 }
4492 
4493 FloatRelation float128_compare(float128 a, float128 b, float_status *s)
4494 {
4495     return float128_do_compare(a, b, s, false);
4496 }
4497 
4498 FloatRelation float128_compare_quiet(float128 a, float128 b, float_status *s)
4499 {
4500     return float128_do_compare(a, b, s, true);
4501 }
4502 
4503 static FloatRelation QEMU_FLATTEN
4504 floatx80_do_compare(floatx80 a, floatx80 b, float_status *s, bool is_quiet)
4505 {
4506     FloatParts128 pa, pb;
4507 
4508     if (!floatx80_unpack_canonical(&pa, a, s) ||
4509         !floatx80_unpack_canonical(&pb, b, s)) {
4510         return float_relation_unordered;
4511     }
4512     return parts_compare(&pa, &pb, s, is_quiet);
4513 }
4514 
4515 FloatRelation floatx80_compare(floatx80 a, floatx80 b, float_status *s)
4516 {
4517     return floatx80_do_compare(a, b, s, false);
4518 }
4519 
4520 FloatRelation floatx80_compare_quiet(floatx80 a, floatx80 b, float_status *s)
4521 {
4522     return floatx80_do_compare(a, b, s, true);
4523 }
4524 
4525 /*
4526  * Scale by 2**N
4527  */
4528 
4529 float16 float16_scalbn(float16 a, int n, float_status *status)
4530 {
4531     FloatParts64 p;
4532 
4533     float16_unpack_canonical(&p, a, status);
4534     parts_scalbn(&p, n, status);
4535     return float16_round_pack_canonical(&p, status);
4536 }
4537 
4538 float32 float32_scalbn(float32 a, int n, float_status *status)
4539 {
4540     FloatParts64 p;
4541 
4542     float32_unpack_canonical(&p, a, status);
4543     parts_scalbn(&p, n, status);
4544     return float32_round_pack_canonical(&p, status);
4545 }
4546 
4547 float64 float64_scalbn(float64 a, int n, float_status *status)
4548 {
4549     FloatParts64 p;
4550 
4551     float64_unpack_canonical(&p, a, status);
4552     parts_scalbn(&p, n, status);
4553     return float64_round_pack_canonical(&p, status);
4554 }
4555 
4556 bfloat16 bfloat16_scalbn(bfloat16 a, int n, float_status *status)
4557 {
4558     FloatParts64 p;
4559 
4560     bfloat16_unpack_canonical(&p, a, status);
4561     parts_scalbn(&p, n, status);
4562     return bfloat16_round_pack_canonical(&p, status);
4563 }
4564 
4565 float128 float128_scalbn(float128 a, int n, float_status *status)
4566 {
4567     FloatParts128 p;
4568 
4569     float128_unpack_canonical(&p, a, status);
4570     parts_scalbn(&p, n, status);
4571     return float128_round_pack_canonical(&p, status);
4572 }
4573 
4574 floatx80 floatx80_scalbn(floatx80 a, int n, float_status *status)
4575 {
4576     FloatParts128 p;
4577 
4578     if (!floatx80_unpack_canonical(&p, a, status)) {
4579         return floatx80_default_nan(status);
4580     }
4581     parts_scalbn(&p, n, status);
4582     return floatx80_round_pack_canonical(&p, status);
4583 }
4584 
4585 /*
4586  * Square Root
4587  */
4588 
4589 float16 QEMU_FLATTEN float16_sqrt(float16 a, float_status *status)
4590 {
4591     FloatParts64 p;
4592 
4593     float16_unpack_canonical(&p, a, status);
4594     parts_sqrt(&p, status, &float16_params);
4595     return float16_round_pack_canonical(&p, status);
4596 }
4597 
4598 static float32 QEMU_SOFTFLOAT_ATTR
4599 soft_f32_sqrt(float32 a, float_status *status)
4600 {
4601     FloatParts64 p;
4602 
4603     float32_unpack_canonical(&p, a, status);
4604     parts_sqrt(&p, status, &float32_params);
4605     return float32_round_pack_canonical(&p, status);
4606 }
4607 
4608 static float64 QEMU_SOFTFLOAT_ATTR
4609 soft_f64_sqrt(float64 a, float_status *status)
4610 {
4611     FloatParts64 p;
4612 
4613     float64_unpack_canonical(&p, a, status);
4614     parts_sqrt(&p, status, &float64_params);
4615     return float64_round_pack_canonical(&p, status);
4616 }
4617 
4618 float32 QEMU_FLATTEN float32_sqrt(float32 xa, float_status *s)
4619 {
4620     union_float32 ua, ur;
4621 
4622     ua.s = xa;
4623     if (unlikely(!can_use_fpu(s))) {
4624         goto soft;
4625     }
4626 
4627     float32_input_flush1(&ua.s, s);
4628     if (QEMU_HARDFLOAT_1F32_USE_FP) {
4629         if (unlikely(!(fpclassify(ua.h) == FP_NORMAL ||
4630                        fpclassify(ua.h) == FP_ZERO) ||
4631                      signbit(ua.h))) {
4632             goto soft;
4633         }
4634     } else if (unlikely(!float32_is_zero_or_normal(ua.s) ||
4635                         float32_is_neg(ua.s))) {
4636         goto soft;
4637     }
4638     ur.h = sqrtf(ua.h);
4639     return ur.s;
4640 
4641  soft:
4642     return soft_f32_sqrt(ua.s, s);
4643 }
4644 
4645 float64 QEMU_FLATTEN float64_sqrt(float64 xa, float_status *s)
4646 {
4647     union_float64 ua, ur;
4648 
4649     ua.s = xa;
4650     if (unlikely(!can_use_fpu(s))) {
4651         goto soft;
4652     }
4653 
4654     float64_input_flush1(&ua.s, s);
4655     if (QEMU_HARDFLOAT_1F64_USE_FP) {
4656         if (unlikely(!(fpclassify(ua.h) == FP_NORMAL ||
4657                        fpclassify(ua.h) == FP_ZERO) ||
4658                      signbit(ua.h))) {
4659             goto soft;
4660         }
4661     } else if (unlikely(!float64_is_zero_or_normal(ua.s) ||
4662                         float64_is_neg(ua.s))) {
4663         goto soft;
4664     }
4665     ur.h = sqrt(ua.h);
4666     return ur.s;
4667 
4668  soft:
4669     return soft_f64_sqrt(ua.s, s);
4670 }
4671 
4672 float64 float64r32_sqrt(float64 a, float_status *status)
4673 {
4674     FloatParts64 p;
4675 
4676     float64_unpack_canonical(&p, a, status);
4677     parts_sqrt(&p, status, &float64_params);
4678     return float64r32_round_pack_canonical(&p, status);
4679 }
4680 
4681 bfloat16 QEMU_FLATTEN bfloat16_sqrt(bfloat16 a, float_status *status)
4682 {
4683     FloatParts64 p;
4684 
4685     bfloat16_unpack_canonical(&p, a, status);
4686     parts_sqrt(&p, status, &bfloat16_params);
4687     return bfloat16_round_pack_canonical(&p, status);
4688 }
4689 
4690 float128 QEMU_FLATTEN float128_sqrt(float128 a, float_status *status)
4691 {
4692     FloatParts128 p;
4693 
4694     float128_unpack_canonical(&p, a, status);
4695     parts_sqrt(&p, status, &float128_params);
4696     return float128_round_pack_canonical(&p, status);
4697 }
4698 
4699 floatx80 floatx80_sqrt(floatx80 a, float_status *s)
4700 {
4701     FloatParts128 p;
4702 
4703     if (!floatx80_unpack_canonical(&p, a, s)) {
4704         return floatx80_default_nan(s);
4705     }
4706     parts_sqrt(&p, s, &floatx80_params[s->floatx80_rounding_precision]);
4707     return floatx80_round_pack_canonical(&p, s);
4708 }
4709 
4710 /*
4711  * log2
4712  */
4713 float32 float32_log2(float32 a, float_status *status)
4714 {
4715     FloatParts64 p;
4716 
4717     float32_unpack_canonical(&p, a, status);
4718     parts_log2(&p, status, &float32_params);
4719     return float32_round_pack_canonical(&p, status);
4720 }
4721 
4722 float64 float64_log2(float64 a, float_status *status)
4723 {
4724     FloatParts64 p;
4725 
4726     float64_unpack_canonical(&p, a, status);
4727     parts_log2(&p, status, &float64_params);
4728     return float64_round_pack_canonical(&p, status);
4729 }
4730 
4731 /*----------------------------------------------------------------------------
4732 | The pattern for a default generated NaN.
4733 *----------------------------------------------------------------------------*/
4734 
4735 float16 float16_default_nan(float_status *status)
4736 {
4737     FloatParts64 p;
4738 
4739     parts_default_nan(&p, status);
4740     p.frac >>= float16_params.frac_shift;
4741     return float16_pack_raw(&p);
4742 }
4743 
4744 float32 float32_default_nan(float_status *status)
4745 {
4746     FloatParts64 p;
4747 
4748     parts_default_nan(&p, status);
4749     p.frac >>= float32_params.frac_shift;
4750     return float32_pack_raw(&p);
4751 }
4752 
4753 float64 float64_default_nan(float_status *status)
4754 {
4755     FloatParts64 p;
4756 
4757     parts_default_nan(&p, status);
4758     p.frac >>= float64_params.frac_shift;
4759     return float64_pack_raw(&p);
4760 }
4761 
4762 float128 float128_default_nan(float_status *status)
4763 {
4764     FloatParts128 p;
4765 
4766     parts_default_nan(&p, status);
4767     frac_shr(&p, float128_params.frac_shift);
4768     return float128_pack_raw(&p);
4769 }
4770 
4771 bfloat16 bfloat16_default_nan(float_status *status)
4772 {
4773     FloatParts64 p;
4774 
4775     parts_default_nan(&p, status);
4776     p.frac >>= bfloat16_params.frac_shift;
4777     return bfloat16_pack_raw(&p);
4778 }
4779 
4780 /*----------------------------------------------------------------------------
4781 | Returns a quiet NaN from a signalling NaN for the floating point value `a'.
4782 *----------------------------------------------------------------------------*/
4783 
4784 float16 float16_silence_nan(float16 a, float_status *status)
4785 {
4786     FloatParts64 p;
4787 
4788     float16_unpack_raw(&p, a);
4789     p.frac <<= float16_params.frac_shift;
4790     parts_silence_nan(&p, status);
4791     p.frac >>= float16_params.frac_shift;
4792     return float16_pack_raw(&p);
4793 }
4794 
4795 float32 float32_silence_nan(float32 a, float_status *status)
4796 {
4797     FloatParts64 p;
4798 
4799     float32_unpack_raw(&p, a);
4800     p.frac <<= float32_params.frac_shift;
4801     parts_silence_nan(&p, status);
4802     p.frac >>= float32_params.frac_shift;
4803     return float32_pack_raw(&p);
4804 }
4805 
4806 float64 float64_silence_nan(float64 a, float_status *status)
4807 {
4808     FloatParts64 p;
4809 
4810     float64_unpack_raw(&p, a);
4811     p.frac <<= float64_params.frac_shift;
4812     parts_silence_nan(&p, status);
4813     p.frac >>= float64_params.frac_shift;
4814     return float64_pack_raw(&p);
4815 }
4816 
4817 bfloat16 bfloat16_silence_nan(bfloat16 a, float_status *status)
4818 {
4819     FloatParts64 p;
4820 
4821     bfloat16_unpack_raw(&p, a);
4822     p.frac <<= bfloat16_params.frac_shift;
4823     parts_silence_nan(&p, status);
4824     p.frac >>= bfloat16_params.frac_shift;
4825     return bfloat16_pack_raw(&p);
4826 }
4827 
4828 float128 float128_silence_nan(float128 a, float_status *status)
4829 {
4830     FloatParts128 p;
4831 
4832     float128_unpack_raw(&p, a);
4833     frac_shl(&p, float128_params.frac_shift);
4834     parts_silence_nan(&p, status);
4835     frac_shr(&p, float128_params.frac_shift);
4836     return float128_pack_raw(&p);
4837 }
4838 
4839 /*----------------------------------------------------------------------------
4840 | If `a' is denormal and we are in flush-to-zero mode then set the
4841 | input-denormal exception and return zero. Otherwise just return the value.
4842 *----------------------------------------------------------------------------*/
4843 
4844 static bool parts_squash_denormal(FloatParts64 p, float_status *status)
4845 {
4846     if (p.exp == 0 && p.frac != 0) {
4847         float_raise(float_flag_input_denormal, status);
4848         return true;
4849     }
4850 
4851     return false;
4852 }
4853 
4854 float16 float16_squash_input_denormal(float16 a, float_status *status)
4855 {
4856     if (status->flush_inputs_to_zero) {
4857         FloatParts64 p;
4858 
4859         float16_unpack_raw(&p, a);
4860         if (parts_squash_denormal(p, status)) {
4861             return float16_set_sign(float16_zero, p.sign);
4862         }
4863     }
4864     return a;
4865 }
4866 
4867 float32 float32_squash_input_denormal(float32 a, float_status *status)
4868 {
4869     if (status->flush_inputs_to_zero) {
4870         FloatParts64 p;
4871 
4872         float32_unpack_raw(&p, a);
4873         if (parts_squash_denormal(p, status)) {
4874             return float32_set_sign(float32_zero, p.sign);
4875         }
4876     }
4877     return a;
4878 }
4879 
4880 float64 float64_squash_input_denormal(float64 a, float_status *status)
4881 {
4882     if (status->flush_inputs_to_zero) {
4883         FloatParts64 p;
4884 
4885         float64_unpack_raw(&p, a);
4886         if (parts_squash_denormal(p, status)) {
4887             return float64_set_sign(float64_zero, p.sign);
4888         }
4889     }
4890     return a;
4891 }
4892 
4893 bfloat16 bfloat16_squash_input_denormal(bfloat16 a, float_status *status)
4894 {
4895     if (status->flush_inputs_to_zero) {
4896         FloatParts64 p;
4897 
4898         bfloat16_unpack_raw(&p, a);
4899         if (parts_squash_denormal(p, status)) {
4900             return bfloat16_set_sign(bfloat16_zero, p.sign);
4901         }
4902     }
4903     return a;
4904 }
4905 
4906 /*----------------------------------------------------------------------------
4907 | Normalizes the subnormal extended double-precision floating-point value
4908 | represented by the denormalized significand `aSig'.  The normalized exponent
4909 | and significand are stored at the locations pointed to by `zExpPtr' and
4910 | `zSigPtr', respectively.
4911 *----------------------------------------------------------------------------*/
4912 
4913 void normalizeFloatx80Subnormal(uint64_t aSig, int32_t *zExpPtr,
4914                                 uint64_t *zSigPtr)
4915 {
4916     int8_t shiftCount;
4917 
4918     shiftCount = clz64(aSig);
4919     *zSigPtr = aSig<<shiftCount;
4920     *zExpPtr = 1 - shiftCount;
4921 }
4922 
4923 /*----------------------------------------------------------------------------
4924 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4925 | and extended significand formed by the concatenation of `zSig0' and `zSig1',
4926 | and returns the proper extended double-precision floating-point value
4927 | corresponding to the abstract input.  Ordinarily, the abstract value is
4928 | rounded and packed into the extended double-precision format, with the
4929 | inexact exception raised if the abstract input cannot be represented
4930 | exactly.  However, if the abstract value is too large, the overflow and
4931 | inexact exceptions are raised and an infinity or maximal finite value is
4932 | returned.  If the abstract value is too small, the input value is rounded to
4933 | a subnormal number, and the underflow and inexact exceptions are raised if
4934 | the abstract input cannot be represented exactly as a subnormal extended
4935 | double-precision floating-point number.
4936 |     If `roundingPrecision' is floatx80_precision_s or floatx80_precision_d,
4937 | the result is rounded to the same number of bits as single or double
4938 | precision, respectively.  Otherwise, the result is rounded to the full
4939 | precision of the extended double-precision format.
4940 |     The input significand must be normalized or smaller.  If the input
4941 | significand is not normalized, `zExp' must be 0; in that case, the result
4942 | returned is a subnormal number, and it must not require rounding.  The
4943 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
4944 | Floating-Point Arithmetic.
4945 *----------------------------------------------------------------------------*/
4946 
4947 floatx80 roundAndPackFloatx80(FloatX80RoundPrec roundingPrecision, bool zSign,
4948                               int32_t zExp, uint64_t zSig0, uint64_t zSig1,
4949                               float_status *status)
4950 {
4951     FloatRoundMode roundingMode;
4952     bool roundNearestEven, increment, isTiny;
4953     int64_t roundIncrement, roundMask, roundBits;
4954 
4955     roundingMode = status->float_rounding_mode;
4956     roundNearestEven = ( roundingMode == float_round_nearest_even );
4957     switch (roundingPrecision) {
4958     case floatx80_precision_x:
4959         goto precision80;
4960     case floatx80_precision_d:
4961         roundIncrement = UINT64_C(0x0000000000000400);
4962         roundMask = UINT64_C(0x00000000000007FF);
4963         break;
4964     case floatx80_precision_s:
4965         roundIncrement = UINT64_C(0x0000008000000000);
4966         roundMask = UINT64_C(0x000000FFFFFFFFFF);
4967         break;
4968     default:
4969         g_assert_not_reached();
4970     }
4971     zSig0 |= ( zSig1 != 0 );
4972     switch (roundingMode) {
4973     case float_round_nearest_even:
4974     case float_round_ties_away:
4975         break;
4976     case float_round_to_zero:
4977         roundIncrement = 0;
4978         break;
4979     case float_round_up:
4980         roundIncrement = zSign ? 0 : roundMask;
4981         break;
4982     case float_round_down:
4983         roundIncrement = zSign ? roundMask : 0;
4984         break;
4985     default:
4986         abort();
4987     }
4988     roundBits = zSig0 & roundMask;
4989     if ( 0x7FFD <= (uint32_t) ( zExp - 1 ) ) {
4990         if (    ( 0x7FFE < zExp )
4991              || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
4992            ) {
4993             goto overflow;
4994         }
4995         if ( zExp <= 0 ) {
4996             if (status->flush_to_zero) {
4997                 float_raise(float_flag_output_denormal, status);
4998                 return packFloatx80(zSign, 0, 0);
4999             }
5000             isTiny = status->tininess_before_rounding
5001                   || (zExp < 0 )
5002                   || (zSig0 <= zSig0 + roundIncrement);
5003             shift64RightJamming( zSig0, 1 - zExp, &zSig0 );
5004             zExp = 0;
5005             roundBits = zSig0 & roundMask;
5006             if (isTiny && roundBits) {
5007                 float_raise(float_flag_underflow, status);
5008             }
5009             if (roundBits) {
5010                 float_raise(float_flag_inexact, status);
5011             }
5012             zSig0 += roundIncrement;
5013             if ( (int64_t) zSig0 < 0 ) zExp = 1;
5014             roundIncrement = roundMask + 1;
5015             if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
5016                 roundMask |= roundIncrement;
5017             }
5018             zSig0 &= ~ roundMask;
5019             return packFloatx80( zSign, zExp, zSig0 );
5020         }
5021     }
5022     if (roundBits) {
5023         float_raise(float_flag_inexact, status);
5024     }
5025     zSig0 += roundIncrement;
5026     if ( zSig0 < roundIncrement ) {
5027         ++zExp;
5028         zSig0 = UINT64_C(0x8000000000000000);
5029     }
5030     roundIncrement = roundMask + 1;
5031     if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
5032         roundMask |= roundIncrement;
5033     }
5034     zSig0 &= ~ roundMask;
5035     if ( zSig0 == 0 ) zExp = 0;
5036     return packFloatx80( zSign, zExp, zSig0 );
5037  precision80:
5038     switch (roundingMode) {
5039     case float_round_nearest_even:
5040     case float_round_ties_away:
5041         increment = ((int64_t)zSig1 < 0);
5042         break;
5043     case float_round_to_zero:
5044         increment = 0;
5045         break;
5046     case float_round_up:
5047         increment = !zSign && zSig1;
5048         break;
5049     case float_round_down:
5050         increment = zSign && zSig1;
5051         break;
5052     default:
5053         abort();
5054     }
5055     if ( 0x7FFD <= (uint32_t) ( zExp - 1 ) ) {
5056         if (    ( 0x7FFE < zExp )
5057              || (    ( zExp == 0x7FFE )
5058                   && ( zSig0 == UINT64_C(0xFFFFFFFFFFFFFFFF) )
5059                   && increment
5060                 )
5061            ) {
5062             roundMask = 0;
5063  overflow:
5064             float_raise(float_flag_overflow | float_flag_inexact, status);
5065             if (    ( roundingMode == float_round_to_zero )
5066                  || ( zSign && ( roundingMode == float_round_up ) )
5067                  || ( ! zSign && ( roundingMode == float_round_down ) )
5068                ) {
5069                 return packFloatx80( zSign, 0x7FFE, ~ roundMask );
5070             }
5071             return packFloatx80(zSign,
5072                                 floatx80_infinity_high,
5073                                 floatx80_infinity_low);
5074         }
5075         if ( zExp <= 0 ) {
5076             isTiny = status->tininess_before_rounding
5077                   || (zExp < 0)
5078                   || !increment
5079                   || (zSig0 < UINT64_C(0xFFFFFFFFFFFFFFFF));
5080             shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );
5081             zExp = 0;
5082             if (isTiny && zSig1) {
5083                 float_raise(float_flag_underflow, status);
5084             }
5085             if (zSig1) {
5086                 float_raise(float_flag_inexact, status);
5087             }
5088             switch (roundingMode) {
5089             case float_round_nearest_even:
5090             case float_round_ties_away:
5091                 increment = ((int64_t)zSig1 < 0);
5092                 break;
5093             case float_round_to_zero:
5094                 increment = 0;
5095                 break;
5096             case float_round_up:
5097                 increment = !zSign && zSig1;
5098                 break;
5099             case float_round_down:
5100                 increment = zSign && zSig1;
5101                 break;
5102             default:
5103                 abort();
5104             }
5105             if ( increment ) {
5106                 ++zSig0;
5107                 if (!(zSig1 << 1) && roundNearestEven) {
5108                     zSig0 &= ~1;
5109                 }
5110                 if ( (int64_t) zSig0 < 0 ) zExp = 1;
5111             }
5112             return packFloatx80( zSign, zExp, zSig0 );
5113         }
5114     }
5115     if (zSig1) {
5116         float_raise(float_flag_inexact, status);
5117     }
5118     if ( increment ) {
5119         ++zSig0;
5120         if ( zSig0 == 0 ) {
5121             ++zExp;
5122             zSig0 = UINT64_C(0x8000000000000000);
5123         }
5124         else {
5125             if (!(zSig1 << 1) && roundNearestEven) {
5126                 zSig0 &= ~1;
5127             }
5128         }
5129     }
5130     else {
5131         if ( zSig0 == 0 ) zExp = 0;
5132     }
5133     return packFloatx80( zSign, zExp, zSig0 );
5134 
5135 }
5136 
5137 /*----------------------------------------------------------------------------
5138 | Takes an abstract floating-point value having sign `zSign', exponent
5139 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
5140 | and returns the proper extended double-precision floating-point value
5141 | corresponding to the abstract input.  This routine is just like
5142 | `roundAndPackFloatx80' except that the input significand does not have to be
5143 | normalized.
5144 *----------------------------------------------------------------------------*/
5145 
5146 floatx80 normalizeRoundAndPackFloatx80(FloatX80RoundPrec roundingPrecision,
5147                                        bool zSign, int32_t zExp,
5148                                        uint64_t zSig0, uint64_t zSig1,
5149                                        float_status *status)
5150 {
5151     int8_t shiftCount;
5152 
5153     if ( zSig0 == 0 ) {
5154         zSig0 = zSig1;
5155         zSig1 = 0;
5156         zExp -= 64;
5157     }
5158     shiftCount = clz64(zSig0);
5159     shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
5160     zExp -= shiftCount;
5161     return roundAndPackFloatx80(roundingPrecision, zSign, zExp,
5162                                 zSig0, zSig1, status);
5163 
5164 }
5165 
5166 /*----------------------------------------------------------------------------
5167 | Returns the binary exponential of the single-precision floating-point value
5168 | `a'. The operation is performed according to the IEC/IEEE Standard for
5169 | Binary Floating-Point Arithmetic.
5170 |
5171 | Uses the following identities:
5172 |
5173 | 1. -------------------------------------------------------------------------
5174 |      x    x*ln(2)
5175 |     2  = e
5176 |
5177 | 2. -------------------------------------------------------------------------
5178 |                      2     3     4     5           n
5179 |      x        x     x     x     x     x           x
5180 |     e  = 1 + --- + --- + --- + --- + --- + ... + --- + ...
5181 |               1!    2!    3!    4!    5!          n!
5182 *----------------------------------------------------------------------------*/
5183 
5184 static const float64 float32_exp2_coefficients[15] =
5185 {
5186     const_float64( 0x3ff0000000000000ll ), /*  1 */
5187     const_float64( 0x3fe0000000000000ll ), /*  2 */
5188     const_float64( 0x3fc5555555555555ll ), /*  3 */
5189     const_float64( 0x3fa5555555555555ll ), /*  4 */
5190     const_float64( 0x3f81111111111111ll ), /*  5 */
5191     const_float64( 0x3f56c16c16c16c17ll ), /*  6 */
5192     const_float64( 0x3f2a01a01a01a01all ), /*  7 */
5193     const_float64( 0x3efa01a01a01a01all ), /*  8 */
5194     const_float64( 0x3ec71de3a556c734ll ), /*  9 */
5195     const_float64( 0x3e927e4fb7789f5cll ), /* 10 */
5196     const_float64( 0x3e5ae64567f544e4ll ), /* 11 */
5197     const_float64( 0x3e21eed8eff8d898ll ), /* 12 */
5198     const_float64( 0x3de6124613a86d09ll ), /* 13 */
5199     const_float64( 0x3da93974a8c07c9dll ), /* 14 */
5200     const_float64( 0x3d6ae7f3e733b81fll ), /* 15 */
5201 };
5202 
5203 float32 float32_exp2(float32 a, float_status *status)
5204 {
5205     FloatParts64 xp, xnp, tp, rp;
5206     int i;
5207 
5208     float32_unpack_canonical(&xp, a, status);
5209     if (unlikely(xp.cls != float_class_normal)) {
5210         switch (xp.cls) {
5211         case float_class_snan:
5212         case float_class_qnan:
5213             parts_return_nan(&xp, status);
5214             return float32_round_pack_canonical(&xp, status);
5215         case float_class_inf:
5216             return xp.sign ? float32_zero : a;
5217         case float_class_zero:
5218             return float32_one;
5219         default:
5220             break;
5221         }
5222         g_assert_not_reached();
5223     }
5224 
5225     float_raise(float_flag_inexact, status);
5226 
5227     float64_unpack_canonical(&tp, float64_ln2, status);
5228     xp = *parts_mul(&xp, &tp, status);
5229     xnp = xp;
5230 
5231     float64_unpack_canonical(&rp, float64_one, status);
5232     for (i = 0 ; i < 15 ; i++) {
5233         float64_unpack_canonical(&tp, float32_exp2_coefficients[i], status);
5234         rp = *parts_muladd(&tp, &xnp, &rp, 0, status);
5235         xnp = *parts_mul(&xnp, &xp, status);
5236     }
5237 
5238     return float32_round_pack_canonical(&rp, status);
5239 }
5240 
5241 /*----------------------------------------------------------------------------
5242 | Rounds the extended double-precision floating-point value `a'
5243 | to the precision provided by floatx80_rounding_precision and returns the
5244 | result as an extended double-precision floating-point value.
5245 | The operation is performed according to the IEC/IEEE Standard for Binary
5246 | Floating-Point Arithmetic.
5247 *----------------------------------------------------------------------------*/
5248 
5249 floatx80 floatx80_round(floatx80 a, float_status *status)
5250 {
5251     FloatParts128 p;
5252 
5253     if (!floatx80_unpack_canonical(&p, a, status)) {
5254         return floatx80_default_nan(status);
5255     }
5256     return floatx80_round_pack_canonical(&p, status);
5257 }
5258 
5259 static void __attribute__((constructor)) softfloat_init(void)
5260 {
5261     union_float64 ua, ub, uc, ur;
5262 
5263     if (QEMU_NO_HARDFLOAT) {
5264         return;
5265     }
5266     /*
5267      * Test that the host's FMA is not obviously broken. For example,
5268      * glibc < 2.23 can perform an incorrect FMA on certain hosts; see
5269      *   https://sourceware.org/bugzilla/show_bug.cgi?id=13304
5270      */
5271     ua.s = 0x0020000000000001ULL;
5272     ub.s = 0x3ca0000000000000ULL;
5273     uc.s = 0x0020000000000000ULL;
5274     ur.h = fma(ua.h, ub.h, uc.h);
5275     if (ur.s != 0x0020000000000001ULL) {
5276         force_soft_fma = true;
5277     }
5278 }
5279