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
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 ===============================================================================
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"
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
96 *----------------------------------------------------------------------------*/
97 #include "fpu/softfloat-macros.h"
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)) { \
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)) { \
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)) { \
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
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
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
215 # define QEMU_HARDFLOAT_USE_ISINF 0
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
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 \
228 # define QEMU_NO_HARDFLOAT 1
229 # define QEMU_SOFTFLOAT_ATTR QEMU_FLATTEN
231 # define QEMU_NO_HARDFLOAT 0
232 # define QEMU_SOFTFLOAT_ATTR QEMU_FLATTEN __attribute__((noinline))
235 static inline bool can_use_fpu(const float_status
*s
)
237 if (QEMU_NO_HARDFLOAT
) {
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.
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 */
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
);
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
) {
331 return float32_is_infinity(a
.s
);
334 static inline bool f64_is_inf(union_float64 a
)
336 if (QEMU_HARDFLOAT_USE_ISINF
) {
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
;
352 if (unlikely(!can_use_fpu(s
))) {
356 float32_input_flush2(&ua
.s
, &ub
.s
, s
);
357 if (unlikely(!pre(ua
, ub
))) {
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
)) {
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
;
383 if (unlikely(!can_use_fpu(s
))) {
387 float64_input_flush2(&ua
.s
, &ub
.s
, s
);
388 if (unlikely(!pre(ua
, ub
))) {
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
)) {
401 return soft(ua
.s
, ub
.s
, s
);
404 /*----------------------------------------------------------------------------
405 | Returns the fraction bits of the single-precision floating-point value `a'.
406 *----------------------------------------------------------------------------*/
408 static inline uint32_t extractFloat32Frac(float32 a
)
410 return float32_val(a
) & 0x007FFFFF;
413 /*----------------------------------------------------------------------------
414 | Returns the exponent bits of the single-precision floating-point value `a'.
415 *----------------------------------------------------------------------------*/
417 static inline int extractFloat32Exp(float32 a
)
419 return (float32_val(a
) >> 23) & 0xFF;
422 /*----------------------------------------------------------------------------
423 | Returns the sign bit of the single-precision floating-point value `a'.
424 *----------------------------------------------------------------------------*/
426 static inline bool extractFloat32Sign(float32 a
)
428 return float32_val(a
) >> 31;
431 /*----------------------------------------------------------------------------
432 | Returns the fraction bits of the double-precision floating-point value `a'.
433 *----------------------------------------------------------------------------*/
435 static inline uint64_t extractFloat64Frac(float64 a
)
437 return float64_val(a
) & UINT64_C(0x000FFFFFFFFFFFFF);
440 /*----------------------------------------------------------------------------
441 | Returns the exponent bits of the double-precision floating-point value `a'.
442 *----------------------------------------------------------------------------*/
444 static inline int extractFloat64Exp(float64 a
)
446 return (float64_val(a
) >> 52) & 0x7FF;
449 /*----------------------------------------------------------------------------
450 | Returns the sign bit of the double-precision floating-point value `a'.
451 *----------------------------------------------------------------------------*/
453 static inline bool extractFloat64Sign(float64 a
)
455 return float64_val(a
) >> 63;
459 * Classify a floating point number. Everything above float_class_qnan
460 * is a NaN so cls >= float_class_qnan is any NaN.
463 typedef enum __attribute__ ((__packed__
)) {
464 float_class_unclassified
,
468 float_class_qnan
, /* all NaNs from here */
472 #define float_cmask(bit) (1u << (bit))
475 float_cmask_zero
= float_cmask(float_class_zero
),
476 float_cmask_normal
= float_cmask(float_class_normal
),
477 float_cmask_inf
= float_cmask(float_class_inf
),
478 float_cmask_qnan
= float_cmask(float_class_qnan
),
479 float_cmask_snan
= float_cmask(float_class_snan
),
481 float_cmask_infzero
= float_cmask_zero
| float_cmask_inf
,
482 float_cmask_anynan
= float_cmask_qnan
| float_cmask_snan
,
485 /* Flags for parts_minmax. */
487 /* Set for minimum; clear for maximum. */
489 /* Set for the IEEE 754-2008 minNum() and maxNum() operations. */
491 /* Set for the IEEE 754-2008 minNumMag() and minNumMag() operations. */
495 /* Simple helpers for checking if, or what kind of, NaN we have */
496 static inline __attribute__((unused
)) bool is_nan(FloatClass c
)
498 return unlikely(c
>= float_class_qnan
);
501 static inline __attribute__((unused
)) bool is_snan(FloatClass c
)
503 return c
== float_class_snan
;
506 static inline __attribute__((unused
)) bool is_qnan(FloatClass c
)
508 return c
== float_class_qnan
;
512 * Structure holding all of the decomposed parts of a float.
513 * The exponent is unbiased and the fraction is normalized.
515 * The fraction words are stored in big-endian word ordering,
516 * so that truncation from a larger format to a smaller format
517 * can be done simply by ignoring subsequent elements.
525 /* Routines that know the structure may reference the singular name. */
528 * Routines expanded with multiple structures reference "hi" and "lo"
529 * depending on the operation. In FloatParts64, "hi" and "lo" are
530 * both the same word and aliased here.
550 uint64_t frac_hm
; /* high-middle */
551 uint64_t frac_lm
; /* low-middle */
555 /* These apply to the most significant word of each FloatPartsN. */
556 #define DECOMPOSED_BINARY_POINT 63
557 #define DECOMPOSED_IMPLICIT_BIT (1ull << DECOMPOSED_BINARY_POINT)
559 /* Structure holding all of the relevant parameters for a format.
560 * exp_size: the size of the exponent field
561 * exp_bias: the offset applied to the exponent field
562 * exp_max: the maximum normalised exponent
563 * frac_size: the size of the fraction field
564 * frac_shift: shift to normalise the fraction with DECOMPOSED_BINARY_POINT
565 * The following are computed based the size of fraction
566 * round_mask: bits below lsb which must be rounded
567 * The following optional modifiers are available:
568 * arm_althp: handle ARM Alternative Half Precision
580 /* Expand fields based on the size of exponent and fraction */
581 #define FLOAT_PARAMS_(E) \
583 .exp_bias = ((1 << E) - 1) >> 1, \
584 .exp_max = (1 << E) - 1
586 #define FLOAT_PARAMS(E, F) \
589 .frac_shift = (-F - 1) & 63, \
590 .round_mask = (1ull << ((-F - 1) & 63)) - 1
592 static const FloatFmt float16_params
= {
596 static const FloatFmt float16_params_ahp
= {
601 static const FloatFmt bfloat16_params
= {
605 static const FloatFmt float32_params
= {
609 static const FloatFmt float64_params
= {
613 static const FloatFmt float128_params
= {
614 FLOAT_PARAMS(15, 112)
617 #define FLOATX80_PARAMS(R) \
619 .frac_size = R == 64 ? 63 : R, \
621 .round_mask = R == 64 ? -1 : (1ull << ((-R - 1) & 63)) - 1
623 static const FloatFmt floatx80_params
[3] = {
624 [floatx80_precision_s
] = { FLOATX80_PARAMS(23) },
625 [floatx80_precision_d
] = { FLOATX80_PARAMS(52) },
626 [floatx80_precision_x
] = { FLOATX80_PARAMS(64) },
629 /* Unpack a float to parts, but do not canonicalize. */
630 static void unpack_raw64(FloatParts64
*r
, const FloatFmt
*fmt
, uint64_t raw
)
632 const int f_size
= fmt
->frac_size
;
633 const int e_size
= fmt
->exp_size
;
635 *r
= (FloatParts64
) {
636 .cls
= float_class_unclassified
,
637 .sign
= extract64(raw
, f_size
+ e_size
, 1),
638 .exp
= extract64(raw
, f_size
, e_size
),
639 .frac
= extract64(raw
, 0, f_size
)
643 static inline void float16_unpack_raw(FloatParts64
*p
, float16 f
)
645 unpack_raw64(p
, &float16_params
, f
);
648 static inline void bfloat16_unpack_raw(FloatParts64
*p
, bfloat16 f
)
650 unpack_raw64(p
, &bfloat16_params
, f
);
653 static inline void float32_unpack_raw(FloatParts64
*p
, float32 f
)
655 unpack_raw64(p
, &float32_params
, f
);
658 static inline void float64_unpack_raw(FloatParts64
*p
, float64 f
)
660 unpack_raw64(p
, &float64_params
, f
);
663 static void floatx80_unpack_raw(FloatParts128
*p
, floatx80 f
)
665 *p
= (FloatParts128
) {
666 .cls
= float_class_unclassified
,
667 .sign
= extract32(f
.high
, 15, 1),
668 .exp
= extract32(f
.high
, 0, 15),
673 static void float128_unpack_raw(FloatParts128
*p
, float128 f
)
675 const int f_size
= float128_params
.frac_size
- 64;
676 const int e_size
= float128_params
.exp_size
;
678 *p
= (FloatParts128
) {
679 .cls
= float_class_unclassified
,
680 .sign
= extract64(f
.high
, f_size
+ e_size
, 1),
681 .exp
= extract64(f
.high
, f_size
, e_size
),
682 .frac_hi
= extract64(f
.high
, 0, f_size
),
687 /* Pack a float from parts, but do not canonicalize. */
688 static uint64_t pack_raw64(const FloatParts64
*p
, const FloatFmt
*fmt
)
690 const int f_size
= fmt
->frac_size
;
691 const int e_size
= fmt
->exp_size
;
694 ret
= (uint64_t)p
->sign
<< (f_size
+ e_size
);
695 ret
= deposit64(ret
, f_size
, e_size
, p
->exp
);
696 ret
= deposit64(ret
, 0, f_size
, p
->frac
);
700 static inline float16
float16_pack_raw(const FloatParts64
*p
)
702 return make_float16(pack_raw64(p
, &float16_params
));
705 static inline bfloat16
bfloat16_pack_raw(const FloatParts64
*p
)
707 return pack_raw64(p
, &bfloat16_params
);
710 static inline float32
float32_pack_raw(const FloatParts64
*p
)
712 return make_float32(pack_raw64(p
, &float32_params
));
715 static inline float64
float64_pack_raw(const FloatParts64
*p
)
717 return make_float64(pack_raw64(p
, &float64_params
));
720 static float128
float128_pack_raw(const FloatParts128
*p
)
722 const int f_size
= float128_params
.frac_size
- 64;
723 const int e_size
= float128_params
.exp_size
;
726 hi
= (uint64_t)p
->sign
<< (f_size
+ e_size
);
727 hi
= deposit64(hi
, f_size
, e_size
, p
->exp
);
728 hi
= deposit64(hi
, 0, f_size
, p
->frac_hi
);
729 return make_float128(hi
, p
->frac_lo
);
732 /*----------------------------------------------------------------------------
733 | Functions and definitions to determine: (1) whether tininess for underflow
734 | is detected before or after rounding by default, (2) what (if anything)
735 | happens when exceptions are raised, (3) how signaling NaNs are distinguished
736 | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
737 | are propagated from function inputs to output. These details are target-
739 *----------------------------------------------------------------------------*/
740 #include "softfloat-specialize.c.inc"
742 #define PARTS_GENERIC_64_128(NAME, P) \
743 QEMU_GENERIC(P, (FloatParts128 *, parts128_##NAME), parts64_##NAME)
745 #define PARTS_GENERIC_64_128_256(NAME, P) \
746 QEMU_GENERIC(P, (FloatParts256 *, parts256_##NAME), \
747 (FloatParts128 *, parts128_##NAME), parts64_##NAME)
749 #define parts_default_nan(P, S) PARTS_GENERIC_64_128(default_nan, P)(P, S)
750 #define parts_silence_nan(P, S) PARTS_GENERIC_64_128(silence_nan, P)(P, S)
752 static void parts64_return_nan(FloatParts64
*a
, float_status
*s
);
753 static void parts128_return_nan(FloatParts128
*a
, float_status
*s
);
755 #define parts_return_nan(P, S) PARTS_GENERIC_64_128(return_nan, P)(P, S)
757 static FloatParts64
*parts64_pick_nan(FloatParts64
*a
, FloatParts64
*b
,
759 static FloatParts128
*parts128_pick_nan(FloatParts128
*a
, FloatParts128
*b
,
762 #define parts_pick_nan(A, B, S) PARTS_GENERIC_64_128(pick_nan, A)(A, B, S)
764 static FloatParts64
*parts64_pick_nan_muladd(FloatParts64
*a
, FloatParts64
*b
,
765 FloatParts64
*c
, float_status
*s
,
766 int ab_mask
, int abc_mask
);
767 static FloatParts128
*parts128_pick_nan_muladd(FloatParts128
*a
,
771 int ab_mask
, int abc_mask
);
773 #define parts_pick_nan_muladd(A, B, C, S, ABM, ABCM) \
774 PARTS_GENERIC_64_128(pick_nan_muladd, A)(A, B, C, S, ABM, ABCM)
776 static void parts64_canonicalize(FloatParts64
*p
, float_status
*status
,
777 const FloatFmt
*fmt
);
778 static void parts128_canonicalize(FloatParts128
*p
, float_status
*status
,
779 const FloatFmt
*fmt
);
781 #define parts_canonicalize(A, S, F) \
782 PARTS_GENERIC_64_128(canonicalize, A)(A, S, F)
784 static void parts64_uncanon_normal(FloatParts64
*p
, float_status
*status
,
785 const FloatFmt
*fmt
);
786 static void parts128_uncanon_normal(FloatParts128
*p
, float_status
*status
,
787 const FloatFmt
*fmt
);
789 #define parts_uncanon_normal(A, S, F) \
790 PARTS_GENERIC_64_128(uncanon_normal, A)(A, S, F)
792 static void parts64_uncanon(FloatParts64
*p
, float_status
*status
,
793 const FloatFmt
*fmt
);
794 static void parts128_uncanon(FloatParts128
*p
, float_status
*status
,
795 const FloatFmt
*fmt
);
797 #define parts_uncanon(A, S, F) \
798 PARTS_GENERIC_64_128(uncanon, A)(A, S, F)
800 static void parts64_add_normal(FloatParts64
*a
, FloatParts64
*b
);
801 static void parts128_add_normal(FloatParts128
*a
, FloatParts128
*b
);
802 static void parts256_add_normal(FloatParts256
*a
, FloatParts256
*b
);
804 #define parts_add_normal(A, B) \
805 PARTS_GENERIC_64_128_256(add_normal, A)(A, B)
807 static bool parts64_sub_normal(FloatParts64
*a
, FloatParts64
*b
);
808 static bool parts128_sub_normal(FloatParts128
*a
, FloatParts128
*b
);
809 static bool parts256_sub_normal(FloatParts256
*a
, FloatParts256
*b
);
811 #define parts_sub_normal(A, B) \
812 PARTS_GENERIC_64_128_256(sub_normal, A)(A, B)
814 static FloatParts64
*parts64_addsub(FloatParts64
*a
, FloatParts64
*b
,
815 float_status
*s
, bool subtract
);
816 static FloatParts128
*parts128_addsub(FloatParts128
*a
, FloatParts128
*b
,
817 float_status
*s
, bool subtract
);
819 #define parts_addsub(A, B, S, Z) \
820 PARTS_GENERIC_64_128(addsub, A)(A, B, S, Z)
822 static FloatParts64
*parts64_mul(FloatParts64
*a
, FloatParts64
*b
,
824 static FloatParts128
*parts128_mul(FloatParts128
*a
, FloatParts128
*b
,
827 #define parts_mul(A, B, S) \
828 PARTS_GENERIC_64_128(mul, A)(A, B, S)
830 static FloatParts64
*parts64_muladd(FloatParts64
*a
, FloatParts64
*b
,
831 FloatParts64
*c
, int flags
,
833 static FloatParts128
*parts128_muladd(FloatParts128
*a
, FloatParts128
*b
,
834 FloatParts128
*c
, int flags
,
837 #define parts_muladd(A, B, C, Z, S) \
838 PARTS_GENERIC_64_128(muladd, A)(A, B, C, Z, S)
840 static FloatParts64
*parts64_div(FloatParts64
*a
, FloatParts64
*b
,
842 static FloatParts128
*parts128_div(FloatParts128
*a
, FloatParts128
*b
,
845 #define parts_div(A, B, S) \
846 PARTS_GENERIC_64_128(div, A)(A, B, S)
848 static void parts64_sqrt(FloatParts64
*a
, float_status
*s
, const FloatFmt
*f
);
849 static void parts128_sqrt(FloatParts128
*a
, float_status
*s
, const FloatFmt
*f
);
851 #define parts_sqrt(A, S, F) \
852 PARTS_GENERIC_64_128(sqrt, A)(A, S, F)
854 static bool parts64_round_to_int_normal(FloatParts64
*a
, FloatRoundMode rm
,
855 int scale
, int frac_size
);
856 static bool parts128_round_to_int_normal(FloatParts128
*a
, FloatRoundMode r
,
857 int scale
, int frac_size
);
859 #define parts_round_to_int_normal(A, R, C, F) \
860 PARTS_GENERIC_64_128(round_to_int_normal, A)(A, R, C, F)
862 static void parts64_round_to_int(FloatParts64
*a
, FloatRoundMode rm
,
863 int scale
, float_status
*s
,
864 const FloatFmt
*fmt
);
865 static void parts128_round_to_int(FloatParts128
*a
, FloatRoundMode r
,
866 int scale
, float_status
*s
,
867 const FloatFmt
*fmt
);
869 #define parts_round_to_int(A, R, C, S, F) \
870 PARTS_GENERIC_64_128(round_to_int, A)(A, R, C, S, F)
872 static int64_t parts64_float_to_sint(FloatParts64
*p
, FloatRoundMode rmode
,
873 int scale
, int64_t min
, int64_t max
,
875 static int64_t parts128_float_to_sint(FloatParts128
*p
, FloatRoundMode rmode
,
876 int scale
, int64_t min
, int64_t max
,
879 #define parts_float_to_sint(P, R, Z, MN, MX, S) \
880 PARTS_GENERIC_64_128(float_to_sint, P)(P, R, Z, MN, MX, S)
882 static uint64_t parts64_float_to_uint(FloatParts64
*p
, FloatRoundMode rmode
,
883 int scale
, uint64_t max
,
885 static uint64_t parts128_float_to_uint(FloatParts128
*p
, FloatRoundMode rmode
,
886 int scale
, uint64_t max
,
889 #define parts_float_to_uint(P, R, Z, M, S) \
890 PARTS_GENERIC_64_128(float_to_uint, P)(P, R, Z, M, S)
892 static void parts64_sint_to_float(FloatParts64
*p
, int64_t a
,
893 int scale
, float_status
*s
);
894 static void parts128_sint_to_float(FloatParts128
*p
, int64_t a
,
895 int scale
, float_status
*s
);
897 #define parts_sint_to_float(P, I, Z, S) \
898 PARTS_GENERIC_64_128(sint_to_float, P)(P, I, Z, S)
900 static void parts64_uint_to_float(FloatParts64
*p
, uint64_t a
,
901 int scale
, float_status
*s
);
902 static void parts128_uint_to_float(FloatParts128
*p
, uint64_t a
,
903 int scale
, float_status
*s
);
905 #define parts_uint_to_float(P, I, Z, S) \
906 PARTS_GENERIC_64_128(uint_to_float, P)(P, I, Z, S)
908 static FloatParts64
*parts64_minmax(FloatParts64
*a
, FloatParts64
*b
,
909 float_status
*s
, int flags
);
910 static FloatParts128
*parts128_minmax(FloatParts128
*a
, FloatParts128
*b
,
911 float_status
*s
, int flags
);
913 #define parts_minmax(A, B, S, F) \
914 PARTS_GENERIC_64_128(minmax, A)(A, B, S, F)
916 static int parts64_compare(FloatParts64
*a
, FloatParts64
*b
,
917 float_status
*s
, bool q
);
918 static int parts128_compare(FloatParts128
*a
, FloatParts128
*b
,
919 float_status
*s
, bool q
);
921 #define parts_compare(A, B, S, Q) \
922 PARTS_GENERIC_64_128(compare, A)(A, B, S, Q)
924 static void parts64_scalbn(FloatParts64
*a
, int n
, float_status
*s
);
925 static void parts128_scalbn(FloatParts128
*a
, int n
, float_status
*s
);
927 #define parts_scalbn(A, N, S) \
928 PARTS_GENERIC_64_128(scalbn, A)(A, N, S)
931 * Helper functions for softfloat-parts.c.inc, per-size operations.
934 #define FRAC_GENERIC_64_128(NAME, P) \
935 QEMU_GENERIC(P, (FloatParts128 *, frac128_##NAME), frac64_##NAME)
937 #define FRAC_GENERIC_64_128_256(NAME, P) \
938 QEMU_GENERIC(P, (FloatParts256 *, frac256_##NAME), \
939 (FloatParts128 *, frac128_##NAME), frac64_##NAME)
941 static bool frac64_add(FloatParts64
*r
, FloatParts64
*a
, FloatParts64
*b
)
943 return uadd64_overflow(a
->frac
, b
->frac
, &r
->frac
);
946 static bool frac128_add(FloatParts128
*r
, FloatParts128
*a
, FloatParts128
*b
)
949 r
->frac_lo
= uadd64_carry(a
->frac_lo
, b
->frac_lo
, &c
);
950 r
->frac_hi
= uadd64_carry(a
->frac_hi
, b
->frac_hi
, &c
);
954 static bool frac256_add(FloatParts256
*r
, FloatParts256
*a
, FloatParts256
*b
)
957 r
->frac_lo
= uadd64_carry(a
->frac_lo
, b
->frac_lo
, &c
);
958 r
->frac_lm
= uadd64_carry(a
->frac_lm
, b
->frac_lm
, &c
);
959 r
->frac_hm
= uadd64_carry(a
->frac_hm
, b
->frac_hm
, &c
);
960 r
->frac_hi
= uadd64_carry(a
->frac_hi
, b
->frac_hi
, &c
);
964 #define frac_add(R, A, B) FRAC_GENERIC_64_128_256(add, R)(R, A, B)
966 static bool frac64_addi(FloatParts64
*r
, FloatParts64
*a
, uint64_t c
)
968 return uadd64_overflow(a
->frac
, c
, &r
->frac
);
971 static bool frac128_addi(FloatParts128
*r
, FloatParts128
*a
, uint64_t c
)
973 c
= uadd64_overflow(a
->frac_lo
, c
, &r
->frac_lo
);
974 return uadd64_overflow(a
->frac_hi
, c
, &r
->frac_hi
);
977 #define frac_addi(R, A, C) FRAC_GENERIC_64_128(addi, R)(R, A, C)
979 static void frac64_allones(FloatParts64
*a
)
984 static void frac128_allones(FloatParts128
*a
)
986 a
->frac_hi
= a
->frac_lo
= -1;
989 #define frac_allones(A) FRAC_GENERIC_64_128(allones, A)(A)
991 static int frac64_cmp(FloatParts64
*a
, FloatParts64
*b
)
993 return a
->frac
== b
->frac
? 0 : a
->frac
< b
->frac
? -1 : 1;
996 static int frac128_cmp(FloatParts128
*a
, FloatParts128
*b
)
998 uint64_t ta
= a
->frac_hi
, tb
= b
->frac_hi
;
1000 ta
= a
->frac_lo
, tb
= b
->frac_lo
;
1005 return ta
< tb
? -1 : 1;
1008 #define frac_cmp(A, B) FRAC_GENERIC_64_128(cmp, A)(A, B)
1010 static void frac64_clear(FloatParts64
*a
)
1015 static void frac128_clear(FloatParts128
*a
)
1017 a
->frac_hi
= a
->frac_lo
= 0;
1020 #define frac_clear(A) FRAC_GENERIC_64_128(clear, A)(A)
1022 static bool frac64_div(FloatParts64
*a
, FloatParts64
*b
)
1024 uint64_t n1
, n0
, r
, q
;
1028 * We want a 2*N / N-bit division to produce exactly an N-bit
1029 * result, so that we do not lose any precision and so that we
1030 * do not have to renormalize afterward. If A.frac < B.frac,
1031 * then division would produce an (N-1)-bit result; shift A left
1032 * by one to produce the an N-bit result, and return true to
1033 * decrement the exponent to match.
1035 * The udiv_qrnnd algorithm that we're using requires normalization,
1036 * i.e. the msb of the denominator must be set, which is already true.
1038 ret
= a
->frac
< b
->frac
;
1046 q
= udiv_qrnnd(&r
, n0
, n1
, b
->frac
);
1048 /* Set lsb if there is a remainder, to set inexact. */
1049 a
->frac
= q
| (r
!= 0);
1054 static bool frac128_div(FloatParts128
*a
, FloatParts128
*b
)
1056 uint64_t q0
, q1
, a0
, a1
, b0
, b1
;
1057 uint64_t r0
, r1
, r2
, r3
, t0
, t1
, t2
, t3
;
1060 a0
= a
->frac_hi
, a1
= a
->frac_lo
;
1061 b0
= b
->frac_hi
, b1
= b
->frac_lo
;
1063 ret
= lt128(a0
, a1
, b0
, b1
);
1065 a1
= shr_double(a0
, a1
, 1);
1069 /* Use 128/64 -> 64 division as estimate for 192/128 -> 128 division. */
1070 q0
= estimateDiv128To64(a0
, a1
, b0
);
1073 * Estimate is high because B1 was not included (unless B1 == 0).
1074 * Reduce quotient and increase remainder until remainder is non-negative.
1075 * This loop will execute 0 to 2 times.
1077 mul128By64To192(b0
, b1
, q0
, &t0
, &t1
, &t2
);
1078 sub192(a0
, a1
, 0, t0
, t1
, t2
, &r0
, &r1
, &r2
);
1081 add192(r0
, r1
, r2
, 0, b0
, b1
, &r0
, &r1
, &r2
);
1084 /* Repeat using the remainder, producing a second word of quotient. */
1085 q1
= estimateDiv128To64(r1
, r2
, b0
);
1086 mul128By64To192(b0
, b1
, q1
, &t1
, &t2
, &t3
);
1087 sub192(r1
, r2
, 0, t1
, t2
, t3
, &r1
, &r2
, &r3
);
1090 add192(r1
, r2
, r3
, 0, b0
, b1
, &r1
, &r2
, &r3
);
1093 /* Any remainder indicates inexact; set sticky bit. */
1094 q1
|= (r2
| r3
) != 0;
1101 #define frac_div(A, B) FRAC_GENERIC_64_128(div, A)(A, B)
1103 static bool frac64_eqz(FloatParts64
*a
)
1105 return a
->frac
== 0;
1108 static bool frac128_eqz(FloatParts128
*a
)
1110 return (a
->frac_hi
| a
->frac_lo
) == 0;
1113 #define frac_eqz(A) FRAC_GENERIC_64_128(eqz, A)(A)
1115 static void frac64_mulw(FloatParts128
*r
, FloatParts64
*a
, FloatParts64
*b
)
1117 mulu64(&r
->frac_lo
, &r
->frac_hi
, a
->frac
, b
->frac
);
1120 static void frac128_mulw(FloatParts256
*r
, FloatParts128
*a
, FloatParts128
*b
)
1122 mul128To256(a
->frac_hi
, a
->frac_lo
, b
->frac_hi
, b
->frac_lo
,
1123 &r
->frac_hi
, &r
->frac_hm
, &r
->frac_lm
, &r
->frac_lo
);
1126 #define frac_mulw(R, A, B) FRAC_GENERIC_64_128(mulw, A)(R, A, B)
1128 static void frac64_neg(FloatParts64
*a
)
1133 static void frac128_neg(FloatParts128
*a
)
1136 a
->frac_lo
= usub64_borrow(0, a
->frac_lo
, &c
);
1137 a
->frac_hi
= usub64_borrow(0, a
->frac_hi
, &c
);
1140 static void frac256_neg(FloatParts256
*a
)
1143 a
->frac_lo
= usub64_borrow(0, a
->frac_lo
, &c
);
1144 a
->frac_lm
= usub64_borrow(0, a
->frac_lm
, &c
);
1145 a
->frac_hm
= usub64_borrow(0, a
->frac_hm
, &c
);
1146 a
->frac_hi
= usub64_borrow(0, a
->frac_hi
, &c
);
1149 #define frac_neg(A) FRAC_GENERIC_64_128_256(neg, A)(A)
1151 static int frac64_normalize(FloatParts64
*a
)
1154 int shift
= clz64(a
->frac
);
1161 static int frac128_normalize(FloatParts128
*a
)
1164 int shl
= clz64(a
->frac_hi
);
1165 a
->frac_hi
= shl_double(a
->frac_hi
, a
->frac_lo
, shl
);
1168 } else if (a
->frac_lo
) {
1169 int shl
= clz64(a
->frac_lo
);
1170 a
->frac_hi
= a
->frac_lo
<< shl
;
1177 static int frac256_normalize(FloatParts256
*a
)
1179 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_hm
;
1180 uint64_t a2
= a
->frac_lm
, a3
= a
->frac_lo
;
1192 a0
= a1
, a1
= a2
, a2
= a3
, a3
= 0;
1195 a0
= a2
, a1
= a3
, a2
= 0, a3
= 0;
1198 a0
= a3
, a1
= 0, a2
= 0, a3
= 0;
1201 a0
= 0, a1
= 0, a2
= 0, a3
= 0;
1211 a0
= shl_double(a0
, a1
, shl
);
1212 a1
= shl_double(a1
, a2
, shl
);
1213 a2
= shl_double(a2
, a3
, shl
);
1224 #define frac_normalize(A) FRAC_GENERIC_64_128_256(normalize, A)(A)
1226 static void frac64_shl(FloatParts64
*a
, int c
)
1231 static void frac128_shl(FloatParts128
*a
, int c
)
1233 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_lo
;
1241 a0
= shl_double(a0
, a1
, c
);
1249 #define frac_shl(A, C) FRAC_GENERIC_64_128(shl, A)(A, C)
1251 static void frac64_shr(FloatParts64
*a
, int c
)
1256 static void frac128_shr(FloatParts128
*a
, int c
)
1258 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_lo
;
1266 a1
= shr_double(a0
, a1
, c
);
1274 #define frac_shr(A, C) FRAC_GENERIC_64_128(shr, A)(A, C)
1276 static void frac64_shrjam(FloatParts64
*a
, int c
)
1278 uint64_t a0
= a
->frac
;
1280 if (likely(c
!= 0)) {
1281 if (likely(c
< 64)) {
1282 a0
= (a0
>> c
) | (shr_double(a0
, 0, c
) != 0);
1290 static void frac128_shrjam(FloatParts128
*a
, int c
)
1292 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_lo
;
1293 uint64_t sticky
= 0;
1295 if (unlikely(c
== 0)) {
1297 } else if (likely(c
< 64)) {
1299 } else if (likely(c
< 128)) {
1313 sticky
|= shr_double(a1
, 0, c
);
1314 a1
= shr_double(a0
, a1
, c
);
1318 a
->frac_lo
= a1
| (sticky
!= 0);
1322 static void frac256_shrjam(FloatParts256
*a
, int c
)
1324 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_hm
;
1325 uint64_t a2
= a
->frac_lm
, a3
= a
->frac_lo
;
1326 uint64_t sticky
= 0;
1328 if (unlikely(c
== 0)) {
1330 } else if (likely(c
< 64)) {
1332 } else if (likely(c
< 256)) {
1333 if (unlikely(c
& 128)) {
1335 a3
= a1
, a2
= a0
, a1
= 0, a0
= 0;
1337 if (unlikely(c
& 64)) {
1339 a3
= a2
, a2
= a1
, a1
= a0
, a0
= 0;
1346 sticky
= a0
| a1
| a2
| a3
;
1347 a0
= a1
= a2
= a3
= 0;
1351 sticky
|= shr_double(a3
, 0, c
);
1352 a3
= shr_double(a2
, a3
, c
);
1353 a2
= shr_double(a1
, a2
, c
);
1354 a1
= shr_double(a0
, a1
, c
);
1358 a
->frac_lo
= a3
| (sticky
!= 0);
1364 #define frac_shrjam(A, C) FRAC_GENERIC_64_128_256(shrjam, A)(A, C)
1366 static bool frac64_sub(FloatParts64
*r
, FloatParts64
*a
, FloatParts64
*b
)
1368 return usub64_overflow(a
->frac
, b
->frac
, &r
->frac
);
1371 static bool frac128_sub(FloatParts128
*r
, FloatParts128
*a
, FloatParts128
*b
)
1374 r
->frac_lo
= usub64_borrow(a
->frac_lo
, b
->frac_lo
, &c
);
1375 r
->frac_hi
= usub64_borrow(a
->frac_hi
, b
->frac_hi
, &c
);
1379 static bool frac256_sub(FloatParts256
*r
, FloatParts256
*a
, FloatParts256
*b
)
1382 r
->frac_lo
= usub64_borrow(a
->frac_lo
, b
->frac_lo
, &c
);
1383 r
->frac_lm
= usub64_borrow(a
->frac_lm
, b
->frac_lm
, &c
);
1384 r
->frac_hm
= usub64_borrow(a
->frac_hm
, b
->frac_hm
, &c
);
1385 r
->frac_hi
= usub64_borrow(a
->frac_hi
, b
->frac_hi
, &c
);
1389 #define frac_sub(R, A, B) FRAC_GENERIC_64_128_256(sub, R)(R, A, B)
1391 static void frac64_truncjam(FloatParts64
*r
, FloatParts128
*a
)
1393 r
->frac
= a
->frac_hi
| (a
->frac_lo
!= 0);
1396 static void frac128_truncjam(FloatParts128
*r
, FloatParts256
*a
)
1398 r
->frac_hi
= a
->frac_hi
;
1399 r
->frac_lo
= a
->frac_hm
| ((a
->frac_lm
| a
->frac_lo
) != 0);
1402 #define frac_truncjam(R, A) FRAC_GENERIC_64_128(truncjam, R)(R, A)
1404 static void frac64_widen(FloatParts128
*r
, FloatParts64
*a
)
1406 r
->frac_hi
= a
->frac
;
1410 static void frac128_widen(FloatParts256
*r
, FloatParts128
*a
)
1412 r
->frac_hi
= a
->frac_hi
;
1413 r
->frac_hm
= a
->frac_lo
;
1418 #define frac_widen(A, B) FRAC_GENERIC_64_128(widen, B)(A, B)
1421 * Reciprocal sqrt table. 1 bit of exponent, 6-bits of mantessa.
1422 * From https://git.musl-libc.org/cgit/musl/tree/src/math/sqrt_data.c
1423 * and thus MIT licenced.
1425 static const uint16_t rsqrt_tab
[128] = {
1426 0xb451, 0xb2f0, 0xb196, 0xb044, 0xaef9, 0xadb6, 0xac79, 0xab43,
1427 0xaa14, 0xa8eb, 0xa7c8, 0xa6aa, 0xa592, 0xa480, 0xa373, 0xa26b,
1428 0xa168, 0xa06a, 0x9f70, 0x9e7b, 0x9d8a, 0x9c9d, 0x9bb5, 0x9ad1,
1429 0x99f0, 0x9913, 0x983a, 0x9765, 0x9693, 0x95c4, 0x94f8, 0x9430,
1430 0x936b, 0x92a9, 0x91ea, 0x912e, 0x9075, 0x8fbe, 0x8f0a, 0x8e59,
1431 0x8daa, 0x8cfe, 0x8c54, 0x8bac, 0x8b07, 0x8a64, 0x89c4, 0x8925,
1432 0x8889, 0x87ee, 0x8756, 0x86c0, 0x862b, 0x8599, 0x8508, 0x8479,
1433 0x83ec, 0x8361, 0x82d8, 0x8250, 0x81c9, 0x8145, 0x80c2, 0x8040,
1434 0xff02, 0xfd0e, 0xfb25, 0xf947, 0xf773, 0xf5aa, 0xf3ea, 0xf234,
1435 0xf087, 0xeee3, 0xed47, 0xebb3, 0xea27, 0xe8a3, 0xe727, 0xe5b2,
1436 0xe443, 0xe2dc, 0xe17a, 0xe020, 0xdecb, 0xdd7d, 0xdc34, 0xdaf1,
1437 0xd9b3, 0xd87b, 0xd748, 0xd61a, 0xd4f1, 0xd3cd, 0xd2ad, 0xd192,
1438 0xd07b, 0xcf69, 0xce5b, 0xcd51, 0xcc4a, 0xcb48, 0xca4a, 0xc94f,
1439 0xc858, 0xc764, 0xc674, 0xc587, 0xc49d, 0xc3b7, 0xc2d4, 0xc1f4,
1440 0xc116, 0xc03c, 0xbf65, 0xbe90, 0xbdbe, 0xbcef, 0xbc23, 0xbb59,
1441 0xba91, 0xb9cc, 0xb90a, 0xb84a, 0xb78c, 0xb6d0, 0xb617, 0xb560,
1444 #define partsN(NAME) glue(glue(glue(parts,N),_),NAME)
1445 #define FloatPartsN glue(FloatParts,N)
1446 #define FloatPartsW glue(FloatParts,W)
1451 #include "softfloat-parts-addsub.c.inc"
1452 #include "softfloat-parts.c.inc"
1459 #include "softfloat-parts-addsub.c.inc"
1460 #include "softfloat-parts.c.inc"
1466 #include "softfloat-parts-addsub.c.inc"
1475 * Pack/unpack routines with a specific FloatFmt.
1478 static void float16a_unpack_canonical(FloatParts64
*p
, float16 f
,
1479 float_status
*s
, const FloatFmt
*params
)
1481 float16_unpack_raw(p
, f
);
1482 parts_canonicalize(p
, s
, params
);
1485 static void float16_unpack_canonical(FloatParts64
*p
, float16 f
,
1488 float16a_unpack_canonical(p
, f
, s
, &float16_params
);
1491 static void bfloat16_unpack_canonical(FloatParts64
*p
, bfloat16 f
,
1494 bfloat16_unpack_raw(p
, f
);
1495 parts_canonicalize(p
, s
, &bfloat16_params
);
1498 static float16
float16a_round_pack_canonical(FloatParts64
*p
,
1500 const FloatFmt
*params
)
1502 parts_uncanon(p
, s
, params
);
1503 return float16_pack_raw(p
);
1506 static float16
float16_round_pack_canonical(FloatParts64
*p
,
1509 return float16a_round_pack_canonical(p
, s
, &float16_params
);
1512 static bfloat16
bfloat16_round_pack_canonical(FloatParts64
*p
,
1515 parts_uncanon(p
, s
, &bfloat16_params
);
1516 return bfloat16_pack_raw(p
);
1519 static void float32_unpack_canonical(FloatParts64
*p
, float32 f
,
1522 float32_unpack_raw(p
, f
);
1523 parts_canonicalize(p
, s
, &float32_params
);
1526 static float32
float32_round_pack_canonical(FloatParts64
*p
,
1529 parts_uncanon(p
, s
, &float32_params
);
1530 return float32_pack_raw(p
);
1533 static void float64_unpack_canonical(FloatParts64
*p
, float64 f
,
1536 float64_unpack_raw(p
, f
);
1537 parts_canonicalize(p
, s
, &float64_params
);
1540 static float64
float64_round_pack_canonical(FloatParts64
*p
,
1543 parts_uncanon(p
, s
, &float64_params
);
1544 return float64_pack_raw(p
);
1547 static void float128_unpack_canonical(FloatParts128
*p
, float128 f
,
1550 float128_unpack_raw(p
, f
);
1551 parts_canonicalize(p
, s
, &float128_params
);
1554 static float128
float128_round_pack_canonical(FloatParts128
*p
,
1557 parts_uncanon(p
, s
, &float128_params
);
1558 return float128_pack_raw(p
);
1561 /* Returns false if the encoding is invalid. */
1562 static bool floatx80_unpack_canonical(FloatParts128
*p
, floatx80 f
,
1565 /* Ensure rounding precision is set before beginning. */
1566 switch (s
->floatx80_rounding_precision
) {
1567 case floatx80_precision_x
:
1568 case floatx80_precision_d
:
1569 case floatx80_precision_s
:
1572 g_assert_not_reached();
1575 if (unlikely(floatx80_invalid_encoding(f
))) {
1576 float_raise(float_flag_invalid
, s
);
1580 floatx80_unpack_raw(p
, f
);
1582 if (likely(p
->exp
!= floatx80_params
[floatx80_precision_x
].exp_max
)) {
1583 parts_canonicalize(p
, s
, &floatx80_params
[floatx80_precision_x
]);
1585 /* The explicit integer bit is ignored, after invalid checks. */
1586 p
->frac_hi
&= MAKE_64BIT_MASK(0, 63);
1587 p
->cls
= (p
->frac_hi
== 0 ? float_class_inf
1588 : parts_is_snan_frac(p
->frac_hi
, s
)
1589 ? float_class_snan
: float_class_qnan
);
1594 static floatx80
floatx80_round_pack_canonical(FloatParts128
*p
,
1597 const FloatFmt
*fmt
= &floatx80_params
[s
->floatx80_rounding_precision
];
1602 case float_class_normal
:
1603 if (s
->floatx80_rounding_precision
== floatx80_precision_x
) {
1604 parts_uncanon_normal(p
, s
, fmt
);
1612 frac_truncjam(&p64
, p
);
1613 parts_uncanon_normal(&p64
, s
, fmt
);
1617 if (exp
!= fmt
->exp_max
) {
1620 /* rounded to inf -- fall through to set frac correctly */
1622 case float_class_inf
:
1623 /* x86 and m68k differ in the setting of the integer bit. */
1624 frac
= floatx80_infinity_low
;
1628 case float_class_zero
:
1633 case float_class_snan
:
1634 case float_class_qnan
:
1635 /* NaNs have the integer bit set. */
1636 frac
= p
->frac_hi
| (1ull << 63);
1641 g_assert_not_reached();
1644 return packFloatx80(p
->sign
, exp
, frac
);
1648 * Addition and subtraction
1651 static float16 QEMU_FLATTEN
1652 float16_addsub(float16 a
, float16 b
, float_status
*status
, bool subtract
)
1654 FloatParts64 pa
, pb
, *pr
;
1656 float16_unpack_canonical(&pa
, a
, status
);
1657 float16_unpack_canonical(&pb
, b
, status
);
1658 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1660 return float16_round_pack_canonical(pr
, status
);
1663 float16
float16_add(float16 a
, float16 b
, float_status
*status
)
1665 return float16_addsub(a
, b
, status
, false);
1668 float16
float16_sub(float16 a
, float16 b
, float_status
*status
)
1670 return float16_addsub(a
, b
, status
, true);
1673 static float32 QEMU_SOFTFLOAT_ATTR
1674 soft_f32_addsub(float32 a
, float32 b
, float_status
*status
, bool subtract
)
1676 FloatParts64 pa
, pb
, *pr
;
1678 float32_unpack_canonical(&pa
, a
, status
);
1679 float32_unpack_canonical(&pb
, b
, status
);
1680 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1682 return float32_round_pack_canonical(pr
, status
);
1685 static float32
soft_f32_add(float32 a
, float32 b
, float_status
*status
)
1687 return soft_f32_addsub(a
, b
, status
, false);
1690 static float32
soft_f32_sub(float32 a
, float32 b
, float_status
*status
)
1692 return soft_f32_addsub(a
, b
, status
, true);
1695 static float64 QEMU_SOFTFLOAT_ATTR
1696 soft_f64_addsub(float64 a
, float64 b
, float_status
*status
, bool subtract
)
1698 FloatParts64 pa
, pb
, *pr
;
1700 float64_unpack_canonical(&pa
, a
, status
);
1701 float64_unpack_canonical(&pb
, b
, status
);
1702 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1704 return float64_round_pack_canonical(pr
, status
);
1707 static float64
soft_f64_add(float64 a
, float64 b
, float_status
*status
)
1709 return soft_f64_addsub(a
, b
, status
, false);
1712 static float64
soft_f64_sub(float64 a
, float64 b
, float_status
*status
)
1714 return soft_f64_addsub(a
, b
, status
, true);
1717 static float hard_f32_add(float a
, float b
)
1722 static float hard_f32_sub(float a
, float b
)
1727 static double hard_f64_add(double a
, double b
)
1732 static double hard_f64_sub(double a
, double b
)
1737 static bool f32_addsubmul_post(union_float32 a
, union_float32 b
)
1739 if (QEMU_HARDFLOAT_2F32_USE_FP
) {
1740 return !(fpclassify(a
.h
) == FP_ZERO
&& fpclassify(b
.h
) == FP_ZERO
);
1742 return !(float32_is_zero(a
.s
) && float32_is_zero(b
.s
));
1745 static bool f64_addsubmul_post(union_float64 a
, union_float64 b
)
1747 if (QEMU_HARDFLOAT_2F64_USE_FP
) {
1748 return !(fpclassify(a
.h
) == FP_ZERO
&& fpclassify(b
.h
) == FP_ZERO
);
1750 return !(float64_is_zero(a
.s
) && float64_is_zero(b
.s
));
1754 static float32
float32_addsub(float32 a
, float32 b
, float_status
*s
,
1755 hard_f32_op2_fn hard
, soft_f32_op2_fn soft
)
1757 return float32_gen2(a
, b
, s
, hard
, soft
,
1758 f32_is_zon2
, f32_addsubmul_post
);
1761 static float64
float64_addsub(float64 a
, float64 b
, float_status
*s
,
1762 hard_f64_op2_fn hard
, soft_f64_op2_fn soft
)
1764 return float64_gen2(a
, b
, s
, hard
, soft
,
1765 f64_is_zon2
, f64_addsubmul_post
);
1768 float32 QEMU_FLATTEN
1769 float32_add(float32 a
, float32 b
, float_status
*s
)
1771 return float32_addsub(a
, b
, s
, hard_f32_add
, soft_f32_add
);
1774 float32 QEMU_FLATTEN
1775 float32_sub(float32 a
, float32 b
, float_status
*s
)
1777 return float32_addsub(a
, b
, s
, hard_f32_sub
, soft_f32_sub
);
1780 float64 QEMU_FLATTEN
1781 float64_add(float64 a
, float64 b
, float_status
*s
)
1783 return float64_addsub(a
, b
, s
, hard_f64_add
, soft_f64_add
);
1786 float64 QEMU_FLATTEN
1787 float64_sub(float64 a
, float64 b
, float_status
*s
)
1789 return float64_addsub(a
, b
, s
, hard_f64_sub
, soft_f64_sub
);
1792 static bfloat16 QEMU_FLATTEN
1793 bfloat16_addsub(bfloat16 a
, bfloat16 b
, float_status
*status
, bool subtract
)
1795 FloatParts64 pa
, pb
, *pr
;
1797 bfloat16_unpack_canonical(&pa
, a
, status
);
1798 bfloat16_unpack_canonical(&pb
, b
, status
);
1799 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1801 return bfloat16_round_pack_canonical(pr
, status
);
1804 bfloat16
bfloat16_add(bfloat16 a
, bfloat16 b
, float_status
*status
)
1806 return bfloat16_addsub(a
, b
, status
, false);
1809 bfloat16
bfloat16_sub(bfloat16 a
, bfloat16 b
, float_status
*status
)
1811 return bfloat16_addsub(a
, b
, status
, true);
1814 static float128 QEMU_FLATTEN
1815 float128_addsub(float128 a
, float128 b
, float_status
*status
, bool subtract
)
1817 FloatParts128 pa
, pb
, *pr
;
1819 float128_unpack_canonical(&pa
, a
, status
);
1820 float128_unpack_canonical(&pb
, b
, status
);
1821 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1823 return float128_round_pack_canonical(pr
, status
);
1826 float128
float128_add(float128 a
, float128 b
, float_status
*status
)
1828 return float128_addsub(a
, b
, status
, false);
1831 float128
float128_sub(float128 a
, float128 b
, float_status
*status
)
1833 return float128_addsub(a
, b
, status
, true);
1836 static floatx80 QEMU_FLATTEN
1837 floatx80_addsub(floatx80 a
, floatx80 b
, float_status
*status
, bool subtract
)
1839 FloatParts128 pa
, pb
, *pr
;
1841 if (!floatx80_unpack_canonical(&pa
, a
, status
) ||
1842 !floatx80_unpack_canonical(&pb
, b
, status
)) {
1843 return floatx80_default_nan(status
);
1846 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1847 return floatx80_round_pack_canonical(pr
, status
);
1850 floatx80
floatx80_add(floatx80 a
, floatx80 b
, float_status
*status
)
1852 return floatx80_addsub(a
, b
, status
, false);
1855 floatx80
floatx80_sub(floatx80 a
, floatx80 b
, float_status
*status
)
1857 return floatx80_addsub(a
, b
, status
, true);
1864 float16 QEMU_FLATTEN
float16_mul(float16 a
, float16 b
, float_status
*status
)
1866 FloatParts64 pa
, pb
, *pr
;
1868 float16_unpack_canonical(&pa
, a
, status
);
1869 float16_unpack_canonical(&pb
, b
, status
);
1870 pr
= parts_mul(&pa
, &pb
, status
);
1872 return float16_round_pack_canonical(pr
, status
);
1875 static float32 QEMU_SOFTFLOAT_ATTR
1876 soft_f32_mul(float32 a
, float32 b
, float_status
*status
)
1878 FloatParts64 pa
, pb
, *pr
;
1880 float32_unpack_canonical(&pa
, a
, status
);
1881 float32_unpack_canonical(&pb
, b
, status
);
1882 pr
= parts_mul(&pa
, &pb
, status
);
1884 return float32_round_pack_canonical(pr
, status
);
1887 static float64 QEMU_SOFTFLOAT_ATTR
1888 soft_f64_mul(float64 a
, float64 b
, float_status
*status
)
1890 FloatParts64 pa
, pb
, *pr
;
1892 float64_unpack_canonical(&pa
, a
, status
);
1893 float64_unpack_canonical(&pb
, b
, status
);
1894 pr
= parts_mul(&pa
, &pb
, status
);
1896 return float64_round_pack_canonical(pr
, status
);
1899 static float hard_f32_mul(float a
, float b
)
1904 static double hard_f64_mul(double a
, double b
)
1909 float32 QEMU_FLATTEN
1910 float32_mul(float32 a
, float32 b
, float_status
*s
)
1912 return float32_gen2(a
, b
, s
, hard_f32_mul
, soft_f32_mul
,
1913 f32_is_zon2
, f32_addsubmul_post
);
1916 float64 QEMU_FLATTEN
1917 float64_mul(float64 a
, float64 b
, float_status
*s
)
1919 return float64_gen2(a
, b
, s
, hard_f64_mul
, soft_f64_mul
,
1920 f64_is_zon2
, f64_addsubmul_post
);
1923 bfloat16 QEMU_FLATTEN
1924 bfloat16_mul(bfloat16 a
, bfloat16 b
, float_status
*status
)
1926 FloatParts64 pa
, pb
, *pr
;
1928 bfloat16_unpack_canonical(&pa
, a
, status
);
1929 bfloat16_unpack_canonical(&pb
, b
, status
);
1930 pr
= parts_mul(&pa
, &pb
, status
);
1932 return bfloat16_round_pack_canonical(pr
, status
);
1935 float128 QEMU_FLATTEN
1936 float128_mul(float128 a
, float128 b
, float_status
*status
)
1938 FloatParts128 pa
, pb
, *pr
;
1940 float128_unpack_canonical(&pa
, a
, status
);
1941 float128_unpack_canonical(&pb
, b
, status
);
1942 pr
= parts_mul(&pa
, &pb
, status
);
1944 return float128_round_pack_canonical(pr
, status
);
1947 floatx80 QEMU_FLATTEN
1948 floatx80_mul(floatx80 a
, floatx80 b
, float_status
*status
)
1950 FloatParts128 pa
, pb
, *pr
;
1952 if (!floatx80_unpack_canonical(&pa
, a
, status
) ||
1953 !floatx80_unpack_canonical(&pb
, b
, status
)) {
1954 return floatx80_default_nan(status
);
1957 pr
= parts_mul(&pa
, &pb
, status
);
1958 return floatx80_round_pack_canonical(pr
, status
);
1962 * Fused multiply-add
1965 float16 QEMU_FLATTEN
float16_muladd(float16 a
, float16 b
, float16 c
,
1966 int flags
, float_status
*status
)
1968 FloatParts64 pa
, pb
, pc
, *pr
;
1970 float16_unpack_canonical(&pa
, a
, status
);
1971 float16_unpack_canonical(&pb
, b
, status
);
1972 float16_unpack_canonical(&pc
, c
, status
);
1973 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1975 return float16_round_pack_canonical(pr
, status
);
1978 static float32 QEMU_SOFTFLOAT_ATTR
1979 soft_f32_muladd(float32 a
, float32 b
, float32 c
, int flags
,
1980 float_status
*status
)
1982 FloatParts64 pa
, pb
, pc
, *pr
;
1984 float32_unpack_canonical(&pa
, a
, status
);
1985 float32_unpack_canonical(&pb
, b
, status
);
1986 float32_unpack_canonical(&pc
, c
, status
);
1987 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1989 return float32_round_pack_canonical(pr
, status
);
1992 static float64 QEMU_SOFTFLOAT_ATTR
1993 soft_f64_muladd(float64 a
, float64 b
, float64 c
, int flags
,
1994 float_status
*status
)
1996 FloatParts64 pa
, pb
, pc
, *pr
;
1998 float64_unpack_canonical(&pa
, a
, status
);
1999 float64_unpack_canonical(&pb
, b
, status
);
2000 float64_unpack_canonical(&pc
, c
, status
);
2001 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
2003 return float64_round_pack_canonical(pr
, status
);
2006 static bool force_soft_fma
;
2008 float32 QEMU_FLATTEN
2009 float32_muladd(float32 xa
, float32 xb
, float32 xc
, int flags
, float_status
*s
)
2011 union_float32 ua
, ub
, uc
, ur
;
2017 if (unlikely(!can_use_fpu(s
))) {
2020 if (unlikely(flags
& float_muladd_halve_result
)) {
2024 float32_input_flush3(&ua
.s
, &ub
.s
, &uc
.s
, s
);
2025 if (unlikely(!f32_is_zon3(ua
, ub
, uc
))) {
2029 if (unlikely(force_soft_fma
)) {
2034 * When (a || b) == 0, there's no need to check for under/over flow,
2035 * since we know the addend is (normal || 0) and the product is 0.
2037 if (float32_is_zero(ua
.s
) || float32_is_zero(ub
.s
)) {
2041 prod_sign
= float32_is_neg(ua
.s
) ^ float32_is_neg(ub
.s
);
2042 prod_sign
^= !!(flags
& float_muladd_negate_product
);
2043 up
.s
= float32_set_sign(float32_zero
, prod_sign
);
2045 if (flags
& float_muladd_negate_c
) {
2050 union_float32 ua_orig
= ua
;
2051 union_float32 uc_orig
= uc
;
2053 if (flags
& float_muladd_negate_product
) {
2056 if (flags
& float_muladd_negate_c
) {
2060 ur
.h
= fmaf(ua
.h
, ub
.h
, uc
.h
);
2062 if (unlikely(f32_is_inf(ur
))) {
2063 float_raise(float_flag_overflow
, s
);
2064 } else if (unlikely(fabsf(ur
.h
) <= FLT_MIN
)) {
2070 if (flags
& float_muladd_negate_result
) {
2071 return float32_chs(ur
.s
);
2076 return soft_f32_muladd(ua
.s
, ub
.s
, uc
.s
, flags
, s
);
2079 float64 QEMU_FLATTEN
2080 float64_muladd(float64 xa
, float64 xb
, float64 xc
, int flags
, float_status
*s
)
2082 union_float64 ua
, ub
, uc
, ur
;
2088 if (unlikely(!can_use_fpu(s
))) {
2091 if (unlikely(flags
& float_muladd_halve_result
)) {
2095 float64_input_flush3(&ua
.s
, &ub
.s
, &uc
.s
, s
);
2096 if (unlikely(!f64_is_zon3(ua
, ub
, uc
))) {
2100 if (unlikely(force_soft_fma
)) {
2105 * When (a || b) == 0, there's no need to check for under/over flow,
2106 * since we know the addend is (normal || 0) and the product is 0.
2108 if (float64_is_zero(ua
.s
) || float64_is_zero(ub
.s
)) {
2112 prod_sign
= float64_is_neg(ua
.s
) ^ float64_is_neg(ub
.s
);
2113 prod_sign
^= !!(flags
& float_muladd_negate_product
);
2114 up
.s
= float64_set_sign(float64_zero
, prod_sign
);
2116 if (flags
& float_muladd_negate_c
) {
2121 union_float64 ua_orig
= ua
;
2122 union_float64 uc_orig
= uc
;
2124 if (flags
& float_muladd_negate_product
) {
2127 if (flags
& float_muladd_negate_c
) {
2131 ur
.h
= fma(ua
.h
, ub
.h
, uc
.h
);
2133 if (unlikely(f64_is_inf(ur
))) {
2134 float_raise(float_flag_overflow
, s
);
2135 } else if (unlikely(fabs(ur
.h
) <= FLT_MIN
)) {
2141 if (flags
& float_muladd_negate_result
) {
2142 return float64_chs(ur
.s
);
2147 return soft_f64_muladd(ua
.s
, ub
.s
, uc
.s
, flags
, s
);
2150 bfloat16 QEMU_FLATTEN
bfloat16_muladd(bfloat16 a
, bfloat16 b
, bfloat16 c
,
2151 int flags
, float_status
*status
)
2153 FloatParts64 pa
, pb
, pc
, *pr
;
2155 bfloat16_unpack_canonical(&pa
, a
, status
);
2156 bfloat16_unpack_canonical(&pb
, b
, status
);
2157 bfloat16_unpack_canonical(&pc
, c
, status
);
2158 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
2160 return bfloat16_round_pack_canonical(pr
, status
);
2163 float128 QEMU_FLATTEN
float128_muladd(float128 a
, float128 b
, float128 c
,
2164 int flags
, float_status
*status
)
2166 FloatParts128 pa
, pb
, pc
, *pr
;
2168 float128_unpack_canonical(&pa
, a
, status
);
2169 float128_unpack_canonical(&pb
, b
, status
);
2170 float128_unpack_canonical(&pc
, c
, status
);
2171 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
2173 return float128_round_pack_canonical(pr
, status
);
2180 float16
float16_div(float16 a
, float16 b
, float_status
*status
)
2182 FloatParts64 pa
, pb
, *pr
;
2184 float16_unpack_canonical(&pa
, a
, status
);
2185 float16_unpack_canonical(&pb
, b
, status
);
2186 pr
= parts_div(&pa
, &pb
, status
);
2188 return float16_round_pack_canonical(pr
, status
);
2191 static float32 QEMU_SOFTFLOAT_ATTR
2192 soft_f32_div(float32 a
, float32 b
, float_status
*status
)
2194 FloatParts64 pa
, pb
, *pr
;
2196 float32_unpack_canonical(&pa
, a
, status
);
2197 float32_unpack_canonical(&pb
, b
, status
);
2198 pr
= parts_div(&pa
, &pb
, status
);
2200 return float32_round_pack_canonical(pr
, status
);
2203 static float64 QEMU_SOFTFLOAT_ATTR
2204 soft_f64_div(float64 a
, float64 b
, float_status
*status
)
2206 FloatParts64 pa
, pb
, *pr
;
2208 float64_unpack_canonical(&pa
, a
, status
);
2209 float64_unpack_canonical(&pb
, b
, status
);
2210 pr
= parts_div(&pa
, &pb
, status
);
2212 return float64_round_pack_canonical(pr
, status
);
2215 static float hard_f32_div(float a
, float b
)
2220 static double hard_f64_div(double a
, double b
)
2225 static bool f32_div_pre(union_float32 a
, union_float32 b
)
2227 if (QEMU_HARDFLOAT_2F32_USE_FP
) {
2228 return (fpclassify(a
.h
) == FP_NORMAL
|| fpclassify(a
.h
) == FP_ZERO
) &&
2229 fpclassify(b
.h
) == FP_NORMAL
;
2231 return float32_is_zero_or_normal(a
.s
) && float32_is_normal(b
.s
);
2234 static bool f64_div_pre(union_float64 a
, union_float64 b
)
2236 if (QEMU_HARDFLOAT_2F64_USE_FP
) {
2237 return (fpclassify(a
.h
) == FP_NORMAL
|| fpclassify(a
.h
) == FP_ZERO
) &&
2238 fpclassify(b
.h
) == FP_NORMAL
;
2240 return float64_is_zero_or_normal(a
.s
) && float64_is_normal(b
.s
);
2243 static bool f32_div_post(union_float32 a
, union_float32 b
)
2245 if (QEMU_HARDFLOAT_2F32_USE_FP
) {
2246 return fpclassify(a
.h
) != FP_ZERO
;
2248 return !float32_is_zero(a
.s
);
2251 static bool f64_div_post(union_float64 a
, union_float64 b
)
2253 if (QEMU_HARDFLOAT_2F64_USE_FP
) {
2254 return fpclassify(a
.h
) != FP_ZERO
;
2256 return !float64_is_zero(a
.s
);
2259 float32 QEMU_FLATTEN
2260 float32_div(float32 a
, float32 b
, float_status
*s
)
2262 return float32_gen2(a
, b
, s
, hard_f32_div
, soft_f32_div
,
2263 f32_div_pre
, f32_div_post
);
2266 float64 QEMU_FLATTEN
2267 float64_div(float64 a
, float64 b
, float_status
*s
)
2269 return float64_gen2(a
, b
, s
, hard_f64_div
, soft_f64_div
,
2270 f64_div_pre
, f64_div_post
);
2273 bfloat16 QEMU_FLATTEN
2274 bfloat16_div(bfloat16 a
, bfloat16 b
, float_status
*status
)
2276 FloatParts64 pa
, pb
, *pr
;
2278 bfloat16_unpack_canonical(&pa
, a
, status
);
2279 bfloat16_unpack_canonical(&pb
, b
, status
);
2280 pr
= parts_div(&pa
, &pb
, status
);
2282 return bfloat16_round_pack_canonical(pr
, status
);
2285 float128 QEMU_FLATTEN
2286 float128_div(float128 a
, float128 b
, float_status
*status
)
2288 FloatParts128 pa
, pb
, *pr
;
2290 float128_unpack_canonical(&pa
, a
, status
);
2291 float128_unpack_canonical(&pb
, b
, status
);
2292 pr
= parts_div(&pa
, &pb
, status
);
2294 return float128_round_pack_canonical(pr
, status
);
2297 floatx80
floatx80_div(floatx80 a
, floatx80 b
, float_status
*status
)
2299 FloatParts128 pa
, pb
, *pr
;
2301 if (!floatx80_unpack_canonical(&pa
, a
, status
) ||
2302 !floatx80_unpack_canonical(&pb
, b
, status
)) {
2303 return floatx80_default_nan(status
);
2306 pr
= parts_div(&pa
, &pb
, status
);
2307 return floatx80_round_pack_canonical(pr
, status
);
2311 * Float to Float conversions
2313 * Returns the result of converting one float format to another. The
2314 * conversion is performed according to the IEC/IEEE Standard for
2315 * Binary Floating-Point Arithmetic.
2317 * Usually this only needs to take care of raising invalid exceptions
2318 * and handling the conversion on NaNs.
2321 static void parts_float_to_ahp(FloatParts64
*a
, float_status
*s
)
2324 case float_class_qnan
:
2325 case float_class_snan
:
2327 * There is no NaN in the destination format. Raise Invalid
2328 * and return a zero with the sign of the input NaN.
2330 float_raise(float_flag_invalid
, s
);
2331 a
->cls
= float_class_zero
;
2334 case float_class_inf
:
2336 * There is no Inf in the destination format. Raise Invalid
2337 * and return the maximum normal with the correct sign.
2339 float_raise(float_flag_invalid
, s
);
2340 a
->cls
= float_class_normal
;
2341 a
->exp
= float16_params_ahp
.exp_max
;
2342 a
->frac
= MAKE_64BIT_MASK(float16_params_ahp
.frac_shift
,
2343 float16_params_ahp
.frac_size
+ 1);
2346 case float_class_normal
:
2347 case float_class_zero
:
2351 g_assert_not_reached();
2355 static void parts64_float_to_float(FloatParts64
*a
, float_status
*s
)
2357 if (is_nan(a
->cls
)) {
2358 parts_return_nan(a
, s
);
2362 static void parts128_float_to_float(FloatParts128
*a
, float_status
*s
)
2364 if (is_nan(a
->cls
)) {
2365 parts_return_nan(a
, s
);
2369 #define parts_float_to_float(P, S) \
2370 PARTS_GENERIC_64_128(float_to_float, P)(P, S)
2372 static void parts_float_to_float_narrow(FloatParts64
*a
, FloatParts128
*b
,
2379 if (a
->cls
== float_class_normal
) {
2380 frac_truncjam(a
, b
);
2381 } else if (is_nan(a
->cls
)) {
2382 /* Discard the low bits of the NaN. */
2383 a
->frac
= b
->frac_hi
;
2384 parts_return_nan(a
, s
);
2388 static void parts_float_to_float_widen(FloatParts128
*a
, FloatParts64
*b
,
2396 if (is_nan(a
->cls
)) {
2397 parts_return_nan(a
, s
);
2401 float32
float16_to_float32(float16 a
, bool ieee
, float_status
*s
)
2403 const FloatFmt
*fmt16
= ieee
? &float16_params
: &float16_params_ahp
;
2406 float16a_unpack_canonical(&p
, a
, s
, fmt16
);
2407 parts_float_to_float(&p
, s
);
2408 return float32_round_pack_canonical(&p
, s
);
2411 float64
float16_to_float64(float16 a
, bool ieee
, float_status
*s
)
2413 const FloatFmt
*fmt16
= ieee
? &float16_params
: &float16_params_ahp
;
2416 float16a_unpack_canonical(&p
, a
, s
, fmt16
);
2417 parts_float_to_float(&p
, s
);
2418 return float64_round_pack_canonical(&p
, s
);
2421 float16
float32_to_float16(float32 a
, bool ieee
, float_status
*s
)
2424 const FloatFmt
*fmt
;
2426 float32_unpack_canonical(&p
, a
, s
);
2428 parts_float_to_float(&p
, s
);
2429 fmt
= &float16_params
;
2431 parts_float_to_ahp(&p
, s
);
2432 fmt
= &float16_params_ahp
;
2434 return float16a_round_pack_canonical(&p
, s
, fmt
);
2437 static float64 QEMU_SOFTFLOAT_ATTR
2438 soft_float32_to_float64(float32 a
, float_status
*s
)
2442 float32_unpack_canonical(&p
, a
, s
);
2443 parts_float_to_float(&p
, s
);
2444 return float64_round_pack_canonical(&p
, s
);
2447 float64
float32_to_float64(float32 a
, float_status
*s
)
2449 if (likely(float32_is_normal(a
))) {
2450 /* Widening conversion can never produce inexact results. */
2456 } else if (float32_is_zero(a
)) {
2457 return float64_set_sign(float64_zero
, float32_is_neg(a
));
2459 return soft_float32_to_float64(a
, s
);
2463 float16
float64_to_float16(float64 a
, bool ieee
, float_status
*s
)
2466 const FloatFmt
*fmt
;
2468 float64_unpack_canonical(&p
, a
, s
);
2470 parts_float_to_float(&p
, s
);
2471 fmt
= &float16_params
;
2473 parts_float_to_ahp(&p
, s
);
2474 fmt
= &float16_params_ahp
;
2476 return float16a_round_pack_canonical(&p
, s
, fmt
);
2479 float32
float64_to_float32(float64 a
, float_status
*s
)
2483 float64_unpack_canonical(&p
, a
, s
);
2484 parts_float_to_float(&p
, s
);
2485 return float32_round_pack_canonical(&p
, s
);
2488 float32
bfloat16_to_float32(bfloat16 a
, float_status
*s
)
2492 bfloat16_unpack_canonical(&p
, a
, s
);
2493 parts_float_to_float(&p
, s
);
2494 return float32_round_pack_canonical(&p
, s
);
2497 float64
bfloat16_to_float64(bfloat16 a
, float_status
*s
)
2501 bfloat16_unpack_canonical(&p
, a
, s
);
2502 parts_float_to_float(&p
, s
);
2503 return float64_round_pack_canonical(&p
, s
);
2506 bfloat16
float32_to_bfloat16(float32 a
, float_status
*s
)
2510 float32_unpack_canonical(&p
, a
, s
);
2511 parts_float_to_float(&p
, s
);
2512 return bfloat16_round_pack_canonical(&p
, s
);
2515 bfloat16
float64_to_bfloat16(float64 a
, float_status
*s
)
2519 float64_unpack_canonical(&p
, a
, s
);
2520 parts_float_to_float(&p
, s
);
2521 return bfloat16_round_pack_canonical(&p
, s
);
2524 float32
float128_to_float32(float128 a
, float_status
*s
)
2529 float128_unpack_canonical(&p128
, a
, s
);
2530 parts_float_to_float_narrow(&p64
, &p128
, s
);
2531 return float32_round_pack_canonical(&p64
, s
);
2534 float64
float128_to_float64(float128 a
, float_status
*s
)
2539 float128_unpack_canonical(&p128
, a
, s
);
2540 parts_float_to_float_narrow(&p64
, &p128
, s
);
2541 return float64_round_pack_canonical(&p64
, s
);
2544 float128
float32_to_float128(float32 a
, float_status
*s
)
2549 float32_unpack_canonical(&p64
, a
, s
);
2550 parts_float_to_float_widen(&p128
, &p64
, s
);
2551 return float128_round_pack_canonical(&p128
, s
);
2554 float128
float64_to_float128(float64 a
, float_status
*s
)
2559 float64_unpack_canonical(&p64
, a
, s
);
2560 parts_float_to_float_widen(&p128
, &p64
, s
);
2561 return float128_round_pack_canonical(&p128
, s
);
2564 float32
floatx80_to_float32(floatx80 a
, float_status
*s
)
2569 if (floatx80_unpack_canonical(&p128
, a
, s
)) {
2570 parts_float_to_float_narrow(&p64
, &p128
, s
);
2572 parts_default_nan(&p64
, s
);
2574 return float32_round_pack_canonical(&p64
, s
);
2577 float64
floatx80_to_float64(floatx80 a
, float_status
*s
)
2582 if (floatx80_unpack_canonical(&p128
, a
, s
)) {
2583 parts_float_to_float_narrow(&p64
, &p128
, s
);
2585 parts_default_nan(&p64
, s
);
2587 return float64_round_pack_canonical(&p64
, s
);
2590 float128
floatx80_to_float128(floatx80 a
, float_status
*s
)
2594 if (floatx80_unpack_canonical(&p
, a
, s
)) {
2595 parts_float_to_float(&p
, s
);
2597 parts_default_nan(&p
, s
);
2599 return float128_round_pack_canonical(&p
, s
);
2602 floatx80
float32_to_floatx80(float32 a
, float_status
*s
)
2607 float32_unpack_canonical(&p64
, a
, s
);
2608 parts_float_to_float_widen(&p128
, &p64
, s
);
2609 return floatx80_round_pack_canonical(&p128
, s
);
2612 floatx80
float64_to_floatx80(float64 a
, float_status
*s
)
2617 float64_unpack_canonical(&p64
, a
, s
);
2618 parts_float_to_float_widen(&p128
, &p64
, s
);
2619 return floatx80_round_pack_canonical(&p128
, s
);
2622 floatx80
float128_to_floatx80(float128 a
, float_status
*s
)
2626 float128_unpack_canonical(&p
, a
, s
);
2627 parts_float_to_float(&p
, s
);
2628 return floatx80_round_pack_canonical(&p
, s
);
2632 * Round to integral value
2635 float16
float16_round_to_int(float16 a
, float_status
*s
)
2639 float16_unpack_canonical(&p
, a
, s
);
2640 parts_round_to_int(&p
, s
->float_rounding_mode
, 0, s
, &float16_params
);
2641 return float16_round_pack_canonical(&p
, s
);
2644 float32
float32_round_to_int(float32 a
, float_status
*s
)
2648 float32_unpack_canonical(&p
, a
, s
);
2649 parts_round_to_int(&p
, s
->float_rounding_mode
, 0, s
, &float32_params
);
2650 return float32_round_pack_canonical(&p
, s
);
2653 float64
float64_round_to_int(float64 a
, float_status
*s
)
2657 float64_unpack_canonical(&p
, a
, s
);
2658 parts_round_to_int(&p
, s
->float_rounding_mode
, 0, s
, &float64_params
);
2659 return float64_round_pack_canonical(&p
, s
);
2662 bfloat16
bfloat16_round_to_int(bfloat16 a
, float_status
*s
)
2666 bfloat16_unpack_canonical(&p
, a
, s
);
2667 parts_round_to_int(&p
, s
->float_rounding_mode
, 0, s
, &bfloat16_params
);
2668 return bfloat16_round_pack_canonical(&p
, s
);
2671 float128
float128_round_to_int(float128 a
, float_status
*s
)
2675 float128_unpack_canonical(&p
, a
, s
);
2676 parts_round_to_int(&p
, s
->float_rounding_mode
, 0, s
, &float128_params
);
2677 return float128_round_pack_canonical(&p
, s
);
2680 floatx80
floatx80_round_to_int(floatx80 a
, float_status
*status
)
2684 if (!floatx80_unpack_canonical(&p
, a
, status
)) {
2685 return floatx80_default_nan(status
);
2688 parts_round_to_int(&p
, status
->float_rounding_mode
, 0, status
,
2689 &floatx80_params
[status
->floatx80_rounding_precision
]);
2690 return floatx80_round_pack_canonical(&p
, status
);
2694 * Floating-point to signed integer conversions
2697 int8_t float16_to_int8_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2702 float16_unpack_canonical(&p
, a
, s
);
2703 return parts_float_to_sint(&p
, rmode
, scale
, INT8_MIN
, INT8_MAX
, s
);
2706 int16_t float16_to_int16_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2711 float16_unpack_canonical(&p
, a
, s
);
2712 return parts_float_to_sint(&p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2715 int32_t float16_to_int32_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2720 float16_unpack_canonical(&p
, a
, s
);
2721 return parts_float_to_sint(&p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2724 int64_t float16_to_int64_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2729 float16_unpack_canonical(&p
, a
, s
);
2730 return parts_float_to_sint(&p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2733 int16_t float32_to_int16_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2738 float32_unpack_canonical(&p
, a
, s
);
2739 return parts_float_to_sint(&p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2742 int32_t float32_to_int32_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2747 float32_unpack_canonical(&p
, a
, s
);
2748 return parts_float_to_sint(&p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2751 int64_t float32_to_int64_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2756 float32_unpack_canonical(&p
, a
, s
);
2757 return parts_float_to_sint(&p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2760 int16_t float64_to_int16_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2765 float64_unpack_canonical(&p
, a
, s
);
2766 return parts_float_to_sint(&p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2769 int32_t float64_to_int32_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2774 float64_unpack_canonical(&p
, a
, s
);
2775 return parts_float_to_sint(&p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2778 int64_t float64_to_int64_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2783 float64_unpack_canonical(&p
, a
, s
);
2784 return parts_float_to_sint(&p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2787 int16_t bfloat16_to_int16_scalbn(bfloat16 a
, FloatRoundMode rmode
, int scale
,
2792 bfloat16_unpack_canonical(&p
, a
, s
);
2793 return parts_float_to_sint(&p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2796 int32_t bfloat16_to_int32_scalbn(bfloat16 a
, FloatRoundMode rmode
, int scale
,
2801 bfloat16_unpack_canonical(&p
, a
, s
);
2802 return parts_float_to_sint(&p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2805 int64_t bfloat16_to_int64_scalbn(bfloat16 a
, FloatRoundMode rmode
, int scale
,
2810 bfloat16_unpack_canonical(&p
, a
, s
);
2811 return parts_float_to_sint(&p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2814 static int32_t float128_to_int32_scalbn(float128 a
, FloatRoundMode rmode
,
2815 int scale
, float_status
*s
)
2819 float128_unpack_canonical(&p
, a
, s
);
2820 return parts_float_to_sint(&p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2823 static int64_t float128_to_int64_scalbn(float128 a
, FloatRoundMode rmode
,
2824 int scale
, float_status
*s
)
2828 float128_unpack_canonical(&p
, a
, s
);
2829 return parts_float_to_sint(&p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2832 static int32_t floatx80_to_int32_scalbn(floatx80 a
, FloatRoundMode rmode
,
2833 int scale
, float_status
*s
)
2837 if (!floatx80_unpack_canonical(&p
, a
, s
)) {
2838 parts_default_nan(&p
, s
);
2840 return parts_float_to_sint(&p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2843 static int64_t floatx80_to_int64_scalbn(floatx80 a
, FloatRoundMode rmode
,
2844 int scale
, float_status
*s
)
2848 if (!floatx80_unpack_canonical(&p
, a
, s
)) {
2849 parts_default_nan(&p
, s
);
2851 return parts_float_to_sint(&p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2854 int8_t float16_to_int8(float16 a
, float_status
*s
)
2856 return float16_to_int8_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2859 int16_t float16_to_int16(float16 a
, float_status
*s
)
2861 return float16_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2864 int32_t float16_to_int32(float16 a
, float_status
*s
)
2866 return float16_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2869 int64_t float16_to_int64(float16 a
, float_status
*s
)
2871 return float16_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2874 int16_t float32_to_int16(float32 a
, float_status
*s
)
2876 return float32_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2879 int32_t float32_to_int32(float32 a
, float_status
*s
)
2881 return float32_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2884 int64_t float32_to_int64(float32 a
, float_status
*s
)
2886 return float32_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2889 int16_t float64_to_int16(float64 a
, float_status
*s
)
2891 return float64_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2894 int32_t float64_to_int32(float64 a
, float_status
*s
)
2896 return float64_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2899 int64_t float64_to_int64(float64 a
, float_status
*s
)
2901 return float64_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2904 int32_t float128_to_int32(float128 a
, float_status
*s
)
2906 return float128_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2909 int64_t float128_to_int64(float128 a
, float_status
*s
)
2911 return float128_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2914 int32_t floatx80_to_int32(floatx80 a
, float_status
*s
)
2916 return floatx80_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2919 int64_t floatx80_to_int64(floatx80 a
, float_status
*s
)
2921 return floatx80_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2924 int16_t float16_to_int16_round_to_zero(float16 a
, float_status
*s
)
2926 return float16_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2929 int32_t float16_to_int32_round_to_zero(float16 a
, float_status
*s
)
2931 return float16_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2934 int64_t float16_to_int64_round_to_zero(float16 a
, float_status
*s
)
2936 return float16_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2939 int16_t float32_to_int16_round_to_zero(float32 a
, float_status
*s
)
2941 return float32_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2944 int32_t float32_to_int32_round_to_zero(float32 a
, float_status
*s
)
2946 return float32_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2949 int64_t float32_to_int64_round_to_zero(float32 a
, float_status
*s
)
2951 return float32_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2954 int16_t float64_to_int16_round_to_zero(float64 a
, float_status
*s
)
2956 return float64_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2959 int32_t float64_to_int32_round_to_zero(float64 a
, float_status
*s
)
2961 return float64_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2964 int64_t float64_to_int64_round_to_zero(float64 a
, float_status
*s
)
2966 return float64_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2969 int32_t float128_to_int32_round_to_zero(float128 a
, float_status
*s
)
2971 return float128_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2974 int64_t float128_to_int64_round_to_zero(float128 a
, float_status
*s
)
2976 return float128_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2979 int32_t floatx80_to_int32_round_to_zero(floatx80 a
, float_status
*s
)
2981 return floatx80_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2984 int64_t floatx80_to_int64_round_to_zero(floatx80 a
, float_status
*s
)
2986 return floatx80_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2989 int16_t bfloat16_to_int16(bfloat16 a
, float_status
*s
)
2991 return bfloat16_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2994 int32_t bfloat16_to_int32(bfloat16 a
, float_status
*s
)
2996 return bfloat16_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2999 int64_t bfloat16_to_int64(bfloat16 a
, float_status
*s
)
3001 return bfloat16_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3004 int16_t bfloat16_to_int16_round_to_zero(bfloat16 a
, float_status
*s
)
3006 return bfloat16_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
3009 int32_t bfloat16_to_int32_round_to_zero(bfloat16 a
, float_status
*s
)
3011 return bfloat16_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
3014 int64_t bfloat16_to_int64_round_to_zero(bfloat16 a
, float_status
*s
)
3016 return bfloat16_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
3020 * Floating-point to unsigned integer conversions
3023 uint8_t float16_to_uint8_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
3028 float16_unpack_canonical(&p
, a
, s
);
3029 return parts_float_to_uint(&p
, rmode
, scale
, UINT8_MAX
, s
);
3032 uint16_t float16_to_uint16_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
3037 float16_unpack_canonical(&p
, a
, s
);
3038 return parts_float_to_uint(&p
, rmode
, scale
, UINT16_MAX
, s
);
3041 uint32_t float16_to_uint32_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
3046 float16_unpack_canonical(&p
, a
, s
);
3047 return parts_float_to_uint(&p
, rmode
, scale
, UINT32_MAX
, s
);
3050 uint64_t float16_to_uint64_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
3055 float16_unpack_canonical(&p
, a
, s
);
3056 return parts_float_to_uint(&p
, rmode
, scale
, UINT64_MAX
, s
);
3059 uint16_t float32_to_uint16_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
3064 float32_unpack_canonical(&p
, a
, s
);
3065 return parts_float_to_uint(&p
, rmode
, scale
, UINT16_MAX
, s
);
3068 uint32_t float32_to_uint32_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
3073 float32_unpack_canonical(&p
, a
, s
);
3074 return parts_float_to_uint(&p
, rmode
, scale
, UINT32_MAX
, s
);
3077 uint64_t float32_to_uint64_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
3082 float32_unpack_canonical(&p
, a
, s
);
3083 return parts_float_to_uint(&p
, rmode
, scale
, UINT64_MAX
, s
);
3086 uint16_t float64_to_uint16_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
3091 float64_unpack_canonical(&p
, a
, s
);
3092 return parts_float_to_uint(&p
, rmode
, scale
, UINT16_MAX
, s
);
3095 uint32_t float64_to_uint32_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
3100 float64_unpack_canonical(&p
, a
, s
);
3101 return parts_float_to_uint(&p
, rmode
, scale
, UINT32_MAX
, s
);
3104 uint64_t float64_to_uint64_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
3109 float64_unpack_canonical(&p
, a
, s
);
3110 return parts_float_to_uint(&p
, rmode
, scale
, UINT64_MAX
, s
);
3113 uint16_t bfloat16_to_uint16_scalbn(bfloat16 a
, FloatRoundMode rmode
,
3114 int scale
, float_status
*s
)
3118 bfloat16_unpack_canonical(&p
, a
, s
);
3119 return parts_float_to_uint(&p
, rmode
, scale
, UINT16_MAX
, s
);
3122 uint32_t bfloat16_to_uint32_scalbn(bfloat16 a
, FloatRoundMode rmode
,
3123 int scale
, float_status
*s
)
3127 bfloat16_unpack_canonical(&p
, a
, s
);
3128 return parts_float_to_uint(&p
, rmode
, scale
, UINT32_MAX
, s
);
3131 uint64_t bfloat16_to_uint64_scalbn(bfloat16 a
, FloatRoundMode rmode
,
3132 int scale
, float_status
*s
)
3136 bfloat16_unpack_canonical(&p
, a
, s
);
3137 return parts_float_to_uint(&p
, rmode
, scale
, UINT64_MAX
, s
);
3140 static uint32_t float128_to_uint32_scalbn(float128 a
, FloatRoundMode rmode
,
3141 int scale
, float_status
*s
)
3145 float128_unpack_canonical(&p
, a
, s
);
3146 return parts_float_to_uint(&p
, rmode
, scale
, UINT32_MAX
, s
);
3149 static uint64_t float128_to_uint64_scalbn(float128 a
, FloatRoundMode rmode
,
3150 int scale
, float_status
*s
)
3154 float128_unpack_canonical(&p
, a
, s
);
3155 return parts_float_to_uint(&p
, rmode
, scale
, UINT64_MAX
, s
);
3158 uint8_t float16_to_uint8(float16 a
, float_status
*s
)
3160 return float16_to_uint8_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3163 uint16_t float16_to_uint16(float16 a
, float_status
*s
)
3165 return float16_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3168 uint32_t float16_to_uint32(float16 a
, float_status
*s
)
3170 return float16_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3173 uint64_t float16_to_uint64(float16 a
, float_status
*s
)
3175 return float16_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3178 uint16_t float32_to_uint16(float32 a
, float_status
*s
)
3180 return float32_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3183 uint32_t float32_to_uint32(float32 a
, float_status
*s
)
3185 return float32_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3188 uint64_t float32_to_uint64(float32 a
, float_status
*s
)
3190 return float32_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3193 uint16_t float64_to_uint16(float64 a
, float_status
*s
)
3195 return float64_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3198 uint32_t float64_to_uint32(float64 a
, float_status
*s
)
3200 return float64_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3203 uint64_t float64_to_uint64(float64 a
, float_status
*s
)
3205 return float64_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3208 uint32_t float128_to_uint32(float128 a
, float_status
*s
)
3210 return float128_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3213 uint64_t float128_to_uint64(float128 a
, float_status
*s
)
3215 return float128_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3218 uint16_t float16_to_uint16_round_to_zero(float16 a
, float_status
*s
)
3220 return float16_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
3223 uint32_t float16_to_uint32_round_to_zero(float16 a
, float_status
*s
)
3225 return float16_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
3228 uint64_t float16_to_uint64_round_to_zero(float16 a
, float_status
*s
)
3230 return float16_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
3233 uint16_t float32_to_uint16_round_to_zero(float32 a
, float_status
*s
)
3235 return float32_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
3238 uint32_t float32_to_uint32_round_to_zero(float32 a
, float_status
*s
)
3240 return float32_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
3243 uint64_t float32_to_uint64_round_to_zero(float32 a
, float_status
*s
)
3245 return float32_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
3248 uint16_t float64_to_uint16_round_to_zero(float64 a
, float_status
*s
)
3250 return float64_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
3253 uint32_t float64_to_uint32_round_to_zero(float64 a
, float_status
*s
)
3255 return float64_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
3258 uint64_t float64_to_uint64_round_to_zero(float64 a
, float_status
*s
)
3260 return float64_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
3263 uint32_t float128_to_uint32_round_to_zero(float128 a
, float_status
*s
)
3265 return float128_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
3268 uint64_t float128_to_uint64_round_to_zero(float128 a
, float_status
*s
)
3270 return float128_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
3273 uint16_t bfloat16_to_uint16(bfloat16 a
, float_status
*s
)
3275 return bfloat16_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3278 uint32_t bfloat16_to_uint32(bfloat16 a
, float_status
*s
)
3280 return bfloat16_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3283 uint64_t bfloat16_to_uint64(bfloat16 a
, float_status
*s
)
3285 return bfloat16_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3288 uint16_t bfloat16_to_uint16_round_to_zero(bfloat16 a
, float_status
*s
)
3290 return bfloat16_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
3293 uint32_t bfloat16_to_uint32_round_to_zero(bfloat16 a
, float_status
*s
)
3295 return bfloat16_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
3298 uint64_t bfloat16_to_uint64_round_to_zero(bfloat16 a
, float_status
*s
)
3300 return bfloat16_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
3304 * Signed integer to floating-point conversions
3307 float16
int64_to_float16_scalbn(int64_t a
, int scale
, float_status
*status
)
3311 parts_sint_to_float(&p
, a
, scale
, status
);
3312 return float16_round_pack_canonical(&p
, status
);
3315 float16
int32_to_float16_scalbn(int32_t a
, int scale
, float_status
*status
)
3317 return int64_to_float16_scalbn(a
, scale
, status
);
3320 float16
int16_to_float16_scalbn(int16_t a
, int scale
, float_status
*status
)
3322 return int64_to_float16_scalbn(a
, scale
, status
);
3325 float16
int64_to_float16(int64_t a
, float_status
*status
)
3327 return int64_to_float16_scalbn(a
, 0, status
);
3330 float16
int32_to_float16(int32_t a
, float_status
*status
)
3332 return int64_to_float16_scalbn(a
, 0, status
);
3335 float16
int16_to_float16(int16_t a
, float_status
*status
)
3337 return int64_to_float16_scalbn(a
, 0, status
);
3340 float16
int8_to_float16(int8_t a
, float_status
*status
)
3342 return int64_to_float16_scalbn(a
, 0, status
);
3345 float32
int64_to_float32_scalbn(int64_t a
, int scale
, float_status
*status
)
3349 parts64_sint_to_float(&p
, a
, scale
, status
);
3350 return float32_round_pack_canonical(&p
, status
);
3353 float32
int32_to_float32_scalbn(int32_t a
, int scale
, float_status
*status
)
3355 return int64_to_float32_scalbn(a
, scale
, status
);
3358 float32
int16_to_float32_scalbn(int16_t a
, int scale
, float_status
*status
)
3360 return int64_to_float32_scalbn(a
, scale
, status
);
3363 float32
int64_to_float32(int64_t a
, float_status
*status
)
3365 return int64_to_float32_scalbn(a
, 0, status
);
3368 float32
int32_to_float32(int32_t a
, float_status
*status
)
3370 return int64_to_float32_scalbn(a
, 0, status
);
3373 float32
int16_to_float32(int16_t a
, float_status
*status
)
3375 return int64_to_float32_scalbn(a
, 0, status
);
3378 float64
int64_to_float64_scalbn(int64_t a
, int scale
, float_status
*status
)
3382 parts_sint_to_float(&p
, a
, scale
, status
);
3383 return float64_round_pack_canonical(&p
, status
);
3386 float64
int32_to_float64_scalbn(int32_t a
, int scale
, float_status
*status
)
3388 return int64_to_float64_scalbn(a
, scale
, status
);
3391 float64
int16_to_float64_scalbn(int16_t a
, int scale
, float_status
*status
)
3393 return int64_to_float64_scalbn(a
, scale
, status
);
3396 float64
int64_to_float64(int64_t a
, float_status
*status
)
3398 return int64_to_float64_scalbn(a
, 0, status
);
3401 float64
int32_to_float64(int32_t a
, float_status
*status
)
3403 return int64_to_float64_scalbn(a
, 0, status
);
3406 float64
int16_to_float64(int16_t a
, float_status
*status
)
3408 return int64_to_float64_scalbn(a
, 0, status
);
3411 bfloat16
int64_to_bfloat16_scalbn(int64_t a
, int scale
, float_status
*status
)
3415 parts_sint_to_float(&p
, a
, scale
, status
);
3416 return bfloat16_round_pack_canonical(&p
, status
);
3419 bfloat16
int32_to_bfloat16_scalbn(int32_t a
, int scale
, float_status
*status
)
3421 return int64_to_bfloat16_scalbn(a
, scale
, status
);
3424 bfloat16
int16_to_bfloat16_scalbn(int16_t a
, int scale
, float_status
*status
)
3426 return int64_to_bfloat16_scalbn(a
, scale
, status
);
3429 bfloat16
int64_to_bfloat16(int64_t a
, float_status
*status
)
3431 return int64_to_bfloat16_scalbn(a
, 0, status
);
3434 bfloat16
int32_to_bfloat16(int32_t a
, float_status
*status
)
3436 return int64_to_bfloat16_scalbn(a
, 0, status
);
3439 bfloat16
int16_to_bfloat16(int16_t a
, float_status
*status
)
3441 return int64_to_bfloat16_scalbn(a
, 0, status
);
3444 float128
int64_to_float128(int64_t a
, float_status
*status
)
3448 parts_sint_to_float(&p
, a
, 0, status
);
3449 return float128_round_pack_canonical(&p
, status
);
3452 float128
int32_to_float128(int32_t a
, float_status
*status
)
3454 return int64_to_float128(a
, status
);
3457 floatx80
int64_to_floatx80(int64_t a
, float_status
*status
)
3461 parts_sint_to_float(&p
, a
, 0, status
);
3462 return floatx80_round_pack_canonical(&p
, status
);
3465 floatx80
int32_to_floatx80(int32_t a
, float_status
*status
)
3467 return int64_to_floatx80(a
, status
);
3471 * Unsigned Integer to floating-point conversions
3474 float16
uint64_to_float16_scalbn(uint64_t a
, int scale
, float_status
*status
)
3478 parts_uint_to_float(&p
, a
, scale
, status
);
3479 return float16_round_pack_canonical(&p
, status
);
3482 float16
uint32_to_float16_scalbn(uint32_t a
, int scale
, float_status
*status
)
3484 return uint64_to_float16_scalbn(a
, scale
, status
);
3487 float16
uint16_to_float16_scalbn(uint16_t a
, int scale
, float_status
*status
)
3489 return uint64_to_float16_scalbn(a
, scale
, status
);
3492 float16
uint64_to_float16(uint64_t a
, float_status
*status
)
3494 return uint64_to_float16_scalbn(a
, 0, status
);
3497 float16
uint32_to_float16(uint32_t a
, float_status
*status
)
3499 return uint64_to_float16_scalbn(a
, 0, status
);
3502 float16
uint16_to_float16(uint16_t a
, float_status
*status
)
3504 return uint64_to_float16_scalbn(a
, 0, status
);
3507 float16
uint8_to_float16(uint8_t a
, float_status
*status
)
3509 return uint64_to_float16_scalbn(a
, 0, status
);
3512 float32
uint64_to_float32_scalbn(uint64_t a
, int scale
, float_status
*status
)
3516 parts_uint_to_float(&p
, a
, scale
, status
);
3517 return float32_round_pack_canonical(&p
, status
);
3520 float32
uint32_to_float32_scalbn(uint32_t a
, int scale
, float_status
*status
)
3522 return uint64_to_float32_scalbn(a
, scale
, status
);
3525 float32
uint16_to_float32_scalbn(uint16_t a
, int scale
, float_status
*status
)
3527 return uint64_to_float32_scalbn(a
, scale
, status
);
3530 float32
uint64_to_float32(uint64_t a
, float_status
*status
)
3532 return uint64_to_float32_scalbn(a
, 0, status
);
3535 float32
uint32_to_float32(uint32_t a
, float_status
*status
)
3537 return uint64_to_float32_scalbn(a
, 0, status
);
3540 float32
uint16_to_float32(uint16_t a
, float_status
*status
)
3542 return uint64_to_float32_scalbn(a
, 0, status
);
3545 float64
uint64_to_float64_scalbn(uint64_t a
, int scale
, float_status
*status
)
3549 parts_uint_to_float(&p
, a
, scale
, status
);
3550 return float64_round_pack_canonical(&p
, status
);
3553 float64
uint32_to_float64_scalbn(uint32_t a
, int scale
, float_status
*status
)
3555 return uint64_to_float64_scalbn(a
, scale
, status
);
3558 float64
uint16_to_float64_scalbn(uint16_t a
, int scale
, float_status
*status
)
3560 return uint64_to_float64_scalbn(a
, scale
, status
);
3563 float64
uint64_to_float64(uint64_t a
, float_status
*status
)
3565 return uint64_to_float64_scalbn(a
, 0, status
);
3568 float64
uint32_to_float64(uint32_t a
, float_status
*status
)
3570 return uint64_to_float64_scalbn(a
, 0, status
);
3573 float64
uint16_to_float64(uint16_t a
, float_status
*status
)
3575 return uint64_to_float64_scalbn(a
, 0, status
);
3578 bfloat16
uint64_to_bfloat16_scalbn(uint64_t a
, int scale
, float_status
*status
)
3582 parts_uint_to_float(&p
, a
, scale
, status
);
3583 return bfloat16_round_pack_canonical(&p
, status
);
3586 bfloat16
uint32_to_bfloat16_scalbn(uint32_t a
, int scale
, float_status
*status
)
3588 return uint64_to_bfloat16_scalbn(a
, scale
, status
);
3591 bfloat16
uint16_to_bfloat16_scalbn(uint16_t a
, int scale
, float_status
*status
)
3593 return uint64_to_bfloat16_scalbn(a
, scale
, status
);
3596 bfloat16
uint64_to_bfloat16(uint64_t a
, float_status
*status
)
3598 return uint64_to_bfloat16_scalbn(a
, 0, status
);
3601 bfloat16
uint32_to_bfloat16(uint32_t a
, float_status
*status
)
3603 return uint64_to_bfloat16_scalbn(a
, 0, status
);
3606 bfloat16
uint16_to_bfloat16(uint16_t a
, float_status
*status
)
3608 return uint64_to_bfloat16_scalbn(a
, 0, status
);
3611 float128
uint64_to_float128(uint64_t a
, float_status
*status
)
3615 parts_uint_to_float(&p
, a
, 0, status
);
3616 return float128_round_pack_canonical(&p
, status
);
3620 * Minimum and maximum
3623 static float16
float16_minmax(float16 a
, float16 b
, float_status
*s
, int flags
)
3625 FloatParts64 pa
, pb
, *pr
;
3627 float16_unpack_canonical(&pa
, a
, s
);
3628 float16_unpack_canonical(&pb
, b
, s
);
3629 pr
= parts_minmax(&pa
, &pb
, s
, flags
);
3631 return float16_round_pack_canonical(pr
, s
);
3634 static bfloat16
bfloat16_minmax(bfloat16 a
, bfloat16 b
,
3635 float_status
*s
, int flags
)
3637 FloatParts64 pa
, pb
, *pr
;
3639 bfloat16_unpack_canonical(&pa
, a
, s
);
3640 bfloat16_unpack_canonical(&pb
, b
, s
);
3641 pr
= parts_minmax(&pa
, &pb
, s
, flags
);
3643 return bfloat16_round_pack_canonical(pr
, s
);
3646 static float32
float32_minmax(float32 a
, float32 b
, float_status
*s
, int flags
)
3648 FloatParts64 pa
, pb
, *pr
;
3650 float32_unpack_canonical(&pa
, a
, s
);
3651 float32_unpack_canonical(&pb
, b
, s
);
3652 pr
= parts_minmax(&pa
, &pb
, s
, flags
);
3654 return float32_round_pack_canonical(pr
, s
);
3657 static float64
float64_minmax(float64 a
, float64 b
, float_status
*s
, int flags
)
3659 FloatParts64 pa
, pb
, *pr
;
3661 float64_unpack_canonical(&pa
, a
, s
);
3662 float64_unpack_canonical(&pb
, b
, s
);
3663 pr
= parts_minmax(&pa
, &pb
, s
, flags
);
3665 return float64_round_pack_canonical(pr
, s
);
3668 static float128
float128_minmax(float128 a
, float128 b
,
3669 float_status
*s
, int flags
)
3671 FloatParts128 pa
, pb
, *pr
;
3673 float128_unpack_canonical(&pa
, a
, s
);
3674 float128_unpack_canonical(&pb
, b
, s
);
3675 pr
= parts_minmax(&pa
, &pb
, s
, flags
);
3677 return float128_round_pack_canonical(pr
, s
);
3680 #define MINMAX_1(type, name, flags) \
3681 type type##_##name(type a, type b, float_status *s) \
3682 { return type##_minmax(a, b, s, flags); }
3684 #define MINMAX_2(type) \
3685 MINMAX_1(type, max, 0) \
3686 MINMAX_1(type, maxnum, minmax_isnum) \
3687 MINMAX_1(type, maxnummag, minmax_isnum | minmax_ismag) \
3688 MINMAX_1(type, min, minmax_ismin) \
3689 MINMAX_1(type, minnum, minmax_ismin | minmax_isnum) \
3690 MINMAX_1(type, minnummag, minmax_ismin | minmax_isnum | minmax_ismag)
3702 * Floating point compare
3705 static FloatRelation QEMU_FLATTEN
3706 float16_do_compare(float16 a
, float16 b
, float_status
*s
, bool is_quiet
)
3708 FloatParts64 pa
, pb
;
3710 float16_unpack_canonical(&pa
, a
, s
);
3711 float16_unpack_canonical(&pb
, b
, s
);
3712 return parts_compare(&pa
, &pb
, s
, is_quiet
);
3715 FloatRelation
float16_compare(float16 a
, float16 b
, float_status
*s
)
3717 return float16_do_compare(a
, b
, s
, false);
3720 FloatRelation
float16_compare_quiet(float16 a
, float16 b
, float_status
*s
)
3722 return float16_do_compare(a
, b
, s
, true);
3725 static FloatRelation QEMU_SOFTFLOAT_ATTR
3726 float32_do_compare(float32 a
, float32 b
, float_status
*s
, bool is_quiet
)
3728 FloatParts64 pa
, pb
;
3730 float32_unpack_canonical(&pa
, a
, s
);
3731 float32_unpack_canonical(&pb
, b
, s
);
3732 return parts_compare(&pa
, &pb
, s
, is_quiet
);
3735 static FloatRelation QEMU_FLATTEN
3736 float32_hs_compare(float32 xa
, float32 xb
, float_status
*s
, bool is_quiet
)
3738 union_float32 ua
, ub
;
3743 if (QEMU_NO_HARDFLOAT
) {
3747 float32_input_flush2(&ua
.s
, &ub
.s
, s
);
3748 if (isgreaterequal(ua
.h
, ub
.h
)) {
3749 if (isgreater(ua
.h
, ub
.h
)) {
3750 return float_relation_greater
;
3752 return float_relation_equal
;
3754 if (likely(isless(ua
.h
, ub
.h
))) {
3755 return float_relation_less
;
3758 * The only condition remaining is unordered.
3759 * Fall through to set flags.
3762 return float32_do_compare(ua
.s
, ub
.s
, s
, is_quiet
);
3765 FloatRelation
float32_compare(float32 a
, float32 b
, float_status
*s
)
3767 return float32_hs_compare(a
, b
, s
, false);
3770 FloatRelation
float32_compare_quiet(float32 a
, float32 b
, float_status
*s
)
3772 return float32_hs_compare(a
, b
, s
, true);
3775 static FloatRelation QEMU_SOFTFLOAT_ATTR
3776 float64_do_compare(float64 a
, float64 b
, float_status
*s
, bool is_quiet
)
3778 FloatParts64 pa
, pb
;
3780 float64_unpack_canonical(&pa
, a
, s
);
3781 float64_unpack_canonical(&pb
, b
, s
);
3782 return parts_compare(&pa
, &pb
, s
, is_quiet
);
3785 static FloatRelation QEMU_FLATTEN
3786 float64_hs_compare(float64 xa
, float64 xb
, float_status
*s
, bool is_quiet
)
3788 union_float64 ua
, ub
;
3793 if (QEMU_NO_HARDFLOAT
) {
3797 float64_input_flush2(&ua
.s
, &ub
.s
, s
);
3798 if (isgreaterequal(ua
.h
, ub
.h
)) {
3799 if (isgreater(ua
.h
, ub
.h
)) {
3800 return float_relation_greater
;
3802 return float_relation_equal
;
3804 if (likely(isless(ua
.h
, ub
.h
))) {
3805 return float_relation_less
;
3808 * The only condition remaining is unordered.
3809 * Fall through to set flags.
3812 return float64_do_compare(ua
.s
, ub
.s
, s
, is_quiet
);
3815 FloatRelation
float64_compare(float64 a
, float64 b
, float_status
*s
)
3817 return float64_hs_compare(a
, b
, s
, false);
3820 FloatRelation
float64_compare_quiet(float64 a
, float64 b
, float_status
*s
)
3822 return float64_hs_compare(a
, b
, s
, true);
3825 static FloatRelation QEMU_FLATTEN
3826 bfloat16_do_compare(bfloat16 a
, bfloat16 b
, float_status
*s
, bool is_quiet
)
3828 FloatParts64 pa
, pb
;
3830 bfloat16_unpack_canonical(&pa
, a
, s
);
3831 bfloat16_unpack_canonical(&pb
, b
, s
);
3832 return parts_compare(&pa
, &pb
, s
, is_quiet
);
3835 FloatRelation
bfloat16_compare(bfloat16 a
, bfloat16 b
, float_status
*s
)
3837 return bfloat16_do_compare(a
, b
, s
, false);
3840 FloatRelation
bfloat16_compare_quiet(bfloat16 a
, bfloat16 b
, float_status
*s
)
3842 return bfloat16_do_compare(a
, b
, s
, true);
3845 static FloatRelation QEMU_FLATTEN
3846 float128_do_compare(float128 a
, float128 b
, float_status
*s
, bool is_quiet
)
3848 FloatParts128 pa
, pb
;
3850 float128_unpack_canonical(&pa
, a
, s
);
3851 float128_unpack_canonical(&pb
, b
, s
);
3852 return parts_compare(&pa
, &pb
, s
, is_quiet
);
3855 FloatRelation
float128_compare(float128 a
, float128 b
, float_status
*s
)
3857 return float128_do_compare(a
, b
, s
, false);
3860 FloatRelation
float128_compare_quiet(float128 a
, float128 b
, float_status
*s
)
3862 return float128_do_compare(a
, b
, s
, true);
3865 static FloatRelation QEMU_FLATTEN
3866 floatx80_do_compare(floatx80 a
, floatx80 b
, float_status
*s
, bool is_quiet
)
3868 FloatParts128 pa
, pb
;
3870 if (!floatx80_unpack_canonical(&pa
, a
, s
) ||
3871 !floatx80_unpack_canonical(&pb
, b
, s
)) {
3872 return float_relation_unordered
;
3874 return parts_compare(&pa
, &pb
, s
, is_quiet
);
3877 FloatRelation
floatx80_compare(floatx80 a
, floatx80 b
, float_status
*s
)
3879 return floatx80_do_compare(a
, b
, s
, false);
3882 FloatRelation
floatx80_compare_quiet(floatx80 a
, floatx80 b
, float_status
*s
)
3884 return floatx80_do_compare(a
, b
, s
, true);
3891 float16
float16_scalbn(float16 a
, int n
, float_status
*status
)
3895 float16_unpack_canonical(&p
, a
, status
);
3896 parts_scalbn(&p
, n
, status
);
3897 return float16_round_pack_canonical(&p
, status
);
3900 float32
float32_scalbn(float32 a
, int n
, float_status
*status
)
3904 float32_unpack_canonical(&p
, a
, status
);
3905 parts_scalbn(&p
, n
, status
);
3906 return float32_round_pack_canonical(&p
, status
);
3909 float64
float64_scalbn(float64 a
, int n
, float_status
*status
)
3913 float64_unpack_canonical(&p
, a
, status
);
3914 parts_scalbn(&p
, n
, status
);
3915 return float64_round_pack_canonical(&p
, status
);
3918 bfloat16
bfloat16_scalbn(bfloat16 a
, int n
, float_status
*status
)
3922 bfloat16_unpack_canonical(&p
, a
, status
);
3923 parts_scalbn(&p
, n
, status
);
3924 return bfloat16_round_pack_canonical(&p
, status
);
3927 float128
float128_scalbn(float128 a
, int n
, float_status
*status
)
3931 float128_unpack_canonical(&p
, a
, status
);
3932 parts_scalbn(&p
, n
, status
);
3933 return float128_round_pack_canonical(&p
, status
);
3936 floatx80
floatx80_scalbn(floatx80 a
, int n
, float_status
*status
)
3940 if (!floatx80_unpack_canonical(&p
, a
, status
)) {
3941 return floatx80_default_nan(status
);
3943 parts_scalbn(&p
, n
, status
);
3944 return floatx80_round_pack_canonical(&p
, status
);
3951 float16 QEMU_FLATTEN
float16_sqrt(float16 a
, float_status
*status
)
3955 float16_unpack_canonical(&p
, a
, status
);
3956 parts_sqrt(&p
, status
, &float16_params
);
3957 return float16_round_pack_canonical(&p
, status
);
3960 static float32 QEMU_SOFTFLOAT_ATTR
3961 soft_f32_sqrt(float32 a
, float_status
*status
)
3965 float32_unpack_canonical(&p
, a
, status
);
3966 parts_sqrt(&p
, status
, &float32_params
);
3967 return float32_round_pack_canonical(&p
, status
);
3970 static float64 QEMU_SOFTFLOAT_ATTR
3971 soft_f64_sqrt(float64 a
, float_status
*status
)
3975 float64_unpack_canonical(&p
, a
, status
);
3976 parts_sqrt(&p
, status
, &float64_params
);
3977 return float64_round_pack_canonical(&p
, status
);
3980 float32 QEMU_FLATTEN
float32_sqrt(float32 xa
, float_status
*s
)
3982 union_float32 ua
, ur
;
3985 if (unlikely(!can_use_fpu(s
))) {
3989 float32_input_flush1(&ua
.s
, s
);
3990 if (QEMU_HARDFLOAT_1F32_USE_FP
) {
3991 if (unlikely(!(fpclassify(ua
.h
) == FP_NORMAL
||
3992 fpclassify(ua
.h
) == FP_ZERO
) ||
3996 } else if (unlikely(!float32_is_zero_or_normal(ua
.s
) ||
3997 float32_is_neg(ua
.s
))) {
4004 return soft_f32_sqrt(ua
.s
, s
);
4007 float64 QEMU_FLATTEN
float64_sqrt(float64 xa
, float_status
*s
)
4009 union_float64 ua
, ur
;
4012 if (unlikely(!can_use_fpu(s
))) {
4016 float64_input_flush1(&ua
.s
, s
);
4017 if (QEMU_HARDFLOAT_1F64_USE_FP
) {
4018 if (unlikely(!(fpclassify(ua
.h
) == FP_NORMAL
||
4019 fpclassify(ua
.h
) == FP_ZERO
) ||
4023 } else if (unlikely(!float64_is_zero_or_normal(ua
.s
) ||
4024 float64_is_neg(ua
.s
))) {
4031 return soft_f64_sqrt(ua
.s
, s
);
4034 bfloat16 QEMU_FLATTEN
bfloat16_sqrt(bfloat16 a
, float_status
*status
)
4038 bfloat16_unpack_canonical(&p
, a
, status
);
4039 parts_sqrt(&p
, status
, &bfloat16_params
);
4040 return bfloat16_round_pack_canonical(&p
, status
);
4043 float128 QEMU_FLATTEN
float128_sqrt(float128 a
, float_status
*status
)
4047 float128_unpack_canonical(&p
, a
, status
);
4048 parts_sqrt(&p
, status
, &float128_params
);
4049 return float128_round_pack_canonical(&p
, status
);
4052 floatx80
floatx80_sqrt(floatx80 a
, float_status
*s
)
4056 if (!floatx80_unpack_canonical(&p
, a
, s
)) {
4057 return floatx80_default_nan(s
);
4059 parts_sqrt(&p
, s
, &floatx80_params
[s
->floatx80_rounding_precision
]);
4060 return floatx80_round_pack_canonical(&p
, s
);
4063 /*----------------------------------------------------------------------------
4064 | The pattern for a default generated NaN.
4065 *----------------------------------------------------------------------------*/
4067 float16
float16_default_nan(float_status
*status
)
4071 parts_default_nan(&p
, status
);
4072 p
.frac
>>= float16_params
.frac_shift
;
4073 return float16_pack_raw(&p
);
4076 float32
float32_default_nan(float_status
*status
)
4080 parts_default_nan(&p
, status
);
4081 p
.frac
>>= float32_params
.frac_shift
;
4082 return float32_pack_raw(&p
);
4085 float64
float64_default_nan(float_status
*status
)
4089 parts_default_nan(&p
, status
);
4090 p
.frac
>>= float64_params
.frac_shift
;
4091 return float64_pack_raw(&p
);
4094 float128
float128_default_nan(float_status
*status
)
4098 parts_default_nan(&p
, status
);
4099 frac_shr(&p
, float128_params
.frac_shift
);
4100 return float128_pack_raw(&p
);
4103 bfloat16
bfloat16_default_nan(float_status
*status
)
4107 parts_default_nan(&p
, status
);
4108 p
.frac
>>= bfloat16_params
.frac_shift
;
4109 return bfloat16_pack_raw(&p
);
4112 /*----------------------------------------------------------------------------
4113 | Returns a quiet NaN from a signalling NaN for the floating point value `a'.
4114 *----------------------------------------------------------------------------*/
4116 float16
float16_silence_nan(float16 a
, float_status
*status
)
4120 float16_unpack_raw(&p
, a
);
4121 p
.frac
<<= float16_params
.frac_shift
;
4122 parts_silence_nan(&p
, status
);
4123 p
.frac
>>= float16_params
.frac_shift
;
4124 return float16_pack_raw(&p
);
4127 float32
float32_silence_nan(float32 a
, float_status
*status
)
4131 float32_unpack_raw(&p
, a
);
4132 p
.frac
<<= float32_params
.frac_shift
;
4133 parts_silence_nan(&p
, status
);
4134 p
.frac
>>= float32_params
.frac_shift
;
4135 return float32_pack_raw(&p
);
4138 float64
float64_silence_nan(float64 a
, float_status
*status
)
4142 float64_unpack_raw(&p
, a
);
4143 p
.frac
<<= float64_params
.frac_shift
;
4144 parts_silence_nan(&p
, status
);
4145 p
.frac
>>= float64_params
.frac_shift
;
4146 return float64_pack_raw(&p
);
4149 bfloat16
bfloat16_silence_nan(bfloat16 a
, float_status
*status
)
4153 bfloat16_unpack_raw(&p
, a
);
4154 p
.frac
<<= bfloat16_params
.frac_shift
;
4155 parts_silence_nan(&p
, status
);
4156 p
.frac
>>= bfloat16_params
.frac_shift
;
4157 return bfloat16_pack_raw(&p
);
4160 float128
float128_silence_nan(float128 a
, float_status
*status
)
4164 float128_unpack_raw(&p
, a
);
4165 frac_shl(&p
, float128_params
.frac_shift
);
4166 parts_silence_nan(&p
, status
);
4167 frac_shr(&p
, float128_params
.frac_shift
);
4168 return float128_pack_raw(&p
);
4171 /*----------------------------------------------------------------------------
4172 | If `a' is denormal and we are in flush-to-zero mode then set the
4173 | input-denormal exception and return zero. Otherwise just return the value.
4174 *----------------------------------------------------------------------------*/
4176 static bool parts_squash_denormal(FloatParts64 p
, float_status
*status
)
4178 if (p
.exp
== 0 && p
.frac
!= 0) {
4179 float_raise(float_flag_input_denormal
, status
);
4186 float16
float16_squash_input_denormal(float16 a
, float_status
*status
)
4188 if (status
->flush_inputs_to_zero
) {
4191 float16_unpack_raw(&p
, a
);
4192 if (parts_squash_denormal(p
, status
)) {
4193 return float16_set_sign(float16_zero
, p
.sign
);
4199 float32
float32_squash_input_denormal(float32 a
, float_status
*status
)
4201 if (status
->flush_inputs_to_zero
) {
4204 float32_unpack_raw(&p
, a
);
4205 if (parts_squash_denormal(p
, status
)) {
4206 return float32_set_sign(float32_zero
, p
.sign
);
4212 float64
float64_squash_input_denormal(float64 a
, float_status
*status
)
4214 if (status
->flush_inputs_to_zero
) {
4217 float64_unpack_raw(&p
, a
);
4218 if (parts_squash_denormal(p
, status
)) {
4219 return float64_set_sign(float64_zero
, p
.sign
);
4225 bfloat16
bfloat16_squash_input_denormal(bfloat16 a
, float_status
*status
)
4227 if (status
->flush_inputs_to_zero
) {
4230 bfloat16_unpack_raw(&p
, a
);
4231 if (parts_squash_denormal(p
, status
)) {
4232 return bfloat16_set_sign(bfloat16_zero
, p
.sign
);
4238 /*----------------------------------------------------------------------------
4239 | Normalizes the subnormal single-precision floating-point value represented
4240 | by the denormalized significand `aSig'. The normalized exponent and
4241 | significand are stored at the locations pointed to by `zExpPtr' and
4242 | `zSigPtr', respectively.
4243 *----------------------------------------------------------------------------*/
4246 normalizeFloat32Subnormal(uint32_t aSig
, int *zExpPtr
, uint32_t *zSigPtr
)
4250 shiftCount
= clz32(aSig
) - 8;
4251 *zSigPtr
= aSig
<<shiftCount
;
4252 *zExpPtr
= 1 - shiftCount
;
4256 /*----------------------------------------------------------------------------
4257 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4258 | and significand `zSig', and returns the proper single-precision floating-
4259 | point value corresponding to the abstract input. Ordinarily, the abstract
4260 | value is simply rounded and packed into the single-precision format, with
4261 | the inexact exception raised if the abstract input cannot be represented
4262 | exactly. However, if the abstract value is too large, the overflow and
4263 | inexact exceptions are raised and an infinity or maximal finite value is
4264 | returned. If the abstract value is too small, the input value is rounded to
4265 | a subnormal number, and the underflow and inexact exceptions are raised if
4266 | the abstract input cannot be represented exactly as a subnormal single-
4267 | precision floating-point number.
4268 | The input significand `zSig' has its binary point between bits 30
4269 | and 29, which is 7 bits to the left of the usual location. This shifted
4270 | significand must be normalized or smaller. If `zSig' is not normalized,
4271 | `zExp' must be 0; in that case, the result returned is a subnormal number,
4272 | and it must not require rounding. In the usual case that `zSig' is
4273 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
4274 | The handling of underflow and overflow follows the IEC/IEEE Standard for
4275 | Binary Floating-Point Arithmetic.
4276 *----------------------------------------------------------------------------*/
4278 static float32
roundAndPackFloat32(bool zSign
, int zExp
, uint32_t zSig
,
4279 float_status
*status
)
4281 int8_t roundingMode
;
4282 bool roundNearestEven
;
4283 int8_t roundIncrement
, roundBits
;
4286 roundingMode
= status
->float_rounding_mode
;
4287 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4288 switch (roundingMode
) {
4289 case float_round_nearest_even
:
4290 case float_round_ties_away
:
4291 roundIncrement
= 0x40;
4293 case float_round_to_zero
:
4296 case float_round_up
:
4297 roundIncrement
= zSign
? 0 : 0x7f;
4299 case float_round_down
:
4300 roundIncrement
= zSign
? 0x7f : 0;
4302 case float_round_to_odd
:
4303 roundIncrement
= zSig
& 0x80 ? 0 : 0x7f;
4309 roundBits
= zSig
& 0x7F;
4310 if ( 0xFD <= (uint16_t) zExp
) {
4311 if ( ( 0xFD < zExp
)
4312 || ( ( zExp
== 0xFD )
4313 && ( (int32_t) ( zSig
+ roundIncrement
) < 0 ) )
4315 bool overflow_to_inf
= roundingMode
!= float_round_to_odd
&&
4316 roundIncrement
!= 0;
4317 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4318 return packFloat32(zSign
, 0xFF, -!overflow_to_inf
);
4321 if (status
->flush_to_zero
) {
4322 float_raise(float_flag_output_denormal
, status
);
4323 return packFloat32(zSign
, 0, 0);
4325 isTiny
= status
->tininess_before_rounding
4327 || (zSig
+ roundIncrement
< 0x80000000);
4328 shift32RightJamming( zSig
, - zExp
, &zSig
);
4330 roundBits
= zSig
& 0x7F;
4331 if (isTiny
&& roundBits
) {
4332 float_raise(float_flag_underflow
, status
);
4334 if (roundingMode
== float_round_to_odd
) {
4336 * For round-to-odd case, the roundIncrement depends on
4337 * zSig which just changed.
4339 roundIncrement
= zSig
& 0x80 ? 0 : 0x7f;
4344 float_raise(float_flag_inexact
, status
);
4346 zSig
= ( zSig
+ roundIncrement
)>>7;
4347 if (!(roundBits
^ 0x40) && roundNearestEven
) {
4350 if ( zSig
== 0 ) zExp
= 0;
4351 return packFloat32( zSign
, zExp
, zSig
);
4355 /*----------------------------------------------------------------------------
4356 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4357 | and significand `zSig', and returns the proper single-precision floating-
4358 | point value corresponding to the abstract input. This routine is just like
4359 | `roundAndPackFloat32' except that `zSig' does not have to be normalized.
4360 | Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
4361 | floating-point exponent.
4362 *----------------------------------------------------------------------------*/
4365 normalizeRoundAndPackFloat32(bool zSign
, int zExp
, uint32_t zSig
,
4366 float_status
*status
)
4370 shiftCount
= clz32(zSig
) - 1;
4371 return roundAndPackFloat32(zSign
, zExp
- shiftCount
, zSig
<<shiftCount
,
4376 /*----------------------------------------------------------------------------
4377 | Normalizes the subnormal double-precision floating-point value represented
4378 | by the denormalized significand `aSig'. The normalized exponent and
4379 | significand are stored at the locations pointed to by `zExpPtr' and
4380 | `zSigPtr', respectively.
4381 *----------------------------------------------------------------------------*/
4384 normalizeFloat64Subnormal(uint64_t aSig
, int *zExpPtr
, uint64_t *zSigPtr
)
4388 shiftCount
= clz64(aSig
) - 11;
4389 *zSigPtr
= aSig
<<shiftCount
;
4390 *zExpPtr
= 1 - shiftCount
;
4394 /*----------------------------------------------------------------------------
4395 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
4396 | double-precision floating-point value, returning the result. After being
4397 | shifted into the proper positions, the three fields are simply added
4398 | together to form the result. This means that any integer portion of `zSig'
4399 | will be added into the exponent. Since a properly normalized significand
4400 | will have an integer portion equal to 1, the `zExp' input should be 1 less
4401 | than the desired result exponent whenever `zSig' is a complete, normalized
4403 *----------------------------------------------------------------------------*/
4405 static inline float64
packFloat64(bool zSign
, int zExp
, uint64_t zSig
)
4408 return make_float64(
4409 ( ( (uint64_t) zSign
)<<63 ) + ( ( (uint64_t) zExp
)<<52 ) + zSig
);
4413 /*----------------------------------------------------------------------------
4414 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4415 | and significand `zSig', and returns the proper double-precision floating-
4416 | point value corresponding to the abstract input. Ordinarily, the abstract
4417 | value is simply rounded and packed into the double-precision format, with
4418 | the inexact exception raised if the abstract input cannot be represented
4419 | exactly. However, if the abstract value is too large, the overflow and
4420 | inexact exceptions are raised and an infinity or maximal finite value is
4421 | returned. If the abstract value is too small, the input value is rounded to
4422 | a subnormal number, and the underflow and inexact exceptions are raised if
4423 | the abstract input cannot be represented exactly as a subnormal double-
4424 | precision floating-point number.
4425 | The input significand `zSig' has its binary point between bits 62
4426 | and 61, which is 10 bits to the left of the usual location. This shifted
4427 | significand must be normalized or smaller. If `zSig' is not normalized,
4428 | `zExp' must be 0; in that case, the result returned is a subnormal number,
4429 | and it must not require rounding. In the usual case that `zSig' is
4430 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
4431 | The handling of underflow and overflow follows the IEC/IEEE Standard for
4432 | Binary Floating-Point Arithmetic.
4433 *----------------------------------------------------------------------------*/
4435 static float64
roundAndPackFloat64(bool zSign
, int zExp
, uint64_t zSig
,
4436 float_status
*status
)
4438 int8_t roundingMode
;
4439 bool roundNearestEven
;
4440 int roundIncrement
, roundBits
;
4443 roundingMode
= status
->float_rounding_mode
;
4444 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4445 switch (roundingMode
) {
4446 case float_round_nearest_even
:
4447 case float_round_ties_away
:
4448 roundIncrement
= 0x200;
4450 case float_round_to_zero
:
4453 case float_round_up
:
4454 roundIncrement
= zSign
? 0 : 0x3ff;
4456 case float_round_down
:
4457 roundIncrement
= zSign
? 0x3ff : 0;
4459 case float_round_to_odd
:
4460 roundIncrement
= (zSig
& 0x400) ? 0 : 0x3ff;
4465 roundBits
= zSig
& 0x3FF;
4466 if ( 0x7FD <= (uint16_t) zExp
) {
4467 if ( ( 0x7FD < zExp
)
4468 || ( ( zExp
== 0x7FD )
4469 && ( (int64_t) ( zSig
+ roundIncrement
) < 0 ) )
4471 bool overflow_to_inf
= roundingMode
!= float_round_to_odd
&&
4472 roundIncrement
!= 0;
4473 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4474 return packFloat64(zSign
, 0x7FF, -(!overflow_to_inf
));
4477 if (status
->flush_to_zero
) {
4478 float_raise(float_flag_output_denormal
, status
);
4479 return packFloat64(zSign
, 0, 0);
4481 isTiny
= status
->tininess_before_rounding
4483 || (zSig
+ roundIncrement
< UINT64_C(0x8000000000000000));
4484 shift64RightJamming( zSig
, - zExp
, &zSig
);
4486 roundBits
= zSig
& 0x3FF;
4487 if (isTiny
&& roundBits
) {
4488 float_raise(float_flag_underflow
, status
);
4490 if (roundingMode
== float_round_to_odd
) {
4492 * For round-to-odd case, the roundIncrement depends on
4493 * zSig which just changed.
4495 roundIncrement
= (zSig
& 0x400) ? 0 : 0x3ff;
4500 float_raise(float_flag_inexact
, status
);
4502 zSig
= ( zSig
+ roundIncrement
)>>10;
4503 if (!(roundBits
^ 0x200) && roundNearestEven
) {
4506 if ( zSig
== 0 ) zExp
= 0;
4507 return packFloat64( zSign
, zExp
, zSig
);
4511 /*----------------------------------------------------------------------------
4512 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4513 | and significand `zSig', and returns the proper double-precision floating-
4514 | point value corresponding to the abstract input. This routine is just like
4515 | `roundAndPackFloat64' except that `zSig' does not have to be normalized.
4516 | Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
4517 | floating-point exponent.
4518 *----------------------------------------------------------------------------*/
4521 normalizeRoundAndPackFloat64(bool zSign
, int zExp
, uint64_t zSig
,
4522 float_status
*status
)
4526 shiftCount
= clz64(zSig
) - 1;
4527 return roundAndPackFloat64(zSign
, zExp
- shiftCount
, zSig
<<shiftCount
,
4532 /*----------------------------------------------------------------------------
4533 | Normalizes the subnormal extended double-precision floating-point value
4534 | represented by the denormalized significand `aSig'. The normalized exponent
4535 | and significand are stored at the locations pointed to by `zExpPtr' and
4536 | `zSigPtr', respectively.
4537 *----------------------------------------------------------------------------*/
4539 void normalizeFloatx80Subnormal(uint64_t aSig
, int32_t *zExpPtr
,
4544 shiftCount
= clz64(aSig
);
4545 *zSigPtr
= aSig
<<shiftCount
;
4546 *zExpPtr
= 1 - shiftCount
;
4549 /*----------------------------------------------------------------------------
4550 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4551 | and extended significand formed by the concatenation of `zSig0' and `zSig1',
4552 | and returns the proper extended double-precision floating-point value
4553 | corresponding to the abstract input. Ordinarily, the abstract value is
4554 | rounded and packed into the extended double-precision format, with the
4555 | inexact exception raised if the abstract input cannot be represented
4556 | exactly. However, if the abstract value is too large, the overflow and
4557 | inexact exceptions are raised and an infinity or maximal finite value is
4558 | returned. If the abstract value is too small, the input value is rounded to
4559 | a subnormal number, and the underflow and inexact exceptions are raised if
4560 | the abstract input cannot be represented exactly as a subnormal extended
4561 | double-precision floating-point number.
4562 | If `roundingPrecision' is floatx80_precision_s or floatx80_precision_d,
4563 | the result is rounded to the same number of bits as single or double
4564 | precision, respectively. Otherwise, the result is rounded to the full
4565 | precision of the extended double-precision format.
4566 | The input significand must be normalized or smaller. If the input
4567 | significand is not normalized, `zExp' must be 0; in that case, the result
4568 | returned is a subnormal number, and it must not require rounding. The
4569 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
4570 | Floating-Point Arithmetic.
4571 *----------------------------------------------------------------------------*/
4573 floatx80
roundAndPackFloatx80(FloatX80RoundPrec roundingPrecision
, bool zSign
,
4574 int32_t zExp
, uint64_t zSig0
, uint64_t zSig1
,
4575 float_status
*status
)
4577 FloatRoundMode roundingMode
;
4578 bool roundNearestEven
, increment
, isTiny
;
4579 int64_t roundIncrement
, roundMask
, roundBits
;
4581 roundingMode
= status
->float_rounding_mode
;
4582 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4583 switch (roundingPrecision
) {
4584 case floatx80_precision_x
:
4586 case floatx80_precision_d
:
4587 roundIncrement
= UINT64_C(0x0000000000000400);
4588 roundMask
= UINT64_C(0x00000000000007FF);
4590 case floatx80_precision_s
:
4591 roundIncrement
= UINT64_C(0x0000008000000000);
4592 roundMask
= UINT64_C(0x000000FFFFFFFFFF);
4595 g_assert_not_reached();
4597 zSig0
|= ( zSig1
!= 0 );
4598 switch (roundingMode
) {
4599 case float_round_nearest_even
:
4600 case float_round_ties_away
:
4602 case float_round_to_zero
:
4605 case float_round_up
:
4606 roundIncrement
= zSign
? 0 : roundMask
;
4608 case float_round_down
:
4609 roundIncrement
= zSign
? roundMask
: 0;
4614 roundBits
= zSig0
& roundMask
;
4615 if ( 0x7FFD <= (uint32_t) ( zExp
- 1 ) ) {
4616 if ( ( 0x7FFE < zExp
)
4617 || ( ( zExp
== 0x7FFE ) && ( zSig0
+ roundIncrement
< zSig0
) )
4622 if (status
->flush_to_zero
) {
4623 float_raise(float_flag_output_denormal
, status
);
4624 return packFloatx80(zSign
, 0, 0);
4626 isTiny
= status
->tininess_before_rounding
4628 || (zSig0
<= zSig0
+ roundIncrement
);
4629 shift64RightJamming( zSig0
, 1 - zExp
, &zSig0
);
4631 roundBits
= zSig0
& roundMask
;
4632 if (isTiny
&& roundBits
) {
4633 float_raise(float_flag_underflow
, status
);
4636 float_raise(float_flag_inexact
, status
);
4638 zSig0
+= roundIncrement
;
4639 if ( (int64_t) zSig0
< 0 ) zExp
= 1;
4640 roundIncrement
= roundMask
+ 1;
4641 if ( roundNearestEven
&& ( roundBits
<<1 == roundIncrement
) ) {
4642 roundMask
|= roundIncrement
;
4644 zSig0
&= ~ roundMask
;
4645 return packFloatx80( zSign
, zExp
, zSig0
);
4649 float_raise(float_flag_inexact
, status
);
4651 zSig0
+= roundIncrement
;
4652 if ( zSig0
< roundIncrement
) {
4654 zSig0
= UINT64_C(0x8000000000000000);
4656 roundIncrement
= roundMask
+ 1;
4657 if ( roundNearestEven
&& ( roundBits
<<1 == roundIncrement
) ) {
4658 roundMask
|= roundIncrement
;
4660 zSig0
&= ~ roundMask
;
4661 if ( zSig0
== 0 ) zExp
= 0;
4662 return packFloatx80( zSign
, zExp
, zSig0
);
4664 switch (roundingMode
) {
4665 case float_round_nearest_even
:
4666 case float_round_ties_away
:
4667 increment
= ((int64_t)zSig1
< 0);
4669 case float_round_to_zero
:
4672 case float_round_up
:
4673 increment
= !zSign
&& zSig1
;
4675 case float_round_down
:
4676 increment
= zSign
&& zSig1
;
4681 if ( 0x7FFD <= (uint32_t) ( zExp
- 1 ) ) {
4682 if ( ( 0x7FFE < zExp
)
4683 || ( ( zExp
== 0x7FFE )
4684 && ( zSig0
== UINT64_C(0xFFFFFFFFFFFFFFFF) )
4690 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4691 if ( ( roundingMode
== float_round_to_zero
)
4692 || ( zSign
&& ( roundingMode
== float_round_up
) )
4693 || ( ! zSign
&& ( roundingMode
== float_round_down
) )
4695 return packFloatx80( zSign
, 0x7FFE, ~ roundMask
);
4697 return packFloatx80(zSign
,
4698 floatx80_infinity_high
,
4699 floatx80_infinity_low
);
4702 isTiny
= status
->tininess_before_rounding
4705 || (zSig0
< UINT64_C(0xFFFFFFFFFFFFFFFF));
4706 shift64ExtraRightJamming( zSig0
, zSig1
, 1 - zExp
, &zSig0
, &zSig1
);
4708 if (isTiny
&& zSig1
) {
4709 float_raise(float_flag_underflow
, status
);
4712 float_raise(float_flag_inexact
, status
);
4714 switch (roundingMode
) {
4715 case float_round_nearest_even
:
4716 case float_round_ties_away
:
4717 increment
= ((int64_t)zSig1
< 0);
4719 case float_round_to_zero
:
4722 case float_round_up
:
4723 increment
= !zSign
&& zSig1
;
4725 case float_round_down
:
4726 increment
= zSign
&& zSig1
;
4733 if (!(zSig1
<< 1) && roundNearestEven
) {
4736 if ( (int64_t) zSig0
< 0 ) zExp
= 1;
4738 return packFloatx80( zSign
, zExp
, zSig0
);
4742 float_raise(float_flag_inexact
, status
);
4748 zSig0
= UINT64_C(0x8000000000000000);
4751 if (!(zSig1
<< 1) && roundNearestEven
) {
4757 if ( zSig0
== 0 ) zExp
= 0;
4759 return packFloatx80( zSign
, zExp
, zSig0
);
4763 /*----------------------------------------------------------------------------
4764 | Takes an abstract floating-point value having sign `zSign', exponent
4765 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
4766 | and returns the proper extended double-precision floating-point value
4767 | corresponding to the abstract input. This routine is just like
4768 | `roundAndPackFloatx80' except that the input significand does not have to be
4770 *----------------------------------------------------------------------------*/
4772 floatx80
normalizeRoundAndPackFloatx80(FloatX80RoundPrec roundingPrecision
,
4773 bool zSign
, int32_t zExp
,
4774 uint64_t zSig0
, uint64_t zSig1
,
4775 float_status
*status
)
4784 shiftCount
= clz64(zSig0
);
4785 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
4787 return roundAndPackFloatx80(roundingPrecision
, zSign
, zExp
,
4788 zSig0
, zSig1
, status
);
4792 /*----------------------------------------------------------------------------
4793 | Returns the least-significant 64 fraction bits of the quadruple-precision
4794 | floating-point value `a'.
4795 *----------------------------------------------------------------------------*/
4797 static inline uint64_t extractFloat128Frac1( float128 a
)
4804 /*----------------------------------------------------------------------------
4805 | Returns the most-significant 48 fraction bits of the quadruple-precision
4806 | floating-point value `a'.
4807 *----------------------------------------------------------------------------*/
4809 static inline uint64_t extractFloat128Frac0( float128 a
)
4812 return a
.high
& UINT64_C(0x0000FFFFFFFFFFFF);
4816 /*----------------------------------------------------------------------------
4817 | Returns the exponent bits of the quadruple-precision floating-point value
4819 *----------------------------------------------------------------------------*/
4821 static inline int32_t extractFloat128Exp( float128 a
)
4824 return ( a
.high
>>48 ) & 0x7FFF;
4828 /*----------------------------------------------------------------------------
4829 | Returns the sign bit of the quadruple-precision floating-point value `a'.
4830 *----------------------------------------------------------------------------*/
4832 static inline bool extractFloat128Sign(float128 a
)
4834 return a
.high
>> 63;
4837 /*----------------------------------------------------------------------------
4838 | Normalizes the subnormal quadruple-precision floating-point value
4839 | represented by the denormalized significand formed by the concatenation of
4840 | `aSig0' and `aSig1'. The normalized exponent is stored at the location
4841 | pointed to by `zExpPtr'. The most significant 49 bits of the normalized
4842 | significand are stored at the location pointed to by `zSig0Ptr', and the
4843 | least significant 64 bits of the normalized significand are stored at the
4844 | location pointed to by `zSig1Ptr'.
4845 *----------------------------------------------------------------------------*/
4848 normalizeFloat128Subnormal(
4859 shiftCount
= clz64(aSig1
) - 15;
4860 if ( shiftCount
< 0 ) {
4861 *zSig0Ptr
= aSig1
>>( - shiftCount
);
4862 *zSig1Ptr
= aSig1
<<( shiftCount
& 63 );
4865 *zSig0Ptr
= aSig1
<<shiftCount
;
4868 *zExpPtr
= - shiftCount
- 63;
4871 shiftCount
= clz64(aSig0
) - 15;
4872 shortShift128Left( aSig0
, aSig1
, shiftCount
, zSig0Ptr
, zSig1Ptr
);
4873 *zExpPtr
= 1 - shiftCount
;
4878 /*----------------------------------------------------------------------------
4879 | Packs the sign `zSign', the exponent `zExp', and the significand formed
4880 | by the concatenation of `zSig0' and `zSig1' into a quadruple-precision
4881 | floating-point value, returning the result. After being shifted into the
4882 | proper positions, the three fields `zSign', `zExp', and `zSig0' are simply
4883 | added together to form the most significant 32 bits of the result. This
4884 | means that any integer portion of `zSig0' will be added into the exponent.
4885 | Since a properly normalized significand will have an integer portion equal
4886 | to 1, the `zExp' input should be 1 less than the desired result exponent
4887 | whenever `zSig0' and `zSig1' concatenated form a complete, normalized
4889 *----------------------------------------------------------------------------*/
4891 static inline float128
4892 packFloat128(bool zSign
, int32_t zExp
, uint64_t zSig0
, uint64_t zSig1
)
4897 z
.high
= ((uint64_t)zSign
<< 63) + ((uint64_t)zExp
<< 48) + zSig0
;
4901 /*----------------------------------------------------------------------------
4902 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4903 | and extended significand formed by the concatenation of `zSig0', `zSig1',
4904 | and `zSig2', and returns the proper quadruple-precision floating-point value
4905 | corresponding to the abstract input. Ordinarily, the abstract value is
4906 | simply rounded and packed into the quadruple-precision format, with the
4907 | inexact exception raised if the abstract input cannot be represented
4908 | exactly. However, if the abstract value is too large, the overflow and
4909 | inexact exceptions are raised and an infinity or maximal finite value is
4910 | returned. If the abstract value is too small, the input value is rounded to
4911 | a subnormal number, and the underflow and inexact exceptions are raised if
4912 | the abstract input cannot be represented exactly as a subnormal quadruple-
4913 | precision floating-point number.
4914 | The input significand must be normalized or smaller. If the input
4915 | significand is not normalized, `zExp' must be 0; in that case, the result
4916 | returned is a subnormal number, and it must not require rounding. In the
4917 | usual case that the input significand is normalized, `zExp' must be 1 less
4918 | than the ``true'' floating-point exponent. The handling of underflow and
4919 | overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4920 *----------------------------------------------------------------------------*/
4922 static float128
roundAndPackFloat128(bool zSign
, int32_t zExp
,
4923 uint64_t zSig0
, uint64_t zSig1
,
4924 uint64_t zSig2
, float_status
*status
)
4926 int8_t roundingMode
;
4927 bool roundNearestEven
, increment
, isTiny
;
4929 roundingMode
= status
->float_rounding_mode
;
4930 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4931 switch (roundingMode
) {
4932 case float_round_nearest_even
:
4933 case float_round_ties_away
:
4934 increment
= ((int64_t)zSig2
< 0);
4936 case float_round_to_zero
:
4939 case float_round_up
:
4940 increment
= !zSign
&& zSig2
;
4942 case float_round_down
:
4943 increment
= zSign
&& zSig2
;
4945 case float_round_to_odd
:
4946 increment
= !(zSig1
& 0x1) && zSig2
;
4951 if ( 0x7FFD <= (uint32_t) zExp
) {
4952 if ( ( 0x7FFD < zExp
)
4953 || ( ( zExp
== 0x7FFD )
4955 UINT64_C(0x0001FFFFFFFFFFFF),
4956 UINT64_C(0xFFFFFFFFFFFFFFFF),
4963 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4964 if ( ( roundingMode
== float_round_to_zero
)
4965 || ( zSign
&& ( roundingMode
== float_round_up
) )
4966 || ( ! zSign
&& ( roundingMode
== float_round_down
) )
4967 || (roundingMode
== float_round_to_odd
)
4973 UINT64_C(0x0000FFFFFFFFFFFF),
4974 UINT64_C(0xFFFFFFFFFFFFFFFF)
4977 return packFloat128( zSign
, 0x7FFF, 0, 0 );
4980 if (status
->flush_to_zero
) {
4981 float_raise(float_flag_output_denormal
, status
);
4982 return packFloat128(zSign
, 0, 0, 0);
4984 isTiny
= status
->tininess_before_rounding
4987 || lt128(zSig0
, zSig1
,
4988 UINT64_C(0x0001FFFFFFFFFFFF),
4989 UINT64_C(0xFFFFFFFFFFFFFFFF));
4990 shift128ExtraRightJamming(
4991 zSig0
, zSig1
, zSig2
, - zExp
, &zSig0
, &zSig1
, &zSig2
);
4993 if (isTiny
&& zSig2
) {
4994 float_raise(float_flag_underflow
, status
);
4996 switch (roundingMode
) {
4997 case float_round_nearest_even
:
4998 case float_round_ties_away
:
4999 increment
= ((int64_t)zSig2
< 0);
5001 case float_round_to_zero
:
5004 case float_round_up
:
5005 increment
= !zSign
&& zSig2
;
5007 case float_round_down
:
5008 increment
= zSign
&& zSig2
;
5010 case float_round_to_odd
:
5011 increment
= !(zSig1
& 0x1) && zSig2
;
5019 float_raise(float_flag_inexact
, status
);
5022 add128( zSig0
, zSig1
, 0, 1, &zSig0
, &zSig1
);
5023 if ((zSig2
+ zSig2
== 0) && roundNearestEven
) {
5028 if ( ( zSig0
| zSig1
) == 0 ) zExp
= 0;
5030 return packFloat128( zSign
, zExp
, zSig0
, zSig1
);
5034 /*----------------------------------------------------------------------------
5035 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
5036 | and significand formed by the concatenation of `zSig0' and `zSig1', and
5037 | returns the proper quadruple-precision floating-point value corresponding
5038 | to the abstract input. This routine is just like `roundAndPackFloat128'
5039 | except that the input significand has fewer bits and does not have to be
5040 | normalized. In all cases, `zExp' must be 1 less than the ``true'' floating-
5042 *----------------------------------------------------------------------------*/
5044 static float128
normalizeRoundAndPackFloat128(bool zSign
, int32_t zExp
,
5045 uint64_t zSig0
, uint64_t zSig1
,
5046 float_status
*status
)
5056 shiftCount
= clz64(zSig0
) - 15;
5057 if ( 0 <= shiftCount
) {
5059 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
5062 shift128ExtraRightJamming(
5063 zSig0
, zSig1
, 0, - shiftCount
, &zSig0
, &zSig1
, &zSig2
);
5066 return roundAndPackFloat128(zSign
, zExp
, zSig0
, zSig1
, zSig2
, status
);
5070 /*----------------------------------------------------------------------------
5071 | Returns the remainder of the single-precision floating-point value `a'
5072 | with respect to the corresponding value `b'. The operation is performed
5073 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5074 *----------------------------------------------------------------------------*/
5076 float32
float32_rem(float32 a
, float32 b
, float_status
*status
)
5079 int aExp
, bExp
, expDiff
;
5080 uint32_t aSig
, bSig
;
5082 uint64_t aSig64
, bSig64
, q64
;
5083 uint32_t alternateASig
;
5085 a
= float32_squash_input_denormal(a
, status
);
5086 b
= float32_squash_input_denormal(b
, status
);
5088 aSig
= extractFloat32Frac( a
);
5089 aExp
= extractFloat32Exp( a
);
5090 aSign
= extractFloat32Sign( a
);
5091 bSig
= extractFloat32Frac( b
);
5092 bExp
= extractFloat32Exp( b
);
5093 if ( aExp
== 0xFF ) {
5094 if ( aSig
|| ( ( bExp
== 0xFF ) && bSig
) ) {
5095 return propagateFloat32NaN(a
, b
, status
);
5097 float_raise(float_flag_invalid
, status
);
5098 return float32_default_nan(status
);
5100 if ( bExp
== 0xFF ) {
5102 return propagateFloat32NaN(a
, b
, status
);
5108 float_raise(float_flag_invalid
, status
);
5109 return float32_default_nan(status
);
5111 normalizeFloat32Subnormal( bSig
, &bExp
, &bSig
);
5114 if ( aSig
== 0 ) return a
;
5115 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
5117 expDiff
= aExp
- bExp
;
5120 if ( expDiff
< 32 ) {
5123 if ( expDiff
< 0 ) {
5124 if ( expDiff
< -1 ) return a
;
5127 q
= ( bSig
<= aSig
);
5128 if ( q
) aSig
-= bSig
;
5129 if ( 0 < expDiff
) {
5130 q
= ( ( (uint64_t) aSig
)<<32 ) / bSig
;
5133 aSig
= ( ( aSig
>>1 )<<( expDiff
- 1 ) ) - bSig
* q
;
5141 if ( bSig
<= aSig
) aSig
-= bSig
;
5142 aSig64
= ( (uint64_t) aSig
)<<40;
5143 bSig64
= ( (uint64_t) bSig
)<<40;
5145 while ( 0 < expDiff
) {
5146 q64
= estimateDiv128To64( aSig64
, 0, bSig64
);
5147 q64
= ( 2 < q64
) ? q64
- 2 : 0;
5148 aSig64
= - ( ( bSig
* q64
)<<38 );
5152 q64
= estimateDiv128To64( aSig64
, 0, bSig64
);
5153 q64
= ( 2 < q64
) ? q64
- 2 : 0;
5154 q
= q64
>>( 64 - expDiff
);
5156 aSig
= ( ( aSig64
>>33 )<<( expDiff
- 1 ) ) - bSig
* q
;
5159 alternateASig
= aSig
;
5162 } while ( 0 <= (int32_t) aSig
);
5163 sigMean
= aSig
+ alternateASig
;
5164 if ( ( sigMean
< 0 ) || ( ( sigMean
== 0 ) && ( q
& 1 ) ) ) {
5165 aSig
= alternateASig
;
5167 zSign
= ( (int32_t) aSig
< 0 );
5168 if ( zSign
) aSig
= - aSig
;
5169 return normalizeRoundAndPackFloat32(aSign
^ zSign
, bExp
, aSig
, status
);
5174 /*----------------------------------------------------------------------------
5175 | Returns the binary exponential of the single-precision floating-point value
5176 | `a'. The operation is performed according to the IEC/IEEE Standard for
5177 | Binary Floating-Point Arithmetic.
5179 | Uses the following identities:
5181 | 1. -------------------------------------------------------------------------
5185 | 2. -------------------------------------------------------------------------
5188 | e = 1 + --- + --- + --- + --- + --- + ... + --- + ...
5190 *----------------------------------------------------------------------------*/
5192 static const float64 float32_exp2_coefficients
[15] =
5194 const_float64( 0x3ff0000000000000ll
), /* 1 */
5195 const_float64( 0x3fe0000000000000ll
), /* 2 */
5196 const_float64( 0x3fc5555555555555ll
), /* 3 */
5197 const_float64( 0x3fa5555555555555ll
), /* 4 */
5198 const_float64( 0x3f81111111111111ll
), /* 5 */
5199 const_float64( 0x3f56c16c16c16c17ll
), /* 6 */
5200 const_float64( 0x3f2a01a01a01a01all
), /* 7 */
5201 const_float64( 0x3efa01a01a01a01all
), /* 8 */
5202 const_float64( 0x3ec71de3a556c734ll
), /* 9 */
5203 const_float64( 0x3e927e4fb7789f5cll
), /* 10 */
5204 const_float64( 0x3e5ae64567f544e4ll
), /* 11 */
5205 const_float64( 0x3e21eed8eff8d898ll
), /* 12 */
5206 const_float64( 0x3de6124613a86d09ll
), /* 13 */
5207 const_float64( 0x3da93974a8c07c9dll
), /* 14 */
5208 const_float64( 0x3d6ae7f3e733b81fll
), /* 15 */
5211 float32
float32_exp2(float32 a
, float_status
*status
)
5213 FloatParts64 xp
, xnp
, tp
, rp
;
5216 float32_unpack_canonical(&xp
, a
, status
);
5217 if (unlikely(xp
.cls
!= float_class_normal
)) {
5219 case float_class_snan
:
5220 case float_class_qnan
:
5221 parts_return_nan(&xp
, status
);
5222 return float32_round_pack_canonical(&xp
, status
);
5223 case float_class_inf
:
5224 return xp
.sign
? float32_zero
: a
;
5225 case float_class_zero
:
5230 g_assert_not_reached();
5233 float_raise(float_flag_inexact
, status
);
5235 float64_unpack_canonical(&xnp
, float64_ln2
, status
);
5236 xp
= *parts_mul(&xp
, &tp
, status
);
5239 float64_unpack_canonical(&rp
, float64_one
, status
);
5240 for (i
= 0 ; i
< 15 ; i
++) {
5241 float64_unpack_canonical(&tp
, float32_exp2_coefficients
[i
], status
);
5242 rp
= *parts_muladd(&tp
, &xp
, &rp
, 0, status
);
5243 xnp
= *parts_mul(&xnp
, &xp
, status
);
5246 return float32_round_pack_canonical(&rp
, status
);
5249 /*----------------------------------------------------------------------------
5250 | Returns the binary log of the single-precision floating-point value `a'.
5251 | The operation is performed according to the IEC/IEEE Standard for Binary
5252 | Floating-Point Arithmetic.
5253 *----------------------------------------------------------------------------*/
5254 float32
float32_log2(float32 a
, float_status
*status
)
5258 uint32_t aSig
, zSig
, i
;
5260 a
= float32_squash_input_denormal(a
, status
);
5261 aSig
= extractFloat32Frac( a
);
5262 aExp
= extractFloat32Exp( a
);
5263 aSign
= extractFloat32Sign( a
);
5266 if ( aSig
== 0 ) return packFloat32( 1, 0xFF, 0 );
5267 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
5270 float_raise(float_flag_invalid
, status
);
5271 return float32_default_nan(status
);
5273 if ( aExp
== 0xFF ) {
5275 return propagateFloat32NaN(a
, float32_zero
, status
);
5285 for (i
= 1 << 22; i
> 0; i
>>= 1) {
5286 aSig
= ( (uint64_t)aSig
* aSig
) >> 23;
5287 if ( aSig
& 0x01000000 ) {
5296 return normalizeRoundAndPackFloat32(zSign
, 0x85, zSig
, status
);
5299 /*----------------------------------------------------------------------------
5300 | Returns the remainder of the double-precision floating-point value `a'
5301 | with respect to the corresponding value `b'. The operation is performed
5302 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5303 *----------------------------------------------------------------------------*/
5305 float64
float64_rem(float64 a
, float64 b
, float_status
*status
)
5308 int aExp
, bExp
, expDiff
;
5309 uint64_t aSig
, bSig
;
5310 uint64_t q
, alternateASig
;
5313 a
= float64_squash_input_denormal(a
, status
);
5314 b
= float64_squash_input_denormal(b
, status
);
5315 aSig
= extractFloat64Frac( a
);
5316 aExp
= extractFloat64Exp( a
);
5317 aSign
= extractFloat64Sign( a
);
5318 bSig
= extractFloat64Frac( b
);
5319 bExp
= extractFloat64Exp( b
);
5320 if ( aExp
== 0x7FF ) {
5321 if ( aSig
|| ( ( bExp
== 0x7FF ) && bSig
) ) {
5322 return propagateFloat64NaN(a
, b
, status
);
5324 float_raise(float_flag_invalid
, status
);
5325 return float64_default_nan(status
);
5327 if ( bExp
== 0x7FF ) {
5329 return propagateFloat64NaN(a
, b
, status
);
5335 float_raise(float_flag_invalid
, status
);
5336 return float64_default_nan(status
);
5338 normalizeFloat64Subnormal( bSig
, &bExp
, &bSig
);
5341 if ( aSig
== 0 ) return a
;
5342 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
5344 expDiff
= aExp
- bExp
;
5345 aSig
= (aSig
| UINT64_C(0x0010000000000000)) << 11;
5346 bSig
= (bSig
| UINT64_C(0x0010000000000000)) << 11;
5347 if ( expDiff
< 0 ) {
5348 if ( expDiff
< -1 ) return a
;
5351 q
= ( bSig
<= aSig
);
5352 if ( q
) aSig
-= bSig
;
5354 while ( 0 < expDiff
) {
5355 q
= estimateDiv128To64( aSig
, 0, bSig
);
5356 q
= ( 2 < q
) ? q
- 2 : 0;
5357 aSig
= - ( ( bSig
>>2 ) * q
);
5361 if ( 0 < expDiff
) {
5362 q
= estimateDiv128To64( aSig
, 0, bSig
);
5363 q
= ( 2 < q
) ? q
- 2 : 0;
5366 aSig
= ( ( aSig
>>1 )<<( expDiff
- 1 ) ) - bSig
* q
;
5373 alternateASig
= aSig
;
5376 } while ( 0 <= (int64_t) aSig
);
5377 sigMean
= aSig
+ alternateASig
;
5378 if ( ( sigMean
< 0 ) || ( ( sigMean
== 0 ) && ( q
& 1 ) ) ) {
5379 aSig
= alternateASig
;
5381 zSign
= ( (int64_t) aSig
< 0 );
5382 if ( zSign
) aSig
= - aSig
;
5383 return normalizeRoundAndPackFloat64(aSign
^ zSign
, bExp
, aSig
, status
);
5387 /*----------------------------------------------------------------------------
5388 | Returns the binary log of the double-precision floating-point value `a'.
5389 | The operation is performed according to the IEC/IEEE Standard for Binary
5390 | Floating-Point Arithmetic.
5391 *----------------------------------------------------------------------------*/
5392 float64
float64_log2(float64 a
, float_status
*status
)
5396 uint64_t aSig
, aSig0
, aSig1
, zSig
, i
;
5397 a
= float64_squash_input_denormal(a
, status
);
5399 aSig
= extractFloat64Frac( a
);
5400 aExp
= extractFloat64Exp( a
);
5401 aSign
= extractFloat64Sign( a
);
5404 if ( aSig
== 0 ) return packFloat64( 1, 0x7FF, 0 );
5405 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
5408 float_raise(float_flag_invalid
, status
);
5409 return float64_default_nan(status
);
5411 if ( aExp
== 0x7FF ) {
5413 return propagateFloat64NaN(a
, float64_zero
, status
);
5419 aSig
|= UINT64_C(0x0010000000000000);
5421 zSig
= (uint64_t)aExp
<< 52;
5422 for (i
= 1LL << 51; i
> 0; i
>>= 1) {
5423 mul64To128( aSig
, aSig
, &aSig0
, &aSig1
);
5424 aSig
= ( aSig0
<< 12 ) | ( aSig1
>> 52 );
5425 if ( aSig
& UINT64_C(0x0020000000000000) ) {
5433 return normalizeRoundAndPackFloat64(zSign
, 0x408, zSig
, status
);
5436 /*----------------------------------------------------------------------------
5437 | Rounds the extended double-precision floating-point value `a'
5438 | to the precision provided by floatx80_rounding_precision and returns the
5439 | result as an extended double-precision floating-point value.
5440 | The operation is performed according to the IEC/IEEE Standard for Binary
5441 | Floating-Point Arithmetic.
5442 *----------------------------------------------------------------------------*/
5444 floatx80
floatx80_round(floatx80 a
, float_status
*status
)
5448 if (!floatx80_unpack_canonical(&p
, a
, status
)) {
5449 return floatx80_default_nan(status
);
5451 return floatx80_round_pack_canonical(&p
, status
);
5454 /*----------------------------------------------------------------------------
5455 | Returns the remainder of the extended double-precision floating-point value
5456 | `a' with respect to the corresponding value `b'. The operation is performed
5457 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic,
5458 | if 'mod' is false; if 'mod' is true, return the remainder based on truncating
5459 | the quotient toward zero instead. '*quotient' is set to the low 64 bits of
5460 | the absolute value of the integer quotient.
5461 *----------------------------------------------------------------------------*/
5463 floatx80
floatx80_modrem(floatx80 a
, floatx80 b
, bool mod
, uint64_t *quotient
,
5464 float_status
*status
)
5467 int32_t aExp
, bExp
, expDiff
, aExpOrig
;
5468 uint64_t aSig0
, aSig1
, bSig
;
5469 uint64_t q
, term0
, term1
, alternateASig0
, alternateASig1
;
5472 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
5473 float_raise(float_flag_invalid
, status
);
5474 return floatx80_default_nan(status
);
5476 aSig0
= extractFloatx80Frac( a
);
5477 aExpOrig
= aExp
= extractFloatx80Exp( a
);
5478 aSign
= extractFloatx80Sign( a
);
5479 bSig
= extractFloatx80Frac( b
);
5480 bExp
= extractFloatx80Exp( b
);
5481 if ( aExp
== 0x7FFF ) {
5482 if ( (uint64_t) ( aSig0
<<1 )
5483 || ( ( bExp
== 0x7FFF ) && (uint64_t) ( bSig
<<1 ) ) ) {
5484 return propagateFloatx80NaN(a
, b
, status
);
5488 if ( bExp
== 0x7FFF ) {
5489 if ((uint64_t)(bSig
<< 1)) {
5490 return propagateFloatx80NaN(a
, b
, status
);
5492 if (aExp
== 0 && aSig0
>> 63) {
5494 * Pseudo-denormal argument must be returned in normalized
5497 return packFloatx80(aSign
, 1, aSig0
);
5504 float_raise(float_flag_invalid
, status
);
5505 return floatx80_default_nan(status
);
5507 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
5510 if ( aSig0
== 0 ) return a
;
5511 normalizeFloatx80Subnormal( aSig0
, &aExp
, &aSig0
);
5514 expDiff
= aExp
- bExp
;
5516 if ( expDiff
< 0 ) {
5517 if ( mod
|| expDiff
< -1 ) {
5518 if (aExp
== 1 && aExpOrig
== 0) {
5520 * Pseudo-denormal argument must be returned in
5523 return packFloatx80(aSign
, aExp
, aSig0
);
5527 shift128Right( aSig0
, 0, 1, &aSig0
, &aSig1
);
5530 *quotient
= q
= ( bSig
<= aSig0
);
5531 if ( q
) aSig0
-= bSig
;
5533 while ( 0 < expDiff
) {
5534 q
= estimateDiv128To64( aSig0
, aSig1
, bSig
);
5535 q
= ( 2 < q
) ? q
- 2 : 0;
5536 mul64To128( bSig
, q
, &term0
, &term1
);
5537 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
5538 shortShift128Left( aSig0
, aSig1
, 62, &aSig0
, &aSig1
);
5544 if ( 0 < expDiff
) {
5545 q
= estimateDiv128To64( aSig0
, aSig1
, bSig
);
5546 q
= ( 2 < q
) ? q
- 2 : 0;
5548 mul64To128( bSig
, q
<<( 64 - expDiff
), &term0
, &term1
);
5549 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
5550 shortShift128Left( 0, bSig
, 64 - expDiff
, &term0
, &term1
);
5551 while ( le128( term0
, term1
, aSig0
, aSig1
) ) {
5553 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
5556 *quotient
<<= expDiff
;
5567 sub128( term0
, term1
, aSig0
, aSig1
, &alternateASig0
, &alternateASig1
);
5568 if ( lt128( alternateASig0
, alternateASig1
, aSig0
, aSig1
)
5569 || ( eq128( alternateASig0
, alternateASig1
, aSig0
, aSig1
)
5572 aSig0
= alternateASig0
;
5573 aSig1
= alternateASig1
;
5579 normalizeRoundAndPackFloatx80(
5580 floatx80_precision_x
, zSign
, bExp
+ expDiff
, aSig0
, aSig1
, status
);
5584 /*----------------------------------------------------------------------------
5585 | Returns the remainder of the extended double-precision floating-point value
5586 | `a' with respect to the corresponding value `b'. The operation is performed
5587 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5588 *----------------------------------------------------------------------------*/
5590 floatx80
floatx80_rem(floatx80 a
, floatx80 b
, float_status
*status
)
5593 return floatx80_modrem(a
, b
, false, "ient
, status
);
5596 /*----------------------------------------------------------------------------
5597 | Returns the remainder of the extended double-precision floating-point value
5598 | `a' with respect to the corresponding value `b', with the quotient truncated
5600 *----------------------------------------------------------------------------*/
5602 floatx80
floatx80_mod(floatx80 a
, floatx80 b
, float_status
*status
)
5605 return floatx80_modrem(a
, b
, true, "ient
, status
);
5608 /*----------------------------------------------------------------------------
5609 | Returns the remainder of the quadruple-precision floating-point value `a'
5610 | with respect to the corresponding value `b'. The operation is performed
5611 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5612 *----------------------------------------------------------------------------*/
5614 float128
float128_rem(float128 a
, float128 b
, float_status
*status
)
5617 int32_t aExp
, bExp
, expDiff
;
5618 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, q
, term0
, term1
, term2
;
5619 uint64_t allZero
, alternateASig0
, alternateASig1
, sigMean1
;
5622 aSig1
= extractFloat128Frac1( a
);
5623 aSig0
= extractFloat128Frac0( a
);
5624 aExp
= extractFloat128Exp( a
);
5625 aSign
= extractFloat128Sign( a
);
5626 bSig1
= extractFloat128Frac1( b
);
5627 bSig0
= extractFloat128Frac0( b
);
5628 bExp
= extractFloat128Exp( b
);
5629 if ( aExp
== 0x7FFF ) {
5630 if ( ( aSig0
| aSig1
)
5631 || ( ( bExp
== 0x7FFF ) && ( bSig0
| bSig1
) ) ) {
5632 return propagateFloat128NaN(a
, b
, status
);
5636 if ( bExp
== 0x7FFF ) {
5637 if (bSig0
| bSig1
) {
5638 return propagateFloat128NaN(a
, b
, status
);
5643 if ( ( bSig0
| bSig1
) == 0 ) {
5645 float_raise(float_flag_invalid
, status
);
5646 return float128_default_nan(status
);
5648 normalizeFloat128Subnormal( bSig0
, bSig1
, &bExp
, &bSig0
, &bSig1
);
5651 if ( ( aSig0
| aSig1
) == 0 ) return a
;
5652 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
5654 expDiff
= aExp
- bExp
;
5655 if ( expDiff
< -1 ) return a
;
5657 aSig0
| UINT64_C(0x0001000000000000),
5659 15 - ( expDiff
< 0 ),
5664 bSig0
| UINT64_C(0x0001000000000000), bSig1
, 15, &bSig0
, &bSig1
);
5665 q
= le128( bSig0
, bSig1
, aSig0
, aSig1
);
5666 if ( q
) sub128( aSig0
, aSig1
, bSig0
, bSig1
, &aSig0
, &aSig1
);
5668 while ( 0 < expDiff
) {
5669 q
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
5670 q
= ( 4 < q
) ? q
- 4 : 0;
5671 mul128By64To192( bSig0
, bSig1
, q
, &term0
, &term1
, &term2
);
5672 shortShift192Left( term0
, term1
, term2
, 61, &term1
, &term2
, &allZero
);
5673 shortShift128Left( aSig0
, aSig1
, 61, &aSig0
, &allZero
);
5674 sub128( aSig0
, 0, term1
, term2
, &aSig0
, &aSig1
);
5677 if ( -64 < expDiff
) {
5678 q
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
5679 q
= ( 4 < q
) ? q
- 4 : 0;
5681 shift128Right( bSig0
, bSig1
, 12, &bSig0
, &bSig1
);
5683 if ( expDiff
< 0 ) {
5684 shift128Right( aSig0
, aSig1
, - expDiff
, &aSig0
, &aSig1
);
5687 shortShift128Left( aSig0
, aSig1
, expDiff
, &aSig0
, &aSig1
);
5689 mul128By64To192( bSig0
, bSig1
, q
, &term0
, &term1
, &term2
);
5690 sub128( aSig0
, aSig1
, term1
, term2
, &aSig0
, &aSig1
);
5693 shift128Right( aSig0
, aSig1
, 12, &aSig0
, &aSig1
);
5694 shift128Right( bSig0
, bSig1
, 12, &bSig0
, &bSig1
);
5697 alternateASig0
= aSig0
;
5698 alternateASig1
= aSig1
;
5700 sub128( aSig0
, aSig1
, bSig0
, bSig1
, &aSig0
, &aSig1
);
5701 } while ( 0 <= (int64_t) aSig0
);
5703 aSig0
, aSig1
, alternateASig0
, alternateASig1
, (uint64_t *)&sigMean0
, &sigMean1
);
5704 if ( ( sigMean0
< 0 )
5705 || ( ( ( sigMean0
| sigMean1
) == 0 ) && ( q
& 1 ) ) ) {
5706 aSig0
= alternateASig0
;
5707 aSig1
= alternateASig1
;
5709 zSign
= ( (int64_t) aSig0
< 0 );
5710 if ( zSign
) sub128( 0, 0, aSig0
, aSig1
, &aSig0
, &aSig1
);
5711 return normalizeRoundAndPackFloat128(aSign
^ zSign
, bExp
- 4, aSig0
, aSig1
,
5714 static void __attribute__((constructor
)) softfloat_init(void)
5716 union_float64 ua
, ub
, uc
, ur
;
5718 if (QEMU_NO_HARDFLOAT
) {
5722 * Test that the host's FMA is not obviously broken. For example,
5723 * glibc < 2.23 can perform an incorrect FMA on certain hosts; see
5724 * https://sourceware.org/bugzilla/show_bug.cgi?id=13304
5726 ua
.s
= 0x0020000000000001ULL
;
5727 ub
.s
= 0x3ca0000000000000ULL
;
5728 uc
.s
= 0x0020000000000000ULL
;
5729 ur
.h
= fma(ua
.h
, ub
.h
, uc
.h
);
5730 if (ur
.s
!= 0x0020000000000001ULL
) {
5731 force_soft_fma
= true;