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
,
486 /* Simple helpers for checking if, or what kind of, NaN we have */
487 static inline __attribute__((unused
)) bool is_nan(FloatClass c
)
489 return unlikely(c
>= float_class_qnan
);
492 static inline __attribute__((unused
)) bool is_snan(FloatClass c
)
494 return c
== float_class_snan
;
497 static inline __attribute__((unused
)) bool is_qnan(FloatClass c
)
499 return c
== float_class_qnan
;
503 * Structure holding all of the decomposed parts of a float.
504 * The exponent is unbiased and the fraction is normalized.
506 * The fraction words are stored in big-endian word ordering,
507 * so that truncation from a larger format to a smaller format
508 * can be done simply by ignoring subsequent elements.
516 /* Routines that know the structure may reference the singular name. */
519 * Routines expanded with multiple structures reference "hi" and "lo"
520 * depending on the operation. In FloatParts64, "hi" and "lo" are
521 * both the same word and aliased here.
541 uint64_t frac_hm
; /* high-middle */
542 uint64_t frac_lm
; /* low-middle */
546 /* These apply to the most significant word of each FloatPartsN. */
547 #define DECOMPOSED_BINARY_POINT 63
548 #define DECOMPOSED_IMPLICIT_BIT (1ull << DECOMPOSED_BINARY_POINT)
550 /* Structure holding all of the relevant parameters for a format.
551 * exp_size: the size of the exponent field
552 * exp_bias: the offset applied to the exponent field
553 * exp_max: the maximum normalised exponent
554 * frac_size: the size of the fraction field
555 * frac_shift: shift to normalise the fraction with DECOMPOSED_BINARY_POINT
556 * The following are computed based the size of fraction
557 * frac_lsb: least significant bit of fraction
558 * frac_lsbm1: the bit below the least significant bit (for rounding)
559 * round_mask/roundeven_mask: masks used for rounding
560 * The following optional modifiers are available:
561 * arm_althp: handle ARM Alternative Half Precision
572 uint64_t roundeven_mask
;
576 /* Expand fields based on the size of exponent and fraction */
577 #define FLOAT_PARAMS(E, F) \
579 .exp_bias = ((1 << E) - 1) >> 1, \
580 .exp_max = (1 << E) - 1, \
582 .frac_shift = (-F - 1) & 63, \
583 .frac_lsb = 1ull << ((-F - 1) & 63), \
584 .frac_lsbm1 = 1ull << ((-F - 2) & 63), \
585 .round_mask = (1ull << ((-F - 1) & 63)) - 1, \
586 .roundeven_mask = (2ull << ((-F - 1) & 63)) - 1
588 static const FloatFmt float16_params
= {
592 static const FloatFmt float16_params_ahp
= {
597 static const FloatFmt bfloat16_params
= {
601 static const FloatFmt float32_params
= {
605 static const FloatFmt float64_params
= {
609 static const FloatFmt float128_params
= {
610 FLOAT_PARAMS(15, 112)
613 /* Unpack a float to parts, but do not canonicalize. */
614 static void unpack_raw64(FloatParts64
*r
, const FloatFmt
*fmt
, uint64_t raw
)
616 const int f_size
= fmt
->frac_size
;
617 const int e_size
= fmt
->exp_size
;
619 *r
= (FloatParts64
) {
620 .cls
= float_class_unclassified
,
621 .sign
= extract64(raw
, f_size
+ e_size
, 1),
622 .exp
= extract64(raw
, f_size
, e_size
),
623 .frac
= extract64(raw
, 0, f_size
)
627 static inline void float16_unpack_raw(FloatParts64
*p
, float16 f
)
629 unpack_raw64(p
, &float16_params
, f
);
632 static inline void bfloat16_unpack_raw(FloatParts64
*p
, bfloat16 f
)
634 unpack_raw64(p
, &bfloat16_params
, f
);
637 static inline void float32_unpack_raw(FloatParts64
*p
, float32 f
)
639 unpack_raw64(p
, &float32_params
, f
);
642 static inline void float64_unpack_raw(FloatParts64
*p
, float64 f
)
644 unpack_raw64(p
, &float64_params
, f
);
647 static void float128_unpack_raw(FloatParts128
*p
, float128 f
)
649 const int f_size
= float128_params
.frac_size
- 64;
650 const int e_size
= float128_params
.exp_size
;
652 *p
= (FloatParts128
) {
653 .cls
= float_class_unclassified
,
654 .sign
= extract64(f
.high
, f_size
+ e_size
, 1),
655 .exp
= extract64(f
.high
, f_size
, e_size
),
656 .frac_hi
= extract64(f
.high
, 0, f_size
),
661 /* Pack a float from parts, but do not canonicalize. */
662 static uint64_t pack_raw64(const FloatParts64
*p
, const FloatFmt
*fmt
)
664 const int f_size
= fmt
->frac_size
;
665 const int e_size
= fmt
->exp_size
;
668 ret
= (uint64_t)p
->sign
<< (f_size
+ e_size
);
669 ret
= deposit64(ret
, f_size
, e_size
, p
->exp
);
670 ret
= deposit64(ret
, 0, f_size
, p
->frac
);
674 static inline float16
float16_pack_raw(const FloatParts64
*p
)
676 return make_float16(pack_raw64(p
, &float16_params
));
679 static inline bfloat16
bfloat16_pack_raw(const FloatParts64
*p
)
681 return pack_raw64(p
, &bfloat16_params
);
684 static inline float32
float32_pack_raw(const FloatParts64
*p
)
686 return make_float32(pack_raw64(p
, &float32_params
));
689 static inline float64
float64_pack_raw(const FloatParts64
*p
)
691 return make_float64(pack_raw64(p
, &float64_params
));
694 static float128
float128_pack_raw(const FloatParts128
*p
)
696 const int f_size
= float128_params
.frac_size
- 64;
697 const int e_size
= float128_params
.exp_size
;
700 hi
= (uint64_t)p
->sign
<< (f_size
+ e_size
);
701 hi
= deposit64(hi
, f_size
, e_size
, p
->exp
);
702 hi
= deposit64(hi
, 0, f_size
, p
->frac_hi
);
703 return make_float128(hi
, p
->frac_lo
);
706 /*----------------------------------------------------------------------------
707 | Functions and definitions to determine: (1) whether tininess for underflow
708 | is detected before or after rounding by default, (2) what (if anything)
709 | happens when exceptions are raised, (3) how signaling NaNs are distinguished
710 | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
711 | are propagated from function inputs to output. These details are target-
713 *----------------------------------------------------------------------------*/
714 #include "softfloat-specialize.c.inc"
716 #define PARTS_GENERIC_64_128(NAME, P) \
717 QEMU_GENERIC(P, (FloatParts128 *, parts128_##NAME), parts64_##NAME)
719 #define PARTS_GENERIC_64_128_256(NAME, P) \
720 QEMU_GENERIC(P, (FloatParts256 *, parts256_##NAME), \
721 (FloatParts128 *, parts128_##NAME), parts64_##NAME)
723 #define parts_default_nan(P, S) PARTS_GENERIC_64_128(default_nan, P)(P, S)
724 #define parts_silence_nan(P, S) PARTS_GENERIC_64_128(silence_nan, P)(P, S)
726 static void parts64_return_nan(FloatParts64
*a
, float_status
*s
);
727 static void parts128_return_nan(FloatParts128
*a
, float_status
*s
);
729 #define parts_return_nan(P, S) PARTS_GENERIC_64_128(return_nan, P)(P, S)
731 static FloatParts64
*parts64_pick_nan(FloatParts64
*a
, FloatParts64
*b
,
733 static FloatParts128
*parts128_pick_nan(FloatParts128
*a
, FloatParts128
*b
,
736 #define parts_pick_nan(A, B, S) PARTS_GENERIC_64_128(pick_nan, A)(A, B, S)
738 static FloatParts64
*parts64_pick_nan_muladd(FloatParts64
*a
, FloatParts64
*b
,
739 FloatParts64
*c
, float_status
*s
,
740 int ab_mask
, int abc_mask
);
741 static FloatParts128
*parts128_pick_nan_muladd(FloatParts128
*a
,
745 int ab_mask
, int abc_mask
);
747 #define parts_pick_nan_muladd(A, B, C, S, ABM, ABCM) \
748 PARTS_GENERIC_64_128(pick_nan_muladd, A)(A, B, C, S, ABM, ABCM)
750 static void parts64_canonicalize(FloatParts64
*p
, float_status
*status
,
751 const FloatFmt
*fmt
);
752 static void parts128_canonicalize(FloatParts128
*p
, float_status
*status
,
753 const FloatFmt
*fmt
);
755 #define parts_canonicalize(A, S, F) \
756 PARTS_GENERIC_64_128(canonicalize, A)(A, S, F)
758 static void parts64_uncanon(FloatParts64
*p
, float_status
*status
,
759 const FloatFmt
*fmt
);
760 static void parts128_uncanon(FloatParts128
*p
, float_status
*status
,
761 const FloatFmt
*fmt
);
763 #define parts_uncanon(A, S, F) \
764 PARTS_GENERIC_64_128(uncanon, A)(A, S, F)
766 static void parts64_add_normal(FloatParts64
*a
, FloatParts64
*b
);
767 static void parts128_add_normal(FloatParts128
*a
, FloatParts128
*b
);
768 static void parts256_add_normal(FloatParts256
*a
, FloatParts256
*b
);
770 #define parts_add_normal(A, B) \
771 PARTS_GENERIC_64_128_256(add_normal, A)(A, B)
773 static bool parts64_sub_normal(FloatParts64
*a
, FloatParts64
*b
);
774 static bool parts128_sub_normal(FloatParts128
*a
, FloatParts128
*b
);
775 static bool parts256_sub_normal(FloatParts256
*a
, FloatParts256
*b
);
777 #define parts_sub_normal(A, B) \
778 PARTS_GENERIC_64_128_256(sub_normal, A)(A, B)
780 static FloatParts64
*parts64_addsub(FloatParts64
*a
, FloatParts64
*b
,
781 float_status
*s
, bool subtract
);
782 static FloatParts128
*parts128_addsub(FloatParts128
*a
, FloatParts128
*b
,
783 float_status
*s
, bool subtract
);
785 #define parts_addsub(A, B, S, Z) \
786 PARTS_GENERIC_64_128(addsub, A)(A, B, S, Z)
788 static FloatParts64
*parts64_mul(FloatParts64
*a
, FloatParts64
*b
,
790 static FloatParts128
*parts128_mul(FloatParts128
*a
, FloatParts128
*b
,
793 #define parts_mul(A, B, S) \
794 PARTS_GENERIC_64_128(mul, A)(A, B, S)
796 static FloatParts64
*parts64_muladd(FloatParts64
*a
, FloatParts64
*b
,
797 FloatParts64
*c
, int flags
,
799 static FloatParts128
*parts128_muladd(FloatParts128
*a
, FloatParts128
*b
,
800 FloatParts128
*c
, int flags
,
803 #define parts_muladd(A, B, C, Z, S) \
804 PARTS_GENERIC_64_128(muladd, A)(A, B, C, Z, S)
806 static FloatParts64
*parts64_div(FloatParts64
*a
, FloatParts64
*b
,
808 static FloatParts128
*parts128_div(FloatParts128
*a
, FloatParts128
*b
,
811 #define parts_div(A, B, S) \
812 PARTS_GENERIC_64_128(div, A)(A, B, S)
815 * Helper functions for softfloat-parts.c.inc, per-size operations.
818 #define FRAC_GENERIC_64_128(NAME, P) \
819 QEMU_GENERIC(P, (FloatParts128 *, frac128_##NAME), frac64_##NAME)
821 #define FRAC_GENERIC_64_128_256(NAME, P) \
822 QEMU_GENERIC(P, (FloatParts256 *, frac256_##NAME), \
823 (FloatParts128 *, frac128_##NAME), frac64_##NAME)
825 static bool frac64_add(FloatParts64
*r
, FloatParts64
*a
, FloatParts64
*b
)
827 return uadd64_overflow(a
->frac
, b
->frac
, &r
->frac
);
830 static bool frac128_add(FloatParts128
*r
, FloatParts128
*a
, FloatParts128
*b
)
833 r
->frac_lo
= uadd64_carry(a
->frac_lo
, b
->frac_lo
, &c
);
834 r
->frac_hi
= uadd64_carry(a
->frac_hi
, b
->frac_hi
, &c
);
838 static bool frac256_add(FloatParts256
*r
, FloatParts256
*a
, FloatParts256
*b
)
841 r
->frac_lo
= uadd64_carry(a
->frac_lo
, b
->frac_lo
, &c
);
842 r
->frac_lm
= uadd64_carry(a
->frac_lm
, b
->frac_lm
, &c
);
843 r
->frac_hm
= uadd64_carry(a
->frac_hm
, b
->frac_hm
, &c
);
844 r
->frac_hi
= uadd64_carry(a
->frac_hi
, b
->frac_hi
, &c
);
848 #define frac_add(R, A, B) FRAC_GENERIC_64_128_256(add, R)(R, A, B)
850 static bool frac64_addi(FloatParts64
*r
, FloatParts64
*a
, uint64_t c
)
852 return uadd64_overflow(a
->frac
, c
, &r
->frac
);
855 static bool frac128_addi(FloatParts128
*r
, FloatParts128
*a
, uint64_t c
)
857 c
= uadd64_overflow(a
->frac_lo
, c
, &r
->frac_lo
);
858 return uadd64_overflow(a
->frac_hi
, c
, &r
->frac_hi
);
861 #define frac_addi(R, A, C) FRAC_GENERIC_64_128(addi, R)(R, A, C)
863 static void frac64_allones(FloatParts64
*a
)
868 static void frac128_allones(FloatParts128
*a
)
870 a
->frac_hi
= a
->frac_lo
= -1;
873 #define frac_allones(A) FRAC_GENERIC_64_128(allones, A)(A)
875 static int frac64_cmp(FloatParts64
*a
, FloatParts64
*b
)
877 return a
->frac
== b
->frac
? 0 : a
->frac
< b
->frac
? -1 : 1;
880 static int frac128_cmp(FloatParts128
*a
, FloatParts128
*b
)
882 uint64_t ta
= a
->frac_hi
, tb
= b
->frac_hi
;
884 ta
= a
->frac_lo
, tb
= b
->frac_lo
;
889 return ta
< tb
? -1 : 1;
892 #define frac_cmp(A, B) FRAC_GENERIC_64_128(cmp, A)(A, B)
894 static void frac64_clear(FloatParts64
*a
)
899 static void frac128_clear(FloatParts128
*a
)
901 a
->frac_hi
= a
->frac_lo
= 0;
904 #define frac_clear(A) FRAC_GENERIC_64_128(clear, A)(A)
906 static bool frac64_div(FloatParts64
*a
, FloatParts64
*b
)
908 uint64_t n1
, n0
, r
, q
;
912 * We want a 2*N / N-bit division to produce exactly an N-bit
913 * result, so that we do not lose any precision and so that we
914 * do not have to renormalize afterward. If A.frac < B.frac,
915 * then division would produce an (N-1)-bit result; shift A left
916 * by one to produce the an N-bit result, and return true to
917 * decrement the exponent to match.
919 * The udiv_qrnnd algorithm that we're using requires normalization,
920 * i.e. the msb of the denominator must be set, which is already true.
922 ret
= a
->frac
< b
->frac
;
930 q
= udiv_qrnnd(&r
, n0
, n1
, b
->frac
);
932 /* Set lsb if there is a remainder, to set inexact. */
933 a
->frac
= q
| (r
!= 0);
938 static bool frac128_div(FloatParts128
*a
, FloatParts128
*b
)
940 uint64_t q0
, q1
, a0
, a1
, b0
, b1
;
941 uint64_t r0
, r1
, r2
, r3
, t0
, t1
, t2
, t3
;
944 a0
= a
->frac_hi
, a1
= a
->frac_lo
;
945 b0
= b
->frac_hi
, b1
= b
->frac_lo
;
947 ret
= lt128(a0
, a1
, b0
, b1
);
949 a1
= shr_double(a0
, a1
, 1);
953 /* Use 128/64 -> 64 division as estimate for 192/128 -> 128 division. */
954 q0
= estimateDiv128To64(a0
, a1
, b0
);
957 * Estimate is high because B1 was not included (unless B1 == 0).
958 * Reduce quotient and increase remainder until remainder is non-negative.
959 * This loop will execute 0 to 2 times.
961 mul128By64To192(b0
, b1
, q0
, &t0
, &t1
, &t2
);
962 sub192(a0
, a1
, 0, t0
, t1
, t2
, &r0
, &r1
, &r2
);
965 add192(r0
, r1
, r2
, 0, b0
, b1
, &r0
, &r1
, &r2
);
968 /* Repeat using the remainder, producing a second word of quotient. */
969 q1
= estimateDiv128To64(r1
, r2
, b0
);
970 mul128By64To192(b0
, b1
, q1
, &t1
, &t2
, &t3
);
971 sub192(r1
, r2
, 0, t1
, t2
, t3
, &r1
, &r2
, &r3
);
974 add192(r1
, r2
, r3
, 0, b0
, b1
, &r1
, &r2
, &r3
);
977 /* Any remainder indicates inexact; set sticky bit. */
978 q1
|= (r2
| r3
) != 0;
985 #define frac_div(A, B) FRAC_GENERIC_64_128(div, A)(A, B)
987 static bool frac64_eqz(FloatParts64
*a
)
992 static bool frac128_eqz(FloatParts128
*a
)
994 return (a
->frac_hi
| a
->frac_lo
) == 0;
997 #define frac_eqz(A) FRAC_GENERIC_64_128(eqz, A)(A)
999 static void frac64_mulw(FloatParts128
*r
, FloatParts64
*a
, FloatParts64
*b
)
1001 mulu64(&r
->frac_lo
, &r
->frac_hi
, a
->frac
, b
->frac
);
1004 static void frac128_mulw(FloatParts256
*r
, FloatParts128
*a
, FloatParts128
*b
)
1006 mul128To256(a
->frac_hi
, a
->frac_lo
, b
->frac_hi
, b
->frac_lo
,
1007 &r
->frac_hi
, &r
->frac_hm
, &r
->frac_lm
, &r
->frac_lo
);
1010 #define frac_mulw(R, A, B) FRAC_GENERIC_64_128(mulw, A)(R, A, B)
1012 static void frac64_neg(FloatParts64
*a
)
1017 static void frac128_neg(FloatParts128
*a
)
1020 a
->frac_lo
= usub64_borrow(0, a
->frac_lo
, &c
);
1021 a
->frac_hi
= usub64_borrow(0, a
->frac_hi
, &c
);
1024 static void frac256_neg(FloatParts256
*a
)
1027 a
->frac_lo
= usub64_borrow(0, a
->frac_lo
, &c
);
1028 a
->frac_lm
= usub64_borrow(0, a
->frac_lm
, &c
);
1029 a
->frac_hm
= usub64_borrow(0, a
->frac_hm
, &c
);
1030 a
->frac_hi
= usub64_borrow(0, a
->frac_hi
, &c
);
1033 #define frac_neg(A) FRAC_GENERIC_64_128_256(neg, A)(A)
1035 static int frac64_normalize(FloatParts64
*a
)
1038 int shift
= clz64(a
->frac
);
1045 static int frac128_normalize(FloatParts128
*a
)
1048 int shl
= clz64(a
->frac_hi
);
1049 a
->frac_hi
= shl_double(a
->frac_hi
, a
->frac_lo
, shl
);
1052 } else if (a
->frac_lo
) {
1053 int shl
= clz64(a
->frac_lo
);
1054 a
->frac_hi
= a
->frac_lo
<< shl
;
1061 static int frac256_normalize(FloatParts256
*a
)
1063 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_hm
;
1064 uint64_t a2
= a
->frac_lm
, a3
= a
->frac_lo
;
1076 a0
= a1
, a1
= a2
, a2
= a3
, a3
= 0;
1079 a0
= a2
, a1
= a3
, a2
= 0, a3
= 0;
1082 a0
= a3
, a1
= 0, a2
= 0, a3
= 0;
1085 a0
= 0, a1
= 0, a2
= 0, a3
= 0;
1095 a0
= shl_double(a0
, a1
, shl
);
1096 a1
= shl_double(a1
, a2
, shl
);
1097 a2
= shl_double(a2
, a3
, shl
);
1108 #define frac_normalize(A) FRAC_GENERIC_64_128_256(normalize, A)(A)
1110 static void frac64_shl(FloatParts64
*a
, int c
)
1115 static void frac128_shl(FloatParts128
*a
, int c
)
1117 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_lo
;
1125 a0
= shl_double(a0
, a1
, c
);
1133 #define frac_shl(A, C) FRAC_GENERIC_64_128(shl, A)(A, C)
1135 static void frac64_shr(FloatParts64
*a
, int c
)
1140 static void frac128_shr(FloatParts128
*a
, int c
)
1142 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_lo
;
1150 a1
= shr_double(a0
, a1
, c
);
1158 #define frac_shr(A, C) FRAC_GENERIC_64_128(shr, A)(A, C)
1160 static void frac64_shrjam(FloatParts64
*a
, int c
)
1162 uint64_t a0
= a
->frac
;
1164 if (likely(c
!= 0)) {
1165 if (likely(c
< 64)) {
1166 a0
= (a0
>> c
) | (shr_double(a0
, 0, c
) != 0);
1174 static void frac128_shrjam(FloatParts128
*a
, int c
)
1176 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_lo
;
1177 uint64_t sticky
= 0;
1179 if (unlikely(c
== 0)) {
1181 } else if (likely(c
< 64)) {
1183 } else if (likely(c
< 128)) {
1197 sticky
|= shr_double(a1
, 0, c
);
1198 a1
= shr_double(a0
, a1
, c
);
1202 a
->frac_lo
= a1
| (sticky
!= 0);
1206 static void frac256_shrjam(FloatParts256
*a
, int c
)
1208 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_hm
;
1209 uint64_t a2
= a
->frac_lm
, a3
= a
->frac_lo
;
1210 uint64_t sticky
= 0;
1212 if (unlikely(c
== 0)) {
1214 } else if (likely(c
< 64)) {
1216 } else if (likely(c
< 256)) {
1217 if (unlikely(c
& 128)) {
1219 a3
= a1
, a2
= a0
, a1
= 0, a0
= 0;
1221 if (unlikely(c
& 64)) {
1223 a3
= a2
, a2
= a1
, a1
= a0
, a0
= 0;
1230 sticky
= a0
| a1
| a2
| a3
;
1231 a0
= a1
= a2
= a3
= 0;
1235 sticky
|= shr_double(a3
, 0, c
);
1236 a3
= shr_double(a2
, a3
, c
);
1237 a2
= shr_double(a1
, a2
, c
);
1238 a1
= shr_double(a0
, a1
, c
);
1242 a
->frac_lo
= a3
| (sticky
!= 0);
1248 #define frac_shrjam(A, C) FRAC_GENERIC_64_128_256(shrjam, A)(A, C)
1250 static bool frac64_sub(FloatParts64
*r
, FloatParts64
*a
, FloatParts64
*b
)
1252 return usub64_overflow(a
->frac
, b
->frac
, &r
->frac
);
1255 static bool frac128_sub(FloatParts128
*r
, FloatParts128
*a
, FloatParts128
*b
)
1258 r
->frac_lo
= usub64_borrow(a
->frac_lo
, b
->frac_lo
, &c
);
1259 r
->frac_hi
= usub64_borrow(a
->frac_hi
, b
->frac_hi
, &c
);
1263 static bool frac256_sub(FloatParts256
*r
, FloatParts256
*a
, FloatParts256
*b
)
1266 r
->frac_lo
= usub64_borrow(a
->frac_lo
, b
->frac_lo
, &c
);
1267 r
->frac_lm
= usub64_borrow(a
->frac_lm
, b
->frac_lm
, &c
);
1268 r
->frac_hm
= usub64_borrow(a
->frac_hm
, b
->frac_hm
, &c
);
1269 r
->frac_hi
= usub64_borrow(a
->frac_hi
, b
->frac_hi
, &c
);
1273 #define frac_sub(R, A, B) FRAC_GENERIC_64_128_256(sub, R)(R, A, B)
1275 static void frac64_truncjam(FloatParts64
*r
, FloatParts128
*a
)
1277 r
->frac
= a
->frac_hi
| (a
->frac_lo
!= 0);
1280 static void frac128_truncjam(FloatParts128
*r
, FloatParts256
*a
)
1282 r
->frac_hi
= a
->frac_hi
;
1283 r
->frac_lo
= a
->frac_hm
| ((a
->frac_lm
| a
->frac_lo
) != 0);
1286 #define frac_truncjam(R, A) FRAC_GENERIC_64_128(truncjam, R)(R, A)
1288 static void frac64_widen(FloatParts128
*r
, FloatParts64
*a
)
1290 r
->frac_hi
= a
->frac
;
1294 static void frac128_widen(FloatParts256
*r
, FloatParts128
*a
)
1296 r
->frac_hi
= a
->frac_hi
;
1297 r
->frac_hm
= a
->frac_lo
;
1302 #define frac_widen(A, B) FRAC_GENERIC_64_128(widen, B)(A, B)
1304 #define partsN(NAME) glue(glue(glue(parts,N),_),NAME)
1305 #define FloatPartsN glue(FloatParts,N)
1306 #define FloatPartsW glue(FloatParts,W)
1311 #include "softfloat-parts-addsub.c.inc"
1312 #include "softfloat-parts.c.inc"
1319 #include "softfloat-parts-addsub.c.inc"
1320 #include "softfloat-parts.c.inc"
1326 #include "softfloat-parts-addsub.c.inc"
1335 * Pack/unpack routines with a specific FloatFmt.
1338 static void float16a_unpack_canonical(FloatParts64
*p
, float16 f
,
1339 float_status
*s
, const FloatFmt
*params
)
1341 float16_unpack_raw(p
, f
);
1342 parts_canonicalize(p
, s
, params
);
1345 static void float16_unpack_canonical(FloatParts64
*p
, float16 f
,
1348 float16a_unpack_canonical(p
, f
, s
, &float16_params
);
1351 static void bfloat16_unpack_canonical(FloatParts64
*p
, bfloat16 f
,
1354 bfloat16_unpack_raw(p
, f
);
1355 parts_canonicalize(p
, s
, &bfloat16_params
);
1358 static float16
float16a_round_pack_canonical(FloatParts64
*p
,
1360 const FloatFmt
*params
)
1362 parts_uncanon(p
, s
, params
);
1363 return float16_pack_raw(p
);
1366 static float16
float16_round_pack_canonical(FloatParts64
*p
,
1369 return float16a_round_pack_canonical(p
, s
, &float16_params
);
1372 static bfloat16
bfloat16_round_pack_canonical(FloatParts64
*p
,
1375 parts_uncanon(p
, s
, &bfloat16_params
);
1376 return bfloat16_pack_raw(p
);
1379 static void float32_unpack_canonical(FloatParts64
*p
, float32 f
,
1382 float32_unpack_raw(p
, f
);
1383 parts_canonicalize(p
, s
, &float32_params
);
1386 static float32
float32_round_pack_canonical(FloatParts64
*p
,
1389 parts_uncanon(p
, s
, &float32_params
);
1390 return float32_pack_raw(p
);
1393 static void float64_unpack_canonical(FloatParts64
*p
, float64 f
,
1396 float64_unpack_raw(p
, f
);
1397 parts_canonicalize(p
, s
, &float64_params
);
1400 static float64
float64_round_pack_canonical(FloatParts64
*p
,
1403 parts_uncanon(p
, s
, &float64_params
);
1404 return float64_pack_raw(p
);
1407 static void float128_unpack_canonical(FloatParts128
*p
, float128 f
,
1410 float128_unpack_raw(p
, f
);
1411 parts_canonicalize(p
, s
, &float128_params
);
1414 static float128
float128_round_pack_canonical(FloatParts128
*p
,
1417 parts_uncanon(p
, s
, &float128_params
);
1418 return float128_pack_raw(p
);
1422 * Addition and subtraction
1425 static float16 QEMU_FLATTEN
1426 float16_addsub(float16 a
, float16 b
, float_status
*status
, bool subtract
)
1428 FloatParts64 pa
, pb
, *pr
;
1430 float16_unpack_canonical(&pa
, a
, status
);
1431 float16_unpack_canonical(&pb
, b
, status
);
1432 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1434 return float16_round_pack_canonical(pr
, status
);
1437 float16
float16_add(float16 a
, float16 b
, float_status
*status
)
1439 return float16_addsub(a
, b
, status
, false);
1442 float16
float16_sub(float16 a
, float16 b
, float_status
*status
)
1444 return float16_addsub(a
, b
, status
, true);
1447 static float32 QEMU_SOFTFLOAT_ATTR
1448 soft_f32_addsub(float32 a
, float32 b
, float_status
*status
, bool subtract
)
1450 FloatParts64 pa
, pb
, *pr
;
1452 float32_unpack_canonical(&pa
, a
, status
);
1453 float32_unpack_canonical(&pb
, b
, status
);
1454 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1456 return float32_round_pack_canonical(pr
, status
);
1459 static float32
soft_f32_add(float32 a
, float32 b
, float_status
*status
)
1461 return soft_f32_addsub(a
, b
, status
, false);
1464 static float32
soft_f32_sub(float32 a
, float32 b
, float_status
*status
)
1466 return soft_f32_addsub(a
, b
, status
, true);
1469 static float64 QEMU_SOFTFLOAT_ATTR
1470 soft_f64_addsub(float64 a
, float64 b
, float_status
*status
, bool subtract
)
1472 FloatParts64 pa
, pb
, *pr
;
1474 float64_unpack_canonical(&pa
, a
, status
);
1475 float64_unpack_canonical(&pb
, b
, status
);
1476 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1478 return float64_round_pack_canonical(pr
, status
);
1481 static float64
soft_f64_add(float64 a
, float64 b
, float_status
*status
)
1483 return soft_f64_addsub(a
, b
, status
, false);
1486 static float64
soft_f64_sub(float64 a
, float64 b
, float_status
*status
)
1488 return soft_f64_addsub(a
, b
, status
, true);
1491 static float hard_f32_add(float a
, float b
)
1496 static float hard_f32_sub(float a
, float b
)
1501 static double hard_f64_add(double a
, double b
)
1506 static double hard_f64_sub(double a
, double b
)
1511 static bool f32_addsubmul_post(union_float32 a
, union_float32 b
)
1513 if (QEMU_HARDFLOAT_2F32_USE_FP
) {
1514 return !(fpclassify(a
.h
) == FP_ZERO
&& fpclassify(b
.h
) == FP_ZERO
);
1516 return !(float32_is_zero(a
.s
) && float32_is_zero(b
.s
));
1519 static bool f64_addsubmul_post(union_float64 a
, union_float64 b
)
1521 if (QEMU_HARDFLOAT_2F64_USE_FP
) {
1522 return !(fpclassify(a
.h
) == FP_ZERO
&& fpclassify(b
.h
) == FP_ZERO
);
1524 return !(float64_is_zero(a
.s
) && float64_is_zero(b
.s
));
1528 static float32
float32_addsub(float32 a
, float32 b
, float_status
*s
,
1529 hard_f32_op2_fn hard
, soft_f32_op2_fn soft
)
1531 return float32_gen2(a
, b
, s
, hard
, soft
,
1532 f32_is_zon2
, f32_addsubmul_post
);
1535 static float64
float64_addsub(float64 a
, float64 b
, float_status
*s
,
1536 hard_f64_op2_fn hard
, soft_f64_op2_fn soft
)
1538 return float64_gen2(a
, b
, s
, hard
, soft
,
1539 f64_is_zon2
, f64_addsubmul_post
);
1542 float32 QEMU_FLATTEN
1543 float32_add(float32 a
, float32 b
, float_status
*s
)
1545 return float32_addsub(a
, b
, s
, hard_f32_add
, soft_f32_add
);
1548 float32 QEMU_FLATTEN
1549 float32_sub(float32 a
, float32 b
, float_status
*s
)
1551 return float32_addsub(a
, b
, s
, hard_f32_sub
, soft_f32_sub
);
1554 float64 QEMU_FLATTEN
1555 float64_add(float64 a
, float64 b
, float_status
*s
)
1557 return float64_addsub(a
, b
, s
, hard_f64_add
, soft_f64_add
);
1560 float64 QEMU_FLATTEN
1561 float64_sub(float64 a
, float64 b
, float_status
*s
)
1563 return float64_addsub(a
, b
, s
, hard_f64_sub
, soft_f64_sub
);
1566 static bfloat16 QEMU_FLATTEN
1567 bfloat16_addsub(bfloat16 a
, bfloat16 b
, float_status
*status
, bool subtract
)
1569 FloatParts64 pa
, pb
, *pr
;
1571 bfloat16_unpack_canonical(&pa
, a
, status
);
1572 bfloat16_unpack_canonical(&pb
, b
, status
);
1573 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1575 return bfloat16_round_pack_canonical(pr
, status
);
1578 bfloat16
bfloat16_add(bfloat16 a
, bfloat16 b
, float_status
*status
)
1580 return bfloat16_addsub(a
, b
, status
, false);
1583 bfloat16
bfloat16_sub(bfloat16 a
, bfloat16 b
, float_status
*status
)
1585 return bfloat16_addsub(a
, b
, status
, true);
1588 static float128 QEMU_FLATTEN
1589 float128_addsub(float128 a
, float128 b
, float_status
*status
, bool subtract
)
1591 FloatParts128 pa
, pb
, *pr
;
1593 float128_unpack_canonical(&pa
, a
, status
);
1594 float128_unpack_canonical(&pb
, b
, status
);
1595 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1597 return float128_round_pack_canonical(pr
, status
);
1600 float128
float128_add(float128 a
, float128 b
, float_status
*status
)
1602 return float128_addsub(a
, b
, status
, false);
1605 float128
float128_sub(float128 a
, float128 b
, float_status
*status
)
1607 return float128_addsub(a
, b
, status
, true);
1614 float16 QEMU_FLATTEN
float16_mul(float16 a
, float16 b
, float_status
*status
)
1616 FloatParts64 pa
, pb
, *pr
;
1618 float16_unpack_canonical(&pa
, a
, status
);
1619 float16_unpack_canonical(&pb
, b
, status
);
1620 pr
= parts_mul(&pa
, &pb
, status
);
1622 return float16_round_pack_canonical(pr
, status
);
1625 static float32 QEMU_SOFTFLOAT_ATTR
1626 soft_f32_mul(float32 a
, float32 b
, float_status
*status
)
1628 FloatParts64 pa
, pb
, *pr
;
1630 float32_unpack_canonical(&pa
, a
, status
);
1631 float32_unpack_canonical(&pb
, b
, status
);
1632 pr
= parts_mul(&pa
, &pb
, status
);
1634 return float32_round_pack_canonical(pr
, status
);
1637 static float64 QEMU_SOFTFLOAT_ATTR
1638 soft_f64_mul(float64 a
, float64 b
, float_status
*status
)
1640 FloatParts64 pa
, pb
, *pr
;
1642 float64_unpack_canonical(&pa
, a
, status
);
1643 float64_unpack_canonical(&pb
, b
, status
);
1644 pr
= parts_mul(&pa
, &pb
, status
);
1646 return float64_round_pack_canonical(pr
, status
);
1649 static float hard_f32_mul(float a
, float b
)
1654 static double hard_f64_mul(double a
, double b
)
1659 float32 QEMU_FLATTEN
1660 float32_mul(float32 a
, float32 b
, float_status
*s
)
1662 return float32_gen2(a
, b
, s
, hard_f32_mul
, soft_f32_mul
,
1663 f32_is_zon2
, f32_addsubmul_post
);
1666 float64 QEMU_FLATTEN
1667 float64_mul(float64 a
, float64 b
, float_status
*s
)
1669 return float64_gen2(a
, b
, s
, hard_f64_mul
, soft_f64_mul
,
1670 f64_is_zon2
, f64_addsubmul_post
);
1673 bfloat16 QEMU_FLATTEN
1674 bfloat16_mul(bfloat16 a
, bfloat16 b
, float_status
*status
)
1676 FloatParts64 pa
, pb
, *pr
;
1678 bfloat16_unpack_canonical(&pa
, a
, status
);
1679 bfloat16_unpack_canonical(&pb
, b
, status
);
1680 pr
= parts_mul(&pa
, &pb
, status
);
1682 return bfloat16_round_pack_canonical(pr
, status
);
1685 float128 QEMU_FLATTEN
1686 float128_mul(float128 a
, float128 b
, float_status
*status
)
1688 FloatParts128 pa
, pb
, *pr
;
1690 float128_unpack_canonical(&pa
, a
, status
);
1691 float128_unpack_canonical(&pb
, b
, status
);
1692 pr
= parts_mul(&pa
, &pb
, status
);
1694 return float128_round_pack_canonical(pr
, status
);
1698 * Fused multiply-add
1701 float16 QEMU_FLATTEN
float16_muladd(float16 a
, float16 b
, float16 c
,
1702 int flags
, float_status
*status
)
1704 FloatParts64 pa
, pb
, pc
, *pr
;
1706 float16_unpack_canonical(&pa
, a
, status
);
1707 float16_unpack_canonical(&pb
, b
, status
);
1708 float16_unpack_canonical(&pc
, c
, status
);
1709 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1711 return float16_round_pack_canonical(pr
, status
);
1714 static float32 QEMU_SOFTFLOAT_ATTR
1715 soft_f32_muladd(float32 a
, float32 b
, float32 c
, int flags
,
1716 float_status
*status
)
1718 FloatParts64 pa
, pb
, pc
, *pr
;
1720 float32_unpack_canonical(&pa
, a
, status
);
1721 float32_unpack_canonical(&pb
, b
, status
);
1722 float32_unpack_canonical(&pc
, c
, status
);
1723 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1725 return float32_round_pack_canonical(pr
, status
);
1728 static float64 QEMU_SOFTFLOAT_ATTR
1729 soft_f64_muladd(float64 a
, float64 b
, float64 c
, int flags
,
1730 float_status
*status
)
1732 FloatParts64 pa
, pb
, pc
, *pr
;
1734 float64_unpack_canonical(&pa
, a
, status
);
1735 float64_unpack_canonical(&pb
, b
, status
);
1736 float64_unpack_canonical(&pc
, c
, status
);
1737 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1739 return float64_round_pack_canonical(pr
, status
);
1742 static bool force_soft_fma
;
1744 float32 QEMU_FLATTEN
1745 float32_muladd(float32 xa
, float32 xb
, float32 xc
, int flags
, float_status
*s
)
1747 union_float32 ua
, ub
, uc
, ur
;
1753 if (unlikely(!can_use_fpu(s
))) {
1756 if (unlikely(flags
& float_muladd_halve_result
)) {
1760 float32_input_flush3(&ua
.s
, &ub
.s
, &uc
.s
, s
);
1761 if (unlikely(!f32_is_zon3(ua
, ub
, uc
))) {
1765 if (unlikely(force_soft_fma
)) {
1770 * When (a || b) == 0, there's no need to check for under/over flow,
1771 * since we know the addend is (normal || 0) and the product is 0.
1773 if (float32_is_zero(ua
.s
) || float32_is_zero(ub
.s
)) {
1777 prod_sign
= float32_is_neg(ua
.s
) ^ float32_is_neg(ub
.s
);
1778 prod_sign
^= !!(flags
& float_muladd_negate_product
);
1779 up
.s
= float32_set_sign(float32_zero
, prod_sign
);
1781 if (flags
& float_muladd_negate_c
) {
1786 union_float32 ua_orig
= ua
;
1787 union_float32 uc_orig
= uc
;
1789 if (flags
& float_muladd_negate_product
) {
1792 if (flags
& float_muladd_negate_c
) {
1796 ur
.h
= fmaf(ua
.h
, ub
.h
, uc
.h
);
1798 if (unlikely(f32_is_inf(ur
))) {
1799 float_raise(float_flag_overflow
, s
);
1800 } else if (unlikely(fabsf(ur
.h
) <= FLT_MIN
)) {
1806 if (flags
& float_muladd_negate_result
) {
1807 return float32_chs(ur
.s
);
1812 return soft_f32_muladd(ua
.s
, ub
.s
, uc
.s
, flags
, s
);
1815 float64 QEMU_FLATTEN
1816 float64_muladd(float64 xa
, float64 xb
, float64 xc
, int flags
, float_status
*s
)
1818 union_float64 ua
, ub
, uc
, ur
;
1824 if (unlikely(!can_use_fpu(s
))) {
1827 if (unlikely(flags
& float_muladd_halve_result
)) {
1831 float64_input_flush3(&ua
.s
, &ub
.s
, &uc
.s
, s
);
1832 if (unlikely(!f64_is_zon3(ua
, ub
, uc
))) {
1836 if (unlikely(force_soft_fma
)) {
1841 * When (a || b) == 0, there's no need to check for under/over flow,
1842 * since we know the addend is (normal || 0) and the product is 0.
1844 if (float64_is_zero(ua
.s
) || float64_is_zero(ub
.s
)) {
1848 prod_sign
= float64_is_neg(ua
.s
) ^ float64_is_neg(ub
.s
);
1849 prod_sign
^= !!(flags
& float_muladd_negate_product
);
1850 up
.s
= float64_set_sign(float64_zero
, prod_sign
);
1852 if (flags
& float_muladd_negate_c
) {
1857 union_float64 ua_orig
= ua
;
1858 union_float64 uc_orig
= uc
;
1860 if (flags
& float_muladd_negate_product
) {
1863 if (flags
& float_muladd_negate_c
) {
1867 ur
.h
= fma(ua
.h
, ub
.h
, uc
.h
);
1869 if (unlikely(f64_is_inf(ur
))) {
1870 float_raise(float_flag_overflow
, s
);
1871 } else if (unlikely(fabs(ur
.h
) <= FLT_MIN
)) {
1877 if (flags
& float_muladd_negate_result
) {
1878 return float64_chs(ur
.s
);
1883 return soft_f64_muladd(ua
.s
, ub
.s
, uc
.s
, flags
, s
);
1886 bfloat16 QEMU_FLATTEN
bfloat16_muladd(bfloat16 a
, bfloat16 b
, bfloat16 c
,
1887 int flags
, float_status
*status
)
1889 FloatParts64 pa
, pb
, pc
, *pr
;
1891 bfloat16_unpack_canonical(&pa
, a
, status
);
1892 bfloat16_unpack_canonical(&pb
, b
, status
);
1893 bfloat16_unpack_canonical(&pc
, c
, status
);
1894 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1896 return bfloat16_round_pack_canonical(pr
, status
);
1899 float128 QEMU_FLATTEN
float128_muladd(float128 a
, float128 b
, float128 c
,
1900 int flags
, float_status
*status
)
1902 FloatParts128 pa
, pb
, pc
, *pr
;
1904 float128_unpack_canonical(&pa
, a
, status
);
1905 float128_unpack_canonical(&pb
, b
, status
);
1906 float128_unpack_canonical(&pc
, c
, status
);
1907 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1909 return float128_round_pack_canonical(pr
, status
);
1916 float16
float16_div(float16 a
, float16 b
, float_status
*status
)
1918 FloatParts64 pa
, pb
, *pr
;
1920 float16_unpack_canonical(&pa
, a
, status
);
1921 float16_unpack_canonical(&pb
, b
, status
);
1922 pr
= parts_div(&pa
, &pb
, status
);
1924 return float16_round_pack_canonical(pr
, status
);
1927 static float32 QEMU_SOFTFLOAT_ATTR
1928 soft_f32_div(float32 a
, float32 b
, float_status
*status
)
1930 FloatParts64 pa
, pb
, *pr
;
1932 float32_unpack_canonical(&pa
, a
, status
);
1933 float32_unpack_canonical(&pb
, b
, status
);
1934 pr
= parts_div(&pa
, &pb
, status
);
1936 return float32_round_pack_canonical(pr
, status
);
1939 static float64 QEMU_SOFTFLOAT_ATTR
1940 soft_f64_div(float64 a
, float64 b
, float_status
*status
)
1942 FloatParts64 pa
, pb
, *pr
;
1944 float64_unpack_canonical(&pa
, a
, status
);
1945 float64_unpack_canonical(&pb
, b
, status
);
1946 pr
= parts_div(&pa
, &pb
, status
);
1948 return float64_round_pack_canonical(pr
, status
);
1951 static float hard_f32_div(float a
, float b
)
1956 static double hard_f64_div(double a
, double b
)
1961 static bool f32_div_pre(union_float32 a
, union_float32 b
)
1963 if (QEMU_HARDFLOAT_2F32_USE_FP
) {
1964 return (fpclassify(a
.h
) == FP_NORMAL
|| fpclassify(a
.h
) == FP_ZERO
) &&
1965 fpclassify(b
.h
) == FP_NORMAL
;
1967 return float32_is_zero_or_normal(a
.s
) && float32_is_normal(b
.s
);
1970 static bool f64_div_pre(union_float64 a
, union_float64 b
)
1972 if (QEMU_HARDFLOAT_2F64_USE_FP
) {
1973 return (fpclassify(a
.h
) == FP_NORMAL
|| fpclassify(a
.h
) == FP_ZERO
) &&
1974 fpclassify(b
.h
) == FP_NORMAL
;
1976 return float64_is_zero_or_normal(a
.s
) && float64_is_normal(b
.s
);
1979 static bool f32_div_post(union_float32 a
, union_float32 b
)
1981 if (QEMU_HARDFLOAT_2F32_USE_FP
) {
1982 return fpclassify(a
.h
) != FP_ZERO
;
1984 return !float32_is_zero(a
.s
);
1987 static bool f64_div_post(union_float64 a
, union_float64 b
)
1989 if (QEMU_HARDFLOAT_2F64_USE_FP
) {
1990 return fpclassify(a
.h
) != FP_ZERO
;
1992 return !float64_is_zero(a
.s
);
1995 float32 QEMU_FLATTEN
1996 float32_div(float32 a
, float32 b
, float_status
*s
)
1998 return float32_gen2(a
, b
, s
, hard_f32_div
, soft_f32_div
,
1999 f32_div_pre
, f32_div_post
);
2002 float64 QEMU_FLATTEN
2003 float64_div(float64 a
, float64 b
, float_status
*s
)
2005 return float64_gen2(a
, b
, s
, hard_f64_div
, soft_f64_div
,
2006 f64_div_pre
, f64_div_post
);
2009 bfloat16 QEMU_FLATTEN
2010 bfloat16_div(bfloat16 a
, bfloat16 b
, float_status
*status
)
2012 FloatParts64 pa
, pb
, *pr
;
2014 bfloat16_unpack_canonical(&pa
, a
, status
);
2015 bfloat16_unpack_canonical(&pb
, b
, status
);
2016 pr
= parts_div(&pa
, &pb
, status
);
2018 return bfloat16_round_pack_canonical(pr
, status
);
2021 float128 QEMU_FLATTEN
2022 float128_div(float128 a
, float128 b
, float_status
*status
)
2024 FloatParts128 pa
, pb
, *pr
;
2026 float128_unpack_canonical(&pa
, a
, status
);
2027 float128_unpack_canonical(&pb
, b
, status
);
2028 pr
= parts_div(&pa
, &pb
, status
);
2030 return float128_round_pack_canonical(pr
, status
);
2034 * Float to Float conversions
2036 * Returns the result of converting one float format to another. The
2037 * conversion is performed according to the IEC/IEEE Standard for
2038 * Binary Floating-Point Arithmetic.
2040 * Usually this only needs to take care of raising invalid exceptions
2041 * and handling the conversion on NaNs.
2044 static void parts_float_to_ahp(FloatParts64
*a
, float_status
*s
)
2047 case float_class_qnan
:
2048 case float_class_snan
:
2050 * There is no NaN in the destination format. Raise Invalid
2051 * and return a zero with the sign of the input NaN.
2053 float_raise(float_flag_invalid
, s
);
2054 a
->cls
= float_class_zero
;
2057 case float_class_inf
:
2059 * There is no Inf in the destination format. Raise Invalid
2060 * and return the maximum normal with the correct sign.
2062 float_raise(float_flag_invalid
, s
);
2063 a
->cls
= float_class_normal
;
2064 a
->exp
= float16_params_ahp
.exp_max
;
2065 a
->frac
= MAKE_64BIT_MASK(float16_params_ahp
.frac_shift
,
2066 float16_params_ahp
.frac_size
+ 1);
2069 case float_class_normal
:
2070 case float_class_zero
:
2074 g_assert_not_reached();
2078 static void parts64_float_to_float(FloatParts64
*a
, float_status
*s
)
2080 if (is_nan(a
->cls
)) {
2081 parts_return_nan(a
, s
);
2085 static void parts128_float_to_float(FloatParts128
*a
, float_status
*s
)
2087 if (is_nan(a
->cls
)) {
2088 parts_return_nan(a
, s
);
2092 #define parts_float_to_float(P, S) \
2093 PARTS_GENERIC_64_128(float_to_float, P)(P, S)
2095 static void parts_float_to_float_narrow(FloatParts64
*a
, FloatParts128
*b
,
2102 if (a
->cls
== float_class_normal
) {
2103 frac_truncjam(a
, b
);
2104 } else if (is_nan(a
->cls
)) {
2105 /* Discard the low bits of the NaN. */
2106 a
->frac
= b
->frac_hi
;
2107 parts_return_nan(a
, s
);
2111 static void parts_float_to_float_widen(FloatParts128
*a
, FloatParts64
*b
,
2119 if (is_nan(a
->cls
)) {
2120 parts_return_nan(a
, s
);
2124 float32
float16_to_float32(float16 a
, bool ieee
, float_status
*s
)
2126 const FloatFmt
*fmt16
= ieee
? &float16_params
: &float16_params_ahp
;
2129 float16a_unpack_canonical(&p
, a
, s
, fmt16
);
2130 parts_float_to_float(&p
, s
);
2131 return float32_round_pack_canonical(&p
, s
);
2134 float64
float16_to_float64(float16 a
, bool ieee
, float_status
*s
)
2136 const FloatFmt
*fmt16
= ieee
? &float16_params
: &float16_params_ahp
;
2139 float16a_unpack_canonical(&p
, a
, s
, fmt16
);
2140 parts_float_to_float(&p
, s
);
2141 return float64_round_pack_canonical(&p
, s
);
2144 float16
float32_to_float16(float32 a
, bool ieee
, float_status
*s
)
2147 const FloatFmt
*fmt
;
2149 float32_unpack_canonical(&p
, a
, s
);
2151 parts_float_to_float(&p
, s
);
2152 fmt
= &float16_params
;
2154 parts_float_to_ahp(&p
, s
);
2155 fmt
= &float16_params_ahp
;
2157 return float16a_round_pack_canonical(&p
, s
, fmt
);
2160 static float64 QEMU_SOFTFLOAT_ATTR
2161 soft_float32_to_float64(float32 a
, float_status
*s
)
2165 float32_unpack_canonical(&p
, a
, s
);
2166 parts_float_to_float(&p
, s
);
2167 return float64_round_pack_canonical(&p
, s
);
2170 float64
float32_to_float64(float32 a
, float_status
*s
)
2172 if (likely(float32_is_normal(a
))) {
2173 /* Widening conversion can never produce inexact results. */
2179 } else if (float32_is_zero(a
)) {
2180 return float64_set_sign(float64_zero
, float32_is_neg(a
));
2182 return soft_float32_to_float64(a
, s
);
2186 float16
float64_to_float16(float64 a
, bool ieee
, float_status
*s
)
2189 const FloatFmt
*fmt
;
2191 float64_unpack_canonical(&p
, a
, s
);
2193 parts_float_to_float(&p
, s
);
2194 fmt
= &float16_params
;
2196 parts_float_to_ahp(&p
, s
);
2197 fmt
= &float16_params_ahp
;
2199 return float16a_round_pack_canonical(&p
, s
, fmt
);
2202 float32
float64_to_float32(float64 a
, float_status
*s
)
2206 float64_unpack_canonical(&p
, a
, s
);
2207 parts_float_to_float(&p
, s
);
2208 return float32_round_pack_canonical(&p
, s
);
2211 float32
bfloat16_to_float32(bfloat16 a
, float_status
*s
)
2215 bfloat16_unpack_canonical(&p
, a
, s
);
2216 parts_float_to_float(&p
, s
);
2217 return float32_round_pack_canonical(&p
, s
);
2220 float64
bfloat16_to_float64(bfloat16 a
, float_status
*s
)
2224 bfloat16_unpack_canonical(&p
, a
, s
);
2225 parts_float_to_float(&p
, s
);
2226 return float64_round_pack_canonical(&p
, s
);
2229 bfloat16
float32_to_bfloat16(float32 a
, float_status
*s
)
2233 float32_unpack_canonical(&p
, a
, s
);
2234 parts_float_to_float(&p
, s
);
2235 return bfloat16_round_pack_canonical(&p
, s
);
2238 bfloat16
float64_to_bfloat16(float64 a
, float_status
*s
)
2242 float64_unpack_canonical(&p
, a
, s
);
2243 parts_float_to_float(&p
, s
);
2244 return bfloat16_round_pack_canonical(&p
, s
);
2247 float32
float128_to_float32(float128 a
, float_status
*s
)
2252 float128_unpack_canonical(&p128
, a
, s
);
2253 parts_float_to_float_narrow(&p64
, &p128
, s
);
2254 return float32_round_pack_canonical(&p64
, s
);
2257 float64
float128_to_float64(float128 a
, float_status
*s
)
2262 float128_unpack_canonical(&p128
, a
, s
);
2263 parts_float_to_float_narrow(&p64
, &p128
, s
);
2264 return float64_round_pack_canonical(&p64
, s
);
2267 float128
float32_to_float128(float32 a
, float_status
*s
)
2272 float32_unpack_canonical(&p64
, a
, s
);
2273 parts_float_to_float_widen(&p128
, &p64
, s
);
2274 return float128_round_pack_canonical(&p128
, s
);
2277 float128
float64_to_float128(float64 a
, float_status
*s
)
2282 float64_unpack_canonical(&p64
, a
, s
);
2283 parts_float_to_float_widen(&p128
, &p64
, s
);
2284 return float128_round_pack_canonical(&p128
, s
);
2288 * Rounds the floating-point value `a' to an integer, and returns the
2289 * result as a floating-point value. The operation is performed
2290 * according to the IEC/IEEE Standard for Binary Floating-Point
2294 static FloatParts64
round_to_int(FloatParts64 a
, FloatRoundMode rmode
,
2295 int scale
, float_status
*s
)
2298 case float_class_qnan
:
2299 case float_class_snan
:
2300 parts_return_nan(&a
, s
);
2303 case float_class_zero
:
2304 case float_class_inf
:
2305 /* already "integral" */
2308 case float_class_normal
:
2309 scale
= MIN(MAX(scale
, -0x10000), 0x10000);
2312 if (a
.exp
>= DECOMPOSED_BINARY_POINT
) {
2313 /* already integral */
2318 /* all fractional */
2319 float_raise(float_flag_inexact
, s
);
2321 case float_round_nearest_even
:
2322 one
= a
.exp
== -1 && a
.frac
> DECOMPOSED_IMPLICIT_BIT
;
2324 case float_round_ties_away
:
2325 one
= a
.exp
== -1 && a
.frac
>= DECOMPOSED_IMPLICIT_BIT
;
2327 case float_round_to_zero
:
2330 case float_round_up
:
2333 case float_round_down
:
2336 case float_round_to_odd
:
2340 g_assert_not_reached();
2344 a
.frac
= DECOMPOSED_IMPLICIT_BIT
;
2347 a
.cls
= float_class_zero
;
2350 uint64_t frac_lsb
= DECOMPOSED_IMPLICIT_BIT
>> a
.exp
;
2351 uint64_t frac_lsbm1
= frac_lsb
>> 1;
2352 uint64_t rnd_even_mask
= (frac_lsb
- 1) | frac_lsb
;
2353 uint64_t rnd_mask
= rnd_even_mask
>> 1;
2357 case float_round_nearest_even
:
2358 inc
= ((a
.frac
& rnd_even_mask
) != frac_lsbm1
? frac_lsbm1
: 0);
2360 case float_round_ties_away
:
2363 case float_round_to_zero
:
2366 case float_round_up
:
2367 inc
= a
.sign
? 0 : rnd_mask
;
2369 case float_round_down
:
2370 inc
= a
.sign
? rnd_mask
: 0;
2372 case float_round_to_odd
:
2373 inc
= a
.frac
& frac_lsb
? 0 : rnd_mask
;
2376 g_assert_not_reached();
2379 if (a
.frac
& rnd_mask
) {
2380 float_raise(float_flag_inexact
, s
);
2381 if (uadd64_overflow(a
.frac
, inc
, &a
.frac
)) {
2383 a
.frac
|= DECOMPOSED_IMPLICIT_BIT
;
2386 a
.frac
&= ~rnd_mask
;
2391 g_assert_not_reached();
2396 float16
float16_round_to_int(float16 a
, float_status
*s
)
2398 FloatParts64 pa
, pr
;
2400 float16_unpack_canonical(&pa
, a
, s
);
2401 pr
= round_to_int(pa
, s
->float_rounding_mode
, 0, s
);
2402 return float16_round_pack_canonical(&pr
, s
);
2405 float32
float32_round_to_int(float32 a
, float_status
*s
)
2407 FloatParts64 pa
, pr
;
2409 float32_unpack_canonical(&pa
, a
, s
);
2410 pr
= round_to_int(pa
, s
->float_rounding_mode
, 0, s
);
2411 return float32_round_pack_canonical(&pr
, s
);
2414 float64
float64_round_to_int(float64 a
, float_status
*s
)
2416 FloatParts64 pa
, pr
;
2418 float64_unpack_canonical(&pa
, a
, s
);
2419 pr
= round_to_int(pa
, s
->float_rounding_mode
, 0, s
);
2420 return float64_round_pack_canonical(&pr
, s
);
2424 * Rounds the bfloat16 value `a' to an integer, and returns the
2425 * result as a bfloat16 value.
2428 bfloat16
bfloat16_round_to_int(bfloat16 a
, float_status
*s
)
2430 FloatParts64 pa
, pr
;
2432 bfloat16_unpack_canonical(&pa
, a
, s
);
2433 pr
= round_to_int(pa
, s
->float_rounding_mode
, 0, s
);
2434 return bfloat16_round_pack_canonical(&pr
, s
);
2438 * Returns the result of converting the floating-point value `a' to
2439 * the two's complement integer format. The conversion is performed
2440 * according to the IEC/IEEE Standard for Binary Floating-Point
2441 * Arithmetic---which means in particular that the conversion is
2442 * rounded according to the current rounding mode. If `a' is a NaN,
2443 * the largest positive integer is returned. Otherwise, if the
2444 * conversion overflows, the largest integer with the same sign as `a'
2448 static int64_t round_to_int_and_pack(FloatParts64 in
, FloatRoundMode rmode
,
2449 int scale
, int64_t min
, int64_t max
,
2453 int orig_flags
= get_float_exception_flags(s
);
2454 FloatParts64 p
= round_to_int(in
, rmode
, scale
, s
);
2457 case float_class_snan
:
2458 case float_class_qnan
:
2459 s
->float_exception_flags
= orig_flags
| float_flag_invalid
;
2461 case float_class_inf
:
2462 s
->float_exception_flags
= orig_flags
| float_flag_invalid
;
2463 return p
.sign
? min
: max
;
2464 case float_class_zero
:
2466 case float_class_normal
:
2467 if (p
.exp
<= DECOMPOSED_BINARY_POINT
) {
2468 r
= p
.frac
>> (DECOMPOSED_BINARY_POINT
- p
.exp
);
2473 if (r
<= -(uint64_t) min
) {
2476 s
->float_exception_flags
= orig_flags
| float_flag_invalid
;
2483 s
->float_exception_flags
= orig_flags
| float_flag_invalid
;
2488 g_assert_not_reached();
2492 int8_t float16_to_int8_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2497 float16_unpack_canonical(&p
, a
, s
);
2498 return round_to_int_and_pack(p
, rmode
, scale
, INT8_MIN
, INT8_MAX
, s
);
2501 int16_t float16_to_int16_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2506 float16_unpack_canonical(&p
, a
, s
);
2507 return round_to_int_and_pack(p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2510 int32_t float16_to_int32_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2515 float16_unpack_canonical(&p
, a
, s
);
2516 return round_to_int_and_pack(p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2519 int64_t float16_to_int64_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2524 float16_unpack_canonical(&p
, a
, s
);
2525 return round_to_int_and_pack(p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2528 int16_t float32_to_int16_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2533 float32_unpack_canonical(&p
, a
, s
);
2534 return round_to_int_and_pack(p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2537 int32_t float32_to_int32_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2542 float32_unpack_canonical(&p
, a
, s
);
2543 return round_to_int_and_pack(p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2546 int64_t float32_to_int64_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2551 float32_unpack_canonical(&p
, a
, s
);
2552 return round_to_int_and_pack(p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2555 int16_t float64_to_int16_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2560 float64_unpack_canonical(&p
, a
, s
);
2561 return round_to_int_and_pack(p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2564 int32_t float64_to_int32_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2569 float64_unpack_canonical(&p
, a
, s
);
2570 return round_to_int_and_pack(p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2573 int64_t float64_to_int64_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2578 float64_unpack_canonical(&p
, a
, s
);
2579 return round_to_int_and_pack(p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2582 int8_t float16_to_int8(float16 a
, float_status
*s
)
2584 return float16_to_int8_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2587 int16_t float16_to_int16(float16 a
, float_status
*s
)
2589 return float16_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2592 int32_t float16_to_int32(float16 a
, float_status
*s
)
2594 return float16_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2597 int64_t float16_to_int64(float16 a
, float_status
*s
)
2599 return float16_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2602 int16_t float32_to_int16(float32 a
, float_status
*s
)
2604 return float32_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2607 int32_t float32_to_int32(float32 a
, float_status
*s
)
2609 return float32_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2612 int64_t float32_to_int64(float32 a
, float_status
*s
)
2614 return float32_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2617 int16_t float64_to_int16(float64 a
, float_status
*s
)
2619 return float64_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2622 int32_t float64_to_int32(float64 a
, float_status
*s
)
2624 return float64_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2627 int64_t float64_to_int64(float64 a
, float_status
*s
)
2629 return float64_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2632 int16_t float16_to_int16_round_to_zero(float16 a
, float_status
*s
)
2634 return float16_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2637 int32_t float16_to_int32_round_to_zero(float16 a
, float_status
*s
)
2639 return float16_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2642 int64_t float16_to_int64_round_to_zero(float16 a
, float_status
*s
)
2644 return float16_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2647 int16_t float32_to_int16_round_to_zero(float32 a
, float_status
*s
)
2649 return float32_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2652 int32_t float32_to_int32_round_to_zero(float32 a
, float_status
*s
)
2654 return float32_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2657 int64_t float32_to_int64_round_to_zero(float32 a
, float_status
*s
)
2659 return float32_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2662 int16_t float64_to_int16_round_to_zero(float64 a
, float_status
*s
)
2664 return float64_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2667 int32_t float64_to_int32_round_to_zero(float64 a
, float_status
*s
)
2669 return float64_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2672 int64_t float64_to_int64_round_to_zero(float64 a
, float_status
*s
)
2674 return float64_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2678 * Returns the result of converting the floating-point value `a' to
2679 * the two's complement integer format.
2682 int16_t bfloat16_to_int16_scalbn(bfloat16 a
, FloatRoundMode rmode
, int scale
,
2687 bfloat16_unpack_canonical(&p
, a
, s
);
2688 return round_to_int_and_pack(p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2691 int32_t bfloat16_to_int32_scalbn(bfloat16 a
, FloatRoundMode rmode
, int scale
,
2696 bfloat16_unpack_canonical(&p
, a
, s
);
2697 return round_to_int_and_pack(p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2700 int64_t bfloat16_to_int64_scalbn(bfloat16 a
, FloatRoundMode rmode
, int scale
,
2705 bfloat16_unpack_canonical(&p
, a
, s
);
2706 return round_to_int_and_pack(p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2709 int16_t bfloat16_to_int16(bfloat16 a
, float_status
*s
)
2711 return bfloat16_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2714 int32_t bfloat16_to_int32(bfloat16 a
, float_status
*s
)
2716 return bfloat16_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2719 int64_t bfloat16_to_int64(bfloat16 a
, float_status
*s
)
2721 return bfloat16_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2724 int16_t bfloat16_to_int16_round_to_zero(bfloat16 a
, float_status
*s
)
2726 return bfloat16_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2729 int32_t bfloat16_to_int32_round_to_zero(bfloat16 a
, float_status
*s
)
2731 return bfloat16_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2734 int64_t bfloat16_to_int64_round_to_zero(bfloat16 a
, float_status
*s
)
2736 return bfloat16_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2740 * Returns the result of converting the floating-point value `a' to
2741 * the unsigned integer format. The conversion is performed according
2742 * to the IEC/IEEE Standard for Binary Floating-Point
2743 * Arithmetic---which means in particular that the conversion is
2744 * rounded according to the current rounding mode. If `a' is a NaN,
2745 * the largest unsigned integer is returned. Otherwise, if the
2746 * conversion overflows, the largest unsigned integer is returned. If
2747 * the 'a' is negative, the result is rounded and zero is returned;
2748 * values that do not round to zero will raise the inexact exception
2752 static uint64_t round_to_uint_and_pack(FloatParts64 in
, FloatRoundMode rmode
,
2753 int scale
, uint64_t max
,
2756 int orig_flags
= get_float_exception_flags(s
);
2757 FloatParts64 p
= round_to_int(in
, rmode
, scale
, s
);
2761 case float_class_snan
:
2762 case float_class_qnan
:
2763 s
->float_exception_flags
= orig_flags
| float_flag_invalid
;
2765 case float_class_inf
:
2766 s
->float_exception_flags
= orig_flags
| float_flag_invalid
;
2767 return p
.sign
? 0 : max
;
2768 case float_class_zero
:
2770 case float_class_normal
:
2772 s
->float_exception_flags
= orig_flags
| float_flag_invalid
;
2776 if (p
.exp
<= DECOMPOSED_BINARY_POINT
) {
2777 r
= p
.frac
>> (DECOMPOSED_BINARY_POINT
- p
.exp
);
2779 s
->float_exception_flags
= orig_flags
| float_flag_invalid
;
2783 /* For uint64 this will never trip, but if p.exp is too large
2784 * to shift a decomposed fraction we shall have exited via the
2788 s
->float_exception_flags
= orig_flags
| float_flag_invalid
;
2793 g_assert_not_reached();
2797 uint8_t float16_to_uint8_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2802 float16_unpack_canonical(&p
, a
, s
);
2803 return round_to_uint_and_pack(p
, rmode
, scale
, UINT8_MAX
, s
);
2806 uint16_t float16_to_uint16_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2811 float16_unpack_canonical(&p
, a
, s
);
2812 return round_to_uint_and_pack(p
, rmode
, scale
, UINT16_MAX
, s
);
2815 uint32_t float16_to_uint32_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2820 float16_unpack_canonical(&p
, a
, s
);
2821 return round_to_uint_and_pack(p
, rmode
, scale
, UINT32_MAX
, s
);
2824 uint64_t float16_to_uint64_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2829 float16_unpack_canonical(&p
, a
, s
);
2830 return round_to_uint_and_pack(p
, rmode
, scale
, UINT64_MAX
, s
);
2833 uint16_t float32_to_uint16_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2838 float32_unpack_canonical(&p
, a
, s
);
2839 return round_to_uint_and_pack(p
, rmode
, scale
, UINT16_MAX
, s
);
2842 uint32_t float32_to_uint32_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2847 float32_unpack_canonical(&p
, a
, s
);
2848 return round_to_uint_and_pack(p
, rmode
, scale
, UINT32_MAX
, s
);
2851 uint64_t float32_to_uint64_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2856 float32_unpack_canonical(&p
, a
, s
);
2857 return round_to_uint_and_pack(p
, rmode
, scale
, UINT64_MAX
, s
);
2860 uint16_t float64_to_uint16_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2865 float64_unpack_canonical(&p
, a
, s
);
2866 return round_to_uint_and_pack(p
, rmode
, scale
, UINT16_MAX
, s
);
2869 uint32_t float64_to_uint32_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2874 float64_unpack_canonical(&p
, a
, s
);
2875 return round_to_uint_and_pack(p
, rmode
, scale
, UINT32_MAX
, s
);
2878 uint64_t float64_to_uint64_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2883 float64_unpack_canonical(&p
, a
, s
);
2884 return round_to_uint_and_pack(p
, rmode
, scale
, UINT64_MAX
, s
);
2887 uint8_t float16_to_uint8(float16 a
, float_status
*s
)
2889 return float16_to_uint8_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2892 uint16_t float16_to_uint16(float16 a
, float_status
*s
)
2894 return float16_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2897 uint32_t float16_to_uint32(float16 a
, float_status
*s
)
2899 return float16_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2902 uint64_t float16_to_uint64(float16 a
, float_status
*s
)
2904 return float16_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2907 uint16_t float32_to_uint16(float32 a
, float_status
*s
)
2909 return float32_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2912 uint32_t float32_to_uint32(float32 a
, float_status
*s
)
2914 return float32_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2917 uint64_t float32_to_uint64(float32 a
, float_status
*s
)
2919 return float32_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2922 uint16_t float64_to_uint16(float64 a
, float_status
*s
)
2924 return float64_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2927 uint32_t float64_to_uint32(float64 a
, float_status
*s
)
2929 return float64_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2932 uint64_t float64_to_uint64(float64 a
, float_status
*s
)
2934 return float64_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2937 uint16_t float16_to_uint16_round_to_zero(float16 a
, float_status
*s
)
2939 return float16_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
2942 uint32_t float16_to_uint32_round_to_zero(float16 a
, float_status
*s
)
2944 return float16_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
2947 uint64_t float16_to_uint64_round_to_zero(float16 a
, float_status
*s
)
2949 return float16_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
2952 uint16_t float32_to_uint16_round_to_zero(float32 a
, float_status
*s
)
2954 return float32_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
2957 uint32_t float32_to_uint32_round_to_zero(float32 a
, float_status
*s
)
2959 return float32_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
2962 uint64_t float32_to_uint64_round_to_zero(float32 a
, float_status
*s
)
2964 return float32_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
2967 uint16_t float64_to_uint16_round_to_zero(float64 a
, float_status
*s
)
2969 return float64_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
2972 uint32_t float64_to_uint32_round_to_zero(float64 a
, float_status
*s
)
2974 return float64_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
2977 uint64_t float64_to_uint64_round_to_zero(float64 a
, float_status
*s
)
2979 return float64_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
2983 * Returns the result of converting the bfloat16 value `a' to
2984 * the unsigned integer format.
2987 uint16_t bfloat16_to_uint16_scalbn(bfloat16 a
, FloatRoundMode rmode
,
2988 int scale
, float_status
*s
)
2992 bfloat16_unpack_canonical(&p
, a
, s
);
2993 return round_to_uint_and_pack(p
, rmode
, scale
, UINT16_MAX
, s
);
2996 uint32_t bfloat16_to_uint32_scalbn(bfloat16 a
, FloatRoundMode rmode
,
2997 int scale
, float_status
*s
)
3001 bfloat16_unpack_canonical(&p
, a
, s
);
3002 return round_to_uint_and_pack(p
, rmode
, scale
, UINT32_MAX
, s
);
3005 uint64_t bfloat16_to_uint64_scalbn(bfloat16 a
, FloatRoundMode rmode
,
3006 int scale
, float_status
*s
)
3010 bfloat16_unpack_canonical(&p
, a
, s
);
3011 return round_to_uint_and_pack(p
, rmode
, scale
, UINT64_MAX
, s
);
3014 uint16_t bfloat16_to_uint16(bfloat16 a
, float_status
*s
)
3016 return bfloat16_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3019 uint32_t bfloat16_to_uint32(bfloat16 a
, float_status
*s
)
3021 return bfloat16_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3024 uint64_t bfloat16_to_uint64(bfloat16 a
, float_status
*s
)
3026 return bfloat16_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
3029 uint16_t bfloat16_to_uint16_round_to_zero(bfloat16 a
, float_status
*s
)
3031 return bfloat16_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
3034 uint32_t bfloat16_to_uint32_round_to_zero(bfloat16 a
, float_status
*s
)
3036 return bfloat16_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
3039 uint64_t bfloat16_to_uint64_round_to_zero(bfloat16 a
, float_status
*s
)
3041 return bfloat16_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
3045 * Integer to float conversions
3047 * Returns the result of converting the two's complement integer `a'
3048 * to the floating-point format. The conversion is performed according
3049 * to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3052 static FloatParts64
int_to_float(int64_t a
, int scale
, float_status
*status
)
3054 FloatParts64 r
= { .sign
= false };
3057 r
.cls
= float_class_zero
;
3062 r
.cls
= float_class_normal
;
3068 scale
= MIN(MAX(scale
, -0x10000), 0x10000);
3070 r
.exp
= DECOMPOSED_BINARY_POINT
- shift
+ scale
;
3071 r
.frac
= f
<< shift
;
3077 float16
int64_to_float16_scalbn(int64_t a
, int scale
, float_status
*status
)
3079 FloatParts64 pa
= int_to_float(a
, scale
, status
);
3080 return float16_round_pack_canonical(&pa
, status
);
3083 float16
int32_to_float16_scalbn(int32_t a
, int scale
, float_status
*status
)
3085 return int64_to_float16_scalbn(a
, scale
, status
);
3088 float16
int16_to_float16_scalbn(int16_t a
, int scale
, float_status
*status
)
3090 return int64_to_float16_scalbn(a
, scale
, status
);
3093 float16
int64_to_float16(int64_t a
, float_status
*status
)
3095 return int64_to_float16_scalbn(a
, 0, status
);
3098 float16
int32_to_float16(int32_t a
, float_status
*status
)
3100 return int64_to_float16_scalbn(a
, 0, status
);
3103 float16
int16_to_float16(int16_t a
, float_status
*status
)
3105 return int64_to_float16_scalbn(a
, 0, status
);
3108 float16
int8_to_float16(int8_t a
, float_status
*status
)
3110 return int64_to_float16_scalbn(a
, 0, status
);
3113 float32
int64_to_float32_scalbn(int64_t a
, int scale
, float_status
*status
)
3115 FloatParts64 pa
= int_to_float(a
, scale
, status
);
3116 return float32_round_pack_canonical(&pa
, status
);
3119 float32
int32_to_float32_scalbn(int32_t a
, int scale
, float_status
*status
)
3121 return int64_to_float32_scalbn(a
, scale
, status
);
3124 float32
int16_to_float32_scalbn(int16_t a
, int scale
, float_status
*status
)
3126 return int64_to_float32_scalbn(a
, scale
, status
);
3129 float32
int64_to_float32(int64_t a
, float_status
*status
)
3131 return int64_to_float32_scalbn(a
, 0, status
);
3134 float32
int32_to_float32(int32_t a
, float_status
*status
)
3136 return int64_to_float32_scalbn(a
, 0, status
);
3139 float32
int16_to_float32(int16_t a
, float_status
*status
)
3141 return int64_to_float32_scalbn(a
, 0, status
);
3144 float64
int64_to_float64_scalbn(int64_t a
, int scale
, float_status
*status
)
3146 FloatParts64 pa
= int_to_float(a
, scale
, status
);
3147 return float64_round_pack_canonical(&pa
, status
);
3150 float64
int32_to_float64_scalbn(int32_t a
, int scale
, float_status
*status
)
3152 return int64_to_float64_scalbn(a
, scale
, status
);
3155 float64
int16_to_float64_scalbn(int16_t a
, int scale
, float_status
*status
)
3157 return int64_to_float64_scalbn(a
, scale
, status
);
3160 float64
int64_to_float64(int64_t a
, float_status
*status
)
3162 return int64_to_float64_scalbn(a
, 0, status
);
3165 float64
int32_to_float64(int32_t a
, float_status
*status
)
3167 return int64_to_float64_scalbn(a
, 0, status
);
3170 float64
int16_to_float64(int16_t a
, float_status
*status
)
3172 return int64_to_float64_scalbn(a
, 0, status
);
3176 * Returns the result of converting the two's complement integer `a'
3177 * to the bfloat16 format.
3180 bfloat16
int64_to_bfloat16_scalbn(int64_t a
, int scale
, float_status
*status
)
3182 FloatParts64 pa
= int_to_float(a
, scale
, status
);
3183 return bfloat16_round_pack_canonical(&pa
, status
);
3186 bfloat16
int32_to_bfloat16_scalbn(int32_t a
, int scale
, float_status
*status
)
3188 return int64_to_bfloat16_scalbn(a
, scale
, status
);
3191 bfloat16
int16_to_bfloat16_scalbn(int16_t a
, int scale
, float_status
*status
)
3193 return int64_to_bfloat16_scalbn(a
, scale
, status
);
3196 bfloat16
int64_to_bfloat16(int64_t a
, float_status
*status
)
3198 return int64_to_bfloat16_scalbn(a
, 0, status
);
3201 bfloat16
int32_to_bfloat16(int32_t a
, float_status
*status
)
3203 return int64_to_bfloat16_scalbn(a
, 0, status
);
3206 bfloat16
int16_to_bfloat16(int16_t a
, float_status
*status
)
3208 return int64_to_bfloat16_scalbn(a
, 0, status
);
3212 * Unsigned Integer to float conversions
3214 * Returns the result of converting the unsigned integer `a' to the
3215 * floating-point format. The conversion is performed according to the
3216 * IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3219 static FloatParts64
uint_to_float(uint64_t a
, int scale
, float_status
*status
)
3221 FloatParts64 r
= { .sign
= false };
3225 r
.cls
= float_class_zero
;
3227 scale
= MIN(MAX(scale
, -0x10000), 0x10000);
3229 r
.cls
= float_class_normal
;
3230 r
.exp
= DECOMPOSED_BINARY_POINT
- shift
+ scale
;
3231 r
.frac
= a
<< shift
;
3237 float16
uint64_to_float16_scalbn(uint64_t a
, int scale
, float_status
*status
)
3239 FloatParts64 pa
= uint_to_float(a
, scale
, status
);
3240 return float16_round_pack_canonical(&pa
, status
);
3243 float16
uint32_to_float16_scalbn(uint32_t a
, int scale
, float_status
*status
)
3245 return uint64_to_float16_scalbn(a
, scale
, status
);
3248 float16
uint16_to_float16_scalbn(uint16_t a
, int scale
, float_status
*status
)
3250 return uint64_to_float16_scalbn(a
, scale
, status
);
3253 float16
uint64_to_float16(uint64_t a
, float_status
*status
)
3255 return uint64_to_float16_scalbn(a
, 0, status
);
3258 float16
uint32_to_float16(uint32_t a
, float_status
*status
)
3260 return uint64_to_float16_scalbn(a
, 0, status
);
3263 float16
uint16_to_float16(uint16_t a
, float_status
*status
)
3265 return uint64_to_float16_scalbn(a
, 0, status
);
3268 float16
uint8_to_float16(uint8_t a
, float_status
*status
)
3270 return uint64_to_float16_scalbn(a
, 0, status
);
3273 float32
uint64_to_float32_scalbn(uint64_t a
, int scale
, float_status
*status
)
3275 FloatParts64 pa
= uint_to_float(a
, scale
, status
);
3276 return float32_round_pack_canonical(&pa
, status
);
3279 float32
uint32_to_float32_scalbn(uint32_t a
, int scale
, float_status
*status
)
3281 return uint64_to_float32_scalbn(a
, scale
, status
);
3284 float32
uint16_to_float32_scalbn(uint16_t a
, int scale
, float_status
*status
)
3286 return uint64_to_float32_scalbn(a
, scale
, status
);
3289 float32
uint64_to_float32(uint64_t a
, float_status
*status
)
3291 return uint64_to_float32_scalbn(a
, 0, status
);
3294 float32
uint32_to_float32(uint32_t a
, float_status
*status
)
3296 return uint64_to_float32_scalbn(a
, 0, status
);
3299 float32
uint16_to_float32(uint16_t a
, float_status
*status
)
3301 return uint64_to_float32_scalbn(a
, 0, status
);
3304 float64
uint64_to_float64_scalbn(uint64_t a
, int scale
, float_status
*status
)
3306 FloatParts64 pa
= uint_to_float(a
, scale
, status
);
3307 return float64_round_pack_canonical(&pa
, status
);
3310 float64
uint32_to_float64_scalbn(uint32_t a
, int scale
, float_status
*status
)
3312 return uint64_to_float64_scalbn(a
, scale
, status
);
3315 float64
uint16_to_float64_scalbn(uint16_t a
, int scale
, float_status
*status
)
3317 return uint64_to_float64_scalbn(a
, scale
, status
);
3320 float64
uint64_to_float64(uint64_t a
, float_status
*status
)
3322 return uint64_to_float64_scalbn(a
, 0, status
);
3325 float64
uint32_to_float64(uint32_t a
, float_status
*status
)
3327 return uint64_to_float64_scalbn(a
, 0, status
);
3330 float64
uint16_to_float64(uint16_t a
, float_status
*status
)
3332 return uint64_to_float64_scalbn(a
, 0, status
);
3336 * Returns the result of converting the unsigned integer `a' to the
3340 bfloat16
uint64_to_bfloat16_scalbn(uint64_t a
, int scale
, float_status
*status
)
3342 FloatParts64 pa
= uint_to_float(a
, scale
, status
);
3343 return bfloat16_round_pack_canonical(&pa
, status
);
3346 bfloat16
uint32_to_bfloat16_scalbn(uint32_t a
, int scale
, float_status
*status
)
3348 return uint64_to_bfloat16_scalbn(a
, scale
, status
);
3351 bfloat16
uint16_to_bfloat16_scalbn(uint16_t a
, int scale
, float_status
*status
)
3353 return uint64_to_bfloat16_scalbn(a
, scale
, status
);
3356 bfloat16
uint64_to_bfloat16(uint64_t a
, float_status
*status
)
3358 return uint64_to_bfloat16_scalbn(a
, 0, status
);
3361 bfloat16
uint32_to_bfloat16(uint32_t a
, float_status
*status
)
3363 return uint64_to_bfloat16_scalbn(a
, 0, status
);
3366 bfloat16
uint16_to_bfloat16(uint16_t a
, float_status
*status
)
3368 return uint64_to_bfloat16_scalbn(a
, 0, status
);
3372 /* min() and max() functions. These can't be implemented as
3373 * 'compare and pick one input' because that would mishandle
3374 * NaNs and +0 vs -0.
3376 * minnum() and maxnum() functions. These are similar to the min()
3377 * and max() functions but if one of the arguments is a QNaN and
3378 * the other is numerical then the numerical argument is returned.
3379 * SNaNs will get quietened before being returned.
3380 * minnum() and maxnum correspond to the IEEE 754-2008 minNum()
3381 * and maxNum() operations. min() and max() are the typical min/max
3382 * semantics provided by many CPUs which predate that specification.
3384 * minnummag() and maxnummag() functions correspond to minNumMag()
3385 * and minNumMag() from the IEEE-754 2008.
3387 static FloatParts64
minmax_floats(FloatParts64 a
, FloatParts64 b
, bool ismin
,
3388 bool ieee
, bool ismag
, float_status
*s
)
3390 if (unlikely(is_nan(a
.cls
) || is_nan(b
.cls
))) {
3392 /* Takes two floating-point values `a' and `b', one of
3393 * which is a NaN, and returns the appropriate NaN
3394 * result. If either `a' or `b' is a signaling NaN,
3395 * the invalid exception is raised.
3397 if (is_snan(a
.cls
) || is_snan(b
.cls
)) {
3398 return *parts_pick_nan(&a
, &b
, s
);
3399 } else if (is_nan(a
.cls
) && !is_nan(b
.cls
)) {
3401 } else if (is_nan(b
.cls
) && !is_nan(a
.cls
)) {
3405 return *parts_pick_nan(&a
, &b
, s
);
3410 case float_class_normal
:
3413 case float_class_inf
:
3416 case float_class_zero
:
3420 g_assert_not_reached();
3424 case float_class_normal
:
3427 case float_class_inf
:
3430 case float_class_zero
:
3434 g_assert_not_reached();
3438 if (ismag
&& (a_exp
!= b_exp
|| a
.frac
!= b
.frac
)) {
3439 bool a_less
= a_exp
< b_exp
;
3440 if (a_exp
== b_exp
) {
3441 a_less
= a
.frac
< b
.frac
;
3443 return a_less
^ ismin
? b
: a
;
3446 if (a
.sign
== b
.sign
) {
3447 bool a_less
= a_exp
< b_exp
;
3448 if (a_exp
== b_exp
) {
3449 a_less
= a
.frac
< b
.frac
;
3451 return a
.sign
^ a_less
^ ismin
? b
: a
;
3453 return a
.sign
^ ismin
? b
: a
;
3458 #define MINMAX(sz, name, ismin, isiee, ismag) \
3459 float ## sz float ## sz ## _ ## name(float ## sz a, float ## sz b, \
3462 FloatParts64 pa, pb, pr; \
3463 float ## sz ## _unpack_canonical(&pa, a, s); \
3464 float ## sz ## _unpack_canonical(&pb, b, s); \
3465 pr = minmax_floats(pa, pb, ismin, isiee, ismag, s); \
3466 return float ## sz ## _round_pack_canonical(&pr, s); \
3469 MINMAX(16, min
, true, false, false)
3470 MINMAX(16, minnum
, true, true, false)
3471 MINMAX(16, minnummag
, true, true, true)
3472 MINMAX(16, max
, false, false, false)
3473 MINMAX(16, maxnum
, false, true, false)
3474 MINMAX(16, maxnummag
, false, true, true)
3476 MINMAX(32, min
, true, false, false)
3477 MINMAX(32, minnum
, true, true, false)
3478 MINMAX(32, minnummag
, true, true, true)
3479 MINMAX(32, max
, false, false, false)
3480 MINMAX(32, maxnum
, false, true, false)
3481 MINMAX(32, maxnummag
, false, true, true)
3483 MINMAX(64, min
, true, false, false)
3484 MINMAX(64, minnum
, true, true, false)
3485 MINMAX(64, minnummag
, true, true, true)
3486 MINMAX(64, max
, false, false, false)
3487 MINMAX(64, maxnum
, false, true, false)
3488 MINMAX(64, maxnummag
, false, true, true)
3492 #define BF16_MINMAX(name, ismin, isiee, ismag) \
3493 bfloat16 bfloat16_ ## name(bfloat16 a, bfloat16 b, float_status *s) \
3495 FloatParts64 pa, pb, pr; \
3496 bfloat16_unpack_canonical(&pa, a, s); \
3497 bfloat16_unpack_canonical(&pb, b, s); \
3498 pr = minmax_floats(pa, pb, ismin, isiee, ismag, s); \
3499 return bfloat16_round_pack_canonical(&pr, s); \
3502 BF16_MINMAX(min
, true, false, false)
3503 BF16_MINMAX(minnum
, true, true, false)
3504 BF16_MINMAX(minnummag
, true, true, true)
3505 BF16_MINMAX(max
, false, false, false)
3506 BF16_MINMAX(maxnum
, false, true, false)
3507 BF16_MINMAX(maxnummag
, false, true, true)
3511 /* Floating point compare */
3512 static FloatRelation
compare_floats(FloatParts64 a
, FloatParts64 b
, bool is_quiet
,
3515 if (is_nan(a
.cls
) || is_nan(b
.cls
)) {
3517 a
.cls
== float_class_snan
||
3518 b
.cls
== float_class_snan
) {
3519 float_raise(float_flag_invalid
, s
);
3521 return float_relation_unordered
;
3524 if (a
.cls
== float_class_zero
) {
3525 if (b
.cls
== float_class_zero
) {
3526 return float_relation_equal
;
3528 return b
.sign
? float_relation_greater
: float_relation_less
;
3529 } else if (b
.cls
== float_class_zero
) {
3530 return a
.sign
? float_relation_less
: float_relation_greater
;
3533 /* The only really important thing about infinity is its sign. If
3534 * both are infinities the sign marks the smallest of the two.
3536 if (a
.cls
== float_class_inf
) {
3537 if ((b
.cls
== float_class_inf
) && (a
.sign
== b
.sign
)) {
3538 return float_relation_equal
;
3540 return a
.sign
? float_relation_less
: float_relation_greater
;
3541 } else if (b
.cls
== float_class_inf
) {
3542 return b
.sign
? float_relation_greater
: float_relation_less
;
3545 if (a
.sign
!= b
.sign
) {
3546 return a
.sign
? float_relation_less
: float_relation_greater
;
3549 if (a
.exp
== b
.exp
) {
3550 if (a
.frac
== b
.frac
) {
3551 return float_relation_equal
;
3554 return a
.frac
> b
.frac
?
3555 float_relation_less
: float_relation_greater
;
3557 return a
.frac
> b
.frac
?
3558 float_relation_greater
: float_relation_less
;
3562 return a
.exp
> b
.exp
? float_relation_less
: float_relation_greater
;
3564 return a
.exp
> b
.exp
? float_relation_greater
: float_relation_less
;
3569 #define COMPARE(name, attr, sz) \
3571 name(float ## sz a, float ## sz b, bool is_quiet, float_status *s) \
3573 FloatParts64 pa, pb; \
3574 float ## sz ## _unpack_canonical(&pa, a, s); \
3575 float ## sz ## _unpack_canonical(&pb, b, s); \
3576 return compare_floats(pa, pb, is_quiet, s); \
3579 COMPARE(soft_f16_compare
, QEMU_FLATTEN
, 16)
3580 COMPARE(soft_f32_compare
, QEMU_SOFTFLOAT_ATTR
, 32)
3581 COMPARE(soft_f64_compare
, QEMU_SOFTFLOAT_ATTR
, 64)
3585 FloatRelation
float16_compare(float16 a
, float16 b
, float_status
*s
)
3587 return soft_f16_compare(a
, b
, false, s
);
3590 FloatRelation
float16_compare_quiet(float16 a
, float16 b
, float_status
*s
)
3592 return soft_f16_compare(a
, b
, true, s
);
3595 static FloatRelation QEMU_FLATTEN
3596 f32_compare(float32 xa
, float32 xb
, bool is_quiet
, float_status
*s
)
3598 union_float32 ua
, ub
;
3603 if (QEMU_NO_HARDFLOAT
) {
3607 float32_input_flush2(&ua
.s
, &ub
.s
, s
);
3608 if (isgreaterequal(ua
.h
, ub
.h
)) {
3609 if (isgreater(ua
.h
, ub
.h
)) {
3610 return float_relation_greater
;
3612 return float_relation_equal
;
3614 if (likely(isless(ua
.h
, ub
.h
))) {
3615 return float_relation_less
;
3617 /* The only condition remaining is unordered.
3618 * Fall through to set flags.
3621 return soft_f32_compare(ua
.s
, ub
.s
, is_quiet
, s
);
3624 FloatRelation
float32_compare(float32 a
, float32 b
, float_status
*s
)
3626 return f32_compare(a
, b
, false, s
);
3629 FloatRelation
float32_compare_quiet(float32 a
, float32 b
, float_status
*s
)
3631 return f32_compare(a
, b
, true, s
);
3634 static FloatRelation QEMU_FLATTEN
3635 f64_compare(float64 xa
, float64 xb
, bool is_quiet
, float_status
*s
)
3637 union_float64 ua
, ub
;
3642 if (QEMU_NO_HARDFLOAT
) {
3646 float64_input_flush2(&ua
.s
, &ub
.s
, s
);
3647 if (isgreaterequal(ua
.h
, ub
.h
)) {
3648 if (isgreater(ua
.h
, ub
.h
)) {
3649 return float_relation_greater
;
3651 return float_relation_equal
;
3653 if (likely(isless(ua
.h
, ub
.h
))) {
3654 return float_relation_less
;
3656 /* The only condition remaining is unordered.
3657 * Fall through to set flags.
3660 return soft_f64_compare(ua
.s
, ub
.s
, is_quiet
, s
);
3663 FloatRelation
float64_compare(float64 a
, float64 b
, float_status
*s
)
3665 return f64_compare(a
, b
, false, s
);
3668 FloatRelation
float64_compare_quiet(float64 a
, float64 b
, float_status
*s
)
3670 return f64_compare(a
, b
, true, s
);
3673 static FloatRelation QEMU_FLATTEN
3674 soft_bf16_compare(bfloat16 a
, bfloat16 b
, bool is_quiet
, float_status
*s
)
3676 FloatParts64 pa
, pb
;
3678 bfloat16_unpack_canonical(&pa
, a
, s
);
3679 bfloat16_unpack_canonical(&pb
, b
, s
);
3680 return compare_floats(pa
, pb
, is_quiet
, s
);
3683 FloatRelation
bfloat16_compare(bfloat16 a
, bfloat16 b
, float_status
*s
)
3685 return soft_bf16_compare(a
, b
, false, s
);
3688 FloatRelation
bfloat16_compare_quiet(bfloat16 a
, bfloat16 b
, float_status
*s
)
3690 return soft_bf16_compare(a
, b
, true, s
);
3693 /* Multiply A by 2 raised to the power N. */
3694 static FloatParts64
scalbn_decomposed(FloatParts64 a
, int n
, float_status
*s
)
3696 if (unlikely(is_nan(a
.cls
))) {
3697 parts_return_nan(&a
, s
);
3699 if (a
.cls
== float_class_normal
) {
3700 /* The largest float type (even though not supported by FloatParts64)
3701 * is float128, which has a 15 bit exponent. Bounding N to 16 bits
3702 * still allows rounding to infinity, without allowing overflow
3703 * within the int32_t that backs FloatParts64.exp.
3705 n
= MIN(MAX(n
, -0x10000), 0x10000);
3711 float16
float16_scalbn(float16 a
, int n
, float_status
*status
)
3713 FloatParts64 pa
, pr
;
3715 float16_unpack_canonical(&pa
, a
, status
);
3716 pr
= scalbn_decomposed(pa
, n
, status
);
3717 return float16_round_pack_canonical(&pr
, status
);
3720 float32
float32_scalbn(float32 a
, int n
, float_status
*status
)
3722 FloatParts64 pa
, pr
;
3724 float32_unpack_canonical(&pa
, a
, status
);
3725 pr
= scalbn_decomposed(pa
, n
, status
);
3726 return float32_round_pack_canonical(&pr
, status
);
3729 float64
float64_scalbn(float64 a
, int n
, float_status
*status
)
3731 FloatParts64 pa
, pr
;
3733 float64_unpack_canonical(&pa
, a
, status
);
3734 pr
= scalbn_decomposed(pa
, n
, status
);
3735 return float64_round_pack_canonical(&pr
, status
);
3738 bfloat16
bfloat16_scalbn(bfloat16 a
, int n
, float_status
*status
)
3740 FloatParts64 pa
, pr
;
3742 bfloat16_unpack_canonical(&pa
, a
, status
);
3743 pr
= scalbn_decomposed(pa
, n
, status
);
3744 return bfloat16_round_pack_canonical(&pr
, status
);
3750 * The old softfloat code did an approximation step before zeroing in
3751 * on the final result. However for simpleness we just compute the
3752 * square root by iterating down from the implicit bit to enough extra
3753 * bits to ensure we get a correctly rounded result.
3755 * This does mean however the calculation is slower than before,
3756 * especially for 64 bit floats.
3759 static FloatParts64
sqrt_float(FloatParts64 a
, float_status
*s
, const FloatFmt
*p
)
3761 uint64_t a_frac
, r_frac
, s_frac
;
3764 if (is_nan(a
.cls
)) {
3765 parts_return_nan(&a
, s
);
3768 if (a
.cls
== float_class_zero
) {
3769 return a
; /* sqrt(+-0) = +-0 */
3772 float_raise(float_flag_invalid
, s
);
3773 parts_default_nan(&a
, s
);
3776 if (a
.cls
== float_class_inf
) {
3777 return a
; /* sqrt(+inf) = +inf */
3780 assert(a
.cls
== float_class_normal
);
3782 /* We need two overflow bits at the top. Adding room for that is a
3783 * right shift. If the exponent is odd, we can discard the low bit
3784 * by multiplying the fraction by 2; that's a left shift. Combine
3785 * those and we shift right by 1 if the exponent is odd, otherwise 2.
3787 a_frac
= a
.frac
>> (2 - (a
.exp
& 1));
3790 /* Bit-by-bit computation of sqrt. */
3794 /* Iterate from implicit bit down to the 3 extra bits to compute a
3795 * properly rounded result. Remember we've inserted two more bits
3796 * at the top, so these positions are two less.
3798 bit
= DECOMPOSED_BINARY_POINT
- 2;
3799 last_bit
= MAX(p
->frac_shift
- 4, 0);
3801 uint64_t q
= 1ULL << bit
;
3802 uint64_t t_frac
= s_frac
+ q
;
3803 if (t_frac
<= a_frac
) {
3804 s_frac
= t_frac
+ q
;
3809 } while (--bit
>= last_bit
);
3811 /* Undo the right shift done above. If there is any remaining
3812 * fraction, the result is inexact. Set the sticky bit.
3814 a
.frac
= (r_frac
<< 2) + (a_frac
!= 0);
3819 float16 QEMU_FLATTEN
float16_sqrt(float16 a
, float_status
*status
)
3821 FloatParts64 pa
, pr
;
3823 float16_unpack_canonical(&pa
, a
, status
);
3824 pr
= sqrt_float(pa
, status
, &float16_params
);
3825 return float16_round_pack_canonical(&pr
, status
);
3828 static float32 QEMU_SOFTFLOAT_ATTR
3829 soft_f32_sqrt(float32 a
, float_status
*status
)
3831 FloatParts64 pa
, pr
;
3833 float32_unpack_canonical(&pa
, a
, status
);
3834 pr
= sqrt_float(pa
, status
, &float32_params
);
3835 return float32_round_pack_canonical(&pr
, status
);
3838 static float64 QEMU_SOFTFLOAT_ATTR
3839 soft_f64_sqrt(float64 a
, float_status
*status
)
3841 FloatParts64 pa
, pr
;
3843 float64_unpack_canonical(&pa
, a
, status
);
3844 pr
= sqrt_float(pa
, status
, &float64_params
);
3845 return float64_round_pack_canonical(&pr
, status
);
3848 float32 QEMU_FLATTEN
float32_sqrt(float32 xa
, float_status
*s
)
3850 union_float32 ua
, ur
;
3853 if (unlikely(!can_use_fpu(s
))) {
3857 float32_input_flush1(&ua
.s
, s
);
3858 if (QEMU_HARDFLOAT_1F32_USE_FP
) {
3859 if (unlikely(!(fpclassify(ua
.h
) == FP_NORMAL
||
3860 fpclassify(ua
.h
) == FP_ZERO
) ||
3864 } else if (unlikely(!float32_is_zero_or_normal(ua
.s
) ||
3865 float32_is_neg(ua
.s
))) {
3872 return soft_f32_sqrt(ua
.s
, s
);
3875 float64 QEMU_FLATTEN
float64_sqrt(float64 xa
, float_status
*s
)
3877 union_float64 ua
, ur
;
3880 if (unlikely(!can_use_fpu(s
))) {
3884 float64_input_flush1(&ua
.s
, s
);
3885 if (QEMU_HARDFLOAT_1F64_USE_FP
) {
3886 if (unlikely(!(fpclassify(ua
.h
) == FP_NORMAL
||
3887 fpclassify(ua
.h
) == FP_ZERO
) ||
3891 } else if (unlikely(!float64_is_zero_or_normal(ua
.s
) ||
3892 float64_is_neg(ua
.s
))) {
3899 return soft_f64_sqrt(ua
.s
, s
);
3902 bfloat16 QEMU_FLATTEN
bfloat16_sqrt(bfloat16 a
, float_status
*status
)
3904 FloatParts64 pa
, pr
;
3906 bfloat16_unpack_canonical(&pa
, a
, status
);
3907 pr
= sqrt_float(pa
, status
, &bfloat16_params
);
3908 return bfloat16_round_pack_canonical(&pr
, status
);
3911 /*----------------------------------------------------------------------------
3912 | The pattern for a default generated NaN.
3913 *----------------------------------------------------------------------------*/
3915 float16
float16_default_nan(float_status
*status
)
3919 parts_default_nan(&p
, status
);
3920 p
.frac
>>= float16_params
.frac_shift
;
3921 return float16_pack_raw(&p
);
3924 float32
float32_default_nan(float_status
*status
)
3928 parts_default_nan(&p
, status
);
3929 p
.frac
>>= float32_params
.frac_shift
;
3930 return float32_pack_raw(&p
);
3933 float64
float64_default_nan(float_status
*status
)
3937 parts_default_nan(&p
, status
);
3938 p
.frac
>>= float64_params
.frac_shift
;
3939 return float64_pack_raw(&p
);
3942 float128
float128_default_nan(float_status
*status
)
3946 parts_default_nan(&p
, status
);
3947 frac_shr(&p
, float128_params
.frac_shift
);
3948 return float128_pack_raw(&p
);
3951 bfloat16
bfloat16_default_nan(float_status
*status
)
3955 parts_default_nan(&p
, status
);
3956 p
.frac
>>= bfloat16_params
.frac_shift
;
3957 return bfloat16_pack_raw(&p
);
3960 /*----------------------------------------------------------------------------
3961 | Returns a quiet NaN from a signalling NaN for the floating point value `a'.
3962 *----------------------------------------------------------------------------*/
3964 float16
float16_silence_nan(float16 a
, float_status
*status
)
3968 float16_unpack_raw(&p
, a
);
3969 p
.frac
<<= float16_params
.frac_shift
;
3970 parts_silence_nan(&p
, status
);
3971 p
.frac
>>= float16_params
.frac_shift
;
3972 return float16_pack_raw(&p
);
3975 float32
float32_silence_nan(float32 a
, float_status
*status
)
3979 float32_unpack_raw(&p
, a
);
3980 p
.frac
<<= float32_params
.frac_shift
;
3981 parts_silence_nan(&p
, status
);
3982 p
.frac
>>= float32_params
.frac_shift
;
3983 return float32_pack_raw(&p
);
3986 float64
float64_silence_nan(float64 a
, float_status
*status
)
3990 float64_unpack_raw(&p
, a
);
3991 p
.frac
<<= float64_params
.frac_shift
;
3992 parts_silence_nan(&p
, status
);
3993 p
.frac
>>= float64_params
.frac_shift
;
3994 return float64_pack_raw(&p
);
3997 bfloat16
bfloat16_silence_nan(bfloat16 a
, float_status
*status
)
4001 bfloat16_unpack_raw(&p
, a
);
4002 p
.frac
<<= bfloat16_params
.frac_shift
;
4003 parts_silence_nan(&p
, status
);
4004 p
.frac
>>= bfloat16_params
.frac_shift
;
4005 return bfloat16_pack_raw(&p
);
4008 float128
float128_silence_nan(float128 a
, float_status
*status
)
4012 float128_unpack_raw(&p
, a
);
4013 frac_shl(&p
, float128_params
.frac_shift
);
4014 parts_silence_nan(&p
, status
);
4015 frac_shr(&p
, float128_params
.frac_shift
);
4016 return float128_pack_raw(&p
);
4019 /*----------------------------------------------------------------------------
4020 | If `a' is denormal and we are in flush-to-zero mode then set the
4021 | input-denormal exception and return zero. Otherwise just return the value.
4022 *----------------------------------------------------------------------------*/
4024 static bool parts_squash_denormal(FloatParts64 p
, float_status
*status
)
4026 if (p
.exp
== 0 && p
.frac
!= 0) {
4027 float_raise(float_flag_input_denormal
, status
);
4034 float16
float16_squash_input_denormal(float16 a
, float_status
*status
)
4036 if (status
->flush_inputs_to_zero
) {
4039 float16_unpack_raw(&p
, a
);
4040 if (parts_squash_denormal(p
, status
)) {
4041 return float16_set_sign(float16_zero
, p
.sign
);
4047 float32
float32_squash_input_denormal(float32 a
, float_status
*status
)
4049 if (status
->flush_inputs_to_zero
) {
4052 float32_unpack_raw(&p
, a
);
4053 if (parts_squash_denormal(p
, status
)) {
4054 return float32_set_sign(float32_zero
, p
.sign
);
4060 float64
float64_squash_input_denormal(float64 a
, float_status
*status
)
4062 if (status
->flush_inputs_to_zero
) {
4065 float64_unpack_raw(&p
, a
);
4066 if (parts_squash_denormal(p
, status
)) {
4067 return float64_set_sign(float64_zero
, p
.sign
);
4073 bfloat16
bfloat16_squash_input_denormal(bfloat16 a
, float_status
*status
)
4075 if (status
->flush_inputs_to_zero
) {
4078 bfloat16_unpack_raw(&p
, a
);
4079 if (parts_squash_denormal(p
, status
)) {
4080 return bfloat16_set_sign(bfloat16_zero
, p
.sign
);
4086 /*----------------------------------------------------------------------------
4087 | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
4088 | and 7, and returns the properly rounded 32-bit integer corresponding to the
4089 | input. If `zSign' is 1, the input is negated before being converted to an
4090 | integer. Bit 63 of `absZ' must be zero. Ordinarily, the fixed-point input
4091 | is simply rounded to an integer, with the inexact exception raised if the
4092 | input cannot be represented exactly as an integer. However, if the fixed-
4093 | point input is too large, the invalid exception is raised and the largest
4094 | positive or negative integer is returned.
4095 *----------------------------------------------------------------------------*/
4097 static int32_t roundAndPackInt32(bool zSign
, uint64_t absZ
,
4098 float_status
*status
)
4100 int8_t roundingMode
;
4101 bool roundNearestEven
;
4102 int8_t roundIncrement
, roundBits
;
4105 roundingMode
= status
->float_rounding_mode
;
4106 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4107 switch (roundingMode
) {
4108 case float_round_nearest_even
:
4109 case float_round_ties_away
:
4110 roundIncrement
= 0x40;
4112 case float_round_to_zero
:
4115 case float_round_up
:
4116 roundIncrement
= zSign
? 0 : 0x7f;
4118 case float_round_down
:
4119 roundIncrement
= zSign
? 0x7f : 0;
4121 case float_round_to_odd
:
4122 roundIncrement
= absZ
& 0x80 ? 0 : 0x7f;
4127 roundBits
= absZ
& 0x7F;
4128 absZ
= ( absZ
+ roundIncrement
)>>7;
4129 if (!(roundBits
^ 0x40) && roundNearestEven
) {
4133 if ( zSign
) z
= - z
;
4134 if ( ( absZ
>>32 ) || ( z
&& ( ( z
< 0 ) ^ zSign
) ) ) {
4135 float_raise(float_flag_invalid
, status
);
4136 return zSign
? INT32_MIN
: INT32_MAX
;
4139 float_raise(float_flag_inexact
, status
);
4145 /*----------------------------------------------------------------------------
4146 | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
4147 | `absZ1', with binary point between bits 63 and 64 (between the input words),
4148 | and returns the properly rounded 64-bit integer corresponding to the input.
4149 | If `zSign' is 1, the input is negated before being converted to an integer.
4150 | Ordinarily, the fixed-point input is simply rounded to an integer, with
4151 | the inexact exception raised if the input cannot be represented exactly as
4152 | an integer. However, if the fixed-point input is too large, the invalid
4153 | exception is raised and the largest positive or negative integer is
4155 *----------------------------------------------------------------------------*/
4157 static int64_t roundAndPackInt64(bool zSign
, uint64_t absZ0
, uint64_t absZ1
,
4158 float_status
*status
)
4160 int8_t roundingMode
;
4161 bool roundNearestEven
, increment
;
4164 roundingMode
= status
->float_rounding_mode
;
4165 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4166 switch (roundingMode
) {
4167 case float_round_nearest_even
:
4168 case float_round_ties_away
:
4169 increment
= ((int64_t) absZ1
< 0);
4171 case float_round_to_zero
:
4174 case float_round_up
:
4175 increment
= !zSign
&& absZ1
;
4177 case float_round_down
:
4178 increment
= zSign
&& absZ1
;
4180 case float_round_to_odd
:
4181 increment
= !(absZ0
& 1) && absZ1
;
4188 if ( absZ0
== 0 ) goto overflow
;
4189 if (!(absZ1
<< 1) && roundNearestEven
) {
4194 if ( zSign
) z
= - z
;
4195 if ( z
&& ( ( z
< 0 ) ^ zSign
) ) {
4197 float_raise(float_flag_invalid
, status
);
4198 return zSign
? INT64_MIN
: INT64_MAX
;
4201 float_raise(float_flag_inexact
, status
);
4207 /*----------------------------------------------------------------------------
4208 | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
4209 | `absZ1', with binary point between bits 63 and 64 (between the input words),
4210 | and returns the properly rounded 64-bit unsigned integer corresponding to the
4211 | input. Ordinarily, the fixed-point input is simply rounded to an integer,
4212 | with the inexact exception raised if the input cannot be represented exactly
4213 | as an integer. However, if the fixed-point input is too large, the invalid
4214 | exception is raised and the largest unsigned integer is returned.
4215 *----------------------------------------------------------------------------*/
4217 static int64_t roundAndPackUint64(bool zSign
, uint64_t absZ0
,
4218 uint64_t absZ1
, float_status
*status
)
4220 int8_t roundingMode
;
4221 bool roundNearestEven
, increment
;
4223 roundingMode
= status
->float_rounding_mode
;
4224 roundNearestEven
= (roundingMode
== float_round_nearest_even
);
4225 switch (roundingMode
) {
4226 case float_round_nearest_even
:
4227 case float_round_ties_away
:
4228 increment
= ((int64_t)absZ1
< 0);
4230 case float_round_to_zero
:
4233 case float_round_up
:
4234 increment
= !zSign
&& absZ1
;
4236 case float_round_down
:
4237 increment
= zSign
&& absZ1
;
4239 case float_round_to_odd
:
4240 increment
= !(absZ0
& 1) && absZ1
;
4248 float_raise(float_flag_invalid
, status
);
4251 if (!(absZ1
<< 1) && roundNearestEven
) {
4256 if (zSign
&& absZ0
) {
4257 float_raise(float_flag_invalid
, status
);
4262 float_raise(float_flag_inexact
, status
);
4267 /*----------------------------------------------------------------------------
4268 | Normalizes the subnormal single-precision floating-point value represented
4269 | by the denormalized significand `aSig'. The normalized exponent and
4270 | significand are stored at the locations pointed to by `zExpPtr' and
4271 | `zSigPtr', respectively.
4272 *----------------------------------------------------------------------------*/
4275 normalizeFloat32Subnormal(uint32_t aSig
, int *zExpPtr
, uint32_t *zSigPtr
)
4279 shiftCount
= clz32(aSig
) - 8;
4280 *zSigPtr
= aSig
<<shiftCount
;
4281 *zExpPtr
= 1 - shiftCount
;
4285 /*----------------------------------------------------------------------------
4286 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4287 | and significand `zSig', and returns the proper single-precision floating-
4288 | point value corresponding to the abstract input. Ordinarily, the abstract
4289 | value is simply rounded and packed into the single-precision format, with
4290 | the inexact exception raised if the abstract input cannot be represented
4291 | exactly. However, if the abstract value is too large, the overflow and
4292 | inexact exceptions are raised and an infinity or maximal finite value is
4293 | returned. If the abstract value is too small, the input value is rounded to
4294 | a subnormal number, and the underflow and inexact exceptions are raised if
4295 | the abstract input cannot be represented exactly as a subnormal single-
4296 | precision floating-point number.
4297 | The input significand `zSig' has its binary point between bits 30
4298 | and 29, which is 7 bits to the left of the usual location. This shifted
4299 | significand must be normalized or smaller. If `zSig' is not normalized,
4300 | `zExp' must be 0; in that case, the result returned is a subnormal number,
4301 | and it must not require rounding. In the usual case that `zSig' is
4302 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
4303 | The handling of underflow and overflow follows the IEC/IEEE Standard for
4304 | Binary Floating-Point Arithmetic.
4305 *----------------------------------------------------------------------------*/
4307 static float32
roundAndPackFloat32(bool zSign
, int zExp
, uint32_t zSig
,
4308 float_status
*status
)
4310 int8_t roundingMode
;
4311 bool roundNearestEven
;
4312 int8_t roundIncrement
, roundBits
;
4315 roundingMode
= status
->float_rounding_mode
;
4316 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4317 switch (roundingMode
) {
4318 case float_round_nearest_even
:
4319 case float_round_ties_away
:
4320 roundIncrement
= 0x40;
4322 case float_round_to_zero
:
4325 case float_round_up
:
4326 roundIncrement
= zSign
? 0 : 0x7f;
4328 case float_round_down
:
4329 roundIncrement
= zSign
? 0x7f : 0;
4331 case float_round_to_odd
:
4332 roundIncrement
= zSig
& 0x80 ? 0 : 0x7f;
4338 roundBits
= zSig
& 0x7F;
4339 if ( 0xFD <= (uint16_t) zExp
) {
4340 if ( ( 0xFD < zExp
)
4341 || ( ( zExp
== 0xFD )
4342 && ( (int32_t) ( zSig
+ roundIncrement
) < 0 ) )
4344 bool overflow_to_inf
= roundingMode
!= float_round_to_odd
&&
4345 roundIncrement
!= 0;
4346 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4347 return packFloat32(zSign
, 0xFF, -!overflow_to_inf
);
4350 if (status
->flush_to_zero
) {
4351 float_raise(float_flag_output_denormal
, status
);
4352 return packFloat32(zSign
, 0, 0);
4354 isTiny
= status
->tininess_before_rounding
4356 || (zSig
+ roundIncrement
< 0x80000000);
4357 shift32RightJamming( zSig
, - zExp
, &zSig
);
4359 roundBits
= zSig
& 0x7F;
4360 if (isTiny
&& roundBits
) {
4361 float_raise(float_flag_underflow
, status
);
4363 if (roundingMode
== float_round_to_odd
) {
4365 * For round-to-odd case, the roundIncrement depends on
4366 * zSig which just changed.
4368 roundIncrement
= zSig
& 0x80 ? 0 : 0x7f;
4373 float_raise(float_flag_inexact
, status
);
4375 zSig
= ( zSig
+ roundIncrement
)>>7;
4376 if (!(roundBits
^ 0x40) && roundNearestEven
) {
4379 if ( zSig
== 0 ) zExp
= 0;
4380 return packFloat32( zSign
, zExp
, zSig
);
4384 /*----------------------------------------------------------------------------
4385 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4386 | and significand `zSig', and returns the proper single-precision floating-
4387 | point value corresponding to the abstract input. This routine is just like
4388 | `roundAndPackFloat32' except that `zSig' does not have to be normalized.
4389 | Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
4390 | floating-point exponent.
4391 *----------------------------------------------------------------------------*/
4394 normalizeRoundAndPackFloat32(bool zSign
, int zExp
, uint32_t zSig
,
4395 float_status
*status
)
4399 shiftCount
= clz32(zSig
) - 1;
4400 return roundAndPackFloat32(zSign
, zExp
- shiftCount
, zSig
<<shiftCount
,
4405 /*----------------------------------------------------------------------------
4406 | Normalizes the subnormal double-precision floating-point value represented
4407 | by the denormalized significand `aSig'. The normalized exponent and
4408 | significand are stored at the locations pointed to by `zExpPtr' and
4409 | `zSigPtr', respectively.
4410 *----------------------------------------------------------------------------*/
4413 normalizeFloat64Subnormal(uint64_t aSig
, int *zExpPtr
, uint64_t *zSigPtr
)
4417 shiftCount
= clz64(aSig
) - 11;
4418 *zSigPtr
= aSig
<<shiftCount
;
4419 *zExpPtr
= 1 - shiftCount
;
4423 /*----------------------------------------------------------------------------
4424 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
4425 | double-precision floating-point value, returning the result. After being
4426 | shifted into the proper positions, the three fields are simply added
4427 | together to form the result. This means that any integer portion of `zSig'
4428 | will be added into the exponent. Since a properly normalized significand
4429 | will have an integer portion equal to 1, the `zExp' input should be 1 less
4430 | than the desired result exponent whenever `zSig' is a complete, normalized
4432 *----------------------------------------------------------------------------*/
4434 static inline float64
packFloat64(bool zSign
, int zExp
, uint64_t zSig
)
4437 return make_float64(
4438 ( ( (uint64_t) zSign
)<<63 ) + ( ( (uint64_t) zExp
)<<52 ) + zSig
);
4442 /*----------------------------------------------------------------------------
4443 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4444 | and significand `zSig', and returns the proper double-precision floating-
4445 | point value corresponding to the abstract input. Ordinarily, the abstract
4446 | value is simply rounded and packed into the double-precision format, with
4447 | the inexact exception raised if the abstract input cannot be represented
4448 | exactly. However, if the abstract value is too large, the overflow and
4449 | inexact exceptions are raised and an infinity or maximal finite value is
4450 | returned. If the abstract value is too small, the input value is rounded to
4451 | a subnormal number, and the underflow and inexact exceptions are raised if
4452 | the abstract input cannot be represented exactly as a subnormal double-
4453 | precision floating-point number.
4454 | The input significand `zSig' has its binary point between bits 62
4455 | and 61, which is 10 bits to the left of the usual location. This shifted
4456 | significand must be normalized or smaller. If `zSig' is not normalized,
4457 | `zExp' must be 0; in that case, the result returned is a subnormal number,
4458 | and it must not require rounding. In the usual case that `zSig' is
4459 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
4460 | The handling of underflow and overflow follows the IEC/IEEE Standard for
4461 | Binary Floating-Point Arithmetic.
4462 *----------------------------------------------------------------------------*/
4464 static float64
roundAndPackFloat64(bool zSign
, int zExp
, uint64_t zSig
,
4465 float_status
*status
)
4467 int8_t roundingMode
;
4468 bool roundNearestEven
;
4469 int roundIncrement
, roundBits
;
4472 roundingMode
= status
->float_rounding_mode
;
4473 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4474 switch (roundingMode
) {
4475 case float_round_nearest_even
:
4476 case float_round_ties_away
:
4477 roundIncrement
= 0x200;
4479 case float_round_to_zero
:
4482 case float_round_up
:
4483 roundIncrement
= zSign
? 0 : 0x3ff;
4485 case float_round_down
:
4486 roundIncrement
= zSign
? 0x3ff : 0;
4488 case float_round_to_odd
:
4489 roundIncrement
= (zSig
& 0x400) ? 0 : 0x3ff;
4494 roundBits
= zSig
& 0x3FF;
4495 if ( 0x7FD <= (uint16_t) zExp
) {
4496 if ( ( 0x7FD < zExp
)
4497 || ( ( zExp
== 0x7FD )
4498 && ( (int64_t) ( zSig
+ roundIncrement
) < 0 ) )
4500 bool overflow_to_inf
= roundingMode
!= float_round_to_odd
&&
4501 roundIncrement
!= 0;
4502 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4503 return packFloat64(zSign
, 0x7FF, -(!overflow_to_inf
));
4506 if (status
->flush_to_zero
) {
4507 float_raise(float_flag_output_denormal
, status
);
4508 return packFloat64(zSign
, 0, 0);
4510 isTiny
= status
->tininess_before_rounding
4512 || (zSig
+ roundIncrement
< UINT64_C(0x8000000000000000));
4513 shift64RightJamming( zSig
, - zExp
, &zSig
);
4515 roundBits
= zSig
& 0x3FF;
4516 if (isTiny
&& roundBits
) {
4517 float_raise(float_flag_underflow
, status
);
4519 if (roundingMode
== float_round_to_odd
) {
4521 * For round-to-odd case, the roundIncrement depends on
4522 * zSig which just changed.
4524 roundIncrement
= (zSig
& 0x400) ? 0 : 0x3ff;
4529 float_raise(float_flag_inexact
, status
);
4531 zSig
= ( zSig
+ roundIncrement
)>>10;
4532 if (!(roundBits
^ 0x200) && roundNearestEven
) {
4535 if ( zSig
== 0 ) zExp
= 0;
4536 return packFloat64( zSign
, zExp
, zSig
);
4540 /*----------------------------------------------------------------------------
4541 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4542 | and significand `zSig', and returns the proper double-precision floating-
4543 | point value corresponding to the abstract input. This routine is just like
4544 | `roundAndPackFloat64' except that `zSig' does not have to be normalized.
4545 | Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
4546 | floating-point exponent.
4547 *----------------------------------------------------------------------------*/
4550 normalizeRoundAndPackFloat64(bool zSign
, int zExp
, uint64_t zSig
,
4551 float_status
*status
)
4555 shiftCount
= clz64(zSig
) - 1;
4556 return roundAndPackFloat64(zSign
, zExp
- shiftCount
, zSig
<<shiftCount
,
4561 /*----------------------------------------------------------------------------
4562 | Normalizes the subnormal extended double-precision floating-point value
4563 | represented by the denormalized significand `aSig'. The normalized exponent
4564 | and significand are stored at the locations pointed to by `zExpPtr' and
4565 | `zSigPtr', respectively.
4566 *----------------------------------------------------------------------------*/
4568 void normalizeFloatx80Subnormal(uint64_t aSig
, int32_t *zExpPtr
,
4573 shiftCount
= clz64(aSig
);
4574 *zSigPtr
= aSig
<<shiftCount
;
4575 *zExpPtr
= 1 - shiftCount
;
4578 /*----------------------------------------------------------------------------
4579 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4580 | and extended significand formed by the concatenation of `zSig0' and `zSig1',
4581 | and returns the proper extended double-precision floating-point value
4582 | corresponding to the abstract input. Ordinarily, the abstract value is
4583 | rounded and packed into the extended double-precision format, with the
4584 | inexact exception raised if the abstract input cannot be represented
4585 | exactly. However, if the abstract value is too large, the overflow and
4586 | inexact exceptions are raised and an infinity or maximal finite value is
4587 | returned. If the abstract value is too small, the input value is rounded to
4588 | a subnormal number, and the underflow and inexact exceptions are raised if
4589 | the abstract input cannot be represented exactly as a subnormal extended
4590 | double-precision floating-point number.
4591 | If `roundingPrecision' is 32 or 64, the result is rounded to the same
4592 | number of bits as single or double precision, respectively. Otherwise, the
4593 | result is rounded to the full precision of the extended double-precision
4595 | The input significand must be normalized or smaller. If the input
4596 | significand is not normalized, `zExp' must be 0; in that case, the result
4597 | returned is a subnormal number, and it must not require rounding. The
4598 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
4599 | Floating-Point Arithmetic.
4600 *----------------------------------------------------------------------------*/
4602 floatx80
roundAndPackFloatx80(int8_t roundingPrecision
, bool zSign
,
4603 int32_t zExp
, uint64_t zSig0
, uint64_t zSig1
,
4604 float_status
*status
)
4606 int8_t roundingMode
;
4607 bool roundNearestEven
, increment
, isTiny
;
4608 int64_t roundIncrement
, roundMask
, roundBits
;
4610 roundingMode
= status
->float_rounding_mode
;
4611 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4612 if ( roundingPrecision
== 80 ) goto precision80
;
4613 if ( roundingPrecision
== 64 ) {
4614 roundIncrement
= UINT64_C(0x0000000000000400);
4615 roundMask
= UINT64_C(0x00000000000007FF);
4617 else if ( roundingPrecision
== 32 ) {
4618 roundIncrement
= UINT64_C(0x0000008000000000);
4619 roundMask
= UINT64_C(0x000000FFFFFFFFFF);
4624 zSig0
|= ( zSig1
!= 0 );
4625 switch (roundingMode
) {
4626 case float_round_nearest_even
:
4627 case float_round_ties_away
:
4629 case float_round_to_zero
:
4632 case float_round_up
:
4633 roundIncrement
= zSign
? 0 : roundMask
;
4635 case float_round_down
:
4636 roundIncrement
= zSign
? roundMask
: 0;
4641 roundBits
= zSig0
& roundMask
;
4642 if ( 0x7FFD <= (uint32_t) ( zExp
- 1 ) ) {
4643 if ( ( 0x7FFE < zExp
)
4644 || ( ( zExp
== 0x7FFE ) && ( zSig0
+ roundIncrement
< zSig0
) )
4649 if (status
->flush_to_zero
) {
4650 float_raise(float_flag_output_denormal
, status
);
4651 return packFloatx80(zSign
, 0, 0);
4653 isTiny
= status
->tininess_before_rounding
4655 || (zSig0
<= zSig0
+ roundIncrement
);
4656 shift64RightJamming( zSig0
, 1 - zExp
, &zSig0
);
4658 roundBits
= zSig0
& roundMask
;
4659 if (isTiny
&& roundBits
) {
4660 float_raise(float_flag_underflow
, status
);
4663 float_raise(float_flag_inexact
, status
);
4665 zSig0
+= roundIncrement
;
4666 if ( (int64_t) zSig0
< 0 ) zExp
= 1;
4667 roundIncrement
= roundMask
+ 1;
4668 if ( roundNearestEven
&& ( roundBits
<<1 == roundIncrement
) ) {
4669 roundMask
|= roundIncrement
;
4671 zSig0
&= ~ roundMask
;
4672 return packFloatx80( zSign
, zExp
, zSig0
);
4676 float_raise(float_flag_inexact
, status
);
4678 zSig0
+= roundIncrement
;
4679 if ( zSig0
< roundIncrement
) {
4681 zSig0
= UINT64_C(0x8000000000000000);
4683 roundIncrement
= roundMask
+ 1;
4684 if ( roundNearestEven
&& ( roundBits
<<1 == roundIncrement
) ) {
4685 roundMask
|= roundIncrement
;
4687 zSig0
&= ~ roundMask
;
4688 if ( zSig0
== 0 ) zExp
= 0;
4689 return packFloatx80( zSign
, zExp
, zSig0
);
4691 switch (roundingMode
) {
4692 case float_round_nearest_even
:
4693 case float_round_ties_away
:
4694 increment
= ((int64_t)zSig1
< 0);
4696 case float_round_to_zero
:
4699 case float_round_up
:
4700 increment
= !zSign
&& zSig1
;
4702 case float_round_down
:
4703 increment
= zSign
&& zSig1
;
4708 if ( 0x7FFD <= (uint32_t) ( zExp
- 1 ) ) {
4709 if ( ( 0x7FFE < zExp
)
4710 || ( ( zExp
== 0x7FFE )
4711 && ( zSig0
== UINT64_C(0xFFFFFFFFFFFFFFFF) )
4717 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4718 if ( ( roundingMode
== float_round_to_zero
)
4719 || ( zSign
&& ( roundingMode
== float_round_up
) )
4720 || ( ! zSign
&& ( roundingMode
== float_round_down
) )
4722 return packFloatx80( zSign
, 0x7FFE, ~ roundMask
);
4724 return packFloatx80(zSign
,
4725 floatx80_infinity_high
,
4726 floatx80_infinity_low
);
4729 isTiny
= status
->tininess_before_rounding
4732 || (zSig0
< UINT64_C(0xFFFFFFFFFFFFFFFF));
4733 shift64ExtraRightJamming( zSig0
, zSig1
, 1 - zExp
, &zSig0
, &zSig1
);
4735 if (isTiny
&& zSig1
) {
4736 float_raise(float_flag_underflow
, status
);
4739 float_raise(float_flag_inexact
, status
);
4741 switch (roundingMode
) {
4742 case float_round_nearest_even
:
4743 case float_round_ties_away
:
4744 increment
= ((int64_t)zSig1
< 0);
4746 case float_round_to_zero
:
4749 case float_round_up
:
4750 increment
= !zSign
&& zSig1
;
4752 case float_round_down
:
4753 increment
= zSign
&& zSig1
;
4760 if (!(zSig1
<< 1) && roundNearestEven
) {
4763 if ( (int64_t) zSig0
< 0 ) zExp
= 1;
4765 return packFloatx80( zSign
, zExp
, zSig0
);
4769 float_raise(float_flag_inexact
, status
);
4775 zSig0
= UINT64_C(0x8000000000000000);
4778 if (!(zSig1
<< 1) && roundNearestEven
) {
4784 if ( zSig0
== 0 ) zExp
= 0;
4786 return packFloatx80( zSign
, zExp
, zSig0
);
4790 /*----------------------------------------------------------------------------
4791 | Takes an abstract floating-point value having sign `zSign', exponent
4792 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
4793 | and returns the proper extended double-precision floating-point value
4794 | corresponding to the abstract input. This routine is just like
4795 | `roundAndPackFloatx80' except that the input significand does not have to be
4797 *----------------------------------------------------------------------------*/
4799 floatx80
normalizeRoundAndPackFloatx80(int8_t roundingPrecision
,
4800 bool zSign
, int32_t zExp
,
4801 uint64_t zSig0
, uint64_t zSig1
,
4802 float_status
*status
)
4811 shiftCount
= clz64(zSig0
);
4812 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
4814 return roundAndPackFloatx80(roundingPrecision
, zSign
, zExp
,
4815 zSig0
, zSig1
, status
);
4819 /*----------------------------------------------------------------------------
4820 | Returns the least-significant 64 fraction bits of the quadruple-precision
4821 | floating-point value `a'.
4822 *----------------------------------------------------------------------------*/
4824 static inline uint64_t extractFloat128Frac1( float128 a
)
4831 /*----------------------------------------------------------------------------
4832 | Returns the most-significant 48 fraction bits of the quadruple-precision
4833 | floating-point value `a'.
4834 *----------------------------------------------------------------------------*/
4836 static inline uint64_t extractFloat128Frac0( float128 a
)
4839 return a
.high
& UINT64_C(0x0000FFFFFFFFFFFF);
4843 /*----------------------------------------------------------------------------
4844 | Returns the exponent bits of the quadruple-precision floating-point value
4846 *----------------------------------------------------------------------------*/
4848 static inline int32_t extractFloat128Exp( float128 a
)
4851 return ( a
.high
>>48 ) & 0x7FFF;
4855 /*----------------------------------------------------------------------------
4856 | Returns the sign bit of the quadruple-precision floating-point value `a'.
4857 *----------------------------------------------------------------------------*/
4859 static inline bool extractFloat128Sign(float128 a
)
4861 return a
.high
>> 63;
4864 /*----------------------------------------------------------------------------
4865 | Normalizes the subnormal quadruple-precision floating-point value
4866 | represented by the denormalized significand formed by the concatenation of
4867 | `aSig0' and `aSig1'. The normalized exponent is stored at the location
4868 | pointed to by `zExpPtr'. The most significant 49 bits of the normalized
4869 | significand are stored at the location pointed to by `zSig0Ptr', and the
4870 | least significant 64 bits of the normalized significand are stored at the
4871 | location pointed to by `zSig1Ptr'.
4872 *----------------------------------------------------------------------------*/
4875 normalizeFloat128Subnormal(
4886 shiftCount
= clz64(aSig1
) - 15;
4887 if ( shiftCount
< 0 ) {
4888 *zSig0Ptr
= aSig1
>>( - shiftCount
);
4889 *zSig1Ptr
= aSig1
<<( shiftCount
& 63 );
4892 *zSig0Ptr
= aSig1
<<shiftCount
;
4895 *zExpPtr
= - shiftCount
- 63;
4898 shiftCount
= clz64(aSig0
) - 15;
4899 shortShift128Left( aSig0
, aSig1
, shiftCount
, zSig0Ptr
, zSig1Ptr
);
4900 *zExpPtr
= 1 - shiftCount
;
4905 /*----------------------------------------------------------------------------
4906 | Packs the sign `zSign', the exponent `zExp', and the significand formed
4907 | by the concatenation of `zSig0' and `zSig1' into a quadruple-precision
4908 | floating-point value, returning the result. After being shifted into the
4909 | proper positions, the three fields `zSign', `zExp', and `zSig0' are simply
4910 | added together to form the most significant 32 bits of the result. This
4911 | means that any integer portion of `zSig0' will be added into the exponent.
4912 | Since a properly normalized significand will have an integer portion equal
4913 | to 1, the `zExp' input should be 1 less than the desired result exponent
4914 | whenever `zSig0' and `zSig1' concatenated form a complete, normalized
4916 *----------------------------------------------------------------------------*/
4918 static inline float128
4919 packFloat128(bool zSign
, int32_t zExp
, uint64_t zSig0
, uint64_t zSig1
)
4924 z
.high
= ((uint64_t)zSign
<< 63) + ((uint64_t)zExp
<< 48) + zSig0
;
4928 /*----------------------------------------------------------------------------
4929 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4930 | and extended significand formed by the concatenation of `zSig0', `zSig1',
4931 | and `zSig2', and returns the proper quadruple-precision floating-point value
4932 | corresponding to the abstract input. Ordinarily, the abstract value is
4933 | simply rounded and packed into the quadruple-precision format, with the
4934 | inexact exception raised if the abstract input cannot be represented
4935 | exactly. However, if the abstract value is too large, the overflow and
4936 | inexact exceptions are raised and an infinity or maximal finite value is
4937 | returned. If the abstract value is too small, the input value is rounded to
4938 | a subnormal number, and the underflow and inexact exceptions are raised if
4939 | the abstract input cannot be represented exactly as a subnormal quadruple-
4940 | precision floating-point number.
4941 | The input significand must be normalized or smaller. If the input
4942 | significand is not normalized, `zExp' must be 0; in that case, the result
4943 | returned is a subnormal number, and it must not require rounding. In the
4944 | usual case that the input significand is normalized, `zExp' must be 1 less
4945 | than the ``true'' floating-point exponent. The handling of underflow and
4946 | overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4947 *----------------------------------------------------------------------------*/
4949 static float128
roundAndPackFloat128(bool zSign
, int32_t zExp
,
4950 uint64_t zSig0
, uint64_t zSig1
,
4951 uint64_t zSig2
, float_status
*status
)
4953 int8_t roundingMode
;
4954 bool roundNearestEven
, increment
, isTiny
;
4956 roundingMode
= status
->float_rounding_mode
;
4957 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4958 switch (roundingMode
) {
4959 case float_round_nearest_even
:
4960 case float_round_ties_away
:
4961 increment
= ((int64_t)zSig2
< 0);
4963 case float_round_to_zero
:
4966 case float_round_up
:
4967 increment
= !zSign
&& zSig2
;
4969 case float_round_down
:
4970 increment
= zSign
&& zSig2
;
4972 case float_round_to_odd
:
4973 increment
= !(zSig1
& 0x1) && zSig2
;
4978 if ( 0x7FFD <= (uint32_t) zExp
) {
4979 if ( ( 0x7FFD < zExp
)
4980 || ( ( zExp
== 0x7FFD )
4982 UINT64_C(0x0001FFFFFFFFFFFF),
4983 UINT64_C(0xFFFFFFFFFFFFFFFF),
4990 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4991 if ( ( roundingMode
== float_round_to_zero
)
4992 || ( zSign
&& ( roundingMode
== float_round_up
) )
4993 || ( ! zSign
&& ( roundingMode
== float_round_down
) )
4994 || (roundingMode
== float_round_to_odd
)
5000 UINT64_C(0x0000FFFFFFFFFFFF),
5001 UINT64_C(0xFFFFFFFFFFFFFFFF)
5004 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5007 if (status
->flush_to_zero
) {
5008 float_raise(float_flag_output_denormal
, status
);
5009 return packFloat128(zSign
, 0, 0, 0);
5011 isTiny
= status
->tininess_before_rounding
5014 || lt128(zSig0
, zSig1
,
5015 UINT64_C(0x0001FFFFFFFFFFFF),
5016 UINT64_C(0xFFFFFFFFFFFFFFFF));
5017 shift128ExtraRightJamming(
5018 zSig0
, zSig1
, zSig2
, - zExp
, &zSig0
, &zSig1
, &zSig2
);
5020 if (isTiny
&& zSig2
) {
5021 float_raise(float_flag_underflow
, status
);
5023 switch (roundingMode
) {
5024 case float_round_nearest_even
:
5025 case float_round_ties_away
:
5026 increment
= ((int64_t)zSig2
< 0);
5028 case float_round_to_zero
:
5031 case float_round_up
:
5032 increment
= !zSign
&& zSig2
;
5034 case float_round_down
:
5035 increment
= zSign
&& zSig2
;
5037 case float_round_to_odd
:
5038 increment
= !(zSig1
& 0x1) && zSig2
;
5046 float_raise(float_flag_inexact
, status
);
5049 add128( zSig0
, zSig1
, 0, 1, &zSig0
, &zSig1
);
5050 if ((zSig2
+ zSig2
== 0) && roundNearestEven
) {
5055 if ( ( zSig0
| zSig1
) == 0 ) zExp
= 0;
5057 return packFloat128( zSign
, zExp
, zSig0
, zSig1
);
5061 /*----------------------------------------------------------------------------
5062 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
5063 | and significand formed by the concatenation of `zSig0' and `zSig1', and
5064 | returns the proper quadruple-precision floating-point value corresponding
5065 | to the abstract input. This routine is just like `roundAndPackFloat128'
5066 | except that the input significand has fewer bits and does not have to be
5067 | normalized. In all cases, `zExp' must be 1 less than the ``true'' floating-
5069 *----------------------------------------------------------------------------*/
5071 static float128
normalizeRoundAndPackFloat128(bool zSign
, int32_t zExp
,
5072 uint64_t zSig0
, uint64_t zSig1
,
5073 float_status
*status
)
5083 shiftCount
= clz64(zSig0
) - 15;
5084 if ( 0 <= shiftCount
) {
5086 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
5089 shift128ExtraRightJamming(
5090 zSig0
, zSig1
, 0, - shiftCount
, &zSig0
, &zSig1
, &zSig2
);
5093 return roundAndPackFloat128(zSign
, zExp
, zSig0
, zSig1
, zSig2
, status
);
5098 /*----------------------------------------------------------------------------
5099 | Returns the result of converting the 32-bit two's complement integer `a'
5100 | to the extended double-precision floating-point format. The conversion
5101 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5103 *----------------------------------------------------------------------------*/
5105 floatx80
int32_to_floatx80(int32_t a
, float_status
*status
)
5112 if ( a
== 0 ) return packFloatx80( 0, 0, 0 );
5114 absA
= zSign
? - a
: a
;
5115 shiftCount
= clz32(absA
) + 32;
5117 return packFloatx80( zSign
, 0x403E - shiftCount
, zSig
<<shiftCount
);
5121 /*----------------------------------------------------------------------------
5122 | Returns the result of converting the 32-bit two's complement integer `a' to
5123 | the quadruple-precision floating-point format. The conversion is performed
5124 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5125 *----------------------------------------------------------------------------*/
5127 float128
int32_to_float128(int32_t a
, float_status
*status
)
5134 if ( a
== 0 ) return packFloat128( 0, 0, 0, 0 );
5136 absA
= zSign
? - a
: a
;
5137 shiftCount
= clz32(absA
) + 17;
5139 return packFloat128( zSign
, 0x402E - shiftCount
, zSig0
<<shiftCount
, 0 );
5143 /*----------------------------------------------------------------------------
5144 | Returns the result of converting the 64-bit two's complement integer `a'
5145 | to the extended double-precision floating-point format. The conversion
5146 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5148 *----------------------------------------------------------------------------*/
5150 floatx80
int64_to_floatx80(int64_t a
, float_status
*status
)
5156 if ( a
== 0 ) return packFloatx80( 0, 0, 0 );
5158 absA
= zSign
? - a
: a
;
5159 shiftCount
= clz64(absA
);
5160 return packFloatx80( zSign
, 0x403E - shiftCount
, absA
<<shiftCount
);
5164 /*----------------------------------------------------------------------------
5165 | Returns the result of converting the 64-bit two's complement integer `a' to
5166 | the quadruple-precision floating-point format. The conversion is performed
5167 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5168 *----------------------------------------------------------------------------*/
5170 float128
int64_to_float128(int64_t a
, float_status
*status
)
5176 uint64_t zSig0
, zSig1
;
5178 if ( a
== 0 ) return packFloat128( 0, 0, 0, 0 );
5180 absA
= zSign
? - a
: a
;
5181 shiftCount
= clz64(absA
) + 49;
5182 zExp
= 0x406E - shiftCount
;
5183 if ( 64 <= shiftCount
) {
5192 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
5193 return packFloat128( zSign
, zExp
, zSig0
, zSig1
);
5197 /*----------------------------------------------------------------------------
5198 | Returns the result of converting the 64-bit unsigned integer `a'
5199 | to the quadruple-precision floating-point format. The conversion is performed
5200 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5201 *----------------------------------------------------------------------------*/
5203 float128
uint64_to_float128(uint64_t a
, float_status
*status
)
5206 return float128_zero
;
5208 return normalizeRoundAndPackFloat128(0, 0x406E, 0, a
, status
);
5211 /*----------------------------------------------------------------------------
5212 | Returns the result of converting the single-precision floating-point value
5213 | `a' to the extended double-precision floating-point format. The conversion
5214 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5216 *----------------------------------------------------------------------------*/
5218 floatx80
float32_to_floatx80(float32 a
, float_status
*status
)
5224 a
= float32_squash_input_denormal(a
, status
);
5225 aSig
= extractFloat32Frac( a
);
5226 aExp
= extractFloat32Exp( a
);
5227 aSign
= extractFloat32Sign( a
);
5228 if ( aExp
== 0xFF ) {
5230 floatx80 res
= commonNaNToFloatx80(float32ToCommonNaN(a
, status
),
5232 return floatx80_silence_nan(res
, status
);
5234 return packFloatx80(aSign
,
5235 floatx80_infinity_high
,
5236 floatx80_infinity_low
);
5239 if ( aSig
== 0 ) return packFloatx80( aSign
, 0, 0 );
5240 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
5243 return packFloatx80( aSign
, aExp
+ 0x3F80, ( (uint64_t) aSig
)<<40 );
5247 /*----------------------------------------------------------------------------
5248 | Returns the remainder of the single-precision floating-point value `a'
5249 | with respect to the corresponding value `b'. The operation is performed
5250 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5251 *----------------------------------------------------------------------------*/
5253 float32
float32_rem(float32 a
, float32 b
, float_status
*status
)
5256 int aExp
, bExp
, expDiff
;
5257 uint32_t aSig
, bSig
;
5259 uint64_t aSig64
, bSig64
, q64
;
5260 uint32_t alternateASig
;
5262 a
= float32_squash_input_denormal(a
, status
);
5263 b
= float32_squash_input_denormal(b
, status
);
5265 aSig
= extractFloat32Frac( a
);
5266 aExp
= extractFloat32Exp( a
);
5267 aSign
= extractFloat32Sign( a
);
5268 bSig
= extractFloat32Frac( b
);
5269 bExp
= extractFloat32Exp( b
);
5270 if ( aExp
== 0xFF ) {
5271 if ( aSig
|| ( ( bExp
== 0xFF ) && bSig
) ) {
5272 return propagateFloat32NaN(a
, b
, status
);
5274 float_raise(float_flag_invalid
, status
);
5275 return float32_default_nan(status
);
5277 if ( bExp
== 0xFF ) {
5279 return propagateFloat32NaN(a
, b
, status
);
5285 float_raise(float_flag_invalid
, status
);
5286 return float32_default_nan(status
);
5288 normalizeFloat32Subnormal( bSig
, &bExp
, &bSig
);
5291 if ( aSig
== 0 ) return a
;
5292 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
5294 expDiff
= aExp
- bExp
;
5297 if ( expDiff
< 32 ) {
5300 if ( expDiff
< 0 ) {
5301 if ( expDiff
< -1 ) return a
;
5304 q
= ( bSig
<= aSig
);
5305 if ( q
) aSig
-= bSig
;
5306 if ( 0 < expDiff
) {
5307 q
= ( ( (uint64_t) aSig
)<<32 ) / bSig
;
5310 aSig
= ( ( aSig
>>1 )<<( expDiff
- 1 ) ) - bSig
* q
;
5318 if ( bSig
<= aSig
) aSig
-= bSig
;
5319 aSig64
= ( (uint64_t) aSig
)<<40;
5320 bSig64
= ( (uint64_t) bSig
)<<40;
5322 while ( 0 < expDiff
) {
5323 q64
= estimateDiv128To64( aSig64
, 0, bSig64
);
5324 q64
= ( 2 < q64
) ? q64
- 2 : 0;
5325 aSig64
= - ( ( bSig
* q64
)<<38 );
5329 q64
= estimateDiv128To64( aSig64
, 0, bSig64
);
5330 q64
= ( 2 < q64
) ? q64
- 2 : 0;
5331 q
= q64
>>( 64 - expDiff
);
5333 aSig
= ( ( aSig64
>>33 )<<( expDiff
- 1 ) ) - bSig
* q
;
5336 alternateASig
= aSig
;
5339 } while ( 0 <= (int32_t) aSig
);
5340 sigMean
= aSig
+ alternateASig
;
5341 if ( ( sigMean
< 0 ) || ( ( sigMean
== 0 ) && ( q
& 1 ) ) ) {
5342 aSig
= alternateASig
;
5344 zSign
= ( (int32_t) aSig
< 0 );
5345 if ( zSign
) aSig
= - aSig
;
5346 return normalizeRoundAndPackFloat32(aSign
^ zSign
, bExp
, aSig
, status
);
5351 /*----------------------------------------------------------------------------
5352 | Returns the binary exponential of the single-precision floating-point value
5353 | `a'. The operation is performed according to the IEC/IEEE Standard for
5354 | Binary Floating-Point Arithmetic.
5356 | Uses the following identities:
5358 | 1. -------------------------------------------------------------------------
5362 | 2. -------------------------------------------------------------------------
5365 | e = 1 + --- + --- + --- + --- + --- + ... + --- + ...
5367 *----------------------------------------------------------------------------*/
5369 static const float64 float32_exp2_coefficients
[15] =
5371 const_float64( 0x3ff0000000000000ll
), /* 1 */
5372 const_float64( 0x3fe0000000000000ll
), /* 2 */
5373 const_float64( 0x3fc5555555555555ll
), /* 3 */
5374 const_float64( 0x3fa5555555555555ll
), /* 4 */
5375 const_float64( 0x3f81111111111111ll
), /* 5 */
5376 const_float64( 0x3f56c16c16c16c17ll
), /* 6 */
5377 const_float64( 0x3f2a01a01a01a01all
), /* 7 */
5378 const_float64( 0x3efa01a01a01a01all
), /* 8 */
5379 const_float64( 0x3ec71de3a556c734ll
), /* 9 */
5380 const_float64( 0x3e927e4fb7789f5cll
), /* 10 */
5381 const_float64( 0x3e5ae64567f544e4ll
), /* 11 */
5382 const_float64( 0x3e21eed8eff8d898ll
), /* 12 */
5383 const_float64( 0x3de6124613a86d09ll
), /* 13 */
5384 const_float64( 0x3da93974a8c07c9dll
), /* 14 */
5385 const_float64( 0x3d6ae7f3e733b81fll
), /* 15 */
5388 float32
float32_exp2(float32 a
, float_status
*status
)
5395 a
= float32_squash_input_denormal(a
, status
);
5397 aSig
= extractFloat32Frac( a
);
5398 aExp
= extractFloat32Exp( a
);
5399 aSign
= extractFloat32Sign( a
);
5401 if ( aExp
== 0xFF) {
5403 return propagateFloat32NaN(a
, float32_zero
, status
);
5405 return (aSign
) ? float32_zero
: a
;
5408 if (aSig
== 0) return float32_one
;
5411 float_raise(float_flag_inexact
, status
);
5413 /* ******************************* */
5414 /* using float64 for approximation */
5415 /* ******************************* */
5416 x
= float32_to_float64(a
, status
);
5417 x
= float64_mul(x
, float64_ln2
, status
);
5421 for (i
= 0 ; i
< 15 ; i
++) {
5424 f
= float64_mul(xn
, float32_exp2_coefficients
[i
], status
);
5425 r
= float64_add(r
, f
, status
);
5427 xn
= float64_mul(xn
, x
, status
);
5430 return float64_to_float32(r
, status
);
5433 /*----------------------------------------------------------------------------
5434 | Returns the binary log of the single-precision floating-point value `a'.
5435 | The operation is performed according to the IEC/IEEE Standard for Binary
5436 | Floating-Point Arithmetic.
5437 *----------------------------------------------------------------------------*/
5438 float32
float32_log2(float32 a
, float_status
*status
)
5442 uint32_t aSig
, zSig
, i
;
5444 a
= float32_squash_input_denormal(a
, status
);
5445 aSig
= extractFloat32Frac( a
);
5446 aExp
= extractFloat32Exp( a
);
5447 aSign
= extractFloat32Sign( a
);
5450 if ( aSig
== 0 ) return packFloat32( 1, 0xFF, 0 );
5451 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
5454 float_raise(float_flag_invalid
, status
);
5455 return float32_default_nan(status
);
5457 if ( aExp
== 0xFF ) {
5459 return propagateFloat32NaN(a
, float32_zero
, status
);
5469 for (i
= 1 << 22; i
> 0; i
>>= 1) {
5470 aSig
= ( (uint64_t)aSig
* aSig
) >> 23;
5471 if ( aSig
& 0x01000000 ) {
5480 return normalizeRoundAndPackFloat32(zSign
, 0x85, zSig
, status
);
5483 /*----------------------------------------------------------------------------
5484 | Returns the result of converting the double-precision floating-point value
5485 | `a' to the extended double-precision floating-point format. The conversion
5486 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5488 *----------------------------------------------------------------------------*/
5490 floatx80
float64_to_floatx80(float64 a
, float_status
*status
)
5496 a
= float64_squash_input_denormal(a
, status
);
5497 aSig
= extractFloat64Frac( a
);
5498 aExp
= extractFloat64Exp( a
);
5499 aSign
= extractFloat64Sign( a
);
5500 if ( aExp
== 0x7FF ) {
5502 floatx80 res
= commonNaNToFloatx80(float64ToCommonNaN(a
, status
),
5504 return floatx80_silence_nan(res
, status
);
5506 return packFloatx80(aSign
,
5507 floatx80_infinity_high
,
5508 floatx80_infinity_low
);
5511 if ( aSig
== 0 ) return packFloatx80( aSign
, 0, 0 );
5512 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
5516 aSign
, aExp
+ 0x3C00, (aSig
| UINT64_C(0x0010000000000000)) << 11);
5520 /*----------------------------------------------------------------------------
5521 | Returns the remainder of the double-precision floating-point value `a'
5522 | with respect to the corresponding value `b'. The operation is performed
5523 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5524 *----------------------------------------------------------------------------*/
5526 float64
float64_rem(float64 a
, float64 b
, float_status
*status
)
5529 int aExp
, bExp
, expDiff
;
5530 uint64_t aSig
, bSig
;
5531 uint64_t q
, alternateASig
;
5534 a
= float64_squash_input_denormal(a
, status
);
5535 b
= float64_squash_input_denormal(b
, status
);
5536 aSig
= extractFloat64Frac( a
);
5537 aExp
= extractFloat64Exp( a
);
5538 aSign
= extractFloat64Sign( a
);
5539 bSig
= extractFloat64Frac( b
);
5540 bExp
= extractFloat64Exp( b
);
5541 if ( aExp
== 0x7FF ) {
5542 if ( aSig
|| ( ( bExp
== 0x7FF ) && bSig
) ) {
5543 return propagateFloat64NaN(a
, b
, status
);
5545 float_raise(float_flag_invalid
, status
);
5546 return float64_default_nan(status
);
5548 if ( bExp
== 0x7FF ) {
5550 return propagateFloat64NaN(a
, b
, status
);
5556 float_raise(float_flag_invalid
, status
);
5557 return float64_default_nan(status
);
5559 normalizeFloat64Subnormal( bSig
, &bExp
, &bSig
);
5562 if ( aSig
== 0 ) return a
;
5563 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
5565 expDiff
= aExp
- bExp
;
5566 aSig
= (aSig
| UINT64_C(0x0010000000000000)) << 11;
5567 bSig
= (bSig
| UINT64_C(0x0010000000000000)) << 11;
5568 if ( expDiff
< 0 ) {
5569 if ( expDiff
< -1 ) return a
;
5572 q
= ( bSig
<= aSig
);
5573 if ( q
) aSig
-= bSig
;
5575 while ( 0 < expDiff
) {
5576 q
= estimateDiv128To64( aSig
, 0, bSig
);
5577 q
= ( 2 < q
) ? q
- 2 : 0;
5578 aSig
= - ( ( bSig
>>2 ) * q
);
5582 if ( 0 < expDiff
) {
5583 q
= estimateDiv128To64( aSig
, 0, bSig
);
5584 q
= ( 2 < q
) ? q
- 2 : 0;
5587 aSig
= ( ( aSig
>>1 )<<( expDiff
- 1 ) ) - bSig
* q
;
5594 alternateASig
= aSig
;
5597 } while ( 0 <= (int64_t) aSig
);
5598 sigMean
= aSig
+ alternateASig
;
5599 if ( ( sigMean
< 0 ) || ( ( sigMean
== 0 ) && ( q
& 1 ) ) ) {
5600 aSig
= alternateASig
;
5602 zSign
= ( (int64_t) aSig
< 0 );
5603 if ( zSign
) aSig
= - aSig
;
5604 return normalizeRoundAndPackFloat64(aSign
^ zSign
, bExp
, aSig
, status
);
5608 /*----------------------------------------------------------------------------
5609 | Returns the binary log of the double-precision floating-point value `a'.
5610 | The operation is performed according to the IEC/IEEE Standard for Binary
5611 | Floating-Point Arithmetic.
5612 *----------------------------------------------------------------------------*/
5613 float64
float64_log2(float64 a
, float_status
*status
)
5617 uint64_t aSig
, aSig0
, aSig1
, zSig
, i
;
5618 a
= float64_squash_input_denormal(a
, status
);
5620 aSig
= extractFloat64Frac( a
);
5621 aExp
= extractFloat64Exp( a
);
5622 aSign
= extractFloat64Sign( a
);
5625 if ( aSig
== 0 ) return packFloat64( 1, 0x7FF, 0 );
5626 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
5629 float_raise(float_flag_invalid
, status
);
5630 return float64_default_nan(status
);
5632 if ( aExp
== 0x7FF ) {
5634 return propagateFloat64NaN(a
, float64_zero
, status
);
5640 aSig
|= UINT64_C(0x0010000000000000);
5642 zSig
= (uint64_t)aExp
<< 52;
5643 for (i
= 1LL << 51; i
> 0; i
>>= 1) {
5644 mul64To128( aSig
, aSig
, &aSig0
, &aSig1
);
5645 aSig
= ( aSig0
<< 12 ) | ( aSig1
>> 52 );
5646 if ( aSig
& UINT64_C(0x0020000000000000) ) {
5654 return normalizeRoundAndPackFloat64(zSign
, 0x408, zSig
, status
);
5657 /*----------------------------------------------------------------------------
5658 | Returns the result of converting the extended double-precision floating-
5659 | point value `a' to the 32-bit two's complement integer format. The
5660 | conversion is performed according to the IEC/IEEE Standard for Binary
5661 | Floating-Point Arithmetic---which means in particular that the conversion
5662 | is rounded according to the current rounding mode. If `a' is a NaN, the
5663 | largest positive integer is returned. Otherwise, if the conversion
5664 | overflows, the largest integer with the same sign as `a' is returned.
5665 *----------------------------------------------------------------------------*/
5667 int32_t floatx80_to_int32(floatx80 a
, float_status
*status
)
5670 int32_t aExp
, shiftCount
;
5673 if (floatx80_invalid_encoding(a
)) {
5674 float_raise(float_flag_invalid
, status
);
5677 aSig
= extractFloatx80Frac( a
);
5678 aExp
= extractFloatx80Exp( a
);
5679 aSign
= extractFloatx80Sign( a
);
5680 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) aSign
= 0;
5681 shiftCount
= 0x4037 - aExp
;
5682 if ( shiftCount
<= 0 ) shiftCount
= 1;
5683 shift64RightJamming( aSig
, shiftCount
, &aSig
);
5684 return roundAndPackInt32(aSign
, aSig
, status
);
5688 /*----------------------------------------------------------------------------
5689 | Returns the result of converting the extended double-precision floating-
5690 | point value `a' to the 32-bit two's complement integer format. The
5691 | conversion is performed according to the IEC/IEEE Standard for Binary
5692 | Floating-Point Arithmetic, except that the conversion is always rounded
5693 | toward zero. If `a' is a NaN, the largest positive integer is returned.
5694 | Otherwise, if the conversion overflows, the largest integer with the same
5695 | sign as `a' is returned.
5696 *----------------------------------------------------------------------------*/
5698 int32_t floatx80_to_int32_round_to_zero(floatx80 a
, float_status
*status
)
5701 int32_t aExp
, shiftCount
;
5702 uint64_t aSig
, savedASig
;
5705 if (floatx80_invalid_encoding(a
)) {
5706 float_raise(float_flag_invalid
, status
);
5709 aSig
= extractFloatx80Frac( a
);
5710 aExp
= extractFloatx80Exp( a
);
5711 aSign
= extractFloatx80Sign( a
);
5712 if ( 0x401E < aExp
) {
5713 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) aSign
= 0;
5716 else if ( aExp
< 0x3FFF ) {
5718 float_raise(float_flag_inexact
, status
);
5722 shiftCount
= 0x403E - aExp
;
5724 aSig
>>= shiftCount
;
5726 if ( aSign
) z
= - z
;
5727 if ( ( z
< 0 ) ^ aSign
) {
5729 float_raise(float_flag_invalid
, status
);
5730 return aSign
? (int32_t) 0x80000000 : 0x7FFFFFFF;
5732 if ( ( aSig
<<shiftCount
) != savedASig
) {
5733 float_raise(float_flag_inexact
, status
);
5739 /*----------------------------------------------------------------------------
5740 | Returns the result of converting the extended double-precision floating-
5741 | point value `a' to the 64-bit two's complement integer format. The
5742 | conversion is performed according to the IEC/IEEE Standard for Binary
5743 | Floating-Point Arithmetic---which means in particular that the conversion
5744 | is rounded according to the current rounding mode. If `a' is a NaN,
5745 | the largest positive integer is returned. Otherwise, if the conversion
5746 | overflows, the largest integer with the same sign as `a' is returned.
5747 *----------------------------------------------------------------------------*/
5749 int64_t floatx80_to_int64(floatx80 a
, float_status
*status
)
5752 int32_t aExp
, shiftCount
;
5753 uint64_t aSig
, aSigExtra
;
5755 if (floatx80_invalid_encoding(a
)) {
5756 float_raise(float_flag_invalid
, status
);
5759 aSig
= extractFloatx80Frac( a
);
5760 aExp
= extractFloatx80Exp( a
);
5761 aSign
= extractFloatx80Sign( a
);
5762 shiftCount
= 0x403E - aExp
;
5763 if ( shiftCount
<= 0 ) {
5765 float_raise(float_flag_invalid
, status
);
5766 if (!aSign
|| floatx80_is_any_nan(a
)) {
5774 shift64ExtraRightJamming( aSig
, 0, shiftCount
, &aSig
, &aSigExtra
);
5776 return roundAndPackInt64(aSign
, aSig
, aSigExtra
, status
);
5780 /*----------------------------------------------------------------------------
5781 | Returns the result of converting the extended double-precision floating-
5782 | point value `a' to the 64-bit two's complement integer format. The
5783 | conversion is performed according to the IEC/IEEE Standard for Binary
5784 | Floating-Point Arithmetic, except that the conversion is always rounded
5785 | toward zero. If `a' is a NaN, the largest positive integer is returned.
5786 | Otherwise, if the conversion overflows, the largest integer with the same
5787 | sign as `a' is returned.
5788 *----------------------------------------------------------------------------*/
5790 int64_t floatx80_to_int64_round_to_zero(floatx80 a
, float_status
*status
)
5793 int32_t aExp
, shiftCount
;
5797 if (floatx80_invalid_encoding(a
)) {
5798 float_raise(float_flag_invalid
, status
);
5801 aSig
= extractFloatx80Frac( a
);
5802 aExp
= extractFloatx80Exp( a
);
5803 aSign
= extractFloatx80Sign( a
);
5804 shiftCount
= aExp
- 0x403E;
5805 if ( 0 <= shiftCount
) {
5806 aSig
&= UINT64_C(0x7FFFFFFFFFFFFFFF);
5807 if ( ( a
.high
!= 0xC03E ) || aSig
) {
5808 float_raise(float_flag_invalid
, status
);
5809 if ( ! aSign
|| ( ( aExp
== 0x7FFF ) && aSig
) ) {
5815 else if ( aExp
< 0x3FFF ) {
5817 float_raise(float_flag_inexact
, status
);
5821 z
= aSig
>>( - shiftCount
);
5822 if ( (uint64_t) ( aSig
<<( shiftCount
& 63 ) ) ) {
5823 float_raise(float_flag_inexact
, status
);
5825 if ( aSign
) z
= - z
;
5830 /*----------------------------------------------------------------------------
5831 | Returns the result of converting the extended double-precision floating-
5832 | point value `a' to the single-precision floating-point format. The
5833 | conversion is performed according to the IEC/IEEE Standard for Binary
5834 | Floating-Point Arithmetic.
5835 *----------------------------------------------------------------------------*/
5837 float32
floatx80_to_float32(floatx80 a
, float_status
*status
)
5843 if (floatx80_invalid_encoding(a
)) {
5844 float_raise(float_flag_invalid
, status
);
5845 return float32_default_nan(status
);
5847 aSig
= extractFloatx80Frac( a
);
5848 aExp
= extractFloatx80Exp( a
);
5849 aSign
= extractFloatx80Sign( a
);
5850 if ( aExp
== 0x7FFF ) {
5851 if ( (uint64_t) ( aSig
<<1 ) ) {
5852 float32 res
= commonNaNToFloat32(floatx80ToCommonNaN(a
, status
),
5854 return float32_silence_nan(res
, status
);
5856 return packFloat32( aSign
, 0xFF, 0 );
5858 shift64RightJamming( aSig
, 33, &aSig
);
5859 if ( aExp
|| aSig
) aExp
-= 0x3F81;
5860 return roundAndPackFloat32(aSign
, aExp
, aSig
, status
);
5864 /*----------------------------------------------------------------------------
5865 | Returns the result of converting the extended double-precision floating-
5866 | point value `a' to the double-precision floating-point format. The
5867 | conversion is performed according to the IEC/IEEE Standard for Binary
5868 | Floating-Point Arithmetic.
5869 *----------------------------------------------------------------------------*/
5871 float64
floatx80_to_float64(floatx80 a
, float_status
*status
)
5875 uint64_t aSig
, zSig
;
5877 if (floatx80_invalid_encoding(a
)) {
5878 float_raise(float_flag_invalid
, status
);
5879 return float64_default_nan(status
);
5881 aSig
= extractFloatx80Frac( a
);
5882 aExp
= extractFloatx80Exp( a
);
5883 aSign
= extractFloatx80Sign( a
);
5884 if ( aExp
== 0x7FFF ) {
5885 if ( (uint64_t) ( aSig
<<1 ) ) {
5886 float64 res
= commonNaNToFloat64(floatx80ToCommonNaN(a
, status
),
5888 return float64_silence_nan(res
, status
);
5890 return packFloat64( aSign
, 0x7FF, 0 );
5892 shift64RightJamming( aSig
, 1, &zSig
);
5893 if ( aExp
|| aSig
) aExp
-= 0x3C01;
5894 return roundAndPackFloat64(aSign
, aExp
, zSig
, status
);
5898 /*----------------------------------------------------------------------------
5899 | Returns the result of converting the extended double-precision floating-
5900 | point value `a' to the quadruple-precision floating-point format. The
5901 | conversion is performed according to the IEC/IEEE Standard for Binary
5902 | Floating-Point Arithmetic.
5903 *----------------------------------------------------------------------------*/
5905 float128
floatx80_to_float128(floatx80 a
, float_status
*status
)
5909 uint64_t aSig
, zSig0
, zSig1
;
5911 if (floatx80_invalid_encoding(a
)) {
5912 float_raise(float_flag_invalid
, status
);
5913 return float128_default_nan(status
);
5915 aSig
= extractFloatx80Frac( a
);
5916 aExp
= extractFloatx80Exp( a
);
5917 aSign
= extractFloatx80Sign( a
);
5918 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) {
5919 float128 res
= commonNaNToFloat128(floatx80ToCommonNaN(a
, status
),
5921 return float128_silence_nan(res
, status
);
5923 shift128Right( aSig
<<1, 0, 16, &zSig0
, &zSig1
);
5924 return packFloat128( aSign
, aExp
, zSig0
, zSig1
);
5928 /*----------------------------------------------------------------------------
5929 | Rounds the extended double-precision floating-point value `a'
5930 | to the precision provided by floatx80_rounding_precision and returns the
5931 | result as an extended double-precision floating-point value.
5932 | The operation is performed according to the IEC/IEEE Standard for Binary
5933 | Floating-Point Arithmetic.
5934 *----------------------------------------------------------------------------*/
5936 floatx80
floatx80_round(floatx80 a
, float_status
*status
)
5938 return roundAndPackFloatx80(status
->floatx80_rounding_precision
,
5939 extractFloatx80Sign(a
),
5940 extractFloatx80Exp(a
),
5941 extractFloatx80Frac(a
), 0, status
);
5944 /*----------------------------------------------------------------------------
5945 | Rounds the extended double-precision floating-point value `a' to an integer,
5946 | and returns the result as an extended quadruple-precision floating-point
5947 | value. The operation is performed according to the IEC/IEEE Standard for
5948 | Binary Floating-Point Arithmetic.
5949 *----------------------------------------------------------------------------*/
5951 floatx80
floatx80_round_to_int(floatx80 a
, float_status
*status
)
5955 uint64_t lastBitMask
, roundBitsMask
;
5958 if (floatx80_invalid_encoding(a
)) {
5959 float_raise(float_flag_invalid
, status
);
5960 return floatx80_default_nan(status
);
5962 aExp
= extractFloatx80Exp( a
);
5963 if ( 0x403E <= aExp
) {
5964 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) ) {
5965 return propagateFloatx80NaN(a
, a
, status
);
5969 if ( aExp
< 0x3FFF ) {
5971 && ( (uint64_t) ( extractFloatx80Frac( a
) ) == 0 ) ) {
5974 float_raise(float_flag_inexact
, status
);
5975 aSign
= extractFloatx80Sign( a
);
5976 switch (status
->float_rounding_mode
) {
5977 case float_round_nearest_even
:
5978 if ( ( aExp
== 0x3FFE ) && (uint64_t) ( extractFloatx80Frac( a
)<<1 )
5981 packFloatx80( aSign
, 0x3FFF, UINT64_C(0x8000000000000000));
5984 case float_round_ties_away
:
5985 if (aExp
== 0x3FFE) {
5986 return packFloatx80(aSign
, 0x3FFF, UINT64_C(0x8000000000000000));
5989 case float_round_down
:
5992 packFloatx80( 1, 0x3FFF, UINT64_C(0x8000000000000000))
5993 : packFloatx80( 0, 0, 0 );
5994 case float_round_up
:
5996 aSign
? packFloatx80( 1, 0, 0 )
5997 : packFloatx80( 0, 0x3FFF, UINT64_C(0x8000000000000000));
5999 case float_round_to_zero
:
6002 g_assert_not_reached();
6004 return packFloatx80( aSign
, 0, 0 );
6007 lastBitMask
<<= 0x403E - aExp
;
6008 roundBitsMask
= lastBitMask
- 1;
6010 switch (status
->float_rounding_mode
) {
6011 case float_round_nearest_even
:
6012 z
.low
+= lastBitMask
>>1;
6013 if ((z
.low
& roundBitsMask
) == 0) {
6014 z
.low
&= ~lastBitMask
;
6017 case float_round_ties_away
:
6018 z
.low
+= lastBitMask
>> 1;
6020 case float_round_to_zero
:
6022 case float_round_up
:
6023 if (!extractFloatx80Sign(z
)) {
6024 z
.low
+= roundBitsMask
;
6027 case float_round_down
:
6028 if (extractFloatx80Sign(z
)) {
6029 z
.low
+= roundBitsMask
;
6035 z
.low
&= ~ roundBitsMask
;
6038 z
.low
= UINT64_C(0x8000000000000000);
6040 if (z
.low
!= a
.low
) {
6041 float_raise(float_flag_inexact
, status
);
6047 /*----------------------------------------------------------------------------
6048 | Returns the result of adding the absolute values of the extended double-
6049 | precision floating-point values `a' and `b'. If `zSign' is 1, the sum is
6050 | negated before being returned. `zSign' is ignored if the result is a NaN.
6051 | The addition is performed according to the IEC/IEEE Standard for Binary
6052 | Floating-Point Arithmetic.
6053 *----------------------------------------------------------------------------*/
6055 static floatx80
addFloatx80Sigs(floatx80 a
, floatx80 b
, bool zSign
,
6056 float_status
*status
)
6058 int32_t aExp
, bExp
, zExp
;
6059 uint64_t aSig
, bSig
, zSig0
, zSig1
;
6062 aSig
= extractFloatx80Frac( a
);
6063 aExp
= extractFloatx80Exp( a
);
6064 bSig
= extractFloatx80Frac( b
);
6065 bExp
= extractFloatx80Exp( b
);
6066 expDiff
= aExp
- bExp
;
6067 if ( 0 < expDiff
) {
6068 if ( aExp
== 0x7FFF ) {
6069 if ((uint64_t)(aSig
<< 1)) {
6070 return propagateFloatx80NaN(a
, b
, status
);
6074 if ( bExp
== 0 ) --expDiff
;
6075 shift64ExtraRightJamming( bSig
, 0, expDiff
, &bSig
, &zSig1
);
6078 else if ( expDiff
< 0 ) {
6079 if ( bExp
== 0x7FFF ) {
6080 if ((uint64_t)(bSig
<< 1)) {
6081 return propagateFloatx80NaN(a
, b
, status
);
6083 return packFloatx80(zSign
,
6084 floatx80_infinity_high
,
6085 floatx80_infinity_low
);
6087 if ( aExp
== 0 ) ++expDiff
;
6088 shift64ExtraRightJamming( aSig
, 0, - expDiff
, &aSig
, &zSig1
);
6092 if ( aExp
== 0x7FFF ) {
6093 if ( (uint64_t) ( ( aSig
| bSig
)<<1 ) ) {
6094 return propagateFloatx80NaN(a
, b
, status
);
6099 zSig0
= aSig
+ bSig
;
6101 if ((aSig
| bSig
) & UINT64_C(0x8000000000000000) && zSig0
< aSig
) {
6102 /* At least one of the values is a pseudo-denormal,
6103 * and there is a carry out of the result. */
6108 return packFloatx80(zSign
, 0, 0);
6110 normalizeFloatx80Subnormal( zSig0
, &zExp
, &zSig0
);
6116 zSig0
= aSig
+ bSig
;
6117 if ( (int64_t) zSig0
< 0 ) goto roundAndPack
;
6119 shift64ExtraRightJamming( zSig0
, zSig1
, 1, &zSig0
, &zSig1
);
6120 zSig0
|= UINT64_C(0x8000000000000000);
6123 return roundAndPackFloatx80(status
->floatx80_rounding_precision
,
6124 zSign
, zExp
, zSig0
, zSig1
, status
);
6127 /*----------------------------------------------------------------------------
6128 | Returns the result of subtracting the absolute values of the extended
6129 | double-precision floating-point values `a' and `b'. If `zSign' is 1, the
6130 | difference is negated before being returned. `zSign' is ignored if the
6131 | result is a NaN. The subtraction is performed according to the IEC/IEEE
6132 | Standard for Binary Floating-Point Arithmetic.
6133 *----------------------------------------------------------------------------*/
6135 static floatx80
subFloatx80Sigs(floatx80 a
, floatx80 b
, bool zSign
,
6136 float_status
*status
)
6138 int32_t aExp
, bExp
, zExp
;
6139 uint64_t aSig
, bSig
, zSig0
, zSig1
;
6142 aSig
= extractFloatx80Frac( a
);
6143 aExp
= extractFloatx80Exp( a
);
6144 bSig
= extractFloatx80Frac( b
);
6145 bExp
= extractFloatx80Exp( b
);
6146 expDiff
= aExp
- bExp
;
6147 if ( 0 < expDiff
) goto aExpBigger
;
6148 if ( expDiff
< 0 ) goto bExpBigger
;
6149 if ( aExp
== 0x7FFF ) {
6150 if ( (uint64_t) ( ( aSig
| bSig
)<<1 ) ) {
6151 return propagateFloatx80NaN(a
, b
, status
);
6153 float_raise(float_flag_invalid
, status
);
6154 return floatx80_default_nan(status
);
6161 if ( bSig
< aSig
) goto aBigger
;
6162 if ( aSig
< bSig
) goto bBigger
;
6163 return packFloatx80(status
->float_rounding_mode
== float_round_down
, 0, 0);
6165 if ( bExp
== 0x7FFF ) {
6166 if ((uint64_t)(bSig
<< 1)) {
6167 return propagateFloatx80NaN(a
, b
, status
);
6169 return packFloatx80(zSign
^ 1, floatx80_infinity_high
,
6170 floatx80_infinity_low
);
6172 if ( aExp
== 0 ) ++expDiff
;
6173 shift128RightJamming( aSig
, 0, - expDiff
, &aSig
, &zSig1
);
6175 sub128( bSig
, 0, aSig
, zSig1
, &zSig0
, &zSig1
);
6178 goto normalizeRoundAndPack
;
6180 if ( aExp
== 0x7FFF ) {
6181 if ((uint64_t)(aSig
<< 1)) {
6182 return propagateFloatx80NaN(a
, b
, status
);
6186 if ( bExp
== 0 ) --expDiff
;
6187 shift128RightJamming( bSig
, 0, expDiff
, &bSig
, &zSig1
);
6189 sub128( aSig
, 0, bSig
, zSig1
, &zSig0
, &zSig1
);
6191 normalizeRoundAndPack
:
6192 return normalizeRoundAndPackFloatx80(status
->floatx80_rounding_precision
,
6193 zSign
, zExp
, zSig0
, zSig1
, status
);
6196 /*----------------------------------------------------------------------------
6197 | Returns the result of adding the extended double-precision floating-point
6198 | values `a' and `b'. The operation is performed according to the IEC/IEEE
6199 | Standard for Binary Floating-Point Arithmetic.
6200 *----------------------------------------------------------------------------*/
6202 floatx80
floatx80_add(floatx80 a
, floatx80 b
, float_status
*status
)
6206 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
6207 float_raise(float_flag_invalid
, status
);
6208 return floatx80_default_nan(status
);
6210 aSign
= extractFloatx80Sign( a
);
6211 bSign
= extractFloatx80Sign( b
);
6212 if ( aSign
== bSign
) {
6213 return addFloatx80Sigs(a
, b
, aSign
, status
);
6216 return subFloatx80Sigs(a
, b
, aSign
, status
);
6221 /*----------------------------------------------------------------------------
6222 | Returns the result of subtracting the extended double-precision floating-
6223 | point values `a' and `b'. The operation is performed according to the
6224 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6225 *----------------------------------------------------------------------------*/
6227 floatx80
floatx80_sub(floatx80 a
, floatx80 b
, float_status
*status
)
6231 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
6232 float_raise(float_flag_invalid
, status
);
6233 return floatx80_default_nan(status
);
6235 aSign
= extractFloatx80Sign( a
);
6236 bSign
= extractFloatx80Sign( b
);
6237 if ( aSign
== bSign
) {
6238 return subFloatx80Sigs(a
, b
, aSign
, status
);
6241 return addFloatx80Sigs(a
, b
, aSign
, status
);
6246 /*----------------------------------------------------------------------------
6247 | Returns the result of multiplying the extended double-precision floating-
6248 | point values `a' and `b'. The operation is performed according to the
6249 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6250 *----------------------------------------------------------------------------*/
6252 floatx80
floatx80_mul(floatx80 a
, floatx80 b
, float_status
*status
)
6254 bool aSign
, bSign
, zSign
;
6255 int32_t aExp
, bExp
, zExp
;
6256 uint64_t aSig
, bSig
, zSig0
, zSig1
;
6258 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
6259 float_raise(float_flag_invalid
, status
);
6260 return floatx80_default_nan(status
);
6262 aSig
= extractFloatx80Frac( a
);
6263 aExp
= extractFloatx80Exp( a
);
6264 aSign
= extractFloatx80Sign( a
);
6265 bSig
= extractFloatx80Frac( b
);
6266 bExp
= extractFloatx80Exp( b
);
6267 bSign
= extractFloatx80Sign( b
);
6268 zSign
= aSign
^ bSign
;
6269 if ( aExp
== 0x7FFF ) {
6270 if ( (uint64_t) ( aSig
<<1 )
6271 || ( ( bExp
== 0x7FFF ) && (uint64_t) ( bSig
<<1 ) ) ) {
6272 return propagateFloatx80NaN(a
, b
, status
);
6274 if ( ( bExp
| bSig
) == 0 ) goto invalid
;
6275 return packFloatx80(zSign
, floatx80_infinity_high
,
6276 floatx80_infinity_low
);
6278 if ( bExp
== 0x7FFF ) {
6279 if ((uint64_t)(bSig
<< 1)) {
6280 return propagateFloatx80NaN(a
, b
, status
);
6282 if ( ( aExp
| aSig
) == 0 ) {
6284 float_raise(float_flag_invalid
, status
);
6285 return floatx80_default_nan(status
);
6287 return packFloatx80(zSign
, floatx80_infinity_high
,
6288 floatx80_infinity_low
);
6291 if ( aSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
6292 normalizeFloatx80Subnormal( aSig
, &aExp
, &aSig
);
6295 if ( bSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
6296 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
6298 zExp
= aExp
+ bExp
- 0x3FFE;
6299 mul64To128( aSig
, bSig
, &zSig0
, &zSig1
);
6300 if ( 0 < (int64_t) zSig0
) {
6301 shortShift128Left( zSig0
, zSig1
, 1, &zSig0
, &zSig1
);
6304 return roundAndPackFloatx80(status
->floatx80_rounding_precision
,
6305 zSign
, zExp
, zSig0
, zSig1
, status
);
6308 /*----------------------------------------------------------------------------
6309 | Returns the result of dividing the extended double-precision floating-point
6310 | value `a' by the corresponding value `b'. The operation is performed
6311 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6312 *----------------------------------------------------------------------------*/
6314 floatx80
floatx80_div(floatx80 a
, floatx80 b
, float_status
*status
)
6316 bool aSign
, bSign
, zSign
;
6317 int32_t aExp
, bExp
, zExp
;
6318 uint64_t aSig
, bSig
, zSig0
, zSig1
;
6319 uint64_t rem0
, rem1
, rem2
, term0
, term1
, term2
;
6321 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
6322 float_raise(float_flag_invalid
, status
);
6323 return floatx80_default_nan(status
);
6325 aSig
= extractFloatx80Frac( a
);
6326 aExp
= extractFloatx80Exp( a
);
6327 aSign
= extractFloatx80Sign( a
);
6328 bSig
= extractFloatx80Frac( b
);
6329 bExp
= extractFloatx80Exp( b
);
6330 bSign
= extractFloatx80Sign( b
);
6331 zSign
= aSign
^ bSign
;
6332 if ( aExp
== 0x7FFF ) {
6333 if ((uint64_t)(aSig
<< 1)) {
6334 return propagateFloatx80NaN(a
, b
, status
);
6336 if ( bExp
== 0x7FFF ) {
6337 if ((uint64_t)(bSig
<< 1)) {
6338 return propagateFloatx80NaN(a
, b
, status
);
6342 return packFloatx80(zSign
, floatx80_infinity_high
,
6343 floatx80_infinity_low
);
6345 if ( bExp
== 0x7FFF ) {
6346 if ((uint64_t)(bSig
<< 1)) {
6347 return propagateFloatx80NaN(a
, b
, status
);
6349 return packFloatx80( zSign
, 0, 0 );
6353 if ( ( aExp
| aSig
) == 0 ) {
6355 float_raise(float_flag_invalid
, status
);
6356 return floatx80_default_nan(status
);
6358 float_raise(float_flag_divbyzero
, status
);
6359 return packFloatx80(zSign
, floatx80_infinity_high
,
6360 floatx80_infinity_low
);
6362 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
6365 if ( aSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
6366 normalizeFloatx80Subnormal( aSig
, &aExp
, &aSig
);
6368 zExp
= aExp
- bExp
+ 0x3FFE;
6370 if ( bSig
<= aSig
) {
6371 shift128Right( aSig
, 0, 1, &aSig
, &rem1
);
6374 zSig0
= estimateDiv128To64( aSig
, rem1
, bSig
);
6375 mul64To128( bSig
, zSig0
, &term0
, &term1
);
6376 sub128( aSig
, rem1
, term0
, term1
, &rem0
, &rem1
);
6377 while ( (int64_t) rem0
< 0 ) {
6379 add128( rem0
, rem1
, 0, bSig
, &rem0
, &rem1
);
6381 zSig1
= estimateDiv128To64( rem1
, 0, bSig
);
6382 if ( (uint64_t) ( zSig1
<<1 ) <= 8 ) {
6383 mul64To128( bSig
, zSig1
, &term1
, &term2
);
6384 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
6385 while ( (int64_t) rem1
< 0 ) {
6387 add128( rem1
, rem2
, 0, bSig
, &rem1
, &rem2
);
6389 zSig1
|= ( ( rem1
| rem2
) != 0 );
6391 return roundAndPackFloatx80(status
->floatx80_rounding_precision
,
6392 zSign
, zExp
, zSig0
, zSig1
, status
);
6395 /*----------------------------------------------------------------------------
6396 | Returns the remainder of the extended double-precision floating-point value
6397 | `a' with respect to the corresponding value `b'. The operation is performed
6398 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic,
6399 | if 'mod' is false; if 'mod' is true, return the remainder based on truncating
6400 | the quotient toward zero instead. '*quotient' is set to the low 64 bits of
6401 | the absolute value of the integer quotient.
6402 *----------------------------------------------------------------------------*/
6404 floatx80
floatx80_modrem(floatx80 a
, floatx80 b
, bool mod
, uint64_t *quotient
,
6405 float_status
*status
)
6408 int32_t aExp
, bExp
, expDiff
, aExpOrig
;
6409 uint64_t aSig0
, aSig1
, bSig
;
6410 uint64_t q
, term0
, term1
, alternateASig0
, alternateASig1
;
6413 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
6414 float_raise(float_flag_invalid
, status
);
6415 return floatx80_default_nan(status
);
6417 aSig0
= extractFloatx80Frac( a
);
6418 aExpOrig
= aExp
= extractFloatx80Exp( a
);
6419 aSign
= extractFloatx80Sign( a
);
6420 bSig
= extractFloatx80Frac( b
);
6421 bExp
= extractFloatx80Exp( b
);
6422 if ( aExp
== 0x7FFF ) {
6423 if ( (uint64_t) ( aSig0
<<1 )
6424 || ( ( bExp
== 0x7FFF ) && (uint64_t) ( bSig
<<1 ) ) ) {
6425 return propagateFloatx80NaN(a
, b
, status
);
6429 if ( bExp
== 0x7FFF ) {
6430 if ((uint64_t)(bSig
<< 1)) {
6431 return propagateFloatx80NaN(a
, b
, status
);
6433 if (aExp
== 0 && aSig0
>> 63) {
6435 * Pseudo-denormal argument must be returned in normalized
6438 return packFloatx80(aSign
, 1, aSig0
);
6445 float_raise(float_flag_invalid
, status
);
6446 return floatx80_default_nan(status
);
6448 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
6451 if ( aSig0
== 0 ) return a
;
6452 normalizeFloatx80Subnormal( aSig0
, &aExp
, &aSig0
);
6455 expDiff
= aExp
- bExp
;
6457 if ( expDiff
< 0 ) {
6458 if ( mod
|| expDiff
< -1 ) {
6459 if (aExp
== 1 && aExpOrig
== 0) {
6461 * Pseudo-denormal argument must be returned in
6464 return packFloatx80(aSign
, aExp
, aSig0
);
6468 shift128Right( aSig0
, 0, 1, &aSig0
, &aSig1
);
6471 *quotient
= q
= ( bSig
<= aSig0
);
6472 if ( q
) aSig0
-= bSig
;
6474 while ( 0 < expDiff
) {
6475 q
= estimateDiv128To64( aSig0
, aSig1
, bSig
);
6476 q
= ( 2 < q
) ? q
- 2 : 0;
6477 mul64To128( bSig
, q
, &term0
, &term1
);
6478 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
6479 shortShift128Left( aSig0
, aSig1
, 62, &aSig0
, &aSig1
);
6485 if ( 0 < expDiff
) {
6486 q
= estimateDiv128To64( aSig0
, aSig1
, bSig
);
6487 q
= ( 2 < q
) ? q
- 2 : 0;
6489 mul64To128( bSig
, q
<<( 64 - expDiff
), &term0
, &term1
);
6490 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
6491 shortShift128Left( 0, bSig
, 64 - expDiff
, &term0
, &term1
);
6492 while ( le128( term0
, term1
, aSig0
, aSig1
) ) {
6494 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
6497 *quotient
<<= expDiff
;
6508 sub128( term0
, term1
, aSig0
, aSig1
, &alternateASig0
, &alternateASig1
);
6509 if ( lt128( alternateASig0
, alternateASig1
, aSig0
, aSig1
)
6510 || ( eq128( alternateASig0
, alternateASig1
, aSig0
, aSig1
)
6513 aSig0
= alternateASig0
;
6514 aSig1
= alternateASig1
;
6520 normalizeRoundAndPackFloatx80(
6521 80, zSign
, bExp
+ expDiff
, aSig0
, aSig1
, status
);
6525 /*----------------------------------------------------------------------------
6526 | Returns the remainder of the extended double-precision floating-point value
6527 | `a' with respect to the corresponding value `b'. The operation is performed
6528 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6529 *----------------------------------------------------------------------------*/
6531 floatx80
floatx80_rem(floatx80 a
, floatx80 b
, float_status
*status
)
6534 return floatx80_modrem(a
, b
, false, "ient
, status
);
6537 /*----------------------------------------------------------------------------
6538 | Returns the remainder of the extended double-precision floating-point value
6539 | `a' with respect to the corresponding value `b', with the quotient truncated
6541 *----------------------------------------------------------------------------*/
6543 floatx80
floatx80_mod(floatx80 a
, floatx80 b
, float_status
*status
)
6546 return floatx80_modrem(a
, b
, true, "ient
, status
);
6549 /*----------------------------------------------------------------------------
6550 | Returns the square root of the extended double-precision floating-point
6551 | value `a'. The operation is performed according to the IEC/IEEE Standard
6552 | for Binary Floating-Point Arithmetic.
6553 *----------------------------------------------------------------------------*/
6555 floatx80
floatx80_sqrt(floatx80 a
, float_status
*status
)
6559 uint64_t aSig0
, aSig1
, zSig0
, zSig1
, doubleZSig0
;
6560 uint64_t rem0
, rem1
, rem2
, rem3
, term0
, term1
, term2
, term3
;
6562 if (floatx80_invalid_encoding(a
)) {
6563 float_raise(float_flag_invalid
, status
);
6564 return floatx80_default_nan(status
);
6566 aSig0
= extractFloatx80Frac( a
);
6567 aExp
= extractFloatx80Exp( a
);
6568 aSign
= extractFloatx80Sign( a
);
6569 if ( aExp
== 0x7FFF ) {
6570 if ((uint64_t)(aSig0
<< 1)) {
6571 return propagateFloatx80NaN(a
, a
, status
);
6573 if ( ! aSign
) return a
;
6577 if ( ( aExp
| aSig0
) == 0 ) return a
;
6579 float_raise(float_flag_invalid
, status
);
6580 return floatx80_default_nan(status
);
6583 if ( aSig0
== 0 ) return packFloatx80( 0, 0, 0 );
6584 normalizeFloatx80Subnormal( aSig0
, &aExp
, &aSig0
);
6586 zExp
= ( ( aExp
- 0x3FFF )>>1 ) + 0x3FFF;
6587 zSig0
= estimateSqrt32( aExp
, aSig0
>>32 );
6588 shift128Right( aSig0
, 0, 2 + ( aExp
& 1 ), &aSig0
, &aSig1
);
6589 zSig0
= estimateDiv128To64( aSig0
, aSig1
, zSig0
<<32 ) + ( zSig0
<<30 );
6590 doubleZSig0
= zSig0
<<1;
6591 mul64To128( zSig0
, zSig0
, &term0
, &term1
);
6592 sub128( aSig0
, aSig1
, term0
, term1
, &rem0
, &rem1
);
6593 while ( (int64_t) rem0
< 0 ) {
6596 add128( rem0
, rem1
, zSig0
>>63, doubleZSig0
| 1, &rem0
, &rem1
);
6598 zSig1
= estimateDiv128To64( rem1
, 0, doubleZSig0
);
6599 if ( ( zSig1
& UINT64_C(0x3FFFFFFFFFFFFFFF) ) <= 5 ) {
6600 if ( zSig1
== 0 ) zSig1
= 1;
6601 mul64To128( doubleZSig0
, zSig1
, &term1
, &term2
);
6602 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
6603 mul64To128( zSig1
, zSig1
, &term2
, &term3
);
6604 sub192( rem1
, rem2
, 0, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
6605 while ( (int64_t) rem1
< 0 ) {
6607 shortShift128Left( 0, zSig1
, 1, &term2
, &term3
);
6609 term2
|= doubleZSig0
;
6610 add192( rem1
, rem2
, rem3
, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
6612 zSig1
|= ( ( rem1
| rem2
| rem3
) != 0 );
6614 shortShift128Left( 0, zSig1
, 1, &zSig0
, &zSig1
);
6615 zSig0
|= doubleZSig0
;
6616 return roundAndPackFloatx80(status
->floatx80_rounding_precision
,
6617 0, zExp
, zSig0
, zSig1
, status
);
6620 /*----------------------------------------------------------------------------
6621 | Returns the result of converting the quadruple-precision floating-point
6622 | value `a' to the 32-bit two's complement integer format. The conversion
6623 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
6624 | Arithmetic---which means in particular that the conversion is rounded
6625 | according to the current rounding mode. If `a' is a NaN, the largest
6626 | positive integer is returned. Otherwise, if the conversion overflows, the
6627 | largest integer with the same sign as `a' is returned.
6628 *----------------------------------------------------------------------------*/
6630 int32_t float128_to_int32(float128 a
, float_status
*status
)
6633 int32_t aExp
, shiftCount
;
6634 uint64_t aSig0
, aSig1
;
6636 aSig1
= extractFloat128Frac1( a
);
6637 aSig0
= extractFloat128Frac0( a
);
6638 aExp
= extractFloat128Exp( a
);
6639 aSign
= extractFloat128Sign( a
);
6640 if ( ( aExp
== 0x7FFF ) && ( aSig0
| aSig1
) ) aSign
= 0;
6641 if ( aExp
) aSig0
|= UINT64_C(0x0001000000000000);
6642 aSig0
|= ( aSig1
!= 0 );
6643 shiftCount
= 0x4028 - aExp
;
6644 if ( 0 < shiftCount
) shift64RightJamming( aSig0
, shiftCount
, &aSig0
);
6645 return roundAndPackInt32(aSign
, aSig0
, status
);
6649 /*----------------------------------------------------------------------------
6650 | Returns the result of converting the quadruple-precision floating-point
6651 | value `a' to the 32-bit two's complement integer format. The conversion
6652 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
6653 | Arithmetic, except that the conversion is always rounded toward zero. If
6654 | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
6655 | conversion overflows, the largest integer with the same sign as `a' is
6657 *----------------------------------------------------------------------------*/
6659 int32_t float128_to_int32_round_to_zero(float128 a
, float_status
*status
)
6662 int32_t aExp
, shiftCount
;
6663 uint64_t aSig0
, aSig1
, savedASig
;
6666 aSig1
= extractFloat128Frac1( a
);
6667 aSig0
= extractFloat128Frac0( a
);
6668 aExp
= extractFloat128Exp( a
);
6669 aSign
= extractFloat128Sign( a
);
6670 aSig0
|= ( aSig1
!= 0 );
6671 if ( 0x401E < aExp
) {
6672 if ( ( aExp
== 0x7FFF ) && aSig0
) aSign
= 0;
6675 else if ( aExp
< 0x3FFF ) {
6676 if (aExp
|| aSig0
) {
6677 float_raise(float_flag_inexact
, status
);
6681 aSig0
|= UINT64_C(0x0001000000000000);
6682 shiftCount
= 0x402F - aExp
;
6684 aSig0
>>= shiftCount
;
6686 if ( aSign
) z
= - z
;
6687 if ( ( z
< 0 ) ^ aSign
) {
6689 float_raise(float_flag_invalid
, status
);
6690 return aSign
? INT32_MIN
: INT32_MAX
;
6692 if ( ( aSig0
<<shiftCount
) != savedASig
) {
6693 float_raise(float_flag_inexact
, status
);
6699 /*----------------------------------------------------------------------------
6700 | Returns the result of converting the quadruple-precision floating-point
6701 | value `a' to the 64-bit two's complement integer format. The conversion
6702 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
6703 | Arithmetic---which means in particular that the conversion is rounded
6704 | according to the current rounding mode. If `a' is a NaN, the largest
6705 | positive integer is returned. Otherwise, if the conversion overflows, the
6706 | largest integer with the same sign as `a' is returned.
6707 *----------------------------------------------------------------------------*/
6709 int64_t float128_to_int64(float128 a
, float_status
*status
)
6712 int32_t aExp
, shiftCount
;
6713 uint64_t aSig0
, aSig1
;
6715 aSig1
= extractFloat128Frac1( a
);
6716 aSig0
= extractFloat128Frac0( a
);
6717 aExp
= extractFloat128Exp( a
);
6718 aSign
= extractFloat128Sign( a
);
6719 if ( aExp
) aSig0
|= UINT64_C(0x0001000000000000);
6720 shiftCount
= 0x402F - aExp
;
6721 if ( shiftCount
<= 0 ) {
6722 if ( 0x403E < aExp
) {
6723 float_raise(float_flag_invalid
, status
);
6725 || ( ( aExp
== 0x7FFF )
6726 && ( aSig1
|| ( aSig0
!= UINT64_C(0x0001000000000000) ) )
6733 shortShift128Left( aSig0
, aSig1
, - shiftCount
, &aSig0
, &aSig1
);
6736 shift64ExtraRightJamming( aSig0
, aSig1
, shiftCount
, &aSig0
, &aSig1
);
6738 return roundAndPackInt64(aSign
, aSig0
, aSig1
, status
);
6742 /*----------------------------------------------------------------------------
6743 | Returns the result of converting the quadruple-precision floating-point
6744 | value `a' to the 64-bit two's complement integer format. The conversion
6745 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
6746 | Arithmetic, except that the conversion is always rounded toward zero.
6747 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
6748 | the conversion overflows, the largest integer with the same sign as `a' is
6750 *----------------------------------------------------------------------------*/
6752 int64_t float128_to_int64_round_to_zero(float128 a
, float_status
*status
)
6755 int32_t aExp
, shiftCount
;
6756 uint64_t aSig0
, aSig1
;
6759 aSig1
= extractFloat128Frac1( a
);
6760 aSig0
= extractFloat128Frac0( a
);
6761 aExp
= extractFloat128Exp( a
);
6762 aSign
= extractFloat128Sign( a
);
6763 if ( aExp
) aSig0
|= UINT64_C(0x0001000000000000);
6764 shiftCount
= aExp
- 0x402F;
6765 if ( 0 < shiftCount
) {
6766 if ( 0x403E <= aExp
) {
6767 aSig0
&= UINT64_C(0x0000FFFFFFFFFFFF);
6768 if ( ( a
.high
== UINT64_C(0xC03E000000000000) )
6769 && ( aSig1
< UINT64_C(0x0002000000000000) ) ) {
6771 float_raise(float_flag_inexact
, status
);
6775 float_raise(float_flag_invalid
, status
);
6776 if ( ! aSign
|| ( ( aExp
== 0x7FFF ) && ( aSig0
| aSig1
) ) ) {
6782 z
= ( aSig0
<<shiftCount
) | ( aSig1
>>( ( - shiftCount
) & 63 ) );
6783 if ( (uint64_t) ( aSig1
<<shiftCount
) ) {
6784 float_raise(float_flag_inexact
, status
);
6788 if ( aExp
< 0x3FFF ) {
6789 if ( aExp
| aSig0
| aSig1
) {
6790 float_raise(float_flag_inexact
, status
);
6794 z
= aSig0
>>( - shiftCount
);
6796 || ( shiftCount
&& (uint64_t) ( aSig0
<<( shiftCount
& 63 ) ) ) ) {
6797 float_raise(float_flag_inexact
, status
);
6800 if ( aSign
) z
= - z
;
6805 /*----------------------------------------------------------------------------
6806 | Returns the result of converting the quadruple-precision floating-point value
6807 | `a' to the 64-bit unsigned integer format. The conversion is
6808 | performed according to the IEC/IEEE Standard for Binary Floating-Point
6809 | Arithmetic---which means in particular that the conversion is rounded
6810 | according to the current rounding mode. If `a' is a NaN, the largest
6811 | positive integer is returned. If the conversion overflows, the
6812 | largest unsigned integer is returned. If 'a' is negative, the value is
6813 | rounded and zero is returned; negative values that do not round to zero
6814 | will raise the inexact exception.
6815 *----------------------------------------------------------------------------*/
6817 uint64_t float128_to_uint64(float128 a
, float_status
*status
)
6822 uint64_t aSig0
, aSig1
;
6824 aSig0
= extractFloat128Frac0(a
);
6825 aSig1
= extractFloat128Frac1(a
);
6826 aExp
= extractFloat128Exp(a
);
6827 aSign
= extractFloat128Sign(a
);
6828 if (aSign
&& (aExp
> 0x3FFE)) {
6829 float_raise(float_flag_invalid
, status
);
6830 if (float128_is_any_nan(a
)) {
6837 aSig0
|= UINT64_C(0x0001000000000000);
6839 shiftCount
= 0x402F - aExp
;
6840 if (shiftCount
<= 0) {
6841 if (0x403E < aExp
) {
6842 float_raise(float_flag_invalid
, status
);
6845 shortShift128Left(aSig0
, aSig1
, -shiftCount
, &aSig0
, &aSig1
);
6847 shift64ExtraRightJamming(aSig0
, aSig1
, shiftCount
, &aSig0
, &aSig1
);
6849 return roundAndPackUint64(aSign
, aSig0
, aSig1
, status
);
6852 uint64_t float128_to_uint64_round_to_zero(float128 a
, float_status
*status
)
6855 signed char current_rounding_mode
= status
->float_rounding_mode
;
6857 set_float_rounding_mode(float_round_to_zero
, status
);
6858 v
= float128_to_uint64(a
, status
);
6859 set_float_rounding_mode(current_rounding_mode
, status
);
6864 /*----------------------------------------------------------------------------
6865 | Returns the result of converting the quadruple-precision floating-point
6866 | value `a' to the 32-bit unsigned integer format. The conversion
6867 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
6868 | Arithmetic except that the conversion is always rounded toward zero.
6869 | If `a' is a NaN, the largest positive integer is returned. Otherwise,
6870 | if the conversion overflows, the largest unsigned integer is returned.
6871 | If 'a' is negative, the value is rounded and zero is returned; negative
6872 | values that do not round to zero will raise the inexact exception.
6873 *----------------------------------------------------------------------------*/
6875 uint32_t float128_to_uint32_round_to_zero(float128 a
, float_status
*status
)
6879 int old_exc_flags
= get_float_exception_flags(status
);
6881 v
= float128_to_uint64_round_to_zero(a
, status
);
6882 if (v
> 0xffffffff) {
6887 set_float_exception_flags(old_exc_flags
, status
);
6888 float_raise(float_flag_invalid
, status
);
6892 /*----------------------------------------------------------------------------
6893 | Returns the result of converting the quadruple-precision floating-point value
6894 | `a' to the 32-bit unsigned integer format. The conversion is
6895 | performed according to the IEC/IEEE Standard for Binary Floating-Point
6896 | Arithmetic---which means in particular that the conversion is rounded
6897 | according to the current rounding mode. If `a' is a NaN, the largest
6898 | positive integer is returned. If the conversion overflows, the
6899 | largest unsigned integer is returned. If 'a' is negative, the value is
6900 | rounded and zero is returned; negative values that do not round to zero
6901 | will raise the inexact exception.
6902 *----------------------------------------------------------------------------*/
6904 uint32_t float128_to_uint32(float128 a
, float_status
*status
)
6908 int old_exc_flags
= get_float_exception_flags(status
);
6910 v
= float128_to_uint64(a
, status
);
6911 if (v
> 0xffffffff) {
6916 set_float_exception_flags(old_exc_flags
, status
);
6917 float_raise(float_flag_invalid
, status
);
6921 /*----------------------------------------------------------------------------
6922 | Returns the result of converting the quadruple-precision floating-point
6923 | value `a' to the extended double-precision floating-point format. The
6924 | conversion is performed according to the IEC/IEEE Standard for Binary
6925 | Floating-Point Arithmetic.
6926 *----------------------------------------------------------------------------*/
6928 floatx80
float128_to_floatx80(float128 a
, float_status
*status
)
6932 uint64_t aSig0
, aSig1
;
6934 aSig1
= extractFloat128Frac1( a
);
6935 aSig0
= extractFloat128Frac0( a
);
6936 aExp
= extractFloat128Exp( a
);
6937 aSign
= extractFloat128Sign( a
);
6938 if ( aExp
== 0x7FFF ) {
6939 if ( aSig0
| aSig1
) {
6940 floatx80 res
= commonNaNToFloatx80(float128ToCommonNaN(a
, status
),
6942 return floatx80_silence_nan(res
, status
);
6944 return packFloatx80(aSign
, floatx80_infinity_high
,
6945 floatx80_infinity_low
);
6948 if ( ( aSig0
| aSig1
) == 0 ) return packFloatx80( aSign
, 0, 0 );
6949 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
6952 aSig0
|= UINT64_C(0x0001000000000000);
6954 shortShift128Left( aSig0
, aSig1
, 15, &aSig0
, &aSig1
);
6955 return roundAndPackFloatx80(80, aSign
, aExp
, aSig0
, aSig1
, status
);
6959 /*----------------------------------------------------------------------------
6960 | Rounds the quadruple-precision floating-point value `a' to an integer, and
6961 | returns the result as a quadruple-precision floating-point value. The
6962 | operation is performed according to the IEC/IEEE Standard for Binary
6963 | Floating-Point Arithmetic.
6964 *----------------------------------------------------------------------------*/
6966 float128
float128_round_to_int(float128 a
, float_status
*status
)
6970 uint64_t lastBitMask
, roundBitsMask
;
6973 aExp
= extractFloat128Exp( a
);
6974 if ( 0x402F <= aExp
) {
6975 if ( 0x406F <= aExp
) {
6976 if ( ( aExp
== 0x7FFF )
6977 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) )
6979 return propagateFloat128NaN(a
, a
, status
);
6984 lastBitMask
= ( lastBitMask
<<( 0x406E - aExp
) )<<1;
6985 roundBitsMask
= lastBitMask
- 1;
6987 switch (status
->float_rounding_mode
) {
6988 case float_round_nearest_even
:
6989 if ( lastBitMask
) {
6990 add128( z
.high
, z
.low
, 0, lastBitMask
>>1, &z
.high
, &z
.low
);
6991 if ( ( z
.low
& roundBitsMask
) == 0 ) z
.low
&= ~ lastBitMask
;
6994 if ( (int64_t) z
.low
< 0 ) {
6996 if ( (uint64_t) ( z
.low
<<1 ) == 0 ) z
.high
&= ~1;
7000 case float_round_ties_away
:
7002 add128(z
.high
, z
.low
, 0, lastBitMask
>> 1, &z
.high
, &z
.low
);
7004 if ((int64_t) z
.low
< 0) {
7009 case float_round_to_zero
:
7011 case float_round_up
:
7012 if (!extractFloat128Sign(z
)) {
7013 add128(z
.high
, z
.low
, 0, roundBitsMask
, &z
.high
, &z
.low
);
7016 case float_round_down
:
7017 if (extractFloat128Sign(z
)) {
7018 add128(z
.high
, z
.low
, 0, roundBitsMask
, &z
.high
, &z
.low
);
7021 case float_round_to_odd
:
7023 * Note that if lastBitMask == 0, the last bit is the lsb
7024 * of high, and roundBitsMask == -1.
7026 if ((lastBitMask
? z
.low
& lastBitMask
: z
.high
& 1) == 0) {
7027 add128(z
.high
, z
.low
, 0, roundBitsMask
, &z
.high
, &z
.low
);
7033 z
.low
&= ~ roundBitsMask
;
7036 if ( aExp
< 0x3FFF ) {
7037 if ( ( ( (uint64_t) ( a
.high
<<1 ) ) | a
.low
) == 0 ) return a
;
7038 float_raise(float_flag_inexact
, status
);
7039 aSign
= extractFloat128Sign( a
);
7040 switch (status
->float_rounding_mode
) {
7041 case float_round_nearest_even
:
7042 if ( ( aExp
== 0x3FFE )
7043 && ( extractFloat128Frac0( a
)
7044 | extractFloat128Frac1( a
) )
7046 return packFloat128( aSign
, 0x3FFF, 0, 0 );
7049 case float_round_ties_away
:
7050 if (aExp
== 0x3FFE) {
7051 return packFloat128(aSign
, 0x3FFF, 0, 0);
7054 case float_round_down
:
7056 aSign
? packFloat128( 1, 0x3FFF, 0, 0 )
7057 : packFloat128( 0, 0, 0, 0 );
7058 case float_round_up
:
7060 aSign
? packFloat128( 1, 0, 0, 0 )
7061 : packFloat128( 0, 0x3FFF, 0, 0 );
7063 case float_round_to_odd
:
7064 return packFloat128(aSign
, 0x3FFF, 0, 0);
7066 case float_round_to_zero
:
7069 return packFloat128( aSign
, 0, 0, 0 );
7072 lastBitMask
<<= 0x402F - aExp
;
7073 roundBitsMask
= lastBitMask
- 1;
7076 switch (status
->float_rounding_mode
) {
7077 case float_round_nearest_even
:
7078 z
.high
+= lastBitMask
>>1;
7079 if ( ( ( z
.high
& roundBitsMask
) | a
.low
) == 0 ) {
7080 z
.high
&= ~ lastBitMask
;
7083 case float_round_ties_away
:
7084 z
.high
+= lastBitMask
>>1;
7086 case float_round_to_zero
:
7088 case float_round_up
:
7089 if (!extractFloat128Sign(z
)) {
7090 z
.high
|= ( a
.low
!= 0 );
7091 z
.high
+= roundBitsMask
;
7094 case float_round_down
:
7095 if (extractFloat128Sign(z
)) {
7096 z
.high
|= (a
.low
!= 0);
7097 z
.high
+= roundBitsMask
;
7100 case float_round_to_odd
:
7101 if ((z
.high
& lastBitMask
) == 0) {
7102 z
.high
|= (a
.low
!= 0);
7103 z
.high
+= roundBitsMask
;
7109 z
.high
&= ~ roundBitsMask
;
7111 if ( ( z
.low
!= a
.low
) || ( z
.high
!= a
.high
) ) {
7112 float_raise(float_flag_inexact
, status
);
7118 /*----------------------------------------------------------------------------
7119 | Returns the remainder of the quadruple-precision floating-point value `a'
7120 | with respect to the corresponding value `b'. The operation is performed
7121 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
7122 *----------------------------------------------------------------------------*/
7124 float128
float128_rem(float128 a
, float128 b
, float_status
*status
)
7127 int32_t aExp
, bExp
, expDiff
;
7128 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, q
, term0
, term1
, term2
;
7129 uint64_t allZero
, alternateASig0
, alternateASig1
, sigMean1
;
7132 aSig1
= extractFloat128Frac1( a
);
7133 aSig0
= extractFloat128Frac0( a
);
7134 aExp
= extractFloat128Exp( a
);
7135 aSign
= extractFloat128Sign( a
);
7136 bSig1
= extractFloat128Frac1( b
);
7137 bSig0
= extractFloat128Frac0( b
);
7138 bExp
= extractFloat128Exp( b
);
7139 if ( aExp
== 0x7FFF ) {
7140 if ( ( aSig0
| aSig1
)
7141 || ( ( bExp
== 0x7FFF ) && ( bSig0
| bSig1
) ) ) {
7142 return propagateFloat128NaN(a
, b
, status
);
7146 if ( bExp
== 0x7FFF ) {
7147 if (bSig0
| bSig1
) {
7148 return propagateFloat128NaN(a
, b
, status
);
7153 if ( ( bSig0
| bSig1
) == 0 ) {
7155 float_raise(float_flag_invalid
, status
);
7156 return float128_default_nan(status
);
7158 normalizeFloat128Subnormal( bSig0
, bSig1
, &bExp
, &bSig0
, &bSig1
);
7161 if ( ( aSig0
| aSig1
) == 0 ) return a
;
7162 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
7164 expDiff
= aExp
- bExp
;
7165 if ( expDiff
< -1 ) return a
;
7167 aSig0
| UINT64_C(0x0001000000000000),
7169 15 - ( expDiff
< 0 ),
7174 bSig0
| UINT64_C(0x0001000000000000), bSig1
, 15, &bSig0
, &bSig1
);
7175 q
= le128( bSig0
, bSig1
, aSig0
, aSig1
);
7176 if ( q
) sub128( aSig0
, aSig1
, bSig0
, bSig1
, &aSig0
, &aSig1
);
7178 while ( 0 < expDiff
) {
7179 q
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
7180 q
= ( 4 < q
) ? q
- 4 : 0;
7181 mul128By64To192( bSig0
, bSig1
, q
, &term0
, &term1
, &term2
);
7182 shortShift192Left( term0
, term1
, term2
, 61, &term1
, &term2
, &allZero
);
7183 shortShift128Left( aSig0
, aSig1
, 61, &aSig0
, &allZero
);
7184 sub128( aSig0
, 0, term1
, term2
, &aSig0
, &aSig1
);
7187 if ( -64 < expDiff
) {
7188 q
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
7189 q
= ( 4 < q
) ? q
- 4 : 0;
7191 shift128Right( bSig0
, bSig1
, 12, &bSig0
, &bSig1
);
7193 if ( expDiff
< 0 ) {
7194 shift128Right( aSig0
, aSig1
, - expDiff
, &aSig0
, &aSig1
);
7197 shortShift128Left( aSig0
, aSig1
, expDiff
, &aSig0
, &aSig1
);
7199 mul128By64To192( bSig0
, bSig1
, q
, &term0
, &term1
, &term2
);
7200 sub128( aSig0
, aSig1
, term1
, term2
, &aSig0
, &aSig1
);
7203 shift128Right( aSig0
, aSig1
, 12, &aSig0
, &aSig1
);
7204 shift128Right( bSig0
, bSig1
, 12, &bSig0
, &bSig1
);
7207 alternateASig0
= aSig0
;
7208 alternateASig1
= aSig1
;
7210 sub128( aSig0
, aSig1
, bSig0
, bSig1
, &aSig0
, &aSig1
);
7211 } while ( 0 <= (int64_t) aSig0
);
7213 aSig0
, aSig1
, alternateASig0
, alternateASig1
, (uint64_t *)&sigMean0
, &sigMean1
);
7214 if ( ( sigMean0
< 0 )
7215 || ( ( ( sigMean0
| sigMean1
) == 0 ) && ( q
& 1 ) ) ) {
7216 aSig0
= alternateASig0
;
7217 aSig1
= alternateASig1
;
7219 zSign
= ( (int64_t) aSig0
< 0 );
7220 if ( zSign
) sub128( 0, 0, aSig0
, aSig1
, &aSig0
, &aSig1
);
7221 return normalizeRoundAndPackFloat128(aSign
^ zSign
, bExp
- 4, aSig0
, aSig1
,
7225 /*----------------------------------------------------------------------------
7226 | Returns the square root of the quadruple-precision floating-point value `a'.
7227 | The operation is performed according to the IEC/IEEE Standard for Binary
7228 | Floating-Point Arithmetic.
7229 *----------------------------------------------------------------------------*/
7231 float128
float128_sqrt(float128 a
, float_status
*status
)
7235 uint64_t aSig0
, aSig1
, zSig0
, zSig1
, zSig2
, doubleZSig0
;
7236 uint64_t rem0
, rem1
, rem2
, rem3
, term0
, term1
, term2
, term3
;
7238 aSig1
= extractFloat128Frac1( a
);
7239 aSig0
= extractFloat128Frac0( a
);
7240 aExp
= extractFloat128Exp( a
);
7241 aSign
= extractFloat128Sign( a
);
7242 if ( aExp
== 0x7FFF ) {
7243 if (aSig0
| aSig1
) {
7244 return propagateFloat128NaN(a
, a
, status
);
7246 if ( ! aSign
) return a
;
7250 if ( ( aExp
| aSig0
| aSig1
) == 0 ) return a
;
7252 float_raise(float_flag_invalid
, status
);
7253 return float128_default_nan(status
);
7256 if ( ( aSig0
| aSig1
) == 0 ) return packFloat128( 0, 0, 0, 0 );
7257 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
7259 zExp
= ( ( aExp
- 0x3FFF )>>1 ) + 0x3FFE;
7260 aSig0
|= UINT64_C(0x0001000000000000);
7261 zSig0
= estimateSqrt32( aExp
, aSig0
>>17 );
7262 shortShift128Left( aSig0
, aSig1
, 13 - ( aExp
& 1 ), &aSig0
, &aSig1
);
7263 zSig0
= estimateDiv128To64( aSig0
, aSig1
, zSig0
<<32 ) + ( zSig0
<<30 );
7264 doubleZSig0
= zSig0
<<1;
7265 mul64To128( zSig0
, zSig0
, &term0
, &term1
);
7266 sub128( aSig0
, aSig1
, term0
, term1
, &rem0
, &rem1
);
7267 while ( (int64_t) rem0
< 0 ) {
7270 add128( rem0
, rem1
, zSig0
>>63, doubleZSig0
| 1, &rem0
, &rem1
);
7272 zSig1
= estimateDiv128To64( rem1
, 0, doubleZSig0
);
7273 if ( ( zSig1
& 0x1FFF ) <= 5 ) {
7274 if ( zSig1
== 0 ) zSig1
= 1;
7275 mul64To128( doubleZSig0
, zSig1
, &term1
, &term2
);
7276 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
7277 mul64To128( zSig1
, zSig1
, &term2
, &term3
);
7278 sub192( rem1
, rem2
, 0, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
7279 while ( (int64_t) rem1
< 0 ) {
7281 shortShift128Left( 0, zSig1
, 1, &term2
, &term3
);
7283 term2
|= doubleZSig0
;
7284 add192( rem1
, rem2
, rem3
, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
7286 zSig1
|= ( ( rem1
| rem2
| rem3
) != 0 );
7288 shift128ExtraRightJamming( zSig0
, zSig1
, 0, 14, &zSig0
, &zSig1
, &zSig2
);
7289 return roundAndPackFloat128(0, zExp
, zSig0
, zSig1
, zSig2
, status
);
7293 static inline FloatRelation
7294 floatx80_compare_internal(floatx80 a
, floatx80 b
, bool is_quiet
,
7295 float_status
*status
)
7299 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
7300 float_raise(float_flag_invalid
, status
);
7301 return float_relation_unordered
;
7303 if (( ( extractFloatx80Exp( a
) == 0x7fff ) &&
7304 ( extractFloatx80Frac( a
)<<1 ) ) ||
7305 ( ( extractFloatx80Exp( b
) == 0x7fff ) &&
7306 ( extractFloatx80Frac( b
)<<1 ) )) {
7308 floatx80_is_signaling_nan(a
, status
) ||
7309 floatx80_is_signaling_nan(b
, status
)) {
7310 float_raise(float_flag_invalid
, status
);
7312 return float_relation_unordered
;
7314 aSign
= extractFloatx80Sign( a
);
7315 bSign
= extractFloatx80Sign( b
);
7316 if ( aSign
!= bSign
) {
7318 if ( ( ( (uint16_t) ( ( a
.high
| b
.high
) << 1 ) ) == 0) &&
7319 ( ( a
.low
| b
.low
) == 0 ) ) {
7321 return float_relation_equal
;
7323 return 1 - (2 * aSign
);
7326 /* Normalize pseudo-denormals before comparison. */
7327 if ((a
.high
& 0x7fff) == 0 && a
.low
& UINT64_C(0x8000000000000000)) {
7330 if ((b
.high
& 0x7fff) == 0 && b
.low
& UINT64_C(0x8000000000000000)) {
7333 if (a
.low
== b
.low
&& a
.high
== b
.high
) {
7334 return float_relation_equal
;
7336 return 1 - 2 * (aSign
^ ( lt128( a
.high
, a
.low
, b
.high
, b
.low
) ));
7341 FloatRelation
floatx80_compare(floatx80 a
, floatx80 b
, float_status
*status
)
7343 return floatx80_compare_internal(a
, b
, 0, status
);
7346 FloatRelation
floatx80_compare_quiet(floatx80 a
, floatx80 b
,
7347 float_status
*status
)
7349 return floatx80_compare_internal(a
, b
, 1, status
);
7352 static inline FloatRelation
7353 float128_compare_internal(float128 a
, float128 b
, bool is_quiet
,
7354 float_status
*status
)
7358 if (( ( extractFloat128Exp( a
) == 0x7fff ) &&
7359 ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) ) ||
7360 ( ( extractFloat128Exp( b
) == 0x7fff ) &&
7361 ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )) {
7363 float128_is_signaling_nan(a
, status
) ||
7364 float128_is_signaling_nan(b
, status
)) {
7365 float_raise(float_flag_invalid
, status
);
7367 return float_relation_unordered
;
7369 aSign
= extractFloat128Sign( a
);
7370 bSign
= extractFloat128Sign( b
);
7371 if ( aSign
!= bSign
) {
7372 if ( ( ( ( a
.high
| b
.high
)<<1 ) | a
.low
| b
.low
) == 0 ) {
7374 return float_relation_equal
;
7376 return 1 - (2 * aSign
);
7379 if (a
.low
== b
.low
&& a
.high
== b
.high
) {
7380 return float_relation_equal
;
7382 return 1 - 2 * (aSign
^ ( lt128( a
.high
, a
.low
, b
.high
, b
.low
) ));
7387 FloatRelation
float128_compare(float128 a
, float128 b
, float_status
*status
)
7389 return float128_compare_internal(a
, b
, 0, status
);
7392 FloatRelation
float128_compare_quiet(float128 a
, float128 b
,
7393 float_status
*status
)
7395 return float128_compare_internal(a
, b
, 1, status
);
7398 floatx80
floatx80_scalbn(floatx80 a
, int n
, float_status
*status
)
7404 if (floatx80_invalid_encoding(a
)) {
7405 float_raise(float_flag_invalid
, status
);
7406 return floatx80_default_nan(status
);
7408 aSig
= extractFloatx80Frac( a
);
7409 aExp
= extractFloatx80Exp( a
);
7410 aSign
= extractFloatx80Sign( a
);
7412 if ( aExp
== 0x7FFF ) {
7414 return propagateFloatx80NaN(a
, a
, status
);
7428 } else if (n
< -0x10000) {
7433 return normalizeRoundAndPackFloatx80(status
->floatx80_rounding_precision
,
7434 aSign
, aExp
, aSig
, 0, status
);
7437 float128
float128_scalbn(float128 a
, int n
, float_status
*status
)
7441 uint64_t aSig0
, aSig1
;
7443 aSig1
= extractFloat128Frac1( a
);
7444 aSig0
= extractFloat128Frac0( a
);
7445 aExp
= extractFloat128Exp( a
);
7446 aSign
= extractFloat128Sign( a
);
7447 if ( aExp
== 0x7FFF ) {
7448 if ( aSig0
| aSig1
) {
7449 return propagateFloat128NaN(a
, a
, status
);
7454 aSig0
|= UINT64_C(0x0001000000000000);
7455 } else if (aSig0
== 0 && aSig1
== 0) {
7463 } else if (n
< -0x10000) {
7468 return normalizeRoundAndPackFloat128( aSign
, aExp
, aSig0
, aSig1
7473 static void __attribute__((constructor
)) softfloat_init(void)
7475 union_float64 ua
, ub
, uc
, ur
;
7477 if (QEMU_NO_HARDFLOAT
) {
7481 * Test that the host's FMA is not obviously broken. For example,
7482 * glibc < 2.23 can perform an incorrect FMA on certain hosts; see
7483 * https://sourceware.org/bugzilla/show_bug.cgi?id=13304
7485 ua
.s
= 0x0020000000000001ULL
;
7486 ub
.s
= 0x3ca0000000000000ULL
;
7487 uc
.s
= 0x0020000000000000ULL
;
7488 ur
.h
= fma(ua
.h
, ub
.h
, uc
.h
);
7489 if (ur
.s
!= 0x0020000000000001ULL
) {
7490 force_soft_fma
= true;