virtio-mem: Implement replay_discarded RamDiscardManager callback
[qemu/kevin.git] / fpu / softfloat.c
blob9a28720d82a5bc87c5fc05677695aad1674854b0
1 /*
2 * QEMU float support
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
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.
19 ===============================================================================
20 This C source file is part of the SoftFloat IEC/IEEE Floating-point
21 Arithmetic Package, Release 2a.
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'.
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.
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.
44 ===============================================================================
47 /* BSD licensing:
48 * Copyright (c) 2006, Fabrice Bellard
49 * All rights reserved.
51 * Redistribution and use in source and binary forms, with or without
52 * modification, are permitted provided that the following conditions are met:
54 * 1. Redistributions of source code must retain the above copyright notice,
55 * this list of conditions and the following disclaimer.
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.
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.
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.
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.
82 /* softfloat (and in particular the code in softfloat-specialize.h) is
83 * target-dependent and needs the TARGET_* macros.
85 #include "qemu/osdep.h"
86 #include <math.h>
87 #include "qemu/bitops.h"
88 #include "fpu/softfloat.h"
90 /* We only need stdlib for abort() */
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"
100 * Hardfloat
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.
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:
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.
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.
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.
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.
129 #define GEN_INPUT_FLUSH__NOCHECK(name, soft_t) \
130 static inline void name(soft_t *a, float_status *s) \
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); \
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
143 #define GEN_INPUT_FLUSH1(name, soft_t) \
144 static inline void name(soft_t *a, float_status *s) \
146 if (likely(!s->flush_inputs_to_zero)) { \
147 return; \
149 soft_t ## _input_flush__nocheck(a, s); \
152 GEN_INPUT_FLUSH1(float32_input_flush1, float32)
153 GEN_INPUT_FLUSH1(float64_input_flush1, float64)
154 #undef GEN_INPUT_FLUSH1
156 #define GEN_INPUT_FLUSH2(name, soft_t) \
157 static inline void name(soft_t *a, soft_t *b, float_status *s) \
159 if (likely(!s->flush_inputs_to_zero)) { \
160 return; \
162 soft_t ## _input_flush__nocheck(a, s); \
163 soft_t ## _input_flush__nocheck(b, s); \
166 GEN_INPUT_FLUSH2(float32_input_flush2, float32)
167 GEN_INPUT_FLUSH2(float64_input_flush2, float64)
168 #undef GEN_INPUT_FLUSH2
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) \
173 if (likely(!s->flush_inputs_to_zero)) { \
174 return; \
176 soft_t ## _input_flush__nocheck(a, s); \
177 soft_t ## _input_flush__nocheck(b, s); \
178 soft_t ## _input_flush__nocheck(c, s); \
181 GEN_INPUT_FLUSH3(float32_input_flush3, float32)
182 GEN_INPUT_FLUSH3(float64_input_flush3, float64)
183 #undef GEN_INPUT_FLUSH3
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.
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
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%.
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
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.
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
235 static inline bool can_use_fpu(const float_status *s)
237 if (QEMU_NO_HARDFLOAT) {
238 return false;
240 return likely(s->float_exception_flags & float_flag_inexact &&
241 s->float_rounding_mode == float_round_nearest_even);
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).
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.
252 * We only generate functions for operations with two inputs, since only
253 * these are common enough to justify consolidating them into common code.
256 typedef union {
257 float32 s;
258 float h;
259 } union_float32;
261 typedef union {
262 float64 s;
263 double h;
264 } union_float64;
266 typedef bool (*f32_check_fn)(union_float32 a, union_float32 b);
267 typedef bool (*f64_check_fn)(union_float64 a, union_float64 b);
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);
274 /* 2-input is-zero-or-normal */
275 static inline bool f32_is_zon2(union_float32 a, union_float32 b)
277 if (QEMU_HARDFLOAT_2F32_USE_FP) {
279 * Not using a temp variable for consecutive fpclassify calls ends up
280 * generating faster code.
282 return (fpclassify(a.h) == FP_NORMAL || fpclassify(a.h) == FP_ZERO) &&
283 (fpclassify(b.h) == FP_NORMAL || fpclassify(b.h) == FP_ZERO);
285 return float32_is_zero_or_normal(a.s) &&
286 float32_is_zero_or_normal(b.s);
289 static inline bool f64_is_zon2(union_float64 a, union_float64 b)
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);
295 return float64_is_zero_or_normal(a.s) &&
296 float64_is_zero_or_normal(b.s);
299 /* 3-input is-zero-or-normal */
300 static inline
301 bool f32_is_zon3(union_float32 a, union_float32 b, union_float32 c)
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);
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);
313 static inline
314 bool f64_is_zon3(union_float64 a, union_float64 b, union_float64 c)
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);
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);
326 static inline bool f32_is_inf(union_float32 a)
328 if (QEMU_HARDFLOAT_USE_ISINF) {
329 return isinf(a.h);
331 return float32_is_infinity(a.s);
334 static inline bool f64_is_inf(union_float64 a)
336 if (QEMU_HARDFLOAT_USE_ISINF) {
337 return isinf(a.h);
339 return float64_is_infinity(a.s);
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)
347 union_float32 ua, ub, ur;
349 ua.s = xa;
350 ub.s = xb;
352 if (unlikely(!can_use_fpu(s))) {
353 goto soft;
356 float32_input_flush2(&ua.s, &ub.s, s);
357 if (unlikely(!pre(ua, ub))) {
358 goto soft;
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;
367 return ur.s;
369 soft:
370 return soft(ua.s, ub.s, s);
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)
378 union_float64 ua, ub, ur;
380 ua.s = xa;
381 ub.s = xb;
383 if (unlikely(!can_use_fpu(s))) {
384 goto soft;
387 float64_input_flush2(&ua.s, &ub.s, s);
388 if (unlikely(!pre(ua, ub))) {
389 goto soft;
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;
398 return ur.s;
400 soft:
401 return soft(ua.s, ub.s, s);
405 * Classify a floating point number. Everything above float_class_qnan
406 * is a NaN so cls >= float_class_qnan is any NaN.
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;
418 #define float_cmask(bit) (1u << (bit))
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),
427 float_cmask_infzero = float_cmask_zero | float_cmask_inf,
428 float_cmask_anynan = float_cmask_qnan | float_cmask_snan,
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,
440 * Set for the IEEE 754-2019 minimumNumber() and maximumNumber()
441 * operations.
443 minmax_isnumber = 8,
446 /* Simple helpers for checking if, or what kind of, NaN we have */
447 static inline __attribute__((unused)) bool is_nan(FloatClass c)
449 return unlikely(c >= float_class_qnan);
452 static inline __attribute__((unused)) bool is_snan(FloatClass c)
454 return c == float_class_snan;
457 static inline __attribute__((unused)) bool is_qnan(FloatClass c)
459 return c == float_class_qnan;
463 * Structure holding all of the decomposed parts of a float.
464 * The exponent is unbiased and the fraction is normalized.
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.
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;
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.
483 uint64_t frac_hi;
484 uint64_t frac_lo;
486 } FloatParts64;
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;
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;
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)
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
521 typedef struct {
522 int exp_size;
523 int exp_bias;
524 int exp_max;
525 int frac_size;
526 int frac_shift;
527 bool arm_althp;
528 uint64_t round_mask;
529 } FloatFmt;
531 /* Expand fields based on the size of exponent and fraction */
532 #define FLOAT_PARAMS_(E) \
533 .exp_size = E, \
534 .exp_bias = ((1 << E) - 1) >> 1, \
535 .exp_max = (1 << E) - 1
537 #define FLOAT_PARAMS(E, F) \
538 FLOAT_PARAMS_(E), \
539 .frac_size = F, \
540 .frac_shift = (-F - 1) & 63, \
541 .round_mask = (1ull << ((-F - 1) & 63)) - 1
543 static const FloatFmt float16_params = {
544 FLOAT_PARAMS(5, 10)
547 static const FloatFmt float16_params_ahp = {
548 FLOAT_PARAMS(5, 10),
549 .arm_althp = true
552 static const FloatFmt bfloat16_params = {
553 FLOAT_PARAMS(8, 7)
556 static const FloatFmt float32_params = {
557 FLOAT_PARAMS(8, 23)
560 static const FloatFmt float64_params = {
561 FLOAT_PARAMS(11, 52)
564 static const FloatFmt float128_params = {
565 FLOAT_PARAMS(15, 112)
568 #define FLOATX80_PARAMS(R) \
569 FLOAT_PARAMS_(15), \
570 .frac_size = R == 64 ? 63 : R, \
571 .frac_shift = 0, \
572 .round_mask = R == 64 ? -1 : (1ull << ((-R - 1) & 63)) - 1
574 static const FloatFmt floatx80_params[3] = {
575 [floatx80_precision_s] = { FLOATX80_PARAMS(23) },
576 [floatx80_precision_d] = { FLOATX80_PARAMS(52) },
577 [floatx80_precision_x] = { FLOATX80_PARAMS(64) },
580 /* Unpack a float to parts, but do not canonicalize. */
581 static void unpack_raw64(FloatParts64 *r, const FloatFmt *fmt, uint64_t raw)
583 const int f_size = fmt->frac_size;
584 const int e_size = fmt->exp_size;
586 *r = (FloatParts64) {
587 .cls = float_class_unclassified,
588 .sign = extract64(raw, f_size + e_size, 1),
589 .exp = extract64(raw, f_size, e_size),
590 .frac = extract64(raw, 0, f_size)
594 static inline void float16_unpack_raw(FloatParts64 *p, float16 f)
596 unpack_raw64(p, &float16_params, f);
599 static inline void bfloat16_unpack_raw(FloatParts64 *p, bfloat16 f)
601 unpack_raw64(p, &bfloat16_params, f);
604 static inline void float32_unpack_raw(FloatParts64 *p, float32 f)
606 unpack_raw64(p, &float32_params, f);
609 static inline void float64_unpack_raw(FloatParts64 *p, float64 f)
611 unpack_raw64(p, &float64_params, f);
614 static void floatx80_unpack_raw(FloatParts128 *p, floatx80 f)
616 *p = (FloatParts128) {
617 .cls = float_class_unclassified,
618 .sign = extract32(f.high, 15, 1),
619 .exp = extract32(f.high, 0, 15),
620 .frac_hi = f.low
624 static void float128_unpack_raw(FloatParts128 *p, float128 f)
626 const int f_size = float128_params.frac_size - 64;
627 const int e_size = float128_params.exp_size;
629 *p = (FloatParts128) {
630 .cls = float_class_unclassified,
631 .sign = extract64(f.high, f_size + e_size, 1),
632 .exp = extract64(f.high, f_size, e_size),
633 .frac_hi = extract64(f.high, 0, f_size),
634 .frac_lo = f.low,
638 /* Pack a float from parts, but do not canonicalize. */
639 static uint64_t pack_raw64(const FloatParts64 *p, const FloatFmt *fmt)
641 const int f_size = fmt->frac_size;
642 const int e_size = fmt->exp_size;
643 uint64_t ret;
645 ret = (uint64_t)p->sign << (f_size + e_size);
646 ret = deposit64(ret, f_size, e_size, p->exp);
647 ret = deposit64(ret, 0, f_size, p->frac);
648 return ret;
651 static inline float16 float16_pack_raw(const FloatParts64 *p)
653 return make_float16(pack_raw64(p, &float16_params));
656 static inline bfloat16 bfloat16_pack_raw(const FloatParts64 *p)
658 return pack_raw64(p, &bfloat16_params);
661 static inline float32 float32_pack_raw(const FloatParts64 *p)
663 return make_float32(pack_raw64(p, &float32_params));
666 static inline float64 float64_pack_raw(const FloatParts64 *p)
668 return make_float64(pack_raw64(p, &float64_params));
671 static float128 float128_pack_raw(const FloatParts128 *p)
673 const int f_size = float128_params.frac_size - 64;
674 const int e_size = float128_params.exp_size;
675 uint64_t hi;
677 hi = (uint64_t)p->sign << (f_size + e_size);
678 hi = deposit64(hi, f_size, e_size, p->exp);
679 hi = deposit64(hi, 0, f_size, p->frac_hi);
680 return make_float128(hi, p->frac_lo);
683 /*----------------------------------------------------------------------------
684 | Functions and definitions to determine: (1) whether tininess for underflow
685 | is detected before or after rounding by default, (2) what (if anything)
686 | happens when exceptions are raised, (3) how signaling NaNs are distinguished
687 | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
688 | are propagated from function inputs to output. These details are target-
689 | specific.
690 *----------------------------------------------------------------------------*/
691 #include "softfloat-specialize.c.inc"
693 #define PARTS_GENERIC_64_128(NAME, P) \
694 _Generic((P), FloatParts64 *: parts64_##NAME, \
695 FloatParts128 *: parts128_##NAME)
697 #define PARTS_GENERIC_64_128_256(NAME, P) \
698 _Generic((P), FloatParts64 *: parts64_##NAME, \
699 FloatParts128 *: parts128_##NAME, \
700 FloatParts256 *: parts256_##NAME)
702 #define parts_default_nan(P, S) PARTS_GENERIC_64_128(default_nan, P)(P, S)
703 #define parts_silence_nan(P, S) PARTS_GENERIC_64_128(silence_nan, P)(P, S)
705 static void parts64_return_nan(FloatParts64 *a, float_status *s);
706 static void parts128_return_nan(FloatParts128 *a, float_status *s);
708 #define parts_return_nan(P, S) PARTS_GENERIC_64_128(return_nan, P)(P, S)
710 static FloatParts64 *parts64_pick_nan(FloatParts64 *a, FloatParts64 *b,
711 float_status *s);
712 static FloatParts128 *parts128_pick_nan(FloatParts128 *a, FloatParts128 *b,
713 float_status *s);
715 #define parts_pick_nan(A, B, S) PARTS_GENERIC_64_128(pick_nan, A)(A, B, S)
717 static FloatParts64 *parts64_pick_nan_muladd(FloatParts64 *a, FloatParts64 *b,
718 FloatParts64 *c, float_status *s,
719 int ab_mask, int abc_mask);
720 static FloatParts128 *parts128_pick_nan_muladd(FloatParts128 *a,
721 FloatParts128 *b,
722 FloatParts128 *c,
723 float_status *s,
724 int ab_mask, int abc_mask);
726 #define parts_pick_nan_muladd(A, B, C, S, ABM, ABCM) \
727 PARTS_GENERIC_64_128(pick_nan_muladd, A)(A, B, C, S, ABM, ABCM)
729 static void parts64_canonicalize(FloatParts64 *p, float_status *status,
730 const FloatFmt *fmt);
731 static void parts128_canonicalize(FloatParts128 *p, float_status *status,
732 const FloatFmt *fmt);
734 #define parts_canonicalize(A, S, F) \
735 PARTS_GENERIC_64_128(canonicalize, A)(A, S, F)
737 static void parts64_uncanon_normal(FloatParts64 *p, float_status *status,
738 const FloatFmt *fmt);
739 static void parts128_uncanon_normal(FloatParts128 *p, float_status *status,
740 const FloatFmt *fmt);
742 #define parts_uncanon_normal(A, S, F) \
743 PARTS_GENERIC_64_128(uncanon_normal, A)(A, S, F)
745 static void parts64_uncanon(FloatParts64 *p, float_status *status,
746 const FloatFmt *fmt);
747 static void parts128_uncanon(FloatParts128 *p, float_status *status,
748 const FloatFmt *fmt);
750 #define parts_uncanon(A, S, F) \
751 PARTS_GENERIC_64_128(uncanon, A)(A, S, F)
753 static void parts64_add_normal(FloatParts64 *a, FloatParts64 *b);
754 static void parts128_add_normal(FloatParts128 *a, FloatParts128 *b);
755 static void parts256_add_normal(FloatParts256 *a, FloatParts256 *b);
757 #define parts_add_normal(A, B) \
758 PARTS_GENERIC_64_128_256(add_normal, A)(A, B)
760 static bool parts64_sub_normal(FloatParts64 *a, FloatParts64 *b);
761 static bool parts128_sub_normal(FloatParts128 *a, FloatParts128 *b);
762 static bool parts256_sub_normal(FloatParts256 *a, FloatParts256 *b);
764 #define parts_sub_normal(A, B) \
765 PARTS_GENERIC_64_128_256(sub_normal, A)(A, B)
767 static FloatParts64 *parts64_addsub(FloatParts64 *a, FloatParts64 *b,
768 float_status *s, bool subtract);
769 static FloatParts128 *parts128_addsub(FloatParts128 *a, FloatParts128 *b,
770 float_status *s, bool subtract);
772 #define parts_addsub(A, B, S, Z) \
773 PARTS_GENERIC_64_128(addsub, A)(A, B, S, Z)
775 static FloatParts64 *parts64_mul(FloatParts64 *a, FloatParts64 *b,
776 float_status *s);
777 static FloatParts128 *parts128_mul(FloatParts128 *a, FloatParts128 *b,
778 float_status *s);
780 #define parts_mul(A, B, S) \
781 PARTS_GENERIC_64_128(mul, A)(A, B, S)
783 static FloatParts64 *parts64_muladd(FloatParts64 *a, FloatParts64 *b,
784 FloatParts64 *c, int flags,
785 float_status *s);
786 static FloatParts128 *parts128_muladd(FloatParts128 *a, FloatParts128 *b,
787 FloatParts128 *c, int flags,
788 float_status *s);
790 #define parts_muladd(A, B, C, Z, S) \
791 PARTS_GENERIC_64_128(muladd, A)(A, B, C, Z, S)
793 static FloatParts64 *parts64_div(FloatParts64 *a, FloatParts64 *b,
794 float_status *s);
795 static FloatParts128 *parts128_div(FloatParts128 *a, FloatParts128 *b,
796 float_status *s);
798 #define parts_div(A, B, S) \
799 PARTS_GENERIC_64_128(div, A)(A, B, S)
801 static FloatParts64 *parts64_modrem(FloatParts64 *a, FloatParts64 *b,
802 uint64_t *mod_quot, float_status *s);
803 static FloatParts128 *parts128_modrem(FloatParts128 *a, FloatParts128 *b,
804 uint64_t *mod_quot, float_status *s);
806 #define parts_modrem(A, B, Q, S) \
807 PARTS_GENERIC_64_128(modrem, A)(A, B, Q, S)
809 static void parts64_sqrt(FloatParts64 *a, float_status *s, const FloatFmt *f);
810 static void parts128_sqrt(FloatParts128 *a, float_status *s, const FloatFmt *f);
812 #define parts_sqrt(A, S, F) \
813 PARTS_GENERIC_64_128(sqrt, A)(A, S, F)
815 static bool parts64_round_to_int_normal(FloatParts64 *a, FloatRoundMode rm,
816 int scale, int frac_size);
817 static bool parts128_round_to_int_normal(FloatParts128 *a, FloatRoundMode r,
818 int scale, int frac_size);
820 #define parts_round_to_int_normal(A, R, C, F) \
821 PARTS_GENERIC_64_128(round_to_int_normal, A)(A, R, C, F)
823 static void parts64_round_to_int(FloatParts64 *a, FloatRoundMode rm,
824 int scale, float_status *s,
825 const FloatFmt *fmt);
826 static void parts128_round_to_int(FloatParts128 *a, FloatRoundMode r,
827 int scale, float_status *s,
828 const FloatFmt *fmt);
830 #define parts_round_to_int(A, R, C, S, F) \
831 PARTS_GENERIC_64_128(round_to_int, A)(A, R, C, S, F)
833 static int64_t parts64_float_to_sint(FloatParts64 *p, FloatRoundMode rmode,
834 int scale, int64_t min, int64_t max,
835 float_status *s);
836 static int64_t parts128_float_to_sint(FloatParts128 *p, FloatRoundMode rmode,
837 int scale, int64_t min, int64_t max,
838 float_status *s);
840 #define parts_float_to_sint(P, R, Z, MN, MX, S) \
841 PARTS_GENERIC_64_128(float_to_sint, P)(P, R, Z, MN, MX, S)
843 static uint64_t parts64_float_to_uint(FloatParts64 *p, FloatRoundMode rmode,
844 int scale, uint64_t max,
845 float_status *s);
846 static uint64_t parts128_float_to_uint(FloatParts128 *p, FloatRoundMode rmode,
847 int scale, uint64_t max,
848 float_status *s);
850 #define parts_float_to_uint(P, R, Z, M, S) \
851 PARTS_GENERIC_64_128(float_to_uint, P)(P, R, Z, M, S)
853 static void parts64_sint_to_float(FloatParts64 *p, int64_t a,
854 int scale, float_status *s);
855 static void parts128_sint_to_float(FloatParts128 *p, int64_t a,
856 int scale, float_status *s);
858 #define parts_sint_to_float(P, I, Z, S) \
859 PARTS_GENERIC_64_128(sint_to_float, P)(P, I, Z, S)
861 static void parts64_uint_to_float(FloatParts64 *p, uint64_t a,
862 int scale, float_status *s);
863 static void parts128_uint_to_float(FloatParts128 *p, uint64_t a,
864 int scale, float_status *s);
866 #define parts_uint_to_float(P, I, Z, S) \
867 PARTS_GENERIC_64_128(uint_to_float, P)(P, I, Z, S)
869 static FloatParts64 *parts64_minmax(FloatParts64 *a, FloatParts64 *b,
870 float_status *s, int flags);
871 static FloatParts128 *parts128_minmax(FloatParts128 *a, FloatParts128 *b,
872 float_status *s, int flags);
874 #define parts_minmax(A, B, S, F) \
875 PARTS_GENERIC_64_128(minmax, A)(A, B, S, F)
877 static int parts64_compare(FloatParts64 *a, FloatParts64 *b,
878 float_status *s, bool q);
879 static int parts128_compare(FloatParts128 *a, FloatParts128 *b,
880 float_status *s, bool q);
882 #define parts_compare(A, B, S, Q) \
883 PARTS_GENERIC_64_128(compare, A)(A, B, S, Q)
885 static void parts64_scalbn(FloatParts64 *a, int n, float_status *s);
886 static void parts128_scalbn(FloatParts128 *a, int n, float_status *s);
888 #define parts_scalbn(A, N, S) \
889 PARTS_GENERIC_64_128(scalbn, A)(A, N, S)
891 static void parts64_log2(FloatParts64 *a, float_status *s, const FloatFmt *f);
892 static void parts128_log2(FloatParts128 *a, float_status *s, const FloatFmt *f);
894 #define parts_log2(A, S, F) \
895 PARTS_GENERIC_64_128(log2, A)(A, S, F)
898 * Helper functions for softfloat-parts.c.inc, per-size operations.
901 #define FRAC_GENERIC_64_128(NAME, P) \
902 _Generic((P), FloatParts64 *: frac64_##NAME, \
903 FloatParts128 *: frac128_##NAME)
905 #define FRAC_GENERIC_64_128_256(NAME, P) \
906 _Generic((P), FloatParts64 *: frac64_##NAME, \
907 FloatParts128 *: frac128_##NAME, \
908 FloatParts256 *: frac256_##NAME)
910 static bool frac64_add(FloatParts64 *r, FloatParts64 *a, FloatParts64 *b)
912 return uadd64_overflow(a->frac, b->frac, &r->frac);
915 static bool frac128_add(FloatParts128 *r, FloatParts128 *a, FloatParts128 *b)
917 bool c = 0;
918 r->frac_lo = uadd64_carry(a->frac_lo, b->frac_lo, &c);
919 r->frac_hi = uadd64_carry(a->frac_hi, b->frac_hi, &c);
920 return c;
923 static bool frac256_add(FloatParts256 *r, FloatParts256 *a, FloatParts256 *b)
925 bool c = 0;
926 r->frac_lo = uadd64_carry(a->frac_lo, b->frac_lo, &c);
927 r->frac_lm = uadd64_carry(a->frac_lm, b->frac_lm, &c);
928 r->frac_hm = uadd64_carry(a->frac_hm, b->frac_hm, &c);
929 r->frac_hi = uadd64_carry(a->frac_hi, b->frac_hi, &c);
930 return c;
933 #define frac_add(R, A, B) FRAC_GENERIC_64_128_256(add, R)(R, A, B)
935 static bool frac64_addi(FloatParts64 *r, FloatParts64 *a, uint64_t c)
937 return uadd64_overflow(a->frac, c, &r->frac);
940 static bool frac128_addi(FloatParts128 *r, FloatParts128 *a, uint64_t c)
942 c = uadd64_overflow(a->frac_lo, c, &r->frac_lo);
943 return uadd64_overflow(a->frac_hi, c, &r->frac_hi);
946 #define frac_addi(R, A, C) FRAC_GENERIC_64_128(addi, R)(R, A, C)
948 static void frac64_allones(FloatParts64 *a)
950 a->frac = -1;
953 static void frac128_allones(FloatParts128 *a)
955 a->frac_hi = a->frac_lo = -1;
958 #define frac_allones(A) FRAC_GENERIC_64_128(allones, A)(A)
960 static int frac64_cmp(FloatParts64 *a, FloatParts64 *b)
962 return a->frac == b->frac ? 0 : a->frac < b->frac ? -1 : 1;
965 static int frac128_cmp(FloatParts128 *a, FloatParts128 *b)
967 uint64_t ta = a->frac_hi, tb = b->frac_hi;
968 if (ta == tb) {
969 ta = a->frac_lo, tb = b->frac_lo;
970 if (ta == tb) {
971 return 0;
974 return ta < tb ? -1 : 1;
977 #define frac_cmp(A, B) FRAC_GENERIC_64_128(cmp, A)(A, B)
979 static void frac64_clear(FloatParts64 *a)
981 a->frac = 0;
984 static void frac128_clear(FloatParts128 *a)
986 a->frac_hi = a->frac_lo = 0;
989 #define frac_clear(A) FRAC_GENERIC_64_128(clear, A)(A)
991 static bool frac64_div(FloatParts64 *a, FloatParts64 *b)
993 uint64_t n1, n0, r, q;
994 bool ret;
997 * We want a 2*N / N-bit division to produce exactly an N-bit
998 * result, so that we do not lose any precision and so that we
999 * do not have to renormalize afterward. If A.frac < B.frac,
1000 * then division would produce an (N-1)-bit result; shift A left
1001 * by one to produce the an N-bit result, and return true to
1002 * decrement the exponent to match.
1004 * The udiv_qrnnd algorithm that we're using requires normalization,
1005 * i.e. the msb of the denominator must be set, which is already true.
1007 ret = a->frac < b->frac;
1008 if (ret) {
1009 n0 = a->frac;
1010 n1 = 0;
1011 } else {
1012 n0 = a->frac >> 1;
1013 n1 = a->frac << 63;
1015 q = udiv_qrnnd(&r, n0, n1, b->frac);
1017 /* Set lsb if there is a remainder, to set inexact. */
1018 a->frac = q | (r != 0);
1020 return ret;
1023 static bool frac128_div(FloatParts128 *a, FloatParts128 *b)
1025 uint64_t q0, q1, a0, a1, b0, b1;
1026 uint64_t r0, r1, r2, r3, t0, t1, t2, t3;
1027 bool ret = false;
1029 a0 = a->frac_hi, a1 = a->frac_lo;
1030 b0 = b->frac_hi, b1 = b->frac_lo;
1032 ret = lt128(a0, a1, b0, b1);
1033 if (!ret) {
1034 a1 = shr_double(a0, a1, 1);
1035 a0 = a0 >> 1;
1038 /* Use 128/64 -> 64 division as estimate for 192/128 -> 128 division. */
1039 q0 = estimateDiv128To64(a0, a1, b0);
1042 * Estimate is high because B1 was not included (unless B1 == 0).
1043 * Reduce quotient and increase remainder until remainder is non-negative.
1044 * This loop will execute 0 to 2 times.
1046 mul128By64To192(b0, b1, q0, &t0, &t1, &t2);
1047 sub192(a0, a1, 0, t0, t1, t2, &r0, &r1, &r2);
1048 while (r0 != 0) {
1049 q0--;
1050 add192(r0, r1, r2, 0, b0, b1, &r0, &r1, &r2);
1053 /* Repeat using the remainder, producing a second word of quotient. */
1054 q1 = estimateDiv128To64(r1, r2, b0);
1055 mul128By64To192(b0, b1, q1, &t1, &t2, &t3);
1056 sub192(r1, r2, 0, t1, t2, t3, &r1, &r2, &r3);
1057 while (r1 != 0) {
1058 q1--;
1059 add192(r1, r2, r3, 0, b0, b1, &r1, &r2, &r3);
1062 /* Any remainder indicates inexact; set sticky bit. */
1063 q1 |= (r2 | r3) != 0;
1065 a->frac_hi = q0;
1066 a->frac_lo = q1;
1067 return ret;
1070 #define frac_div(A, B) FRAC_GENERIC_64_128(div, A)(A, B)
1072 static bool frac64_eqz(FloatParts64 *a)
1074 return a->frac == 0;
1077 static bool frac128_eqz(FloatParts128 *a)
1079 return (a->frac_hi | a->frac_lo) == 0;
1082 #define frac_eqz(A) FRAC_GENERIC_64_128(eqz, A)(A)
1084 static void frac64_mulw(FloatParts128 *r, FloatParts64 *a, FloatParts64 *b)
1086 mulu64(&r->frac_lo, &r->frac_hi, a->frac, b->frac);
1089 static void frac128_mulw(FloatParts256 *r, FloatParts128 *a, FloatParts128 *b)
1091 mul128To256(a->frac_hi, a->frac_lo, b->frac_hi, b->frac_lo,
1092 &r->frac_hi, &r->frac_hm, &r->frac_lm, &r->frac_lo);
1095 #define frac_mulw(R, A, B) FRAC_GENERIC_64_128(mulw, A)(R, A, B)
1097 static void frac64_neg(FloatParts64 *a)
1099 a->frac = -a->frac;
1102 static void frac128_neg(FloatParts128 *a)
1104 bool c = 0;
1105 a->frac_lo = usub64_borrow(0, a->frac_lo, &c);
1106 a->frac_hi = usub64_borrow(0, a->frac_hi, &c);
1109 static void frac256_neg(FloatParts256 *a)
1111 bool c = 0;
1112 a->frac_lo = usub64_borrow(0, a->frac_lo, &c);
1113 a->frac_lm = usub64_borrow(0, a->frac_lm, &c);
1114 a->frac_hm = usub64_borrow(0, a->frac_hm, &c);
1115 a->frac_hi = usub64_borrow(0, a->frac_hi, &c);
1118 #define frac_neg(A) FRAC_GENERIC_64_128_256(neg, A)(A)
1120 static int frac64_normalize(FloatParts64 *a)
1122 if (a->frac) {
1123 int shift = clz64(a->frac);
1124 a->frac <<= shift;
1125 return shift;
1127 return 64;
1130 static int frac128_normalize(FloatParts128 *a)
1132 if (a->frac_hi) {
1133 int shl = clz64(a->frac_hi);
1134 a->frac_hi = shl_double(a->frac_hi, a->frac_lo, shl);
1135 a->frac_lo <<= shl;
1136 return shl;
1137 } else if (a->frac_lo) {
1138 int shl = clz64(a->frac_lo);
1139 a->frac_hi = a->frac_lo << shl;
1140 a->frac_lo = 0;
1141 return shl + 64;
1143 return 128;
1146 static int frac256_normalize(FloatParts256 *a)
1148 uint64_t a0 = a->frac_hi, a1 = a->frac_hm;
1149 uint64_t a2 = a->frac_lm, a3 = a->frac_lo;
1150 int ret, shl;
1152 if (likely(a0)) {
1153 shl = clz64(a0);
1154 if (shl == 0) {
1155 return 0;
1157 ret = shl;
1158 } else {
1159 if (a1) {
1160 ret = 64;
1161 a0 = a1, a1 = a2, a2 = a3, a3 = 0;
1162 } else if (a2) {
1163 ret = 128;
1164 a0 = a2, a1 = a3, a2 = 0, a3 = 0;
1165 } else if (a3) {
1166 ret = 192;
1167 a0 = a3, a1 = 0, a2 = 0, a3 = 0;
1168 } else {
1169 ret = 256;
1170 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
1171 goto done;
1173 shl = clz64(a0);
1174 if (shl == 0) {
1175 goto done;
1177 ret += shl;
1180 a0 = shl_double(a0, a1, shl);
1181 a1 = shl_double(a1, a2, shl);
1182 a2 = shl_double(a2, a3, shl);
1183 a3 <<= shl;
1185 done:
1186 a->frac_hi = a0;
1187 a->frac_hm = a1;
1188 a->frac_lm = a2;
1189 a->frac_lo = a3;
1190 return ret;
1193 #define frac_normalize(A) FRAC_GENERIC_64_128_256(normalize, A)(A)
1195 static void frac64_modrem(FloatParts64 *a, FloatParts64 *b, uint64_t *mod_quot)
1197 uint64_t a0, a1, b0, t0, t1, q, quot;
1198 int exp_diff = a->exp - b->exp;
1199 int shift;
1201 a0 = a->frac;
1202 a1 = 0;
1204 if (exp_diff < -1) {
1205 if (mod_quot) {
1206 *mod_quot = 0;
1208 return;
1210 if (exp_diff == -1) {
1211 a0 >>= 1;
1212 exp_diff = 0;
1215 b0 = b->frac;
1216 quot = q = b0 <= a0;
1217 if (q) {
1218 a0 -= b0;
1221 exp_diff -= 64;
1222 while (exp_diff > 0) {
1223 q = estimateDiv128To64(a0, a1, b0);
1224 q = q > 2 ? q - 2 : 0;
1225 mul64To128(b0, q, &t0, &t1);
1226 sub128(a0, a1, t0, t1, &a0, &a1);
1227 shortShift128Left(a0, a1, 62, &a0, &a1);
1228 exp_diff -= 62;
1229 quot = (quot << 62) + q;
1232 exp_diff += 64;
1233 if (exp_diff > 0) {
1234 q = estimateDiv128To64(a0, a1, b0);
1235 q = q > 2 ? (q - 2) >> (64 - exp_diff) : 0;
1236 mul64To128(b0, q << (64 - exp_diff), &t0, &t1);
1237 sub128(a0, a1, t0, t1, &a0, &a1);
1238 shortShift128Left(0, b0, 64 - exp_diff, &t0, &t1);
1239 while (le128(t0, t1, a0, a1)) {
1240 ++q;
1241 sub128(a0, a1, t0, t1, &a0, &a1);
1243 quot = (exp_diff < 64 ? quot << exp_diff : 0) + q;
1244 } else {
1245 t0 = b0;
1246 t1 = 0;
1249 if (mod_quot) {
1250 *mod_quot = quot;
1251 } else {
1252 sub128(t0, t1, a0, a1, &t0, &t1);
1253 if (lt128(t0, t1, a0, a1) ||
1254 (eq128(t0, t1, a0, a1) && (q & 1))) {
1255 a0 = t0;
1256 a1 = t1;
1257 a->sign = !a->sign;
1261 if (likely(a0)) {
1262 shift = clz64(a0);
1263 shortShift128Left(a0, a1, shift, &a0, &a1);
1264 } else if (likely(a1)) {
1265 shift = clz64(a1);
1266 a0 = a1 << shift;
1267 a1 = 0;
1268 shift += 64;
1269 } else {
1270 a->cls = float_class_zero;
1271 return;
1274 a->exp = b->exp + exp_diff - shift;
1275 a->frac = a0 | (a1 != 0);
1278 static void frac128_modrem(FloatParts128 *a, FloatParts128 *b,
1279 uint64_t *mod_quot)
1281 uint64_t a0, a1, a2, b0, b1, t0, t1, t2, q, quot;
1282 int exp_diff = a->exp - b->exp;
1283 int shift;
1285 a0 = a->frac_hi;
1286 a1 = a->frac_lo;
1287 a2 = 0;
1289 if (exp_diff < -1) {
1290 if (mod_quot) {
1291 *mod_quot = 0;
1293 return;
1295 if (exp_diff == -1) {
1296 shift128Right(a0, a1, 1, &a0, &a1);
1297 exp_diff = 0;
1300 b0 = b->frac_hi;
1301 b1 = b->frac_lo;
1303 quot = q = le128(b0, b1, a0, a1);
1304 if (q) {
1305 sub128(a0, a1, b0, b1, &a0, &a1);
1308 exp_diff -= 64;
1309 while (exp_diff > 0) {
1310 q = estimateDiv128To64(a0, a1, b0);
1311 q = q > 4 ? q - 4 : 0;
1312 mul128By64To192(b0, b1, q, &t0, &t1, &t2);
1313 sub192(a0, a1, a2, t0, t1, t2, &a0, &a1, &a2);
1314 shortShift192Left(a0, a1, a2, 61, &a0, &a1, &a2);
1315 exp_diff -= 61;
1316 quot = (quot << 61) + q;
1319 exp_diff += 64;
1320 if (exp_diff > 0) {
1321 q = estimateDiv128To64(a0, a1, b0);
1322 q = q > 4 ? (q - 4) >> (64 - exp_diff) : 0;
1323 mul128By64To192(b0, b1, q << (64 - exp_diff), &t0, &t1, &t2);
1324 sub192(a0, a1, a2, t0, t1, t2, &a0, &a1, &a2);
1325 shortShift192Left(0, b0, b1, 64 - exp_diff, &t0, &t1, &t2);
1326 while (le192(t0, t1, t2, a0, a1, a2)) {
1327 ++q;
1328 sub192(a0, a1, a2, t0, t1, t2, &a0, &a1, &a2);
1330 quot = (exp_diff < 64 ? quot << exp_diff : 0) + q;
1331 } else {
1332 t0 = b0;
1333 t1 = b1;
1334 t2 = 0;
1337 if (mod_quot) {
1338 *mod_quot = quot;
1339 } else {
1340 sub192(t0, t1, t2, a0, a1, a2, &t0, &t1, &t2);
1341 if (lt192(t0, t1, t2, a0, a1, a2) ||
1342 (eq192(t0, t1, t2, a0, a1, a2) && (q & 1))) {
1343 a0 = t0;
1344 a1 = t1;
1345 a2 = t2;
1346 a->sign = !a->sign;
1350 if (likely(a0)) {
1351 shift = clz64(a0);
1352 shortShift192Left(a0, a1, a2, shift, &a0, &a1, &a2);
1353 } else if (likely(a1)) {
1354 shift = clz64(a1);
1355 shortShift128Left(a1, a2, shift, &a0, &a1);
1356 a2 = 0;
1357 shift += 64;
1358 } else if (likely(a2)) {
1359 shift = clz64(a2);
1360 a0 = a2 << shift;
1361 a1 = a2 = 0;
1362 shift += 128;
1363 } else {
1364 a->cls = float_class_zero;
1365 return;
1368 a->exp = b->exp + exp_diff - shift;
1369 a->frac_hi = a0;
1370 a->frac_lo = a1 | (a2 != 0);
1373 #define frac_modrem(A, B, Q) FRAC_GENERIC_64_128(modrem, A)(A, B, Q)
1375 static void frac64_shl(FloatParts64 *a, int c)
1377 a->frac <<= c;
1380 static void frac128_shl(FloatParts128 *a, int c)
1382 uint64_t a0 = a->frac_hi, a1 = a->frac_lo;
1384 if (c & 64) {
1385 a0 = a1, a1 = 0;
1388 c &= 63;
1389 if (c) {
1390 a0 = shl_double(a0, a1, c);
1391 a1 = a1 << c;
1394 a->frac_hi = a0;
1395 a->frac_lo = a1;
1398 #define frac_shl(A, C) FRAC_GENERIC_64_128(shl, A)(A, C)
1400 static void frac64_shr(FloatParts64 *a, int c)
1402 a->frac >>= c;
1405 static void frac128_shr(FloatParts128 *a, int c)
1407 uint64_t a0 = a->frac_hi, a1 = a->frac_lo;
1409 if (c & 64) {
1410 a1 = a0, a0 = 0;
1413 c &= 63;
1414 if (c) {
1415 a1 = shr_double(a0, a1, c);
1416 a0 = a0 >> c;
1419 a->frac_hi = a0;
1420 a->frac_lo = a1;
1423 #define frac_shr(A, C) FRAC_GENERIC_64_128(shr, A)(A, C)
1425 static void frac64_shrjam(FloatParts64 *a, int c)
1427 uint64_t a0 = a->frac;
1429 if (likely(c != 0)) {
1430 if (likely(c < 64)) {
1431 a0 = (a0 >> c) | (shr_double(a0, 0, c) != 0);
1432 } else {
1433 a0 = a0 != 0;
1435 a->frac = a0;
1439 static void frac128_shrjam(FloatParts128 *a, int c)
1441 uint64_t a0 = a->frac_hi, a1 = a->frac_lo;
1442 uint64_t sticky = 0;
1444 if (unlikely(c == 0)) {
1445 return;
1446 } else if (likely(c < 64)) {
1447 /* nothing */
1448 } else if (likely(c < 128)) {
1449 sticky = a1;
1450 a1 = a0;
1451 a0 = 0;
1452 c &= 63;
1453 if (c == 0) {
1454 goto done;
1456 } else {
1457 sticky = a0 | a1;
1458 a0 = a1 = 0;
1459 goto done;
1462 sticky |= shr_double(a1, 0, c);
1463 a1 = shr_double(a0, a1, c);
1464 a0 = a0 >> c;
1466 done:
1467 a->frac_lo = a1 | (sticky != 0);
1468 a->frac_hi = a0;
1471 static void frac256_shrjam(FloatParts256 *a, int c)
1473 uint64_t a0 = a->frac_hi, a1 = a->frac_hm;
1474 uint64_t a2 = a->frac_lm, a3 = a->frac_lo;
1475 uint64_t sticky = 0;
1477 if (unlikely(c == 0)) {
1478 return;
1479 } else if (likely(c < 64)) {
1480 /* nothing */
1481 } else if (likely(c < 256)) {
1482 if (unlikely(c & 128)) {
1483 sticky |= a2 | a3;
1484 a3 = a1, a2 = a0, a1 = 0, a0 = 0;
1486 if (unlikely(c & 64)) {
1487 sticky |= a3;
1488 a3 = a2, a2 = a1, a1 = a0, a0 = 0;
1490 c &= 63;
1491 if (c == 0) {
1492 goto done;
1494 } else {
1495 sticky = a0 | a1 | a2 | a3;
1496 a0 = a1 = a2 = a3 = 0;
1497 goto done;
1500 sticky |= shr_double(a3, 0, c);
1501 a3 = shr_double(a2, a3, c);
1502 a2 = shr_double(a1, a2, c);
1503 a1 = shr_double(a0, a1, c);
1504 a0 = a0 >> c;
1506 done:
1507 a->frac_lo = a3 | (sticky != 0);
1508 a->frac_lm = a2;
1509 a->frac_hm = a1;
1510 a->frac_hi = a0;
1513 #define frac_shrjam(A, C) FRAC_GENERIC_64_128_256(shrjam, A)(A, C)
1515 static bool frac64_sub(FloatParts64 *r, FloatParts64 *a, FloatParts64 *b)
1517 return usub64_overflow(a->frac, b->frac, &r->frac);
1520 static bool frac128_sub(FloatParts128 *r, FloatParts128 *a, FloatParts128 *b)
1522 bool c = 0;
1523 r->frac_lo = usub64_borrow(a->frac_lo, b->frac_lo, &c);
1524 r->frac_hi = usub64_borrow(a->frac_hi, b->frac_hi, &c);
1525 return c;
1528 static bool frac256_sub(FloatParts256 *r, FloatParts256 *a, FloatParts256 *b)
1530 bool c = 0;
1531 r->frac_lo = usub64_borrow(a->frac_lo, b->frac_lo, &c);
1532 r->frac_lm = usub64_borrow(a->frac_lm, b->frac_lm, &c);
1533 r->frac_hm = usub64_borrow(a->frac_hm, b->frac_hm, &c);
1534 r->frac_hi = usub64_borrow(a->frac_hi, b->frac_hi, &c);
1535 return c;
1538 #define frac_sub(R, A, B) FRAC_GENERIC_64_128_256(sub, R)(R, A, B)
1540 static void frac64_truncjam(FloatParts64 *r, FloatParts128 *a)
1542 r->frac = a->frac_hi | (a->frac_lo != 0);
1545 static void frac128_truncjam(FloatParts128 *r, FloatParts256 *a)
1547 r->frac_hi = a->frac_hi;
1548 r->frac_lo = a->frac_hm | ((a->frac_lm | a->frac_lo) != 0);
1551 #define frac_truncjam(R, A) FRAC_GENERIC_64_128(truncjam, R)(R, A)
1553 static void frac64_widen(FloatParts128 *r, FloatParts64 *a)
1555 r->frac_hi = a->frac;
1556 r->frac_lo = 0;
1559 static void frac128_widen(FloatParts256 *r, FloatParts128 *a)
1561 r->frac_hi = a->frac_hi;
1562 r->frac_hm = a->frac_lo;
1563 r->frac_lm = 0;
1564 r->frac_lo = 0;
1567 #define frac_widen(A, B) FRAC_GENERIC_64_128(widen, B)(A, B)
1570 * Reciprocal sqrt table. 1 bit of exponent, 6-bits of mantessa.
1571 * From https://git.musl-libc.org/cgit/musl/tree/src/math/sqrt_data.c
1572 * and thus MIT licenced.
1574 static const uint16_t rsqrt_tab[128] = {
1575 0xb451, 0xb2f0, 0xb196, 0xb044, 0xaef9, 0xadb6, 0xac79, 0xab43,
1576 0xaa14, 0xa8eb, 0xa7c8, 0xa6aa, 0xa592, 0xa480, 0xa373, 0xa26b,
1577 0xa168, 0xa06a, 0x9f70, 0x9e7b, 0x9d8a, 0x9c9d, 0x9bb5, 0x9ad1,
1578 0x99f0, 0x9913, 0x983a, 0x9765, 0x9693, 0x95c4, 0x94f8, 0x9430,
1579 0x936b, 0x92a9, 0x91ea, 0x912e, 0x9075, 0x8fbe, 0x8f0a, 0x8e59,
1580 0x8daa, 0x8cfe, 0x8c54, 0x8bac, 0x8b07, 0x8a64, 0x89c4, 0x8925,
1581 0x8889, 0x87ee, 0x8756, 0x86c0, 0x862b, 0x8599, 0x8508, 0x8479,
1582 0x83ec, 0x8361, 0x82d8, 0x8250, 0x81c9, 0x8145, 0x80c2, 0x8040,
1583 0xff02, 0xfd0e, 0xfb25, 0xf947, 0xf773, 0xf5aa, 0xf3ea, 0xf234,
1584 0xf087, 0xeee3, 0xed47, 0xebb3, 0xea27, 0xe8a3, 0xe727, 0xe5b2,
1585 0xe443, 0xe2dc, 0xe17a, 0xe020, 0xdecb, 0xdd7d, 0xdc34, 0xdaf1,
1586 0xd9b3, 0xd87b, 0xd748, 0xd61a, 0xd4f1, 0xd3cd, 0xd2ad, 0xd192,
1587 0xd07b, 0xcf69, 0xce5b, 0xcd51, 0xcc4a, 0xcb48, 0xca4a, 0xc94f,
1588 0xc858, 0xc764, 0xc674, 0xc587, 0xc49d, 0xc3b7, 0xc2d4, 0xc1f4,
1589 0xc116, 0xc03c, 0xbf65, 0xbe90, 0xbdbe, 0xbcef, 0xbc23, 0xbb59,
1590 0xba91, 0xb9cc, 0xb90a, 0xb84a, 0xb78c, 0xb6d0, 0xb617, 0xb560,
1593 #define partsN(NAME) glue(glue(glue(parts,N),_),NAME)
1594 #define FloatPartsN glue(FloatParts,N)
1595 #define FloatPartsW glue(FloatParts,W)
1597 #define N 64
1598 #define W 128
1600 #include "softfloat-parts-addsub.c.inc"
1601 #include "softfloat-parts.c.inc"
1603 #undef N
1604 #undef W
1605 #define N 128
1606 #define W 256
1608 #include "softfloat-parts-addsub.c.inc"
1609 #include "softfloat-parts.c.inc"
1611 #undef N
1612 #undef W
1613 #define N 256
1615 #include "softfloat-parts-addsub.c.inc"
1617 #undef N
1618 #undef W
1619 #undef partsN
1620 #undef FloatPartsN
1621 #undef FloatPartsW
1624 * Pack/unpack routines with a specific FloatFmt.
1627 static void float16a_unpack_canonical(FloatParts64 *p, float16 f,
1628 float_status *s, const FloatFmt *params)
1630 float16_unpack_raw(p, f);
1631 parts_canonicalize(p, s, params);
1634 static void float16_unpack_canonical(FloatParts64 *p, float16 f,
1635 float_status *s)
1637 float16a_unpack_canonical(p, f, s, &float16_params);
1640 static void bfloat16_unpack_canonical(FloatParts64 *p, bfloat16 f,
1641 float_status *s)
1643 bfloat16_unpack_raw(p, f);
1644 parts_canonicalize(p, s, &bfloat16_params);
1647 static float16 float16a_round_pack_canonical(FloatParts64 *p,
1648 float_status *s,
1649 const FloatFmt *params)
1651 parts_uncanon(p, s, params);
1652 return float16_pack_raw(p);
1655 static float16 float16_round_pack_canonical(FloatParts64 *p,
1656 float_status *s)
1658 return float16a_round_pack_canonical(p, s, &float16_params);
1661 static bfloat16 bfloat16_round_pack_canonical(FloatParts64 *p,
1662 float_status *s)
1664 parts_uncanon(p, s, &bfloat16_params);
1665 return bfloat16_pack_raw(p);
1668 static void float32_unpack_canonical(FloatParts64 *p, float32 f,
1669 float_status *s)
1671 float32_unpack_raw(p, f);
1672 parts_canonicalize(p, s, &float32_params);
1675 static float32 float32_round_pack_canonical(FloatParts64 *p,
1676 float_status *s)
1678 parts_uncanon(p, s, &float32_params);
1679 return float32_pack_raw(p);
1682 static void float64_unpack_canonical(FloatParts64 *p, float64 f,
1683 float_status *s)
1685 float64_unpack_raw(p, f);
1686 parts_canonicalize(p, s, &float64_params);
1689 static float64 float64_round_pack_canonical(FloatParts64 *p,
1690 float_status *s)
1692 parts_uncanon(p, s, &float64_params);
1693 return float64_pack_raw(p);
1696 static void float128_unpack_canonical(FloatParts128 *p, float128 f,
1697 float_status *s)
1699 float128_unpack_raw(p, f);
1700 parts_canonicalize(p, s, &float128_params);
1703 static float128 float128_round_pack_canonical(FloatParts128 *p,
1704 float_status *s)
1706 parts_uncanon(p, s, &float128_params);
1707 return float128_pack_raw(p);
1710 /* Returns false if the encoding is invalid. */
1711 static bool floatx80_unpack_canonical(FloatParts128 *p, floatx80 f,
1712 float_status *s)
1714 /* Ensure rounding precision is set before beginning. */
1715 switch (s->floatx80_rounding_precision) {
1716 case floatx80_precision_x:
1717 case floatx80_precision_d:
1718 case floatx80_precision_s:
1719 break;
1720 default:
1721 g_assert_not_reached();
1724 if (unlikely(floatx80_invalid_encoding(f))) {
1725 float_raise(float_flag_invalid, s);
1726 return false;
1729 floatx80_unpack_raw(p, f);
1731 if (likely(p->exp != floatx80_params[floatx80_precision_x].exp_max)) {
1732 parts_canonicalize(p, s, &floatx80_params[floatx80_precision_x]);
1733 } else {
1734 /* The explicit integer bit is ignored, after invalid checks. */
1735 p->frac_hi &= MAKE_64BIT_MASK(0, 63);
1736 p->cls = (p->frac_hi == 0 ? float_class_inf
1737 : parts_is_snan_frac(p->frac_hi, s)
1738 ? float_class_snan : float_class_qnan);
1740 return true;
1743 static floatx80 floatx80_round_pack_canonical(FloatParts128 *p,
1744 float_status *s)
1746 const FloatFmt *fmt = &floatx80_params[s->floatx80_rounding_precision];
1747 uint64_t frac;
1748 int exp;
1750 switch (p->cls) {
1751 case float_class_normal:
1752 if (s->floatx80_rounding_precision == floatx80_precision_x) {
1753 parts_uncanon_normal(p, s, fmt);
1754 frac = p->frac_hi;
1755 exp = p->exp;
1756 } else {
1757 FloatParts64 p64;
1759 p64.sign = p->sign;
1760 p64.exp = p->exp;
1761 frac_truncjam(&p64, p);
1762 parts_uncanon_normal(&p64, s, fmt);
1763 frac = p64.frac;
1764 exp = p64.exp;
1766 if (exp != fmt->exp_max) {
1767 break;
1769 /* rounded to inf -- fall through to set frac correctly */
1771 case float_class_inf:
1772 /* x86 and m68k differ in the setting of the integer bit. */
1773 frac = floatx80_infinity_low;
1774 exp = fmt->exp_max;
1775 break;
1777 case float_class_zero:
1778 frac = 0;
1779 exp = 0;
1780 break;
1782 case float_class_snan:
1783 case float_class_qnan:
1784 /* NaNs have the integer bit set. */
1785 frac = p->frac_hi | (1ull << 63);
1786 exp = fmt->exp_max;
1787 break;
1789 default:
1790 g_assert_not_reached();
1793 return packFloatx80(p->sign, exp, frac);
1797 * Addition and subtraction
1800 static float16 QEMU_FLATTEN
1801 float16_addsub(float16 a, float16 b, float_status *status, bool subtract)
1803 FloatParts64 pa, pb, *pr;
1805 float16_unpack_canonical(&pa, a, status);
1806 float16_unpack_canonical(&pb, b, status);
1807 pr = parts_addsub(&pa, &pb, status, subtract);
1809 return float16_round_pack_canonical(pr, status);
1812 float16 float16_add(float16 a, float16 b, float_status *status)
1814 return float16_addsub(a, b, status, false);
1817 float16 float16_sub(float16 a, float16 b, float_status *status)
1819 return float16_addsub(a, b, status, true);
1822 static float32 QEMU_SOFTFLOAT_ATTR
1823 soft_f32_addsub(float32 a, float32 b, float_status *status, bool subtract)
1825 FloatParts64 pa, pb, *pr;
1827 float32_unpack_canonical(&pa, a, status);
1828 float32_unpack_canonical(&pb, b, status);
1829 pr = parts_addsub(&pa, &pb, status, subtract);
1831 return float32_round_pack_canonical(pr, status);
1834 static float32 soft_f32_add(float32 a, float32 b, float_status *status)
1836 return soft_f32_addsub(a, b, status, false);
1839 static float32 soft_f32_sub(float32 a, float32 b, float_status *status)
1841 return soft_f32_addsub(a, b, status, true);
1844 static float64 QEMU_SOFTFLOAT_ATTR
1845 soft_f64_addsub(float64 a, float64 b, float_status *status, bool subtract)
1847 FloatParts64 pa, pb, *pr;
1849 float64_unpack_canonical(&pa, a, status);
1850 float64_unpack_canonical(&pb, b, status);
1851 pr = parts_addsub(&pa, &pb, status, subtract);
1853 return float64_round_pack_canonical(pr, status);
1856 static float64 soft_f64_add(float64 a, float64 b, float_status *status)
1858 return soft_f64_addsub(a, b, status, false);
1861 static float64 soft_f64_sub(float64 a, float64 b, float_status *status)
1863 return soft_f64_addsub(a, b, status, true);
1866 static float hard_f32_add(float a, float b)
1868 return a + b;
1871 static float hard_f32_sub(float a, float b)
1873 return a - b;
1876 static double hard_f64_add(double a, double b)
1878 return a + b;
1881 static double hard_f64_sub(double a, double b)
1883 return a - b;
1886 static bool f32_addsubmul_post(union_float32 a, union_float32 b)
1888 if (QEMU_HARDFLOAT_2F32_USE_FP) {
1889 return !(fpclassify(a.h) == FP_ZERO && fpclassify(b.h) == FP_ZERO);
1891 return !(float32_is_zero(a.s) && float32_is_zero(b.s));
1894 static bool f64_addsubmul_post(union_float64 a, union_float64 b)
1896 if (QEMU_HARDFLOAT_2F64_USE_FP) {
1897 return !(fpclassify(a.h) == FP_ZERO && fpclassify(b.h) == FP_ZERO);
1898 } else {
1899 return !(float64_is_zero(a.s) && float64_is_zero(b.s));
1903 static float32 float32_addsub(float32 a, float32 b, float_status *s,
1904 hard_f32_op2_fn hard, soft_f32_op2_fn soft)
1906 return float32_gen2(a, b, s, hard, soft,
1907 f32_is_zon2, f32_addsubmul_post);
1910 static float64 float64_addsub(float64 a, float64 b, float_status *s,
1911 hard_f64_op2_fn hard, soft_f64_op2_fn soft)
1913 return float64_gen2(a, b, s, hard, soft,
1914 f64_is_zon2, f64_addsubmul_post);
1917 float32 QEMU_FLATTEN
1918 float32_add(float32 a, float32 b, float_status *s)
1920 return float32_addsub(a, b, s, hard_f32_add, soft_f32_add);
1923 float32 QEMU_FLATTEN
1924 float32_sub(float32 a, float32 b, float_status *s)
1926 return float32_addsub(a, b, s, hard_f32_sub, soft_f32_sub);
1929 float64 QEMU_FLATTEN
1930 float64_add(float64 a, float64 b, float_status *s)
1932 return float64_addsub(a, b, s, hard_f64_add, soft_f64_add);
1935 float64 QEMU_FLATTEN
1936 float64_sub(float64 a, float64 b, float_status *s)
1938 return float64_addsub(a, b, s, hard_f64_sub, soft_f64_sub);
1941 static bfloat16 QEMU_FLATTEN
1942 bfloat16_addsub(bfloat16 a, bfloat16 b, float_status *status, bool subtract)
1944 FloatParts64 pa, pb, *pr;
1946 bfloat16_unpack_canonical(&pa, a, status);
1947 bfloat16_unpack_canonical(&pb, b, status);
1948 pr = parts_addsub(&pa, &pb, status, subtract);
1950 return bfloat16_round_pack_canonical(pr, status);
1953 bfloat16 bfloat16_add(bfloat16 a, bfloat16 b, float_status *status)
1955 return bfloat16_addsub(a, b, status, false);
1958 bfloat16 bfloat16_sub(bfloat16 a, bfloat16 b, float_status *status)
1960 return bfloat16_addsub(a, b, status, true);
1963 static float128 QEMU_FLATTEN
1964 float128_addsub(float128 a, float128 b, float_status *status, bool subtract)
1966 FloatParts128 pa, pb, *pr;
1968 float128_unpack_canonical(&pa, a, status);
1969 float128_unpack_canonical(&pb, b, status);
1970 pr = parts_addsub(&pa, &pb, status, subtract);
1972 return float128_round_pack_canonical(pr, status);
1975 float128 float128_add(float128 a, float128 b, float_status *status)
1977 return float128_addsub(a, b, status, false);
1980 float128 float128_sub(float128 a, float128 b, float_status *status)
1982 return float128_addsub(a, b, status, true);
1985 static floatx80 QEMU_FLATTEN
1986 floatx80_addsub(floatx80 a, floatx80 b, float_status *status, bool subtract)
1988 FloatParts128 pa, pb, *pr;
1990 if (!floatx80_unpack_canonical(&pa, a, status) ||
1991 !floatx80_unpack_canonical(&pb, b, status)) {
1992 return floatx80_default_nan(status);
1995 pr = parts_addsub(&pa, &pb, status, subtract);
1996 return floatx80_round_pack_canonical(pr, status);
1999 floatx80 floatx80_add(floatx80 a, floatx80 b, float_status *status)
2001 return floatx80_addsub(a, b, status, false);
2004 floatx80 floatx80_sub(floatx80 a, floatx80 b, float_status *status)
2006 return floatx80_addsub(a, b, status, true);
2010 * Multiplication
2013 float16 QEMU_FLATTEN float16_mul(float16 a, float16 b, float_status *status)
2015 FloatParts64 pa, pb, *pr;
2017 float16_unpack_canonical(&pa, a, status);
2018 float16_unpack_canonical(&pb, b, status);
2019 pr = parts_mul(&pa, &pb, status);
2021 return float16_round_pack_canonical(pr, status);
2024 static float32 QEMU_SOFTFLOAT_ATTR
2025 soft_f32_mul(float32 a, float32 b, float_status *status)
2027 FloatParts64 pa, pb, *pr;
2029 float32_unpack_canonical(&pa, a, status);
2030 float32_unpack_canonical(&pb, b, status);
2031 pr = parts_mul(&pa, &pb, status);
2033 return float32_round_pack_canonical(pr, status);
2036 static float64 QEMU_SOFTFLOAT_ATTR
2037 soft_f64_mul(float64 a, float64 b, float_status *status)
2039 FloatParts64 pa, pb, *pr;
2041 float64_unpack_canonical(&pa, a, status);
2042 float64_unpack_canonical(&pb, b, status);
2043 pr = parts_mul(&pa, &pb, status);
2045 return float64_round_pack_canonical(pr, status);
2048 static float hard_f32_mul(float a, float b)
2050 return a * b;
2053 static double hard_f64_mul(double a, double b)
2055 return a * b;
2058 float32 QEMU_FLATTEN
2059 float32_mul(float32 a, float32 b, float_status *s)
2061 return float32_gen2(a, b, s, hard_f32_mul, soft_f32_mul,
2062 f32_is_zon2, f32_addsubmul_post);
2065 float64 QEMU_FLATTEN
2066 float64_mul(float64 a, float64 b, float_status *s)
2068 return float64_gen2(a, b, s, hard_f64_mul, soft_f64_mul,
2069 f64_is_zon2, f64_addsubmul_post);
2072 bfloat16 QEMU_FLATTEN
2073 bfloat16_mul(bfloat16 a, bfloat16 b, float_status *status)
2075 FloatParts64 pa, pb, *pr;
2077 bfloat16_unpack_canonical(&pa, a, status);
2078 bfloat16_unpack_canonical(&pb, b, status);
2079 pr = parts_mul(&pa, &pb, status);
2081 return bfloat16_round_pack_canonical(pr, status);
2084 float128 QEMU_FLATTEN
2085 float128_mul(float128 a, float128 b, float_status *status)
2087 FloatParts128 pa, pb, *pr;
2089 float128_unpack_canonical(&pa, a, status);
2090 float128_unpack_canonical(&pb, b, status);
2091 pr = parts_mul(&pa, &pb, status);
2093 return float128_round_pack_canonical(pr, status);
2096 floatx80 QEMU_FLATTEN
2097 floatx80_mul(floatx80 a, floatx80 b, float_status *status)
2099 FloatParts128 pa, pb, *pr;
2101 if (!floatx80_unpack_canonical(&pa, a, status) ||
2102 !floatx80_unpack_canonical(&pb, b, status)) {
2103 return floatx80_default_nan(status);
2106 pr = parts_mul(&pa, &pb, status);
2107 return floatx80_round_pack_canonical(pr, status);
2111 * Fused multiply-add
2114 float16 QEMU_FLATTEN float16_muladd(float16 a, float16 b, float16 c,
2115 int flags, float_status *status)
2117 FloatParts64 pa, pb, pc, *pr;
2119 float16_unpack_canonical(&pa, a, status);
2120 float16_unpack_canonical(&pb, b, status);
2121 float16_unpack_canonical(&pc, c, status);
2122 pr = parts_muladd(&pa, &pb, &pc, flags, status);
2124 return float16_round_pack_canonical(pr, status);
2127 static float32 QEMU_SOFTFLOAT_ATTR
2128 soft_f32_muladd(float32 a, float32 b, float32 c, int flags,
2129 float_status *status)
2131 FloatParts64 pa, pb, pc, *pr;
2133 float32_unpack_canonical(&pa, a, status);
2134 float32_unpack_canonical(&pb, b, status);
2135 float32_unpack_canonical(&pc, c, status);
2136 pr = parts_muladd(&pa, &pb, &pc, flags, status);
2138 return float32_round_pack_canonical(pr, status);
2141 static float64 QEMU_SOFTFLOAT_ATTR
2142 soft_f64_muladd(float64 a, float64 b, float64 c, int flags,
2143 float_status *status)
2145 FloatParts64 pa, pb, pc, *pr;
2147 float64_unpack_canonical(&pa, a, status);
2148 float64_unpack_canonical(&pb, b, status);
2149 float64_unpack_canonical(&pc, c, status);
2150 pr = parts_muladd(&pa, &pb, &pc, flags, status);
2152 return float64_round_pack_canonical(pr, status);
2155 static bool force_soft_fma;
2157 float32 QEMU_FLATTEN
2158 float32_muladd(float32 xa, float32 xb, float32 xc, int flags, float_status *s)
2160 union_float32 ua, ub, uc, ur;
2162 ua.s = xa;
2163 ub.s = xb;
2164 uc.s = xc;
2166 if (unlikely(!can_use_fpu(s))) {
2167 goto soft;
2169 if (unlikely(flags & float_muladd_halve_result)) {
2170 goto soft;
2173 float32_input_flush3(&ua.s, &ub.s, &uc.s, s);
2174 if (unlikely(!f32_is_zon3(ua, ub, uc))) {
2175 goto soft;
2178 if (unlikely(force_soft_fma)) {
2179 goto soft;
2183 * When (a || b) == 0, there's no need to check for under/over flow,
2184 * since we know the addend is (normal || 0) and the product is 0.
2186 if (float32_is_zero(ua.s) || float32_is_zero(ub.s)) {
2187 union_float32 up;
2188 bool prod_sign;
2190 prod_sign = float32_is_neg(ua.s) ^ float32_is_neg(ub.s);
2191 prod_sign ^= !!(flags & float_muladd_negate_product);
2192 up.s = float32_set_sign(float32_zero, prod_sign);
2194 if (flags & float_muladd_negate_c) {
2195 uc.h = -uc.h;
2197 ur.h = up.h + uc.h;
2198 } else {
2199 union_float32 ua_orig = ua;
2200 union_float32 uc_orig = uc;
2202 if (flags & float_muladd_negate_product) {
2203 ua.h = -ua.h;
2205 if (flags & float_muladd_negate_c) {
2206 uc.h = -uc.h;
2209 ur.h = fmaf(ua.h, ub.h, uc.h);
2211 if (unlikely(f32_is_inf(ur))) {
2212 float_raise(float_flag_overflow, s);
2213 } else if (unlikely(fabsf(ur.h) <= FLT_MIN)) {
2214 ua = ua_orig;
2215 uc = uc_orig;
2216 goto soft;
2219 if (flags & float_muladd_negate_result) {
2220 return float32_chs(ur.s);
2222 return ur.s;
2224 soft:
2225 return soft_f32_muladd(ua.s, ub.s, uc.s, flags, s);
2228 float64 QEMU_FLATTEN
2229 float64_muladd(float64 xa, float64 xb, float64 xc, int flags, float_status *s)
2231 union_float64 ua, ub, uc, ur;
2233 ua.s = xa;
2234 ub.s = xb;
2235 uc.s = xc;
2237 if (unlikely(!can_use_fpu(s))) {
2238 goto soft;
2240 if (unlikely(flags & float_muladd_halve_result)) {
2241 goto soft;
2244 float64_input_flush3(&ua.s, &ub.s, &uc.s, s);
2245 if (unlikely(!f64_is_zon3(ua, ub, uc))) {
2246 goto soft;
2249 if (unlikely(force_soft_fma)) {
2250 goto soft;
2254 * When (a || b) == 0, there's no need to check for under/over flow,
2255 * since we know the addend is (normal || 0) and the product is 0.
2257 if (float64_is_zero(ua.s) || float64_is_zero(ub.s)) {
2258 union_float64 up;
2259 bool prod_sign;
2261 prod_sign = float64_is_neg(ua.s) ^ float64_is_neg(ub.s);
2262 prod_sign ^= !!(flags & float_muladd_negate_product);
2263 up.s = float64_set_sign(float64_zero, prod_sign);
2265 if (flags & float_muladd_negate_c) {
2266 uc.h = -uc.h;
2268 ur.h = up.h + uc.h;
2269 } else {
2270 union_float64 ua_orig = ua;
2271 union_float64 uc_orig = uc;
2273 if (flags & float_muladd_negate_product) {
2274 ua.h = -ua.h;
2276 if (flags & float_muladd_negate_c) {
2277 uc.h = -uc.h;
2280 ur.h = fma(ua.h, ub.h, uc.h);
2282 if (unlikely(f64_is_inf(ur))) {
2283 float_raise(float_flag_overflow, s);
2284 } else if (unlikely(fabs(ur.h) <= FLT_MIN)) {
2285 ua = ua_orig;
2286 uc = uc_orig;
2287 goto soft;
2290 if (flags & float_muladd_negate_result) {
2291 return float64_chs(ur.s);
2293 return ur.s;
2295 soft:
2296 return soft_f64_muladd(ua.s, ub.s, uc.s, flags, s);
2299 bfloat16 QEMU_FLATTEN bfloat16_muladd(bfloat16 a, bfloat16 b, bfloat16 c,
2300 int flags, float_status *status)
2302 FloatParts64 pa, pb, pc, *pr;
2304 bfloat16_unpack_canonical(&pa, a, status);
2305 bfloat16_unpack_canonical(&pb, b, status);
2306 bfloat16_unpack_canonical(&pc, c, status);
2307 pr = parts_muladd(&pa, &pb, &pc, flags, status);
2309 return bfloat16_round_pack_canonical(pr, status);
2312 float128 QEMU_FLATTEN float128_muladd(float128 a, float128 b, float128 c,
2313 int flags, float_status *status)
2315 FloatParts128 pa, pb, pc, *pr;
2317 float128_unpack_canonical(&pa, a, status);
2318 float128_unpack_canonical(&pb, b, status);
2319 float128_unpack_canonical(&pc, c, status);
2320 pr = parts_muladd(&pa, &pb, &pc, flags, status);
2322 return float128_round_pack_canonical(pr, status);
2326 * Division
2329 float16 float16_div(float16 a, float16 b, float_status *status)
2331 FloatParts64 pa, pb, *pr;
2333 float16_unpack_canonical(&pa, a, status);
2334 float16_unpack_canonical(&pb, b, status);
2335 pr = parts_div(&pa, &pb, status);
2337 return float16_round_pack_canonical(pr, status);
2340 static float32 QEMU_SOFTFLOAT_ATTR
2341 soft_f32_div(float32 a, float32 b, float_status *status)
2343 FloatParts64 pa, pb, *pr;
2345 float32_unpack_canonical(&pa, a, status);
2346 float32_unpack_canonical(&pb, b, status);
2347 pr = parts_div(&pa, &pb, status);
2349 return float32_round_pack_canonical(pr, status);
2352 static float64 QEMU_SOFTFLOAT_ATTR
2353 soft_f64_div(float64 a, float64 b, float_status *status)
2355 FloatParts64 pa, pb, *pr;
2357 float64_unpack_canonical(&pa, a, status);
2358 float64_unpack_canonical(&pb, b, status);
2359 pr = parts_div(&pa, &pb, status);
2361 return float64_round_pack_canonical(pr, status);
2364 static float hard_f32_div(float a, float b)
2366 return a / b;
2369 static double hard_f64_div(double a, double b)
2371 return a / b;
2374 static bool f32_div_pre(union_float32 a, union_float32 b)
2376 if (QEMU_HARDFLOAT_2F32_USE_FP) {
2377 return (fpclassify(a.h) == FP_NORMAL || fpclassify(a.h) == FP_ZERO) &&
2378 fpclassify(b.h) == FP_NORMAL;
2380 return float32_is_zero_or_normal(a.s) && float32_is_normal(b.s);
2383 static bool f64_div_pre(union_float64 a, union_float64 b)
2385 if (QEMU_HARDFLOAT_2F64_USE_FP) {
2386 return (fpclassify(a.h) == FP_NORMAL || fpclassify(a.h) == FP_ZERO) &&
2387 fpclassify(b.h) == FP_NORMAL;
2389 return float64_is_zero_or_normal(a.s) && float64_is_normal(b.s);
2392 static bool f32_div_post(union_float32 a, union_float32 b)
2394 if (QEMU_HARDFLOAT_2F32_USE_FP) {
2395 return fpclassify(a.h) != FP_ZERO;
2397 return !float32_is_zero(a.s);
2400 static bool f64_div_post(union_float64 a, union_float64 b)
2402 if (QEMU_HARDFLOAT_2F64_USE_FP) {
2403 return fpclassify(a.h) != FP_ZERO;
2405 return !float64_is_zero(a.s);
2408 float32 QEMU_FLATTEN
2409 float32_div(float32 a, float32 b, float_status *s)
2411 return float32_gen2(a, b, s, hard_f32_div, soft_f32_div,
2412 f32_div_pre, f32_div_post);
2415 float64 QEMU_FLATTEN
2416 float64_div(float64 a, float64 b, float_status *s)
2418 return float64_gen2(a, b, s, hard_f64_div, soft_f64_div,
2419 f64_div_pre, f64_div_post);
2422 bfloat16 QEMU_FLATTEN
2423 bfloat16_div(bfloat16 a, bfloat16 b, float_status *status)
2425 FloatParts64 pa, pb, *pr;
2427 bfloat16_unpack_canonical(&pa, a, status);
2428 bfloat16_unpack_canonical(&pb, b, status);
2429 pr = parts_div(&pa, &pb, status);
2431 return bfloat16_round_pack_canonical(pr, status);
2434 float128 QEMU_FLATTEN
2435 float128_div(float128 a, float128 b, float_status *status)
2437 FloatParts128 pa, pb, *pr;
2439 float128_unpack_canonical(&pa, a, status);
2440 float128_unpack_canonical(&pb, b, status);
2441 pr = parts_div(&pa, &pb, status);
2443 return float128_round_pack_canonical(pr, status);
2446 floatx80 floatx80_div(floatx80 a, floatx80 b, float_status *status)
2448 FloatParts128 pa, pb, *pr;
2450 if (!floatx80_unpack_canonical(&pa, a, status) ||
2451 !floatx80_unpack_canonical(&pb, b, status)) {
2452 return floatx80_default_nan(status);
2455 pr = parts_div(&pa, &pb, status);
2456 return floatx80_round_pack_canonical(pr, status);
2460 * Remainder
2463 float32 float32_rem(float32 a, float32 b, float_status *status)
2465 FloatParts64 pa, pb, *pr;
2467 float32_unpack_canonical(&pa, a, status);
2468 float32_unpack_canonical(&pb, b, status);
2469 pr = parts_modrem(&pa, &pb, NULL, status);
2471 return float32_round_pack_canonical(pr, status);
2474 float64 float64_rem(float64 a, float64 b, float_status *status)
2476 FloatParts64 pa, pb, *pr;
2478 float64_unpack_canonical(&pa, a, status);
2479 float64_unpack_canonical(&pb, b, status);
2480 pr = parts_modrem(&pa, &pb, NULL, status);
2482 return float64_round_pack_canonical(pr, status);
2485 float128 float128_rem(float128 a, float128 b, float_status *status)
2487 FloatParts128 pa, pb, *pr;
2489 float128_unpack_canonical(&pa, a, status);
2490 float128_unpack_canonical(&pb, b, status);
2491 pr = parts_modrem(&pa, &pb, NULL, status);
2493 return float128_round_pack_canonical(pr, status);
2497 * Returns the remainder of the extended double-precision floating-point value
2498 * `a' with respect to the corresponding value `b'.
2499 * If 'mod' is false, the operation is performed according to the IEC/IEEE
2500 * Standard for Binary Floating-Point Arithmetic. If 'mod' is true, return
2501 * the remainder based on truncating the quotient toward zero instead and
2502 * *quotient is set to the low 64 bits of the absolute value of the integer
2503 * quotient.
2505 floatx80 floatx80_modrem(floatx80 a, floatx80 b, bool mod,
2506 uint64_t *quotient, float_status *status)
2508 FloatParts128 pa, pb, *pr;
2510 *quotient = 0;
2511 if (!floatx80_unpack_canonical(&pa, a, status) ||
2512 !floatx80_unpack_canonical(&pb, b, status)) {
2513 return floatx80_default_nan(status);
2515 pr = parts_modrem(&pa, &pb, mod ? quotient : NULL, status);
2517 return floatx80_round_pack_canonical(pr, status);
2520 floatx80 floatx80_rem(floatx80 a, floatx80 b, float_status *status)
2522 uint64_t quotient;
2523 return floatx80_modrem(a, b, false, &quotient, status);
2526 floatx80 floatx80_mod(floatx80 a, floatx80 b, float_status *status)
2528 uint64_t quotient;
2529 return floatx80_modrem(a, b, true, &quotient, status);
2533 * Float to Float conversions
2535 * Returns the result of converting one float format to another. The
2536 * conversion is performed according to the IEC/IEEE Standard for
2537 * Binary Floating-Point Arithmetic.
2539 * Usually this only needs to take care of raising invalid exceptions
2540 * and handling the conversion on NaNs.
2543 static void parts_float_to_ahp(FloatParts64 *a, float_status *s)
2545 switch (a->cls) {
2546 case float_class_qnan:
2547 case float_class_snan:
2549 * There is no NaN in the destination format. Raise Invalid
2550 * and return a zero with the sign of the input NaN.
2552 float_raise(float_flag_invalid, s);
2553 a->cls = float_class_zero;
2554 break;
2556 case float_class_inf:
2558 * There is no Inf in the destination format. Raise Invalid
2559 * and return the maximum normal with the correct sign.
2561 float_raise(float_flag_invalid, s);
2562 a->cls = float_class_normal;
2563 a->exp = float16_params_ahp.exp_max;
2564 a->frac = MAKE_64BIT_MASK(float16_params_ahp.frac_shift,
2565 float16_params_ahp.frac_size + 1);
2566 break;
2568 case float_class_normal:
2569 case float_class_zero:
2570 break;
2572 default:
2573 g_assert_not_reached();
2577 static void parts64_float_to_float(FloatParts64 *a, float_status *s)
2579 if (is_nan(a->cls)) {
2580 parts_return_nan(a, s);
2584 static void parts128_float_to_float(FloatParts128 *a, float_status *s)
2586 if (is_nan(a->cls)) {
2587 parts_return_nan(a, s);
2591 #define parts_float_to_float(P, S) \
2592 PARTS_GENERIC_64_128(float_to_float, P)(P, S)
2594 static void parts_float_to_float_narrow(FloatParts64 *a, FloatParts128 *b,
2595 float_status *s)
2597 a->cls = b->cls;
2598 a->sign = b->sign;
2599 a->exp = b->exp;
2601 if (a->cls == float_class_normal) {
2602 frac_truncjam(a, b);
2603 } else if (is_nan(a->cls)) {
2604 /* Discard the low bits of the NaN. */
2605 a->frac = b->frac_hi;
2606 parts_return_nan(a, s);
2610 static void parts_float_to_float_widen(FloatParts128 *a, FloatParts64 *b,
2611 float_status *s)
2613 a->cls = b->cls;
2614 a->sign = b->sign;
2615 a->exp = b->exp;
2616 frac_widen(a, b);
2618 if (is_nan(a->cls)) {
2619 parts_return_nan(a, s);
2623 float32 float16_to_float32(float16 a, bool ieee, float_status *s)
2625 const FloatFmt *fmt16 = ieee ? &float16_params : &float16_params_ahp;
2626 FloatParts64 p;
2628 float16a_unpack_canonical(&p, a, s, fmt16);
2629 parts_float_to_float(&p, s);
2630 return float32_round_pack_canonical(&p, s);
2633 float64 float16_to_float64(float16 a, bool ieee, float_status *s)
2635 const FloatFmt *fmt16 = ieee ? &float16_params : &float16_params_ahp;
2636 FloatParts64 p;
2638 float16a_unpack_canonical(&p, a, s, fmt16);
2639 parts_float_to_float(&p, s);
2640 return float64_round_pack_canonical(&p, s);
2643 float16 float32_to_float16(float32 a, bool ieee, float_status *s)
2645 FloatParts64 p;
2646 const FloatFmt *fmt;
2648 float32_unpack_canonical(&p, a, s);
2649 if (ieee) {
2650 parts_float_to_float(&p, s);
2651 fmt = &float16_params;
2652 } else {
2653 parts_float_to_ahp(&p, s);
2654 fmt = &float16_params_ahp;
2656 return float16a_round_pack_canonical(&p, s, fmt);
2659 static float64 QEMU_SOFTFLOAT_ATTR
2660 soft_float32_to_float64(float32 a, float_status *s)
2662 FloatParts64 p;
2664 float32_unpack_canonical(&p, a, s);
2665 parts_float_to_float(&p, s);
2666 return float64_round_pack_canonical(&p, s);
2669 float64 float32_to_float64(float32 a, float_status *s)
2671 if (likely(float32_is_normal(a))) {
2672 /* Widening conversion can never produce inexact results. */
2673 union_float32 uf;
2674 union_float64 ud;
2675 uf.s = a;
2676 ud.h = uf.h;
2677 return ud.s;
2678 } else if (float32_is_zero(a)) {
2679 return float64_set_sign(float64_zero, float32_is_neg(a));
2680 } else {
2681 return soft_float32_to_float64(a, s);
2685 float16 float64_to_float16(float64 a, bool ieee, float_status *s)
2687 FloatParts64 p;
2688 const FloatFmt *fmt;
2690 float64_unpack_canonical(&p, a, s);
2691 if (ieee) {
2692 parts_float_to_float(&p, s);
2693 fmt = &float16_params;
2694 } else {
2695 parts_float_to_ahp(&p, s);
2696 fmt = &float16_params_ahp;
2698 return float16a_round_pack_canonical(&p, s, fmt);
2701 float32 float64_to_float32(float64 a, float_status *s)
2703 FloatParts64 p;
2705 float64_unpack_canonical(&p, a, s);
2706 parts_float_to_float(&p, s);
2707 return float32_round_pack_canonical(&p, s);
2710 float32 bfloat16_to_float32(bfloat16 a, float_status *s)
2712 FloatParts64 p;
2714 bfloat16_unpack_canonical(&p, a, s);
2715 parts_float_to_float(&p, s);
2716 return float32_round_pack_canonical(&p, s);
2719 float64 bfloat16_to_float64(bfloat16 a, float_status *s)
2721 FloatParts64 p;
2723 bfloat16_unpack_canonical(&p, a, s);
2724 parts_float_to_float(&p, s);
2725 return float64_round_pack_canonical(&p, s);
2728 bfloat16 float32_to_bfloat16(float32 a, float_status *s)
2730 FloatParts64 p;
2732 float32_unpack_canonical(&p, a, s);
2733 parts_float_to_float(&p, s);
2734 return bfloat16_round_pack_canonical(&p, s);
2737 bfloat16 float64_to_bfloat16(float64 a, float_status *s)
2739 FloatParts64 p;
2741 float64_unpack_canonical(&p, a, s);
2742 parts_float_to_float(&p, s);
2743 return bfloat16_round_pack_canonical(&p, s);
2746 float32 float128_to_float32(float128 a, float_status *s)
2748 FloatParts64 p64;
2749 FloatParts128 p128;
2751 float128_unpack_canonical(&p128, a, s);
2752 parts_float_to_float_narrow(&p64, &p128, s);
2753 return float32_round_pack_canonical(&p64, s);
2756 float64 float128_to_float64(float128 a, float_status *s)
2758 FloatParts64 p64;
2759 FloatParts128 p128;
2761 float128_unpack_canonical(&p128, a, s);
2762 parts_float_to_float_narrow(&p64, &p128, s);
2763 return float64_round_pack_canonical(&p64, s);
2766 float128 float32_to_float128(float32 a, float_status *s)
2768 FloatParts64 p64;
2769 FloatParts128 p128;
2771 float32_unpack_canonical(&p64, a, s);
2772 parts_float_to_float_widen(&p128, &p64, s);
2773 return float128_round_pack_canonical(&p128, s);
2776 float128 float64_to_float128(float64 a, float_status *s)
2778 FloatParts64 p64;
2779 FloatParts128 p128;
2781 float64_unpack_canonical(&p64, a, s);
2782 parts_float_to_float_widen(&p128, &p64, s);
2783 return float128_round_pack_canonical(&p128, s);
2786 float32 floatx80_to_float32(floatx80 a, float_status *s)
2788 FloatParts64 p64;
2789 FloatParts128 p128;
2791 if (floatx80_unpack_canonical(&p128, a, s)) {
2792 parts_float_to_float_narrow(&p64, &p128, s);
2793 } else {
2794 parts_default_nan(&p64, s);
2796 return float32_round_pack_canonical(&p64, s);
2799 float64 floatx80_to_float64(floatx80 a, float_status *s)
2801 FloatParts64 p64;
2802 FloatParts128 p128;
2804 if (floatx80_unpack_canonical(&p128, a, s)) {
2805 parts_float_to_float_narrow(&p64, &p128, s);
2806 } else {
2807 parts_default_nan(&p64, s);
2809 return float64_round_pack_canonical(&p64, s);
2812 float128 floatx80_to_float128(floatx80 a, float_status *s)
2814 FloatParts128 p;
2816 if (floatx80_unpack_canonical(&p, a, s)) {
2817 parts_float_to_float(&p, s);
2818 } else {
2819 parts_default_nan(&p, s);
2821 return float128_round_pack_canonical(&p, s);
2824 floatx80 float32_to_floatx80(float32 a, float_status *s)
2826 FloatParts64 p64;
2827 FloatParts128 p128;
2829 float32_unpack_canonical(&p64, a, s);
2830 parts_float_to_float_widen(&p128, &p64, s);
2831 return floatx80_round_pack_canonical(&p128, s);
2834 floatx80 float64_to_floatx80(float64 a, float_status *s)
2836 FloatParts64 p64;
2837 FloatParts128 p128;
2839 float64_unpack_canonical(&p64, a, s);
2840 parts_float_to_float_widen(&p128, &p64, s);
2841 return floatx80_round_pack_canonical(&p128, s);
2844 floatx80 float128_to_floatx80(float128 a, float_status *s)
2846 FloatParts128 p;
2848 float128_unpack_canonical(&p, a, s);
2849 parts_float_to_float(&p, s);
2850 return floatx80_round_pack_canonical(&p, s);
2854 * Round to integral value
2857 float16 float16_round_to_int(float16 a, float_status *s)
2859 FloatParts64 p;
2861 float16_unpack_canonical(&p, a, s);
2862 parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float16_params);
2863 return float16_round_pack_canonical(&p, s);
2866 float32 float32_round_to_int(float32 a, float_status *s)
2868 FloatParts64 p;
2870 float32_unpack_canonical(&p, a, s);
2871 parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float32_params);
2872 return float32_round_pack_canonical(&p, s);
2875 float64 float64_round_to_int(float64 a, float_status *s)
2877 FloatParts64 p;
2879 float64_unpack_canonical(&p, a, s);
2880 parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float64_params);
2881 return float64_round_pack_canonical(&p, s);
2884 bfloat16 bfloat16_round_to_int(bfloat16 a, float_status *s)
2886 FloatParts64 p;
2888 bfloat16_unpack_canonical(&p, a, s);
2889 parts_round_to_int(&p, s->float_rounding_mode, 0, s, &bfloat16_params);
2890 return bfloat16_round_pack_canonical(&p, s);
2893 float128 float128_round_to_int(float128 a, float_status *s)
2895 FloatParts128 p;
2897 float128_unpack_canonical(&p, a, s);
2898 parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float128_params);
2899 return float128_round_pack_canonical(&p, s);
2902 floatx80 floatx80_round_to_int(floatx80 a, float_status *status)
2904 FloatParts128 p;
2906 if (!floatx80_unpack_canonical(&p, a, status)) {
2907 return floatx80_default_nan(status);
2910 parts_round_to_int(&p, status->float_rounding_mode, 0, status,
2911 &floatx80_params[status->floatx80_rounding_precision]);
2912 return floatx80_round_pack_canonical(&p, status);
2916 * Floating-point to signed integer conversions
2919 int8_t float16_to_int8_scalbn(float16 a, FloatRoundMode rmode, int scale,
2920 float_status *s)
2922 FloatParts64 p;
2924 float16_unpack_canonical(&p, a, s);
2925 return parts_float_to_sint(&p, rmode, scale, INT8_MIN, INT8_MAX, s);
2928 int16_t float16_to_int16_scalbn(float16 a, FloatRoundMode rmode, int scale,
2929 float_status *s)
2931 FloatParts64 p;
2933 float16_unpack_canonical(&p, a, s);
2934 return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
2937 int32_t float16_to_int32_scalbn(float16 a, FloatRoundMode rmode, int scale,
2938 float_status *s)
2940 FloatParts64 p;
2942 float16_unpack_canonical(&p, a, s);
2943 return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
2946 int64_t float16_to_int64_scalbn(float16 a, FloatRoundMode rmode, int scale,
2947 float_status *s)
2949 FloatParts64 p;
2951 float16_unpack_canonical(&p, a, s);
2952 return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
2955 int16_t float32_to_int16_scalbn(float32 a, FloatRoundMode rmode, int scale,
2956 float_status *s)
2958 FloatParts64 p;
2960 float32_unpack_canonical(&p, a, s);
2961 return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
2964 int32_t float32_to_int32_scalbn(float32 a, FloatRoundMode rmode, int scale,
2965 float_status *s)
2967 FloatParts64 p;
2969 float32_unpack_canonical(&p, a, s);
2970 return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
2973 int64_t float32_to_int64_scalbn(float32 a, FloatRoundMode rmode, int scale,
2974 float_status *s)
2976 FloatParts64 p;
2978 float32_unpack_canonical(&p, a, s);
2979 return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
2982 int16_t float64_to_int16_scalbn(float64 a, FloatRoundMode rmode, int scale,
2983 float_status *s)
2985 FloatParts64 p;
2987 float64_unpack_canonical(&p, a, s);
2988 return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
2991 int32_t float64_to_int32_scalbn(float64 a, FloatRoundMode rmode, int scale,
2992 float_status *s)
2994 FloatParts64 p;
2996 float64_unpack_canonical(&p, a, s);
2997 return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
3000 int64_t float64_to_int64_scalbn(float64 a, FloatRoundMode rmode, int scale,
3001 float_status *s)
3003 FloatParts64 p;
3005 float64_unpack_canonical(&p, a, s);
3006 return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
3009 int16_t bfloat16_to_int16_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
3010 float_status *s)
3012 FloatParts64 p;
3014 bfloat16_unpack_canonical(&p, a, s);
3015 return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
3018 int32_t bfloat16_to_int32_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
3019 float_status *s)
3021 FloatParts64 p;
3023 bfloat16_unpack_canonical(&p, a, s);
3024 return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
3027 int64_t bfloat16_to_int64_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
3028 float_status *s)
3030 FloatParts64 p;
3032 bfloat16_unpack_canonical(&p, a, s);
3033 return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
3036 static int32_t float128_to_int32_scalbn(float128 a, FloatRoundMode rmode,
3037 int scale, float_status *s)
3039 FloatParts128 p;
3041 float128_unpack_canonical(&p, a, s);
3042 return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
3045 static int64_t float128_to_int64_scalbn(float128 a, FloatRoundMode rmode,
3046 int scale, float_status *s)
3048 FloatParts128 p;
3050 float128_unpack_canonical(&p, a, s);
3051 return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
3054 static int32_t floatx80_to_int32_scalbn(floatx80 a, FloatRoundMode rmode,
3055 int scale, float_status *s)
3057 FloatParts128 p;
3059 if (!floatx80_unpack_canonical(&p, a, s)) {
3060 parts_default_nan(&p, s);
3062 return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
3065 static int64_t floatx80_to_int64_scalbn(floatx80 a, FloatRoundMode rmode,
3066 int scale, float_status *s)
3068 FloatParts128 p;
3070 if (!floatx80_unpack_canonical(&p, a, s)) {
3071 parts_default_nan(&p, s);
3073 return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
3076 int8_t float16_to_int8(float16 a, float_status *s)
3078 return float16_to_int8_scalbn(a, s->float_rounding_mode, 0, s);
3081 int16_t float16_to_int16(float16 a, float_status *s)
3083 return float16_to_int16_scalbn(a, s->float_rounding_mode, 0, s);
3086 int32_t float16_to_int32(float16 a, float_status *s)
3088 return float16_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3091 int64_t float16_to_int64(float16 a, float_status *s)
3093 return float16_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3096 int16_t float32_to_int16(float32 a, float_status *s)
3098 return float32_to_int16_scalbn(a, s->float_rounding_mode, 0, s);
3101 int32_t float32_to_int32(float32 a, float_status *s)
3103 return float32_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3106 int64_t float32_to_int64(float32 a, float_status *s)
3108 return float32_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3111 int16_t float64_to_int16(float64 a, float_status *s)
3113 return float64_to_int16_scalbn(a, s->float_rounding_mode, 0, s);
3116 int32_t float64_to_int32(float64 a, float_status *s)
3118 return float64_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3121 int64_t float64_to_int64(float64 a, float_status *s)
3123 return float64_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3126 int32_t float128_to_int32(float128 a, float_status *s)
3128 return float128_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3131 int64_t float128_to_int64(float128 a, float_status *s)
3133 return float128_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3136 int32_t floatx80_to_int32(floatx80 a, float_status *s)
3138 return floatx80_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3141 int64_t floatx80_to_int64(floatx80 a, float_status *s)
3143 return floatx80_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3146 int16_t float16_to_int16_round_to_zero(float16 a, float_status *s)
3148 return float16_to_int16_scalbn(a, float_round_to_zero, 0, s);
3151 int32_t float16_to_int32_round_to_zero(float16 a, float_status *s)
3153 return float16_to_int32_scalbn(a, float_round_to_zero, 0, s);
3156 int64_t float16_to_int64_round_to_zero(float16 a, float_status *s)
3158 return float16_to_int64_scalbn(a, float_round_to_zero, 0, s);
3161 int16_t float32_to_int16_round_to_zero(float32 a, float_status *s)
3163 return float32_to_int16_scalbn(a, float_round_to_zero, 0, s);
3166 int32_t float32_to_int32_round_to_zero(float32 a, float_status *s)
3168 return float32_to_int32_scalbn(a, float_round_to_zero, 0, s);
3171 int64_t float32_to_int64_round_to_zero(float32 a, float_status *s)
3173 return float32_to_int64_scalbn(a, float_round_to_zero, 0, s);
3176 int16_t float64_to_int16_round_to_zero(float64 a, float_status *s)
3178 return float64_to_int16_scalbn(a, float_round_to_zero, 0, s);
3181 int32_t float64_to_int32_round_to_zero(float64 a, float_status *s)
3183 return float64_to_int32_scalbn(a, float_round_to_zero, 0, s);
3186 int64_t float64_to_int64_round_to_zero(float64 a, float_status *s)
3188 return float64_to_int64_scalbn(a, float_round_to_zero, 0, s);
3191 int32_t float128_to_int32_round_to_zero(float128 a, float_status *s)
3193 return float128_to_int32_scalbn(a, float_round_to_zero, 0, s);
3196 int64_t float128_to_int64_round_to_zero(float128 a, float_status *s)
3198 return float128_to_int64_scalbn(a, float_round_to_zero, 0, s);
3201 int32_t floatx80_to_int32_round_to_zero(floatx80 a, float_status *s)
3203 return floatx80_to_int32_scalbn(a, float_round_to_zero, 0, s);
3206 int64_t floatx80_to_int64_round_to_zero(floatx80 a, float_status *s)
3208 return floatx80_to_int64_scalbn(a, float_round_to_zero, 0, s);
3211 int16_t bfloat16_to_int16(bfloat16 a, float_status *s)
3213 return bfloat16_to_int16_scalbn(a, s->float_rounding_mode, 0, s);
3216 int32_t bfloat16_to_int32(bfloat16 a, float_status *s)
3218 return bfloat16_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
3221 int64_t bfloat16_to_int64(bfloat16 a, float_status *s)
3223 return bfloat16_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
3226 int16_t bfloat16_to_int16_round_to_zero(bfloat16 a, float_status *s)
3228 return bfloat16_to_int16_scalbn(a, float_round_to_zero, 0, s);
3231 int32_t bfloat16_to_int32_round_to_zero(bfloat16 a, float_status *s)
3233 return bfloat16_to_int32_scalbn(a, float_round_to_zero, 0, s);
3236 int64_t bfloat16_to_int64_round_to_zero(bfloat16 a, float_status *s)
3238 return bfloat16_to_int64_scalbn(a, float_round_to_zero, 0, s);
3242 * Floating-point to unsigned integer conversions
3245 uint8_t float16_to_uint8_scalbn(float16 a, FloatRoundMode rmode, int scale,
3246 float_status *s)
3248 FloatParts64 p;
3250 float16_unpack_canonical(&p, a, s);
3251 return parts_float_to_uint(&p, rmode, scale, UINT8_MAX, s);
3254 uint16_t float16_to_uint16_scalbn(float16 a, FloatRoundMode rmode, int scale,
3255 float_status *s)
3257 FloatParts64 p;
3259 float16_unpack_canonical(&p, a, s);
3260 return parts_float_to_uint(&p, rmode, scale, UINT16_MAX, s);
3263 uint32_t float16_to_uint32_scalbn(float16 a, FloatRoundMode rmode, int scale,
3264 float_status *s)
3266 FloatParts64 p;
3268 float16_unpack_canonical(&p, a, s);
3269 return parts_float_to_uint(&p, rmode, scale, UINT32_MAX, s);
3272 uint64_t float16_to_uint64_scalbn(float16 a, FloatRoundMode rmode, int scale,
3273 float_status *s)
3275 FloatParts64 p;
3277 float16_unpack_canonical(&p, a, s);
3278 return parts_float_to_uint(&p, rmode, scale, UINT64_MAX, s);
3281 uint16_t float32_to_uint16_scalbn(float32 a, FloatRoundMode rmode, int scale,
3282 float_status *s)
3284 FloatParts64 p;
3286 float32_unpack_canonical(&p, a, s);
3287 return parts_float_to_uint(&p, rmode, scale, UINT16_MAX, s);
3290 uint32_t float32_to_uint32_scalbn(float32 a, FloatRoundMode rmode, int scale,
3291 float_status *s)
3293 FloatParts64 p;
3295 float32_unpack_canonical(&p, a, s);
3296 return parts_float_to_uint(&p, rmode, scale, UINT32_MAX, s);
3299 uint64_t float32_to_uint64_scalbn(float32 a, FloatRoundMode rmode, int scale,
3300 float_status *s)
3302 FloatParts64 p;
3304 float32_unpack_canonical(&p, a, s);
3305 return parts_float_to_uint(&p, rmode, scale, UINT64_MAX, s);
3308 uint16_t float64_to_uint16_scalbn(float64 a, FloatRoundMode rmode, int scale,
3309 float_status *s)
3311 FloatParts64 p;
3313 float64_unpack_canonical(&p, a, s);
3314 return parts_float_to_uint(&p, rmode, scale, UINT16_MAX, s);
3317 uint32_t float64_to_uint32_scalbn(float64 a, FloatRoundMode rmode, int scale,
3318 float_status *s)
3320 FloatParts64 p;
3322 float64_unpack_canonical(&p, a, s);
3323 return parts_float_to_uint(&p, rmode, scale, UINT32_MAX, s);
3326 uint64_t float64_to_uint64_scalbn(float64 a, FloatRoundMode rmode, int scale,
3327 float_status *s)
3329 FloatParts64 p;
3331 float64_unpack_canonical(&p, a, s);
3332 return parts_float_to_uint(&p, rmode, scale, UINT64_MAX, s);
3335 uint16_t bfloat16_to_uint16_scalbn(bfloat16 a, FloatRoundMode rmode,
3336 int scale, float_status *s)
3338 FloatParts64 p;
3340 bfloat16_unpack_canonical(&p, a, s);
3341 return parts_float_to_uint(&p, rmode, scale, UINT16_MAX, s);
3344 uint32_t bfloat16_to_uint32_scalbn(bfloat16 a, FloatRoundMode rmode,
3345 int scale, float_status *s)
3347 FloatParts64 p;
3349 bfloat16_unpack_canonical(&p, a, s);
3350 return parts_float_to_uint(&p, rmode, scale, UINT32_MAX, s);
3353 uint64_t bfloat16_to_uint64_scalbn(bfloat16 a, FloatRoundMode rmode,
3354 int scale, float_status *s)
3356 FloatParts64 p;
3358 bfloat16_unpack_canonical(&p, a, s);
3359 return parts_float_to_uint(&p, rmode, scale, UINT64_MAX, s);
3362 static uint32_t float128_to_uint32_scalbn(float128 a, FloatRoundMode rmode,
3363 int scale, float_status *s)
3365 FloatParts128 p;
3367 float128_unpack_canonical(&p, a, s);
3368 return parts_float_to_uint(&p, rmode, scale, UINT32_MAX, s);
3371 static uint64_t float128_to_uint64_scalbn(float128 a, FloatRoundMode rmode,
3372 int scale, float_status *s)
3374 FloatParts128 p;
3376 float128_unpack_canonical(&p, a, s);
3377 return parts_float_to_uint(&p, rmode, scale, UINT64_MAX, s);
3380 uint8_t float16_to_uint8(float16 a, float_status *s)
3382 return float16_to_uint8_scalbn(a, s->float_rounding_mode, 0, s);
3385 uint16_t float16_to_uint16(float16 a, float_status *s)
3387 return float16_to_uint16_scalbn(a, s->float_rounding_mode, 0, s);
3390 uint32_t float16_to_uint32(float16 a, float_status *s)
3392 return float16_to_uint32_scalbn(a, s->float_rounding_mode, 0, s);
3395 uint64_t float16_to_uint64(float16 a, float_status *s)
3397 return float16_to_uint64_scalbn(a, s->float_rounding_mode, 0, s);
3400 uint16_t float32_to_uint16(float32 a, float_status *s)
3402 return float32_to_uint16_scalbn(a, s->float_rounding_mode, 0, s);
3405 uint32_t float32_to_uint32(float32 a, float_status *s)
3407 return float32_to_uint32_scalbn(a, s->float_rounding_mode, 0, s);
3410 uint64_t float32_to_uint64(float32 a, float_status *s)
3412 return float32_to_uint64_scalbn(a, s->float_rounding_mode, 0, s);
3415 uint16_t float64_to_uint16(float64 a, float_status *s)
3417 return float64_to_uint16_scalbn(a, s->float_rounding_mode, 0, s);
3420 uint32_t float64_to_uint32(float64 a, float_status *s)
3422 return float64_to_uint32_scalbn(a, s->float_rounding_mode, 0, s);
3425 uint64_t float64_to_uint64(float64 a, float_status *s)
3427 return float64_to_uint64_scalbn(a, s->float_rounding_mode, 0, s);
3430 uint32_t float128_to_uint32(float128 a, float_status *s)
3432 return float128_to_uint32_scalbn(a, s->float_rounding_mode, 0, s);
3435 uint64_t float128_to_uint64(float128 a, float_status *s)
3437 return float128_to_uint64_scalbn(a, s->float_rounding_mode, 0, s);
3440 uint16_t float16_to_uint16_round_to_zero(float16 a, float_status *s)
3442 return float16_to_uint16_scalbn(a, float_round_to_zero, 0, s);
3445 uint32_t float16_to_uint32_round_to_zero(float16 a, float_status *s)
3447 return float16_to_uint32_scalbn(a, float_round_to_zero, 0, s);
3450 uint64_t float16_to_uint64_round_to_zero(float16 a, float_status *s)
3452 return float16_to_uint64_scalbn(a, float_round_to_zero, 0, s);
3455 uint16_t float32_to_uint16_round_to_zero(float32 a, float_status *s)
3457 return float32_to_uint16_scalbn(a, float_round_to_zero, 0, s);
3460 uint32_t float32_to_uint32_round_to_zero(float32 a, float_status *s)
3462 return float32_to_uint32_scalbn(a, float_round_to_zero, 0, s);
3465 uint64_t float32_to_uint64_round_to_zero(float32 a, float_status *s)
3467 return float32_to_uint64_scalbn(a, float_round_to_zero, 0, s);
3470 uint16_t float64_to_uint16_round_to_zero(float64 a, float_status *s)
3472 return float64_to_uint16_scalbn(a, float_round_to_zero, 0, s);
3475 uint32_t float64_to_uint32_round_to_zero(float64 a, float_status *s)
3477 return float64_to_uint32_scalbn(a, float_round_to_zero, 0, s);
3480 uint64_t float64_to_uint64_round_to_zero(float64 a, float_status *s)
3482 return float64_to_uint64_scalbn(a, float_round_to_zero, 0, s);
3485 uint32_t float128_to_uint32_round_to_zero(float128 a, float_status *s)
3487 return float128_to_uint32_scalbn(a, float_round_to_zero, 0, s);
3490 uint64_t float128_to_uint64_round_to_zero(float128 a, float_status *s)
3492 return float128_to_uint64_scalbn(a, float_round_to_zero, 0, s);
3495 uint16_t bfloat16_to_uint16(bfloat16 a, float_status *s)
3497 return bfloat16_to_uint16_scalbn(a, s->float_rounding_mode, 0, s);
3500 uint32_t bfloat16_to_uint32(bfloat16 a, float_status *s)
3502 return bfloat16_to_uint32_scalbn(a, s->float_rounding_mode, 0, s);
3505 uint64_t bfloat16_to_uint64(bfloat16 a, float_status *s)
3507 return bfloat16_to_uint64_scalbn(a, s->float_rounding_mode, 0, s);
3510 uint16_t bfloat16_to_uint16_round_to_zero(bfloat16 a, float_status *s)
3512 return bfloat16_to_uint16_scalbn(a, float_round_to_zero, 0, s);
3515 uint32_t bfloat16_to_uint32_round_to_zero(bfloat16 a, float_status *s)
3517 return bfloat16_to_uint32_scalbn(a, float_round_to_zero, 0, s);
3520 uint64_t bfloat16_to_uint64_round_to_zero(bfloat16 a, float_status *s)
3522 return bfloat16_to_uint64_scalbn(a, float_round_to_zero, 0, s);
3526 * Signed integer to floating-point conversions
3529 float16 int64_to_float16_scalbn(int64_t a, int scale, float_status *status)
3531 FloatParts64 p;
3533 parts_sint_to_float(&p, a, scale, status);
3534 return float16_round_pack_canonical(&p, status);
3537 float16 int32_to_float16_scalbn(int32_t a, int scale, float_status *status)
3539 return int64_to_float16_scalbn(a, scale, status);
3542 float16 int16_to_float16_scalbn(int16_t a, int scale, float_status *status)
3544 return int64_to_float16_scalbn(a, scale, status);
3547 float16 int64_to_float16(int64_t a, float_status *status)
3549 return int64_to_float16_scalbn(a, 0, status);
3552 float16 int32_to_float16(int32_t a, float_status *status)
3554 return int64_to_float16_scalbn(a, 0, status);
3557 float16 int16_to_float16(int16_t a, float_status *status)
3559 return int64_to_float16_scalbn(a, 0, status);
3562 float16 int8_to_float16(int8_t a, float_status *status)
3564 return int64_to_float16_scalbn(a, 0, status);
3567 float32 int64_to_float32_scalbn(int64_t a, int scale, float_status *status)
3569 FloatParts64 p;
3571 /* Without scaling, there are no overflow concerns. */
3572 if (likely(scale == 0) && can_use_fpu(status)) {
3573 union_float32 ur;
3574 ur.h = a;
3575 return ur.s;
3578 parts64_sint_to_float(&p, a, scale, status);
3579 return float32_round_pack_canonical(&p, status);
3582 float32 int32_to_float32_scalbn(int32_t a, int scale, float_status *status)
3584 return int64_to_float32_scalbn(a, scale, status);
3587 float32 int16_to_float32_scalbn(int16_t a, int scale, float_status *status)
3589 return int64_to_float32_scalbn(a, scale, status);
3592 float32 int64_to_float32(int64_t a, float_status *status)
3594 return int64_to_float32_scalbn(a, 0, status);
3597 float32 int32_to_float32(int32_t a, float_status *status)
3599 return int64_to_float32_scalbn(a, 0, status);
3602 float32 int16_to_float32(int16_t a, float_status *status)
3604 return int64_to_float32_scalbn(a, 0, status);
3607 float64 int64_to_float64_scalbn(int64_t a, int scale, float_status *status)
3609 FloatParts64 p;
3611 /* Without scaling, there are no overflow concerns. */
3612 if (likely(scale == 0) && can_use_fpu(status)) {
3613 union_float64 ur;
3614 ur.h = a;
3615 return ur.s;
3618 parts_sint_to_float(&p, a, scale, status);
3619 return float64_round_pack_canonical(&p, status);
3622 float64 int32_to_float64_scalbn(int32_t a, int scale, float_status *status)
3624 return int64_to_float64_scalbn(a, scale, status);
3627 float64 int16_to_float64_scalbn(int16_t a, int scale, float_status *status)
3629 return int64_to_float64_scalbn(a, scale, status);
3632 float64 int64_to_float64(int64_t a, float_status *status)
3634 return int64_to_float64_scalbn(a, 0, status);
3637 float64 int32_to_float64(int32_t a, float_status *status)
3639 return int64_to_float64_scalbn(a, 0, status);
3642 float64 int16_to_float64(int16_t a, float_status *status)
3644 return int64_to_float64_scalbn(a, 0, status);
3647 bfloat16 int64_to_bfloat16_scalbn(int64_t a, int scale, float_status *status)
3649 FloatParts64 p;
3651 parts_sint_to_float(&p, a, scale, status);
3652 return bfloat16_round_pack_canonical(&p, status);
3655 bfloat16 int32_to_bfloat16_scalbn(int32_t a, int scale, float_status *status)
3657 return int64_to_bfloat16_scalbn(a, scale, status);
3660 bfloat16 int16_to_bfloat16_scalbn(int16_t a, int scale, float_status *status)
3662 return int64_to_bfloat16_scalbn(a, scale, status);
3665 bfloat16 int64_to_bfloat16(int64_t a, float_status *status)
3667 return int64_to_bfloat16_scalbn(a, 0, status);
3670 bfloat16 int32_to_bfloat16(int32_t a, float_status *status)
3672 return int64_to_bfloat16_scalbn(a, 0, status);
3675 bfloat16 int16_to_bfloat16(int16_t a, float_status *status)
3677 return int64_to_bfloat16_scalbn(a, 0, status);
3680 float128 int64_to_float128(int64_t a, float_status *status)
3682 FloatParts128 p;
3684 parts_sint_to_float(&p, a, 0, status);
3685 return float128_round_pack_canonical(&p, status);
3688 float128 int32_to_float128(int32_t a, float_status *status)
3690 return int64_to_float128(a, status);
3693 floatx80 int64_to_floatx80(int64_t a, float_status *status)
3695 FloatParts128 p;
3697 parts_sint_to_float(&p, a, 0, status);
3698 return floatx80_round_pack_canonical(&p, status);
3701 floatx80 int32_to_floatx80(int32_t a, float_status *status)
3703 return int64_to_floatx80(a, status);
3707 * Unsigned Integer to floating-point conversions
3710 float16 uint64_to_float16_scalbn(uint64_t a, int scale, float_status *status)
3712 FloatParts64 p;
3714 parts_uint_to_float(&p, a, scale, status);
3715 return float16_round_pack_canonical(&p, status);
3718 float16 uint32_to_float16_scalbn(uint32_t a, int scale, float_status *status)
3720 return uint64_to_float16_scalbn(a, scale, status);
3723 float16 uint16_to_float16_scalbn(uint16_t a, int scale, float_status *status)
3725 return uint64_to_float16_scalbn(a, scale, status);
3728 float16 uint64_to_float16(uint64_t a, float_status *status)
3730 return uint64_to_float16_scalbn(a, 0, status);
3733 float16 uint32_to_float16(uint32_t a, float_status *status)
3735 return uint64_to_float16_scalbn(a, 0, status);
3738 float16 uint16_to_float16(uint16_t a, float_status *status)
3740 return uint64_to_float16_scalbn(a, 0, status);
3743 float16 uint8_to_float16(uint8_t a, float_status *status)
3745 return uint64_to_float16_scalbn(a, 0, status);
3748 float32 uint64_to_float32_scalbn(uint64_t a, int scale, float_status *status)
3750 FloatParts64 p;
3752 /* Without scaling, there are no overflow concerns. */
3753 if (likely(scale == 0) && can_use_fpu(status)) {
3754 union_float32 ur;
3755 ur.h = a;
3756 return ur.s;
3759 parts_uint_to_float(&p, a, scale, status);
3760 return float32_round_pack_canonical(&p, status);
3763 float32 uint32_to_float32_scalbn(uint32_t a, int scale, float_status *status)
3765 return uint64_to_float32_scalbn(a, scale, status);
3768 float32 uint16_to_float32_scalbn(uint16_t a, int scale, float_status *status)
3770 return uint64_to_float32_scalbn(a, scale, status);
3773 float32 uint64_to_float32(uint64_t a, float_status *status)
3775 return uint64_to_float32_scalbn(a, 0, status);
3778 float32 uint32_to_float32(uint32_t a, float_status *status)
3780 return uint64_to_float32_scalbn(a, 0, status);
3783 float32 uint16_to_float32(uint16_t a, float_status *status)
3785 return uint64_to_float32_scalbn(a, 0, status);
3788 float64 uint64_to_float64_scalbn(uint64_t a, int scale, float_status *status)
3790 FloatParts64 p;
3792 /* Without scaling, there are no overflow concerns. */
3793 if (likely(scale == 0) && can_use_fpu(status)) {
3794 union_float64 ur;
3795 ur.h = a;
3796 return ur.s;
3799 parts_uint_to_float(&p, a, scale, status);
3800 return float64_round_pack_canonical(&p, status);
3803 float64 uint32_to_float64_scalbn(uint32_t a, int scale, float_status *status)
3805 return uint64_to_float64_scalbn(a, scale, status);
3808 float64 uint16_to_float64_scalbn(uint16_t a, int scale, float_status *status)
3810 return uint64_to_float64_scalbn(a, scale, status);
3813 float64 uint64_to_float64(uint64_t a, float_status *status)
3815 return uint64_to_float64_scalbn(a, 0, status);
3818 float64 uint32_to_float64(uint32_t a, float_status *status)
3820 return uint64_to_float64_scalbn(a, 0, status);
3823 float64 uint16_to_float64(uint16_t a, float_status *status)
3825 return uint64_to_float64_scalbn(a, 0, status);
3828 bfloat16 uint64_to_bfloat16_scalbn(uint64_t a, int scale, float_status *status)
3830 FloatParts64 p;
3832 parts_uint_to_float(&p, a, scale, status);
3833 return bfloat16_round_pack_canonical(&p, status);
3836 bfloat16 uint32_to_bfloat16_scalbn(uint32_t a, int scale, float_status *status)
3838 return uint64_to_bfloat16_scalbn(a, scale, status);
3841 bfloat16 uint16_to_bfloat16_scalbn(uint16_t a, int scale, float_status *status)
3843 return uint64_to_bfloat16_scalbn(a, scale, status);
3846 bfloat16 uint64_to_bfloat16(uint64_t a, float_status *status)
3848 return uint64_to_bfloat16_scalbn(a, 0, status);
3851 bfloat16 uint32_to_bfloat16(uint32_t a, float_status *status)
3853 return uint64_to_bfloat16_scalbn(a, 0, status);
3856 bfloat16 uint16_to_bfloat16(uint16_t a, float_status *status)
3858 return uint64_to_bfloat16_scalbn(a, 0, status);
3861 float128 uint64_to_float128(uint64_t a, float_status *status)
3863 FloatParts128 p;
3865 parts_uint_to_float(&p, a, 0, status);
3866 return float128_round_pack_canonical(&p, status);
3870 * Minimum and maximum
3873 static float16 float16_minmax(float16 a, float16 b, float_status *s, int flags)
3875 FloatParts64 pa, pb, *pr;
3877 float16_unpack_canonical(&pa, a, s);
3878 float16_unpack_canonical(&pb, b, s);
3879 pr = parts_minmax(&pa, &pb, s, flags);
3881 return float16_round_pack_canonical(pr, s);
3884 static bfloat16 bfloat16_minmax(bfloat16 a, bfloat16 b,
3885 float_status *s, int flags)
3887 FloatParts64 pa, pb, *pr;
3889 bfloat16_unpack_canonical(&pa, a, s);
3890 bfloat16_unpack_canonical(&pb, b, s);
3891 pr = parts_minmax(&pa, &pb, s, flags);
3893 return bfloat16_round_pack_canonical(pr, s);
3896 static float32 float32_minmax(float32 a, float32 b, float_status *s, int flags)
3898 FloatParts64 pa, pb, *pr;
3900 float32_unpack_canonical(&pa, a, s);
3901 float32_unpack_canonical(&pb, b, s);
3902 pr = parts_minmax(&pa, &pb, s, flags);
3904 return float32_round_pack_canonical(pr, s);
3907 static float64 float64_minmax(float64 a, float64 b, float_status *s, int flags)
3909 FloatParts64 pa, pb, *pr;
3911 float64_unpack_canonical(&pa, a, s);
3912 float64_unpack_canonical(&pb, b, s);
3913 pr = parts_minmax(&pa, &pb, s, flags);
3915 return float64_round_pack_canonical(pr, s);
3918 static float128 float128_minmax(float128 a, float128 b,
3919 float_status *s, int flags)
3921 FloatParts128 pa, pb, *pr;
3923 float128_unpack_canonical(&pa, a, s);
3924 float128_unpack_canonical(&pb, b, s);
3925 pr = parts_minmax(&pa, &pb, s, flags);
3927 return float128_round_pack_canonical(pr, s);
3930 #define MINMAX_1(type, name, flags) \
3931 type type##_##name(type a, type b, float_status *s) \
3932 { return type##_minmax(a, b, s, flags); }
3934 #define MINMAX_2(type) \
3935 MINMAX_1(type, max, 0) \
3936 MINMAX_1(type, maxnum, minmax_isnum) \
3937 MINMAX_1(type, maxnummag, minmax_isnum | minmax_ismag) \
3938 MINMAX_1(type, maximum_number, minmax_isnumber) \
3939 MINMAX_1(type, min, minmax_ismin) \
3940 MINMAX_1(type, minnum, minmax_ismin | minmax_isnum) \
3941 MINMAX_1(type, minnummag, minmax_ismin | minmax_isnum | minmax_ismag) \
3942 MINMAX_1(type, minimum_number, minmax_ismin | minmax_isnumber) \
3944 MINMAX_2(float16)
3945 MINMAX_2(bfloat16)
3946 MINMAX_2(float32)
3947 MINMAX_2(float64)
3948 MINMAX_2(float128)
3950 #undef MINMAX_1
3951 #undef MINMAX_2
3954 * Floating point compare
3957 static FloatRelation QEMU_FLATTEN
3958 float16_do_compare(float16 a, float16 b, float_status *s, bool is_quiet)
3960 FloatParts64 pa, pb;
3962 float16_unpack_canonical(&pa, a, s);
3963 float16_unpack_canonical(&pb, b, s);
3964 return parts_compare(&pa, &pb, s, is_quiet);
3967 FloatRelation float16_compare(float16 a, float16 b, float_status *s)
3969 return float16_do_compare(a, b, s, false);
3972 FloatRelation float16_compare_quiet(float16 a, float16 b, float_status *s)
3974 return float16_do_compare(a, b, s, true);
3977 static FloatRelation QEMU_SOFTFLOAT_ATTR
3978 float32_do_compare(float32 a, float32 b, float_status *s, bool is_quiet)
3980 FloatParts64 pa, pb;
3982 float32_unpack_canonical(&pa, a, s);
3983 float32_unpack_canonical(&pb, b, s);
3984 return parts_compare(&pa, &pb, s, is_quiet);
3987 static FloatRelation QEMU_FLATTEN
3988 float32_hs_compare(float32 xa, float32 xb, float_status *s, bool is_quiet)
3990 union_float32 ua, ub;
3992 ua.s = xa;
3993 ub.s = xb;
3995 if (QEMU_NO_HARDFLOAT) {
3996 goto soft;
3999 float32_input_flush2(&ua.s, &ub.s, s);
4000 if (isgreaterequal(ua.h, ub.h)) {
4001 if (isgreater(ua.h, ub.h)) {
4002 return float_relation_greater;
4004 return float_relation_equal;
4006 if (likely(isless(ua.h, ub.h))) {
4007 return float_relation_less;
4010 * The only condition remaining is unordered.
4011 * Fall through to set flags.
4013 soft:
4014 return float32_do_compare(ua.s, ub.s, s, is_quiet);
4017 FloatRelation float32_compare(float32 a, float32 b, float_status *s)
4019 return float32_hs_compare(a, b, s, false);
4022 FloatRelation float32_compare_quiet(float32 a, float32 b, float_status *s)
4024 return float32_hs_compare(a, b, s, true);
4027 static FloatRelation QEMU_SOFTFLOAT_ATTR
4028 float64_do_compare(float64 a, float64 b, float_status *s, bool is_quiet)
4030 FloatParts64 pa, pb;
4032 float64_unpack_canonical(&pa, a, s);
4033 float64_unpack_canonical(&pb, b, s);
4034 return parts_compare(&pa, &pb, s, is_quiet);
4037 static FloatRelation QEMU_FLATTEN
4038 float64_hs_compare(float64 xa, float64 xb, float_status *s, bool is_quiet)
4040 union_float64 ua, ub;
4042 ua.s = xa;
4043 ub.s = xb;
4045 if (QEMU_NO_HARDFLOAT) {
4046 goto soft;
4049 float64_input_flush2(&ua.s, &ub.s, s);
4050 if (isgreaterequal(ua.h, ub.h)) {
4051 if (isgreater(ua.h, ub.h)) {
4052 return float_relation_greater;
4054 return float_relation_equal;
4056 if (likely(isless(ua.h, ub.h))) {
4057 return float_relation_less;
4060 * The only condition remaining is unordered.
4061 * Fall through to set flags.
4063 soft:
4064 return float64_do_compare(ua.s, ub.s, s, is_quiet);
4067 FloatRelation float64_compare(float64 a, float64 b, float_status *s)
4069 return float64_hs_compare(a, b, s, false);
4072 FloatRelation float64_compare_quiet(float64 a, float64 b, float_status *s)
4074 return float64_hs_compare(a, b, s, true);
4077 static FloatRelation QEMU_FLATTEN
4078 bfloat16_do_compare(bfloat16 a, bfloat16 b, float_status *s, bool is_quiet)
4080 FloatParts64 pa, pb;
4082 bfloat16_unpack_canonical(&pa, a, s);
4083 bfloat16_unpack_canonical(&pb, b, s);
4084 return parts_compare(&pa, &pb, s, is_quiet);
4087 FloatRelation bfloat16_compare(bfloat16 a, bfloat16 b, float_status *s)
4089 return bfloat16_do_compare(a, b, s, false);
4092 FloatRelation bfloat16_compare_quiet(bfloat16 a, bfloat16 b, float_status *s)
4094 return bfloat16_do_compare(a, b, s, true);
4097 static FloatRelation QEMU_FLATTEN
4098 float128_do_compare(float128 a, float128 b, float_status *s, bool is_quiet)
4100 FloatParts128 pa, pb;
4102 float128_unpack_canonical(&pa, a, s);
4103 float128_unpack_canonical(&pb, b, s);
4104 return parts_compare(&pa, &pb, s, is_quiet);
4107 FloatRelation float128_compare(float128 a, float128 b, float_status *s)
4109 return float128_do_compare(a, b, s, false);
4112 FloatRelation float128_compare_quiet(float128 a, float128 b, float_status *s)
4114 return float128_do_compare(a, b, s, true);
4117 static FloatRelation QEMU_FLATTEN
4118 floatx80_do_compare(floatx80 a, floatx80 b, float_status *s, bool is_quiet)
4120 FloatParts128 pa, pb;
4122 if (!floatx80_unpack_canonical(&pa, a, s) ||
4123 !floatx80_unpack_canonical(&pb, b, s)) {
4124 return float_relation_unordered;
4126 return parts_compare(&pa, &pb, s, is_quiet);
4129 FloatRelation floatx80_compare(floatx80 a, floatx80 b, float_status *s)
4131 return floatx80_do_compare(a, b, s, false);
4134 FloatRelation floatx80_compare_quiet(floatx80 a, floatx80 b, float_status *s)
4136 return floatx80_do_compare(a, b, s, true);
4140 * Scale by 2**N
4143 float16 float16_scalbn(float16 a, int n, float_status *status)
4145 FloatParts64 p;
4147 float16_unpack_canonical(&p, a, status);
4148 parts_scalbn(&p, n, status);
4149 return float16_round_pack_canonical(&p, status);
4152 float32 float32_scalbn(float32 a, int n, float_status *status)
4154 FloatParts64 p;
4156 float32_unpack_canonical(&p, a, status);
4157 parts_scalbn(&p, n, status);
4158 return float32_round_pack_canonical(&p, status);
4161 float64 float64_scalbn(float64 a, int n, float_status *status)
4163 FloatParts64 p;
4165 float64_unpack_canonical(&p, a, status);
4166 parts_scalbn(&p, n, status);
4167 return float64_round_pack_canonical(&p, status);
4170 bfloat16 bfloat16_scalbn(bfloat16 a, int n, float_status *status)
4172 FloatParts64 p;
4174 bfloat16_unpack_canonical(&p, a, status);
4175 parts_scalbn(&p, n, status);
4176 return bfloat16_round_pack_canonical(&p, status);
4179 float128 float128_scalbn(float128 a, int n, float_status *status)
4181 FloatParts128 p;
4183 float128_unpack_canonical(&p, a, status);
4184 parts_scalbn(&p, n, status);
4185 return float128_round_pack_canonical(&p, status);
4188 floatx80 floatx80_scalbn(floatx80 a, int n, float_status *status)
4190 FloatParts128 p;
4192 if (!floatx80_unpack_canonical(&p, a, status)) {
4193 return floatx80_default_nan(status);
4195 parts_scalbn(&p, n, status);
4196 return floatx80_round_pack_canonical(&p, status);
4200 * Square Root
4203 float16 QEMU_FLATTEN float16_sqrt(float16 a, float_status *status)
4205 FloatParts64 p;
4207 float16_unpack_canonical(&p, a, status);
4208 parts_sqrt(&p, status, &float16_params);
4209 return float16_round_pack_canonical(&p, status);
4212 static float32 QEMU_SOFTFLOAT_ATTR
4213 soft_f32_sqrt(float32 a, float_status *status)
4215 FloatParts64 p;
4217 float32_unpack_canonical(&p, a, status);
4218 parts_sqrt(&p, status, &float32_params);
4219 return float32_round_pack_canonical(&p, status);
4222 static float64 QEMU_SOFTFLOAT_ATTR
4223 soft_f64_sqrt(float64 a, float_status *status)
4225 FloatParts64 p;
4227 float64_unpack_canonical(&p, a, status);
4228 parts_sqrt(&p, status, &float64_params);
4229 return float64_round_pack_canonical(&p, status);
4232 float32 QEMU_FLATTEN float32_sqrt(float32 xa, float_status *s)
4234 union_float32 ua, ur;
4236 ua.s = xa;
4237 if (unlikely(!can_use_fpu(s))) {
4238 goto soft;
4241 float32_input_flush1(&ua.s, s);
4242 if (QEMU_HARDFLOAT_1F32_USE_FP) {
4243 if (unlikely(!(fpclassify(ua.h) == FP_NORMAL ||
4244 fpclassify(ua.h) == FP_ZERO) ||
4245 signbit(ua.h))) {
4246 goto soft;
4248 } else if (unlikely(!float32_is_zero_or_normal(ua.s) ||
4249 float32_is_neg(ua.s))) {
4250 goto soft;
4252 ur.h = sqrtf(ua.h);
4253 return ur.s;
4255 soft:
4256 return soft_f32_sqrt(ua.s, s);
4259 float64 QEMU_FLATTEN float64_sqrt(float64 xa, float_status *s)
4261 union_float64 ua, ur;
4263 ua.s = xa;
4264 if (unlikely(!can_use_fpu(s))) {
4265 goto soft;
4268 float64_input_flush1(&ua.s, s);
4269 if (QEMU_HARDFLOAT_1F64_USE_FP) {
4270 if (unlikely(!(fpclassify(ua.h) == FP_NORMAL ||
4271 fpclassify(ua.h) == FP_ZERO) ||
4272 signbit(ua.h))) {
4273 goto soft;
4275 } else if (unlikely(!float64_is_zero_or_normal(ua.s) ||
4276 float64_is_neg(ua.s))) {
4277 goto soft;
4279 ur.h = sqrt(ua.h);
4280 return ur.s;
4282 soft:
4283 return soft_f64_sqrt(ua.s, s);
4286 bfloat16 QEMU_FLATTEN bfloat16_sqrt(bfloat16 a, float_status *status)
4288 FloatParts64 p;
4290 bfloat16_unpack_canonical(&p, a, status);
4291 parts_sqrt(&p, status, &bfloat16_params);
4292 return bfloat16_round_pack_canonical(&p, status);
4295 float128 QEMU_FLATTEN float128_sqrt(float128 a, float_status *status)
4297 FloatParts128 p;
4299 float128_unpack_canonical(&p, a, status);
4300 parts_sqrt(&p, status, &float128_params);
4301 return float128_round_pack_canonical(&p, status);
4304 floatx80 floatx80_sqrt(floatx80 a, float_status *s)
4306 FloatParts128 p;
4308 if (!floatx80_unpack_canonical(&p, a, s)) {
4309 return floatx80_default_nan(s);
4311 parts_sqrt(&p, s, &floatx80_params[s->floatx80_rounding_precision]);
4312 return floatx80_round_pack_canonical(&p, s);
4316 * log2
4318 float32 float32_log2(float32 a, float_status *status)
4320 FloatParts64 p;
4322 float32_unpack_canonical(&p, a, status);
4323 parts_log2(&p, status, &float32_params);
4324 return float32_round_pack_canonical(&p, status);
4327 float64 float64_log2(float64 a, float_status *status)
4329 FloatParts64 p;
4331 float64_unpack_canonical(&p, a, status);
4332 parts_log2(&p, status, &float64_params);
4333 return float64_round_pack_canonical(&p, status);
4336 /*----------------------------------------------------------------------------
4337 | The pattern for a default generated NaN.
4338 *----------------------------------------------------------------------------*/
4340 float16 float16_default_nan(float_status *status)
4342 FloatParts64 p;
4344 parts_default_nan(&p, status);
4345 p.frac >>= float16_params.frac_shift;
4346 return float16_pack_raw(&p);
4349 float32 float32_default_nan(float_status *status)
4351 FloatParts64 p;
4353 parts_default_nan(&p, status);
4354 p.frac >>= float32_params.frac_shift;
4355 return float32_pack_raw(&p);
4358 float64 float64_default_nan(float_status *status)
4360 FloatParts64 p;
4362 parts_default_nan(&p, status);
4363 p.frac >>= float64_params.frac_shift;
4364 return float64_pack_raw(&p);
4367 float128 float128_default_nan(float_status *status)
4369 FloatParts128 p;
4371 parts_default_nan(&p, status);
4372 frac_shr(&p, float128_params.frac_shift);
4373 return float128_pack_raw(&p);
4376 bfloat16 bfloat16_default_nan(float_status *status)
4378 FloatParts64 p;
4380 parts_default_nan(&p, status);
4381 p.frac >>= bfloat16_params.frac_shift;
4382 return bfloat16_pack_raw(&p);
4385 /*----------------------------------------------------------------------------
4386 | Returns a quiet NaN from a signalling NaN for the floating point value `a'.
4387 *----------------------------------------------------------------------------*/
4389 float16 float16_silence_nan(float16 a, float_status *status)
4391 FloatParts64 p;
4393 float16_unpack_raw(&p, a);
4394 p.frac <<= float16_params.frac_shift;
4395 parts_silence_nan(&p, status);
4396 p.frac >>= float16_params.frac_shift;
4397 return float16_pack_raw(&p);
4400 float32 float32_silence_nan(float32 a, float_status *status)
4402 FloatParts64 p;
4404 float32_unpack_raw(&p, a);
4405 p.frac <<= float32_params.frac_shift;
4406 parts_silence_nan(&p, status);
4407 p.frac >>= float32_params.frac_shift;
4408 return float32_pack_raw(&p);
4411 float64 float64_silence_nan(float64 a, float_status *status)
4413 FloatParts64 p;
4415 float64_unpack_raw(&p, a);
4416 p.frac <<= float64_params.frac_shift;
4417 parts_silence_nan(&p, status);
4418 p.frac >>= float64_params.frac_shift;
4419 return float64_pack_raw(&p);
4422 bfloat16 bfloat16_silence_nan(bfloat16 a, float_status *status)
4424 FloatParts64 p;
4426 bfloat16_unpack_raw(&p, a);
4427 p.frac <<= bfloat16_params.frac_shift;
4428 parts_silence_nan(&p, status);
4429 p.frac >>= bfloat16_params.frac_shift;
4430 return bfloat16_pack_raw(&p);
4433 float128 float128_silence_nan(float128 a, float_status *status)
4435 FloatParts128 p;
4437 float128_unpack_raw(&p, a);
4438 frac_shl(&p, float128_params.frac_shift);
4439 parts_silence_nan(&p, status);
4440 frac_shr(&p, float128_params.frac_shift);
4441 return float128_pack_raw(&p);
4444 /*----------------------------------------------------------------------------
4445 | If `a' is denormal and we are in flush-to-zero mode then set the
4446 | input-denormal exception and return zero. Otherwise just return the value.
4447 *----------------------------------------------------------------------------*/
4449 static bool parts_squash_denormal(FloatParts64 p, float_status *status)
4451 if (p.exp == 0 && p.frac != 0) {
4452 float_raise(float_flag_input_denormal, status);
4453 return true;
4456 return false;
4459 float16 float16_squash_input_denormal(float16 a, float_status *status)
4461 if (status->flush_inputs_to_zero) {
4462 FloatParts64 p;
4464 float16_unpack_raw(&p, a);
4465 if (parts_squash_denormal(p, status)) {
4466 return float16_set_sign(float16_zero, p.sign);
4469 return a;
4472 float32 float32_squash_input_denormal(float32 a, float_status *status)
4474 if (status->flush_inputs_to_zero) {
4475 FloatParts64 p;
4477 float32_unpack_raw(&p, a);
4478 if (parts_squash_denormal(p, status)) {
4479 return float32_set_sign(float32_zero, p.sign);
4482 return a;
4485 float64 float64_squash_input_denormal(float64 a, float_status *status)
4487 if (status->flush_inputs_to_zero) {
4488 FloatParts64 p;
4490 float64_unpack_raw(&p, a);
4491 if (parts_squash_denormal(p, status)) {
4492 return float64_set_sign(float64_zero, p.sign);
4495 return a;
4498 bfloat16 bfloat16_squash_input_denormal(bfloat16 a, float_status *status)
4500 if (status->flush_inputs_to_zero) {
4501 FloatParts64 p;
4503 bfloat16_unpack_raw(&p, a);
4504 if (parts_squash_denormal(p, status)) {
4505 return bfloat16_set_sign(bfloat16_zero, p.sign);
4508 return a;
4511 /*----------------------------------------------------------------------------
4512 | Normalizes the subnormal extended double-precision floating-point value
4513 | represented by the denormalized significand `aSig'. The normalized exponent
4514 | and significand are stored at the locations pointed to by `zExpPtr' and
4515 | `zSigPtr', respectively.
4516 *----------------------------------------------------------------------------*/
4518 void normalizeFloatx80Subnormal(uint64_t aSig, int32_t *zExpPtr,
4519 uint64_t *zSigPtr)
4521 int8_t shiftCount;
4523 shiftCount = clz64(aSig);
4524 *zSigPtr = aSig<<shiftCount;
4525 *zExpPtr = 1 - shiftCount;
4528 /*----------------------------------------------------------------------------
4529 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4530 | and extended significand formed by the concatenation of `zSig0' and `zSig1',
4531 | and returns the proper extended double-precision floating-point value
4532 | corresponding to the abstract input. Ordinarily, the abstract value is
4533 | rounded and packed into the extended double-precision format, with the
4534 | inexact exception raised if the abstract input cannot be represented
4535 | exactly. However, if the abstract value is too large, the overflow and
4536 | inexact exceptions are raised and an infinity or maximal finite value is
4537 | returned. If the abstract value is too small, the input value is rounded to
4538 | a subnormal number, and the underflow and inexact exceptions are raised if
4539 | the abstract input cannot be represented exactly as a subnormal extended
4540 | double-precision floating-point number.
4541 | If `roundingPrecision' is floatx80_precision_s or floatx80_precision_d,
4542 | the result is rounded to the same number of bits as single or double
4543 | precision, respectively. Otherwise, the result is rounded to the full
4544 | precision of the extended double-precision format.
4545 | The input significand must be normalized or smaller. If the input
4546 | significand is not normalized, `zExp' must be 0; in that case, the result
4547 | returned is a subnormal number, and it must not require rounding. The
4548 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
4549 | Floating-Point Arithmetic.
4550 *----------------------------------------------------------------------------*/
4552 floatx80 roundAndPackFloatx80(FloatX80RoundPrec roundingPrecision, bool zSign,
4553 int32_t zExp, uint64_t zSig0, uint64_t zSig1,
4554 float_status *status)
4556 FloatRoundMode roundingMode;
4557 bool roundNearestEven, increment, isTiny;
4558 int64_t roundIncrement, roundMask, roundBits;
4560 roundingMode = status->float_rounding_mode;
4561 roundNearestEven = ( roundingMode == float_round_nearest_even );
4562 switch (roundingPrecision) {
4563 case floatx80_precision_x:
4564 goto precision80;
4565 case floatx80_precision_d:
4566 roundIncrement = UINT64_C(0x0000000000000400);
4567 roundMask = UINT64_C(0x00000000000007FF);
4568 break;
4569 case floatx80_precision_s:
4570 roundIncrement = UINT64_C(0x0000008000000000);
4571 roundMask = UINT64_C(0x000000FFFFFFFFFF);
4572 break;
4573 default:
4574 g_assert_not_reached();
4576 zSig0 |= ( zSig1 != 0 );
4577 switch (roundingMode) {
4578 case float_round_nearest_even:
4579 case float_round_ties_away:
4580 break;
4581 case float_round_to_zero:
4582 roundIncrement = 0;
4583 break;
4584 case float_round_up:
4585 roundIncrement = zSign ? 0 : roundMask;
4586 break;
4587 case float_round_down:
4588 roundIncrement = zSign ? roundMask : 0;
4589 break;
4590 default:
4591 abort();
4593 roundBits = zSig0 & roundMask;
4594 if ( 0x7FFD <= (uint32_t) ( zExp - 1 ) ) {
4595 if ( ( 0x7FFE < zExp )
4596 || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
4598 goto overflow;
4600 if ( zExp <= 0 ) {
4601 if (status->flush_to_zero) {
4602 float_raise(float_flag_output_denormal, status);
4603 return packFloatx80(zSign, 0, 0);
4605 isTiny = status->tininess_before_rounding
4606 || (zExp < 0 )
4607 || (zSig0 <= zSig0 + roundIncrement);
4608 shift64RightJamming( zSig0, 1 - zExp, &zSig0 );
4609 zExp = 0;
4610 roundBits = zSig0 & roundMask;
4611 if (isTiny && roundBits) {
4612 float_raise(float_flag_underflow, status);
4614 if (roundBits) {
4615 float_raise(float_flag_inexact, status);
4617 zSig0 += roundIncrement;
4618 if ( (int64_t) zSig0 < 0 ) zExp = 1;
4619 roundIncrement = roundMask + 1;
4620 if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
4621 roundMask |= roundIncrement;
4623 zSig0 &= ~ roundMask;
4624 return packFloatx80( zSign, zExp, zSig0 );
4627 if (roundBits) {
4628 float_raise(float_flag_inexact, status);
4630 zSig0 += roundIncrement;
4631 if ( zSig0 < roundIncrement ) {
4632 ++zExp;
4633 zSig0 = UINT64_C(0x8000000000000000);
4635 roundIncrement = roundMask + 1;
4636 if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
4637 roundMask |= roundIncrement;
4639 zSig0 &= ~ roundMask;
4640 if ( zSig0 == 0 ) zExp = 0;
4641 return packFloatx80( zSign, zExp, zSig0 );
4642 precision80:
4643 switch (roundingMode) {
4644 case float_round_nearest_even:
4645 case float_round_ties_away:
4646 increment = ((int64_t)zSig1 < 0);
4647 break;
4648 case float_round_to_zero:
4649 increment = 0;
4650 break;
4651 case float_round_up:
4652 increment = !zSign && zSig1;
4653 break;
4654 case float_round_down:
4655 increment = zSign && zSig1;
4656 break;
4657 default:
4658 abort();
4660 if ( 0x7FFD <= (uint32_t) ( zExp - 1 ) ) {
4661 if ( ( 0x7FFE < zExp )
4662 || ( ( zExp == 0x7FFE )
4663 && ( zSig0 == UINT64_C(0xFFFFFFFFFFFFFFFF) )
4664 && increment
4667 roundMask = 0;
4668 overflow:
4669 float_raise(float_flag_overflow | float_flag_inexact, status);
4670 if ( ( roundingMode == float_round_to_zero )
4671 || ( zSign && ( roundingMode == float_round_up ) )
4672 || ( ! zSign && ( roundingMode == float_round_down ) )
4674 return packFloatx80( zSign, 0x7FFE, ~ roundMask );
4676 return packFloatx80(zSign,
4677 floatx80_infinity_high,
4678 floatx80_infinity_low);
4680 if ( zExp <= 0 ) {
4681 isTiny = status->tininess_before_rounding
4682 || (zExp < 0)
4683 || !increment
4684 || (zSig0 < UINT64_C(0xFFFFFFFFFFFFFFFF));
4685 shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );
4686 zExp = 0;
4687 if (isTiny && zSig1) {
4688 float_raise(float_flag_underflow, status);
4690 if (zSig1) {
4691 float_raise(float_flag_inexact, status);
4693 switch (roundingMode) {
4694 case float_round_nearest_even:
4695 case float_round_ties_away:
4696 increment = ((int64_t)zSig1 < 0);
4697 break;
4698 case float_round_to_zero:
4699 increment = 0;
4700 break;
4701 case float_round_up:
4702 increment = !zSign && zSig1;
4703 break;
4704 case float_round_down:
4705 increment = zSign && zSig1;
4706 break;
4707 default:
4708 abort();
4710 if ( increment ) {
4711 ++zSig0;
4712 if (!(zSig1 << 1) && roundNearestEven) {
4713 zSig0 &= ~1;
4715 if ( (int64_t) zSig0 < 0 ) zExp = 1;
4717 return packFloatx80( zSign, zExp, zSig0 );
4720 if (zSig1) {
4721 float_raise(float_flag_inexact, status);
4723 if ( increment ) {
4724 ++zSig0;
4725 if ( zSig0 == 0 ) {
4726 ++zExp;
4727 zSig0 = UINT64_C(0x8000000000000000);
4729 else {
4730 if (!(zSig1 << 1) && roundNearestEven) {
4731 zSig0 &= ~1;
4735 else {
4736 if ( zSig0 == 0 ) zExp = 0;
4738 return packFloatx80( zSign, zExp, zSig0 );
4742 /*----------------------------------------------------------------------------
4743 | Takes an abstract floating-point value having sign `zSign', exponent
4744 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
4745 | and returns the proper extended double-precision floating-point value
4746 | corresponding to the abstract input. This routine is just like
4747 | `roundAndPackFloatx80' except that the input significand does not have to be
4748 | normalized.
4749 *----------------------------------------------------------------------------*/
4751 floatx80 normalizeRoundAndPackFloatx80(FloatX80RoundPrec roundingPrecision,
4752 bool zSign, int32_t zExp,
4753 uint64_t zSig0, uint64_t zSig1,
4754 float_status *status)
4756 int8_t shiftCount;
4758 if ( zSig0 == 0 ) {
4759 zSig0 = zSig1;
4760 zSig1 = 0;
4761 zExp -= 64;
4763 shiftCount = clz64(zSig0);
4764 shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
4765 zExp -= shiftCount;
4766 return roundAndPackFloatx80(roundingPrecision, zSign, zExp,
4767 zSig0, zSig1, status);
4771 /*----------------------------------------------------------------------------
4772 | Returns the binary exponential of the single-precision floating-point value
4773 | `a'. The operation is performed according to the IEC/IEEE Standard for
4774 | Binary Floating-Point Arithmetic.
4776 | Uses the following identities:
4778 | 1. -------------------------------------------------------------------------
4779 | x x*ln(2)
4780 | 2 = e
4782 | 2. -------------------------------------------------------------------------
4783 | 2 3 4 5 n
4784 | x x x x x x x
4785 | e = 1 + --- + --- + --- + --- + --- + ... + --- + ...
4786 | 1! 2! 3! 4! 5! n!
4787 *----------------------------------------------------------------------------*/
4789 static const float64 float32_exp2_coefficients[15] =
4791 const_float64( 0x3ff0000000000000ll ), /* 1 */
4792 const_float64( 0x3fe0000000000000ll ), /* 2 */
4793 const_float64( 0x3fc5555555555555ll ), /* 3 */
4794 const_float64( 0x3fa5555555555555ll ), /* 4 */
4795 const_float64( 0x3f81111111111111ll ), /* 5 */
4796 const_float64( 0x3f56c16c16c16c17ll ), /* 6 */
4797 const_float64( 0x3f2a01a01a01a01all ), /* 7 */
4798 const_float64( 0x3efa01a01a01a01all ), /* 8 */
4799 const_float64( 0x3ec71de3a556c734ll ), /* 9 */
4800 const_float64( 0x3e927e4fb7789f5cll ), /* 10 */
4801 const_float64( 0x3e5ae64567f544e4ll ), /* 11 */
4802 const_float64( 0x3e21eed8eff8d898ll ), /* 12 */
4803 const_float64( 0x3de6124613a86d09ll ), /* 13 */
4804 const_float64( 0x3da93974a8c07c9dll ), /* 14 */
4805 const_float64( 0x3d6ae7f3e733b81fll ), /* 15 */
4808 float32 float32_exp2(float32 a, float_status *status)
4810 FloatParts64 xp, xnp, tp, rp;
4811 int i;
4813 float32_unpack_canonical(&xp, a, status);
4814 if (unlikely(xp.cls != float_class_normal)) {
4815 switch (xp.cls) {
4816 case float_class_snan:
4817 case float_class_qnan:
4818 parts_return_nan(&xp, status);
4819 return float32_round_pack_canonical(&xp, status);
4820 case float_class_inf:
4821 return xp.sign ? float32_zero : a;
4822 case float_class_zero:
4823 return float32_one;
4824 default:
4825 break;
4827 g_assert_not_reached();
4830 float_raise(float_flag_inexact, status);
4832 float64_unpack_canonical(&tp, float64_ln2, status);
4833 xp = *parts_mul(&xp, &tp, status);
4834 xnp = xp;
4836 float64_unpack_canonical(&rp, float64_one, status);
4837 for (i = 0 ; i < 15 ; i++) {
4838 float64_unpack_canonical(&tp, float32_exp2_coefficients[i], status);
4839 rp = *parts_muladd(&tp, &xp, &rp, 0, status);
4840 xnp = *parts_mul(&xnp, &xp, status);
4843 return float32_round_pack_canonical(&rp, status);
4846 /*----------------------------------------------------------------------------
4847 | Rounds the extended double-precision floating-point value `a'
4848 | to the precision provided by floatx80_rounding_precision and returns the
4849 | result as an extended double-precision floating-point value.
4850 | The operation is performed according to the IEC/IEEE Standard for Binary
4851 | Floating-Point Arithmetic.
4852 *----------------------------------------------------------------------------*/
4854 floatx80 floatx80_round(floatx80 a, float_status *status)
4856 FloatParts128 p;
4858 if (!floatx80_unpack_canonical(&p, a, status)) {
4859 return floatx80_default_nan(status);
4861 return floatx80_round_pack_canonical(&p, status);
4864 static void __attribute__((constructor)) softfloat_init(void)
4866 union_float64 ua, ub, uc, ur;
4868 if (QEMU_NO_HARDFLOAT) {
4869 return;
4872 * Test that the host's FMA is not obviously broken. For example,
4873 * glibc < 2.23 can perform an incorrect FMA on certain hosts; see
4874 * https://sourceware.org/bugzilla/show_bug.cgi?id=13304
4876 ua.s = 0x0020000000000001ULL;
4877 ub.s = 0x3ca0000000000000ULL;
4878 uc.s = 0x0020000000000000ULL;
4879 ur.h = fma(ua.h, ub.h, uc.h);
4880 if (ur.s != 0x0020000000000001ULL) {
4881 force_soft_fma = true;