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)
814 static bool parts64_round_to_int_normal(FloatParts64
*a
, FloatRoundMode rm
,
815 int scale
, int frac_size
);
816 static bool parts128_round_to_int_normal(FloatParts128
*a
, FloatRoundMode r
,
817 int scale
, int frac_size
);
819 #define parts_round_to_int_normal(A, R, C, F) \
820 PARTS_GENERIC_64_128(round_to_int_normal, A)(A, R, C, F)
822 static void parts64_round_to_int(FloatParts64
*a
, FloatRoundMode rm
,
823 int scale
, float_status
*s
,
824 const FloatFmt
*fmt
);
825 static void parts128_round_to_int(FloatParts128
*a
, FloatRoundMode r
,
826 int scale
, float_status
*s
,
827 const FloatFmt
*fmt
);
829 #define parts_round_to_int(A, R, C, S, F) \
830 PARTS_GENERIC_64_128(round_to_int, A)(A, R, C, S, F)
833 * Helper functions for softfloat-parts.c.inc, per-size operations.
836 #define FRAC_GENERIC_64_128(NAME, P) \
837 QEMU_GENERIC(P, (FloatParts128 *, frac128_##NAME), frac64_##NAME)
839 #define FRAC_GENERIC_64_128_256(NAME, P) \
840 QEMU_GENERIC(P, (FloatParts256 *, frac256_##NAME), \
841 (FloatParts128 *, frac128_##NAME), frac64_##NAME)
843 static bool frac64_add(FloatParts64
*r
, FloatParts64
*a
, FloatParts64
*b
)
845 return uadd64_overflow(a
->frac
, b
->frac
, &r
->frac
);
848 static bool frac128_add(FloatParts128
*r
, FloatParts128
*a
, FloatParts128
*b
)
851 r
->frac_lo
= uadd64_carry(a
->frac_lo
, b
->frac_lo
, &c
);
852 r
->frac_hi
= uadd64_carry(a
->frac_hi
, b
->frac_hi
, &c
);
856 static bool frac256_add(FloatParts256
*r
, FloatParts256
*a
, FloatParts256
*b
)
859 r
->frac_lo
= uadd64_carry(a
->frac_lo
, b
->frac_lo
, &c
);
860 r
->frac_lm
= uadd64_carry(a
->frac_lm
, b
->frac_lm
, &c
);
861 r
->frac_hm
= uadd64_carry(a
->frac_hm
, b
->frac_hm
, &c
);
862 r
->frac_hi
= uadd64_carry(a
->frac_hi
, b
->frac_hi
, &c
);
866 #define frac_add(R, A, B) FRAC_GENERIC_64_128_256(add, R)(R, A, B)
868 static bool frac64_addi(FloatParts64
*r
, FloatParts64
*a
, uint64_t c
)
870 return uadd64_overflow(a
->frac
, c
, &r
->frac
);
873 static bool frac128_addi(FloatParts128
*r
, FloatParts128
*a
, uint64_t c
)
875 c
= uadd64_overflow(a
->frac_lo
, c
, &r
->frac_lo
);
876 return uadd64_overflow(a
->frac_hi
, c
, &r
->frac_hi
);
879 #define frac_addi(R, A, C) FRAC_GENERIC_64_128(addi, R)(R, A, C)
881 static void frac64_allones(FloatParts64
*a
)
886 static void frac128_allones(FloatParts128
*a
)
888 a
->frac_hi
= a
->frac_lo
= -1;
891 #define frac_allones(A) FRAC_GENERIC_64_128(allones, A)(A)
893 static int frac64_cmp(FloatParts64
*a
, FloatParts64
*b
)
895 return a
->frac
== b
->frac
? 0 : a
->frac
< b
->frac
? -1 : 1;
898 static int frac128_cmp(FloatParts128
*a
, FloatParts128
*b
)
900 uint64_t ta
= a
->frac_hi
, tb
= b
->frac_hi
;
902 ta
= a
->frac_lo
, tb
= b
->frac_lo
;
907 return ta
< tb
? -1 : 1;
910 #define frac_cmp(A, B) FRAC_GENERIC_64_128(cmp, A)(A, B)
912 static void frac64_clear(FloatParts64
*a
)
917 static void frac128_clear(FloatParts128
*a
)
919 a
->frac_hi
= a
->frac_lo
= 0;
922 #define frac_clear(A) FRAC_GENERIC_64_128(clear, A)(A)
924 static bool frac64_div(FloatParts64
*a
, FloatParts64
*b
)
926 uint64_t n1
, n0
, r
, q
;
930 * We want a 2*N / N-bit division to produce exactly an N-bit
931 * result, so that we do not lose any precision and so that we
932 * do not have to renormalize afterward. If A.frac < B.frac,
933 * then division would produce an (N-1)-bit result; shift A left
934 * by one to produce the an N-bit result, and return true to
935 * decrement the exponent to match.
937 * The udiv_qrnnd algorithm that we're using requires normalization,
938 * i.e. the msb of the denominator must be set, which is already true.
940 ret
= a
->frac
< b
->frac
;
948 q
= udiv_qrnnd(&r
, n0
, n1
, b
->frac
);
950 /* Set lsb if there is a remainder, to set inexact. */
951 a
->frac
= q
| (r
!= 0);
956 static bool frac128_div(FloatParts128
*a
, FloatParts128
*b
)
958 uint64_t q0
, q1
, a0
, a1
, b0
, b1
;
959 uint64_t r0
, r1
, r2
, r3
, t0
, t1
, t2
, t3
;
962 a0
= a
->frac_hi
, a1
= a
->frac_lo
;
963 b0
= b
->frac_hi
, b1
= b
->frac_lo
;
965 ret
= lt128(a0
, a1
, b0
, b1
);
967 a1
= shr_double(a0
, a1
, 1);
971 /* Use 128/64 -> 64 division as estimate for 192/128 -> 128 division. */
972 q0
= estimateDiv128To64(a0
, a1
, b0
);
975 * Estimate is high because B1 was not included (unless B1 == 0).
976 * Reduce quotient and increase remainder until remainder is non-negative.
977 * This loop will execute 0 to 2 times.
979 mul128By64To192(b0
, b1
, q0
, &t0
, &t1
, &t2
);
980 sub192(a0
, a1
, 0, t0
, t1
, t2
, &r0
, &r1
, &r2
);
983 add192(r0
, r1
, r2
, 0, b0
, b1
, &r0
, &r1
, &r2
);
986 /* Repeat using the remainder, producing a second word of quotient. */
987 q1
= estimateDiv128To64(r1
, r2
, b0
);
988 mul128By64To192(b0
, b1
, q1
, &t1
, &t2
, &t3
);
989 sub192(r1
, r2
, 0, t1
, t2
, t3
, &r1
, &r2
, &r3
);
992 add192(r1
, r2
, r3
, 0, b0
, b1
, &r1
, &r2
, &r3
);
995 /* Any remainder indicates inexact; set sticky bit. */
996 q1
|= (r2
| r3
) != 0;
1003 #define frac_div(A, B) FRAC_GENERIC_64_128(div, A)(A, B)
1005 static bool frac64_eqz(FloatParts64
*a
)
1007 return a
->frac
== 0;
1010 static bool frac128_eqz(FloatParts128
*a
)
1012 return (a
->frac_hi
| a
->frac_lo
) == 0;
1015 #define frac_eqz(A) FRAC_GENERIC_64_128(eqz, A)(A)
1017 static void frac64_mulw(FloatParts128
*r
, FloatParts64
*a
, FloatParts64
*b
)
1019 mulu64(&r
->frac_lo
, &r
->frac_hi
, a
->frac
, b
->frac
);
1022 static void frac128_mulw(FloatParts256
*r
, FloatParts128
*a
, FloatParts128
*b
)
1024 mul128To256(a
->frac_hi
, a
->frac_lo
, b
->frac_hi
, b
->frac_lo
,
1025 &r
->frac_hi
, &r
->frac_hm
, &r
->frac_lm
, &r
->frac_lo
);
1028 #define frac_mulw(R, A, B) FRAC_GENERIC_64_128(mulw, A)(R, A, B)
1030 static void frac64_neg(FloatParts64
*a
)
1035 static void frac128_neg(FloatParts128
*a
)
1038 a
->frac_lo
= usub64_borrow(0, a
->frac_lo
, &c
);
1039 a
->frac_hi
= usub64_borrow(0, a
->frac_hi
, &c
);
1042 static void frac256_neg(FloatParts256
*a
)
1045 a
->frac_lo
= usub64_borrow(0, a
->frac_lo
, &c
);
1046 a
->frac_lm
= usub64_borrow(0, a
->frac_lm
, &c
);
1047 a
->frac_hm
= usub64_borrow(0, a
->frac_hm
, &c
);
1048 a
->frac_hi
= usub64_borrow(0, a
->frac_hi
, &c
);
1051 #define frac_neg(A) FRAC_GENERIC_64_128_256(neg, A)(A)
1053 static int frac64_normalize(FloatParts64
*a
)
1056 int shift
= clz64(a
->frac
);
1063 static int frac128_normalize(FloatParts128
*a
)
1066 int shl
= clz64(a
->frac_hi
);
1067 a
->frac_hi
= shl_double(a
->frac_hi
, a
->frac_lo
, shl
);
1070 } else if (a
->frac_lo
) {
1071 int shl
= clz64(a
->frac_lo
);
1072 a
->frac_hi
= a
->frac_lo
<< shl
;
1079 static int frac256_normalize(FloatParts256
*a
)
1081 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_hm
;
1082 uint64_t a2
= a
->frac_lm
, a3
= a
->frac_lo
;
1094 a0
= a1
, a1
= a2
, a2
= a3
, a3
= 0;
1097 a0
= a2
, a1
= a3
, a2
= 0, a3
= 0;
1100 a0
= a3
, a1
= 0, a2
= 0, a3
= 0;
1103 a0
= 0, a1
= 0, a2
= 0, a3
= 0;
1113 a0
= shl_double(a0
, a1
, shl
);
1114 a1
= shl_double(a1
, a2
, shl
);
1115 a2
= shl_double(a2
, a3
, shl
);
1126 #define frac_normalize(A) FRAC_GENERIC_64_128_256(normalize, A)(A)
1128 static void frac64_shl(FloatParts64
*a
, int c
)
1133 static void frac128_shl(FloatParts128
*a
, int c
)
1135 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_lo
;
1143 a0
= shl_double(a0
, a1
, c
);
1151 #define frac_shl(A, C) FRAC_GENERIC_64_128(shl, A)(A, C)
1153 static void frac64_shr(FloatParts64
*a
, int c
)
1158 static void frac128_shr(FloatParts128
*a
, int c
)
1160 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_lo
;
1168 a1
= shr_double(a0
, a1
, c
);
1176 #define frac_shr(A, C) FRAC_GENERIC_64_128(shr, A)(A, C)
1178 static void frac64_shrjam(FloatParts64
*a
, int c
)
1180 uint64_t a0
= a
->frac
;
1182 if (likely(c
!= 0)) {
1183 if (likely(c
< 64)) {
1184 a0
= (a0
>> c
) | (shr_double(a0
, 0, c
) != 0);
1192 static void frac128_shrjam(FloatParts128
*a
, int c
)
1194 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_lo
;
1195 uint64_t sticky
= 0;
1197 if (unlikely(c
== 0)) {
1199 } else if (likely(c
< 64)) {
1201 } else if (likely(c
< 128)) {
1215 sticky
|= shr_double(a1
, 0, c
);
1216 a1
= shr_double(a0
, a1
, c
);
1220 a
->frac_lo
= a1
| (sticky
!= 0);
1224 static void frac256_shrjam(FloatParts256
*a
, int c
)
1226 uint64_t a0
= a
->frac_hi
, a1
= a
->frac_hm
;
1227 uint64_t a2
= a
->frac_lm
, a3
= a
->frac_lo
;
1228 uint64_t sticky
= 0;
1230 if (unlikely(c
== 0)) {
1232 } else if (likely(c
< 64)) {
1234 } else if (likely(c
< 256)) {
1235 if (unlikely(c
& 128)) {
1237 a3
= a1
, a2
= a0
, a1
= 0, a0
= 0;
1239 if (unlikely(c
& 64)) {
1241 a3
= a2
, a2
= a1
, a1
= a0
, a0
= 0;
1248 sticky
= a0
| a1
| a2
| a3
;
1249 a0
= a1
= a2
= a3
= 0;
1253 sticky
|= shr_double(a3
, 0, c
);
1254 a3
= shr_double(a2
, a3
, c
);
1255 a2
= shr_double(a1
, a2
, c
);
1256 a1
= shr_double(a0
, a1
, c
);
1260 a
->frac_lo
= a3
| (sticky
!= 0);
1266 #define frac_shrjam(A, C) FRAC_GENERIC_64_128_256(shrjam, A)(A, C)
1268 static bool frac64_sub(FloatParts64
*r
, FloatParts64
*a
, FloatParts64
*b
)
1270 return usub64_overflow(a
->frac
, b
->frac
, &r
->frac
);
1273 static bool frac128_sub(FloatParts128
*r
, FloatParts128
*a
, FloatParts128
*b
)
1276 r
->frac_lo
= usub64_borrow(a
->frac_lo
, b
->frac_lo
, &c
);
1277 r
->frac_hi
= usub64_borrow(a
->frac_hi
, b
->frac_hi
, &c
);
1281 static bool frac256_sub(FloatParts256
*r
, FloatParts256
*a
, FloatParts256
*b
)
1284 r
->frac_lo
= usub64_borrow(a
->frac_lo
, b
->frac_lo
, &c
);
1285 r
->frac_lm
= usub64_borrow(a
->frac_lm
, b
->frac_lm
, &c
);
1286 r
->frac_hm
= usub64_borrow(a
->frac_hm
, b
->frac_hm
, &c
);
1287 r
->frac_hi
= usub64_borrow(a
->frac_hi
, b
->frac_hi
, &c
);
1291 #define frac_sub(R, A, B) FRAC_GENERIC_64_128_256(sub, R)(R, A, B)
1293 static void frac64_truncjam(FloatParts64
*r
, FloatParts128
*a
)
1295 r
->frac
= a
->frac_hi
| (a
->frac_lo
!= 0);
1298 static void frac128_truncjam(FloatParts128
*r
, FloatParts256
*a
)
1300 r
->frac_hi
= a
->frac_hi
;
1301 r
->frac_lo
= a
->frac_hm
| ((a
->frac_lm
| a
->frac_lo
) != 0);
1304 #define frac_truncjam(R, A) FRAC_GENERIC_64_128(truncjam, R)(R, A)
1306 static void frac64_widen(FloatParts128
*r
, FloatParts64
*a
)
1308 r
->frac_hi
= a
->frac
;
1312 static void frac128_widen(FloatParts256
*r
, FloatParts128
*a
)
1314 r
->frac_hi
= a
->frac_hi
;
1315 r
->frac_hm
= a
->frac_lo
;
1320 #define frac_widen(A, B) FRAC_GENERIC_64_128(widen, B)(A, B)
1322 #define partsN(NAME) glue(glue(glue(parts,N),_),NAME)
1323 #define FloatPartsN glue(FloatParts,N)
1324 #define FloatPartsW glue(FloatParts,W)
1329 #include "softfloat-parts-addsub.c.inc"
1330 #include "softfloat-parts.c.inc"
1337 #include "softfloat-parts-addsub.c.inc"
1338 #include "softfloat-parts.c.inc"
1344 #include "softfloat-parts-addsub.c.inc"
1353 * Pack/unpack routines with a specific FloatFmt.
1356 static void float16a_unpack_canonical(FloatParts64
*p
, float16 f
,
1357 float_status
*s
, const FloatFmt
*params
)
1359 float16_unpack_raw(p
, f
);
1360 parts_canonicalize(p
, s
, params
);
1363 static void float16_unpack_canonical(FloatParts64
*p
, float16 f
,
1366 float16a_unpack_canonical(p
, f
, s
, &float16_params
);
1369 static void bfloat16_unpack_canonical(FloatParts64
*p
, bfloat16 f
,
1372 bfloat16_unpack_raw(p
, f
);
1373 parts_canonicalize(p
, s
, &bfloat16_params
);
1376 static float16
float16a_round_pack_canonical(FloatParts64
*p
,
1378 const FloatFmt
*params
)
1380 parts_uncanon(p
, s
, params
);
1381 return float16_pack_raw(p
);
1384 static float16
float16_round_pack_canonical(FloatParts64
*p
,
1387 return float16a_round_pack_canonical(p
, s
, &float16_params
);
1390 static bfloat16
bfloat16_round_pack_canonical(FloatParts64
*p
,
1393 parts_uncanon(p
, s
, &bfloat16_params
);
1394 return bfloat16_pack_raw(p
);
1397 static void float32_unpack_canonical(FloatParts64
*p
, float32 f
,
1400 float32_unpack_raw(p
, f
);
1401 parts_canonicalize(p
, s
, &float32_params
);
1404 static float32
float32_round_pack_canonical(FloatParts64
*p
,
1407 parts_uncanon(p
, s
, &float32_params
);
1408 return float32_pack_raw(p
);
1411 static void float64_unpack_canonical(FloatParts64
*p
, float64 f
,
1414 float64_unpack_raw(p
, f
);
1415 parts_canonicalize(p
, s
, &float64_params
);
1418 static float64
float64_round_pack_canonical(FloatParts64
*p
,
1421 parts_uncanon(p
, s
, &float64_params
);
1422 return float64_pack_raw(p
);
1425 static void float128_unpack_canonical(FloatParts128
*p
, float128 f
,
1428 float128_unpack_raw(p
, f
);
1429 parts_canonicalize(p
, s
, &float128_params
);
1432 static float128
float128_round_pack_canonical(FloatParts128
*p
,
1435 parts_uncanon(p
, s
, &float128_params
);
1436 return float128_pack_raw(p
);
1440 * Addition and subtraction
1443 static float16 QEMU_FLATTEN
1444 float16_addsub(float16 a
, float16 b
, float_status
*status
, bool subtract
)
1446 FloatParts64 pa
, pb
, *pr
;
1448 float16_unpack_canonical(&pa
, a
, status
);
1449 float16_unpack_canonical(&pb
, b
, status
);
1450 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1452 return float16_round_pack_canonical(pr
, status
);
1455 float16
float16_add(float16 a
, float16 b
, float_status
*status
)
1457 return float16_addsub(a
, b
, status
, false);
1460 float16
float16_sub(float16 a
, float16 b
, float_status
*status
)
1462 return float16_addsub(a
, b
, status
, true);
1465 static float32 QEMU_SOFTFLOAT_ATTR
1466 soft_f32_addsub(float32 a
, float32 b
, float_status
*status
, bool subtract
)
1468 FloatParts64 pa
, pb
, *pr
;
1470 float32_unpack_canonical(&pa
, a
, status
);
1471 float32_unpack_canonical(&pb
, b
, status
);
1472 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1474 return float32_round_pack_canonical(pr
, status
);
1477 static float32
soft_f32_add(float32 a
, float32 b
, float_status
*status
)
1479 return soft_f32_addsub(a
, b
, status
, false);
1482 static float32
soft_f32_sub(float32 a
, float32 b
, float_status
*status
)
1484 return soft_f32_addsub(a
, b
, status
, true);
1487 static float64 QEMU_SOFTFLOAT_ATTR
1488 soft_f64_addsub(float64 a
, float64 b
, float_status
*status
, bool subtract
)
1490 FloatParts64 pa
, pb
, *pr
;
1492 float64_unpack_canonical(&pa
, a
, status
);
1493 float64_unpack_canonical(&pb
, b
, status
);
1494 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1496 return float64_round_pack_canonical(pr
, status
);
1499 static float64
soft_f64_add(float64 a
, float64 b
, float_status
*status
)
1501 return soft_f64_addsub(a
, b
, status
, false);
1504 static float64
soft_f64_sub(float64 a
, float64 b
, float_status
*status
)
1506 return soft_f64_addsub(a
, b
, status
, true);
1509 static float hard_f32_add(float a
, float b
)
1514 static float hard_f32_sub(float a
, float b
)
1519 static double hard_f64_add(double a
, double b
)
1524 static double hard_f64_sub(double a
, double b
)
1529 static bool f32_addsubmul_post(union_float32 a
, union_float32 b
)
1531 if (QEMU_HARDFLOAT_2F32_USE_FP
) {
1532 return !(fpclassify(a
.h
) == FP_ZERO
&& fpclassify(b
.h
) == FP_ZERO
);
1534 return !(float32_is_zero(a
.s
) && float32_is_zero(b
.s
));
1537 static bool f64_addsubmul_post(union_float64 a
, union_float64 b
)
1539 if (QEMU_HARDFLOAT_2F64_USE_FP
) {
1540 return !(fpclassify(a
.h
) == FP_ZERO
&& fpclassify(b
.h
) == FP_ZERO
);
1542 return !(float64_is_zero(a
.s
) && float64_is_zero(b
.s
));
1546 static float32
float32_addsub(float32 a
, float32 b
, float_status
*s
,
1547 hard_f32_op2_fn hard
, soft_f32_op2_fn soft
)
1549 return float32_gen2(a
, b
, s
, hard
, soft
,
1550 f32_is_zon2
, f32_addsubmul_post
);
1553 static float64
float64_addsub(float64 a
, float64 b
, float_status
*s
,
1554 hard_f64_op2_fn hard
, soft_f64_op2_fn soft
)
1556 return float64_gen2(a
, b
, s
, hard
, soft
,
1557 f64_is_zon2
, f64_addsubmul_post
);
1560 float32 QEMU_FLATTEN
1561 float32_add(float32 a
, float32 b
, float_status
*s
)
1563 return float32_addsub(a
, b
, s
, hard_f32_add
, soft_f32_add
);
1566 float32 QEMU_FLATTEN
1567 float32_sub(float32 a
, float32 b
, float_status
*s
)
1569 return float32_addsub(a
, b
, s
, hard_f32_sub
, soft_f32_sub
);
1572 float64 QEMU_FLATTEN
1573 float64_add(float64 a
, float64 b
, float_status
*s
)
1575 return float64_addsub(a
, b
, s
, hard_f64_add
, soft_f64_add
);
1578 float64 QEMU_FLATTEN
1579 float64_sub(float64 a
, float64 b
, float_status
*s
)
1581 return float64_addsub(a
, b
, s
, hard_f64_sub
, soft_f64_sub
);
1584 static bfloat16 QEMU_FLATTEN
1585 bfloat16_addsub(bfloat16 a
, bfloat16 b
, float_status
*status
, bool subtract
)
1587 FloatParts64 pa
, pb
, *pr
;
1589 bfloat16_unpack_canonical(&pa
, a
, status
);
1590 bfloat16_unpack_canonical(&pb
, b
, status
);
1591 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1593 return bfloat16_round_pack_canonical(pr
, status
);
1596 bfloat16
bfloat16_add(bfloat16 a
, bfloat16 b
, float_status
*status
)
1598 return bfloat16_addsub(a
, b
, status
, false);
1601 bfloat16
bfloat16_sub(bfloat16 a
, bfloat16 b
, float_status
*status
)
1603 return bfloat16_addsub(a
, b
, status
, true);
1606 static float128 QEMU_FLATTEN
1607 float128_addsub(float128 a
, float128 b
, float_status
*status
, bool subtract
)
1609 FloatParts128 pa
, pb
, *pr
;
1611 float128_unpack_canonical(&pa
, a
, status
);
1612 float128_unpack_canonical(&pb
, b
, status
);
1613 pr
= parts_addsub(&pa
, &pb
, status
, subtract
);
1615 return float128_round_pack_canonical(pr
, status
);
1618 float128
float128_add(float128 a
, float128 b
, float_status
*status
)
1620 return float128_addsub(a
, b
, status
, false);
1623 float128
float128_sub(float128 a
, float128 b
, float_status
*status
)
1625 return float128_addsub(a
, b
, status
, true);
1632 float16 QEMU_FLATTEN
float16_mul(float16 a
, float16 b
, float_status
*status
)
1634 FloatParts64 pa
, pb
, *pr
;
1636 float16_unpack_canonical(&pa
, a
, status
);
1637 float16_unpack_canonical(&pb
, b
, status
);
1638 pr
= parts_mul(&pa
, &pb
, status
);
1640 return float16_round_pack_canonical(pr
, status
);
1643 static float32 QEMU_SOFTFLOAT_ATTR
1644 soft_f32_mul(float32 a
, float32 b
, float_status
*status
)
1646 FloatParts64 pa
, pb
, *pr
;
1648 float32_unpack_canonical(&pa
, a
, status
);
1649 float32_unpack_canonical(&pb
, b
, status
);
1650 pr
= parts_mul(&pa
, &pb
, status
);
1652 return float32_round_pack_canonical(pr
, status
);
1655 static float64 QEMU_SOFTFLOAT_ATTR
1656 soft_f64_mul(float64 a
, float64 b
, float_status
*status
)
1658 FloatParts64 pa
, pb
, *pr
;
1660 float64_unpack_canonical(&pa
, a
, status
);
1661 float64_unpack_canonical(&pb
, b
, status
);
1662 pr
= parts_mul(&pa
, &pb
, status
);
1664 return float64_round_pack_canonical(pr
, status
);
1667 static float hard_f32_mul(float a
, float b
)
1672 static double hard_f64_mul(double a
, double b
)
1677 float32 QEMU_FLATTEN
1678 float32_mul(float32 a
, float32 b
, float_status
*s
)
1680 return float32_gen2(a
, b
, s
, hard_f32_mul
, soft_f32_mul
,
1681 f32_is_zon2
, f32_addsubmul_post
);
1684 float64 QEMU_FLATTEN
1685 float64_mul(float64 a
, float64 b
, float_status
*s
)
1687 return float64_gen2(a
, b
, s
, hard_f64_mul
, soft_f64_mul
,
1688 f64_is_zon2
, f64_addsubmul_post
);
1691 bfloat16 QEMU_FLATTEN
1692 bfloat16_mul(bfloat16 a
, bfloat16 b
, float_status
*status
)
1694 FloatParts64 pa
, pb
, *pr
;
1696 bfloat16_unpack_canonical(&pa
, a
, status
);
1697 bfloat16_unpack_canonical(&pb
, b
, status
);
1698 pr
= parts_mul(&pa
, &pb
, status
);
1700 return bfloat16_round_pack_canonical(pr
, status
);
1703 float128 QEMU_FLATTEN
1704 float128_mul(float128 a
, float128 b
, float_status
*status
)
1706 FloatParts128 pa
, pb
, *pr
;
1708 float128_unpack_canonical(&pa
, a
, status
);
1709 float128_unpack_canonical(&pb
, b
, status
);
1710 pr
= parts_mul(&pa
, &pb
, status
);
1712 return float128_round_pack_canonical(pr
, status
);
1716 * Fused multiply-add
1719 float16 QEMU_FLATTEN
float16_muladd(float16 a
, float16 b
, float16 c
,
1720 int flags
, float_status
*status
)
1722 FloatParts64 pa
, pb
, pc
, *pr
;
1724 float16_unpack_canonical(&pa
, a
, status
);
1725 float16_unpack_canonical(&pb
, b
, status
);
1726 float16_unpack_canonical(&pc
, c
, status
);
1727 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1729 return float16_round_pack_canonical(pr
, status
);
1732 static float32 QEMU_SOFTFLOAT_ATTR
1733 soft_f32_muladd(float32 a
, float32 b
, float32 c
, int flags
,
1734 float_status
*status
)
1736 FloatParts64 pa
, pb
, pc
, *pr
;
1738 float32_unpack_canonical(&pa
, a
, status
);
1739 float32_unpack_canonical(&pb
, b
, status
);
1740 float32_unpack_canonical(&pc
, c
, status
);
1741 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1743 return float32_round_pack_canonical(pr
, status
);
1746 static float64 QEMU_SOFTFLOAT_ATTR
1747 soft_f64_muladd(float64 a
, float64 b
, float64 c
, int flags
,
1748 float_status
*status
)
1750 FloatParts64 pa
, pb
, pc
, *pr
;
1752 float64_unpack_canonical(&pa
, a
, status
);
1753 float64_unpack_canonical(&pb
, b
, status
);
1754 float64_unpack_canonical(&pc
, c
, status
);
1755 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1757 return float64_round_pack_canonical(pr
, status
);
1760 static bool force_soft_fma
;
1762 float32 QEMU_FLATTEN
1763 float32_muladd(float32 xa
, float32 xb
, float32 xc
, int flags
, float_status
*s
)
1765 union_float32 ua
, ub
, uc
, ur
;
1771 if (unlikely(!can_use_fpu(s
))) {
1774 if (unlikely(flags
& float_muladd_halve_result
)) {
1778 float32_input_flush3(&ua
.s
, &ub
.s
, &uc
.s
, s
);
1779 if (unlikely(!f32_is_zon3(ua
, ub
, uc
))) {
1783 if (unlikely(force_soft_fma
)) {
1788 * When (a || b) == 0, there's no need to check for under/over flow,
1789 * since we know the addend is (normal || 0) and the product is 0.
1791 if (float32_is_zero(ua
.s
) || float32_is_zero(ub
.s
)) {
1795 prod_sign
= float32_is_neg(ua
.s
) ^ float32_is_neg(ub
.s
);
1796 prod_sign
^= !!(flags
& float_muladd_negate_product
);
1797 up
.s
= float32_set_sign(float32_zero
, prod_sign
);
1799 if (flags
& float_muladd_negate_c
) {
1804 union_float32 ua_orig
= ua
;
1805 union_float32 uc_orig
= uc
;
1807 if (flags
& float_muladd_negate_product
) {
1810 if (flags
& float_muladd_negate_c
) {
1814 ur
.h
= fmaf(ua
.h
, ub
.h
, uc
.h
);
1816 if (unlikely(f32_is_inf(ur
))) {
1817 float_raise(float_flag_overflow
, s
);
1818 } else if (unlikely(fabsf(ur
.h
) <= FLT_MIN
)) {
1824 if (flags
& float_muladd_negate_result
) {
1825 return float32_chs(ur
.s
);
1830 return soft_f32_muladd(ua
.s
, ub
.s
, uc
.s
, flags
, s
);
1833 float64 QEMU_FLATTEN
1834 float64_muladd(float64 xa
, float64 xb
, float64 xc
, int flags
, float_status
*s
)
1836 union_float64 ua
, ub
, uc
, ur
;
1842 if (unlikely(!can_use_fpu(s
))) {
1845 if (unlikely(flags
& float_muladd_halve_result
)) {
1849 float64_input_flush3(&ua
.s
, &ub
.s
, &uc
.s
, s
);
1850 if (unlikely(!f64_is_zon3(ua
, ub
, uc
))) {
1854 if (unlikely(force_soft_fma
)) {
1859 * When (a || b) == 0, there's no need to check for under/over flow,
1860 * since we know the addend is (normal || 0) and the product is 0.
1862 if (float64_is_zero(ua
.s
) || float64_is_zero(ub
.s
)) {
1866 prod_sign
= float64_is_neg(ua
.s
) ^ float64_is_neg(ub
.s
);
1867 prod_sign
^= !!(flags
& float_muladd_negate_product
);
1868 up
.s
= float64_set_sign(float64_zero
, prod_sign
);
1870 if (flags
& float_muladd_negate_c
) {
1875 union_float64 ua_orig
= ua
;
1876 union_float64 uc_orig
= uc
;
1878 if (flags
& float_muladd_negate_product
) {
1881 if (flags
& float_muladd_negate_c
) {
1885 ur
.h
= fma(ua
.h
, ub
.h
, uc
.h
);
1887 if (unlikely(f64_is_inf(ur
))) {
1888 float_raise(float_flag_overflow
, s
);
1889 } else if (unlikely(fabs(ur
.h
) <= FLT_MIN
)) {
1895 if (flags
& float_muladd_negate_result
) {
1896 return float64_chs(ur
.s
);
1901 return soft_f64_muladd(ua
.s
, ub
.s
, uc
.s
, flags
, s
);
1904 bfloat16 QEMU_FLATTEN
bfloat16_muladd(bfloat16 a
, bfloat16 b
, bfloat16 c
,
1905 int flags
, float_status
*status
)
1907 FloatParts64 pa
, pb
, pc
, *pr
;
1909 bfloat16_unpack_canonical(&pa
, a
, status
);
1910 bfloat16_unpack_canonical(&pb
, b
, status
);
1911 bfloat16_unpack_canonical(&pc
, c
, status
);
1912 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1914 return bfloat16_round_pack_canonical(pr
, status
);
1917 float128 QEMU_FLATTEN
float128_muladd(float128 a
, float128 b
, float128 c
,
1918 int flags
, float_status
*status
)
1920 FloatParts128 pa
, pb
, pc
, *pr
;
1922 float128_unpack_canonical(&pa
, a
, status
);
1923 float128_unpack_canonical(&pb
, b
, status
);
1924 float128_unpack_canonical(&pc
, c
, status
);
1925 pr
= parts_muladd(&pa
, &pb
, &pc
, flags
, status
);
1927 return float128_round_pack_canonical(pr
, status
);
1934 float16
float16_div(float16 a
, float16 b
, float_status
*status
)
1936 FloatParts64 pa
, pb
, *pr
;
1938 float16_unpack_canonical(&pa
, a
, status
);
1939 float16_unpack_canonical(&pb
, b
, status
);
1940 pr
= parts_div(&pa
, &pb
, status
);
1942 return float16_round_pack_canonical(pr
, status
);
1945 static float32 QEMU_SOFTFLOAT_ATTR
1946 soft_f32_div(float32 a
, float32 b
, float_status
*status
)
1948 FloatParts64 pa
, pb
, *pr
;
1950 float32_unpack_canonical(&pa
, a
, status
);
1951 float32_unpack_canonical(&pb
, b
, status
);
1952 pr
= parts_div(&pa
, &pb
, status
);
1954 return float32_round_pack_canonical(pr
, status
);
1957 static float64 QEMU_SOFTFLOAT_ATTR
1958 soft_f64_div(float64 a
, float64 b
, float_status
*status
)
1960 FloatParts64 pa
, pb
, *pr
;
1962 float64_unpack_canonical(&pa
, a
, status
);
1963 float64_unpack_canonical(&pb
, b
, status
);
1964 pr
= parts_div(&pa
, &pb
, status
);
1966 return float64_round_pack_canonical(pr
, status
);
1969 static float hard_f32_div(float a
, float b
)
1974 static double hard_f64_div(double a
, double b
)
1979 static bool f32_div_pre(union_float32 a
, union_float32 b
)
1981 if (QEMU_HARDFLOAT_2F32_USE_FP
) {
1982 return (fpclassify(a
.h
) == FP_NORMAL
|| fpclassify(a
.h
) == FP_ZERO
) &&
1983 fpclassify(b
.h
) == FP_NORMAL
;
1985 return float32_is_zero_or_normal(a
.s
) && float32_is_normal(b
.s
);
1988 static bool f64_div_pre(union_float64 a
, union_float64 b
)
1990 if (QEMU_HARDFLOAT_2F64_USE_FP
) {
1991 return (fpclassify(a
.h
) == FP_NORMAL
|| fpclassify(a
.h
) == FP_ZERO
) &&
1992 fpclassify(b
.h
) == FP_NORMAL
;
1994 return float64_is_zero_or_normal(a
.s
) && float64_is_normal(b
.s
);
1997 static bool f32_div_post(union_float32 a
, union_float32 b
)
1999 if (QEMU_HARDFLOAT_2F32_USE_FP
) {
2000 return fpclassify(a
.h
) != FP_ZERO
;
2002 return !float32_is_zero(a
.s
);
2005 static bool f64_div_post(union_float64 a
, union_float64 b
)
2007 if (QEMU_HARDFLOAT_2F64_USE_FP
) {
2008 return fpclassify(a
.h
) != FP_ZERO
;
2010 return !float64_is_zero(a
.s
);
2013 float32 QEMU_FLATTEN
2014 float32_div(float32 a
, float32 b
, float_status
*s
)
2016 return float32_gen2(a
, b
, s
, hard_f32_div
, soft_f32_div
,
2017 f32_div_pre
, f32_div_post
);
2020 float64 QEMU_FLATTEN
2021 float64_div(float64 a
, float64 b
, float_status
*s
)
2023 return float64_gen2(a
, b
, s
, hard_f64_div
, soft_f64_div
,
2024 f64_div_pre
, f64_div_post
);
2027 bfloat16 QEMU_FLATTEN
2028 bfloat16_div(bfloat16 a
, bfloat16 b
, float_status
*status
)
2030 FloatParts64 pa
, pb
, *pr
;
2032 bfloat16_unpack_canonical(&pa
, a
, status
);
2033 bfloat16_unpack_canonical(&pb
, b
, status
);
2034 pr
= parts_div(&pa
, &pb
, status
);
2036 return bfloat16_round_pack_canonical(pr
, status
);
2039 float128 QEMU_FLATTEN
2040 float128_div(float128 a
, float128 b
, float_status
*status
)
2042 FloatParts128 pa
, pb
, *pr
;
2044 float128_unpack_canonical(&pa
, a
, status
);
2045 float128_unpack_canonical(&pb
, b
, status
);
2046 pr
= parts_div(&pa
, &pb
, status
);
2048 return float128_round_pack_canonical(pr
, status
);
2052 * Float to Float conversions
2054 * Returns the result of converting one float format to another. The
2055 * conversion is performed according to the IEC/IEEE Standard for
2056 * Binary Floating-Point Arithmetic.
2058 * Usually this only needs to take care of raising invalid exceptions
2059 * and handling the conversion on NaNs.
2062 static void parts_float_to_ahp(FloatParts64
*a
, float_status
*s
)
2065 case float_class_qnan
:
2066 case float_class_snan
:
2068 * There is no NaN in the destination format. Raise Invalid
2069 * and return a zero with the sign of the input NaN.
2071 float_raise(float_flag_invalid
, s
);
2072 a
->cls
= float_class_zero
;
2075 case float_class_inf
:
2077 * There is no Inf in the destination format. Raise Invalid
2078 * and return the maximum normal with the correct sign.
2080 float_raise(float_flag_invalid
, s
);
2081 a
->cls
= float_class_normal
;
2082 a
->exp
= float16_params_ahp
.exp_max
;
2083 a
->frac
= MAKE_64BIT_MASK(float16_params_ahp
.frac_shift
,
2084 float16_params_ahp
.frac_size
+ 1);
2087 case float_class_normal
:
2088 case float_class_zero
:
2092 g_assert_not_reached();
2096 static void parts64_float_to_float(FloatParts64
*a
, float_status
*s
)
2098 if (is_nan(a
->cls
)) {
2099 parts_return_nan(a
, s
);
2103 static void parts128_float_to_float(FloatParts128
*a
, float_status
*s
)
2105 if (is_nan(a
->cls
)) {
2106 parts_return_nan(a
, s
);
2110 #define parts_float_to_float(P, S) \
2111 PARTS_GENERIC_64_128(float_to_float, P)(P, S)
2113 static void parts_float_to_float_narrow(FloatParts64
*a
, FloatParts128
*b
,
2120 if (a
->cls
== float_class_normal
) {
2121 frac_truncjam(a
, b
);
2122 } else if (is_nan(a
->cls
)) {
2123 /* Discard the low bits of the NaN. */
2124 a
->frac
= b
->frac_hi
;
2125 parts_return_nan(a
, s
);
2129 static void parts_float_to_float_widen(FloatParts128
*a
, FloatParts64
*b
,
2137 if (is_nan(a
->cls
)) {
2138 parts_return_nan(a
, s
);
2142 float32
float16_to_float32(float16 a
, bool ieee
, float_status
*s
)
2144 const FloatFmt
*fmt16
= ieee
? &float16_params
: &float16_params_ahp
;
2147 float16a_unpack_canonical(&p
, a
, s
, fmt16
);
2148 parts_float_to_float(&p
, s
);
2149 return float32_round_pack_canonical(&p
, s
);
2152 float64
float16_to_float64(float16 a
, bool ieee
, float_status
*s
)
2154 const FloatFmt
*fmt16
= ieee
? &float16_params
: &float16_params_ahp
;
2157 float16a_unpack_canonical(&p
, a
, s
, fmt16
);
2158 parts_float_to_float(&p
, s
);
2159 return float64_round_pack_canonical(&p
, s
);
2162 float16
float32_to_float16(float32 a
, bool ieee
, float_status
*s
)
2165 const FloatFmt
*fmt
;
2167 float32_unpack_canonical(&p
, a
, s
);
2169 parts_float_to_float(&p
, s
);
2170 fmt
= &float16_params
;
2172 parts_float_to_ahp(&p
, s
);
2173 fmt
= &float16_params_ahp
;
2175 return float16a_round_pack_canonical(&p
, s
, fmt
);
2178 static float64 QEMU_SOFTFLOAT_ATTR
2179 soft_float32_to_float64(float32 a
, float_status
*s
)
2183 float32_unpack_canonical(&p
, a
, s
);
2184 parts_float_to_float(&p
, s
);
2185 return float64_round_pack_canonical(&p
, s
);
2188 float64
float32_to_float64(float32 a
, float_status
*s
)
2190 if (likely(float32_is_normal(a
))) {
2191 /* Widening conversion can never produce inexact results. */
2197 } else if (float32_is_zero(a
)) {
2198 return float64_set_sign(float64_zero
, float32_is_neg(a
));
2200 return soft_float32_to_float64(a
, s
);
2204 float16
float64_to_float16(float64 a
, bool ieee
, float_status
*s
)
2207 const FloatFmt
*fmt
;
2209 float64_unpack_canonical(&p
, a
, s
);
2211 parts_float_to_float(&p
, s
);
2212 fmt
= &float16_params
;
2214 parts_float_to_ahp(&p
, s
);
2215 fmt
= &float16_params_ahp
;
2217 return float16a_round_pack_canonical(&p
, s
, fmt
);
2220 float32
float64_to_float32(float64 a
, float_status
*s
)
2224 float64_unpack_canonical(&p
, a
, s
);
2225 parts_float_to_float(&p
, s
);
2226 return float32_round_pack_canonical(&p
, s
);
2229 float32
bfloat16_to_float32(bfloat16 a
, float_status
*s
)
2233 bfloat16_unpack_canonical(&p
, a
, s
);
2234 parts_float_to_float(&p
, s
);
2235 return float32_round_pack_canonical(&p
, s
);
2238 float64
bfloat16_to_float64(bfloat16 a
, float_status
*s
)
2242 bfloat16_unpack_canonical(&p
, a
, s
);
2243 parts_float_to_float(&p
, s
);
2244 return float64_round_pack_canonical(&p
, s
);
2247 bfloat16
float32_to_bfloat16(float32 a
, float_status
*s
)
2251 float32_unpack_canonical(&p
, a
, s
);
2252 parts_float_to_float(&p
, s
);
2253 return bfloat16_round_pack_canonical(&p
, s
);
2256 bfloat16
float64_to_bfloat16(float64 a
, float_status
*s
)
2260 float64_unpack_canonical(&p
, a
, s
);
2261 parts_float_to_float(&p
, s
);
2262 return bfloat16_round_pack_canonical(&p
, s
);
2265 float32
float128_to_float32(float128 a
, float_status
*s
)
2270 float128_unpack_canonical(&p128
, a
, s
);
2271 parts_float_to_float_narrow(&p64
, &p128
, s
);
2272 return float32_round_pack_canonical(&p64
, s
);
2275 float64
float128_to_float64(float128 a
, float_status
*s
)
2280 float128_unpack_canonical(&p128
, a
, s
);
2281 parts_float_to_float_narrow(&p64
, &p128
, s
);
2282 return float64_round_pack_canonical(&p64
, s
);
2285 float128
float32_to_float128(float32 a
, float_status
*s
)
2290 float32_unpack_canonical(&p64
, a
, s
);
2291 parts_float_to_float_widen(&p128
, &p64
, s
);
2292 return float128_round_pack_canonical(&p128
, s
);
2295 float128
float64_to_float128(float64 a
, float_status
*s
)
2300 float64_unpack_canonical(&p64
, a
, s
);
2301 parts_float_to_float_widen(&p128
, &p64
, s
);
2302 return float128_round_pack_canonical(&p128
, s
);
2306 * Round to integral value
2309 float16
float16_round_to_int(float16 a
, float_status
*s
)
2313 float16_unpack_canonical(&p
, a
, s
);
2314 parts_round_to_int(&p
, s
->float_rounding_mode
, 0, s
, &float16_params
);
2315 return float16_round_pack_canonical(&p
, s
);
2318 float32
float32_round_to_int(float32 a
, float_status
*s
)
2322 float32_unpack_canonical(&p
, a
, s
);
2323 parts_round_to_int(&p
, s
->float_rounding_mode
, 0, s
, &float32_params
);
2324 return float32_round_pack_canonical(&p
, s
);
2327 float64
float64_round_to_int(float64 a
, float_status
*s
)
2331 float64_unpack_canonical(&p
, a
, s
);
2332 parts_round_to_int(&p
, s
->float_rounding_mode
, 0, s
, &float64_params
);
2333 return float64_round_pack_canonical(&p
, s
);
2336 bfloat16
bfloat16_round_to_int(bfloat16 a
, float_status
*s
)
2340 bfloat16_unpack_canonical(&p
, a
, s
);
2341 parts_round_to_int(&p
, s
->float_rounding_mode
, 0, s
, &bfloat16_params
);
2342 return bfloat16_round_pack_canonical(&p
, s
);
2345 float128
float128_round_to_int(float128 a
, float_status
*s
)
2349 float128_unpack_canonical(&p
, a
, s
);
2350 parts_round_to_int(&p
, s
->float_rounding_mode
, 0, s
, &float128_params
);
2351 return float128_round_pack_canonical(&p
, s
);
2355 * Returns the result of converting the floating-point value `a' to
2356 * the two's complement integer format. The conversion is performed
2357 * according to the IEC/IEEE Standard for Binary Floating-Point
2358 * Arithmetic---which means in particular that the conversion is
2359 * rounded according to the current rounding mode. If `a' is a NaN,
2360 * the largest positive integer is returned. Otherwise, if the
2361 * conversion overflows, the largest integer with the same sign as `a'
2365 static int64_t round_to_int_and_pack(FloatParts64 p
, FloatRoundMode rmode
,
2366 int scale
, int64_t min
, int64_t max
,
2373 case float_class_snan
:
2374 case float_class_qnan
:
2375 flags
= float_flag_invalid
;
2379 case float_class_inf
:
2380 flags
= float_flag_invalid
;
2381 r
= p
.sign
? min
: max
;
2384 case float_class_zero
:
2387 case float_class_normal
:
2388 /* TODO: 62 = N - 2, frac_size for rounding */
2389 if (parts_round_to_int_normal(&p
, rmode
, scale
, 62)) {
2390 flags
= float_flag_inexact
;
2393 if (p
.exp
<= DECOMPOSED_BINARY_POINT
) {
2394 r
= p
.frac
>> (DECOMPOSED_BINARY_POINT
- p
.exp
);
2399 if (r
<= -(uint64_t)min
) {
2402 flags
= float_flag_invalid
;
2405 } else if (r
> max
) {
2406 flags
= float_flag_invalid
;
2412 g_assert_not_reached();
2415 float_raise(flags
, s
);
2419 int8_t float16_to_int8_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2424 float16_unpack_canonical(&p
, a
, s
);
2425 return round_to_int_and_pack(p
, rmode
, scale
, INT8_MIN
, INT8_MAX
, s
);
2428 int16_t float16_to_int16_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2433 float16_unpack_canonical(&p
, a
, s
);
2434 return round_to_int_and_pack(p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2437 int32_t float16_to_int32_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2442 float16_unpack_canonical(&p
, a
, s
);
2443 return round_to_int_and_pack(p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2446 int64_t float16_to_int64_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2451 float16_unpack_canonical(&p
, a
, s
);
2452 return round_to_int_and_pack(p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2455 int16_t float32_to_int16_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2460 float32_unpack_canonical(&p
, a
, s
);
2461 return round_to_int_and_pack(p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2464 int32_t float32_to_int32_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2469 float32_unpack_canonical(&p
, a
, s
);
2470 return round_to_int_and_pack(p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2473 int64_t float32_to_int64_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2478 float32_unpack_canonical(&p
, a
, s
);
2479 return round_to_int_and_pack(p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2482 int16_t float64_to_int16_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2487 float64_unpack_canonical(&p
, a
, s
);
2488 return round_to_int_and_pack(p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2491 int32_t float64_to_int32_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2496 float64_unpack_canonical(&p
, a
, s
);
2497 return round_to_int_and_pack(p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2500 int64_t float64_to_int64_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2505 float64_unpack_canonical(&p
, a
, s
);
2506 return round_to_int_and_pack(p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2509 int8_t float16_to_int8(float16 a
, float_status
*s
)
2511 return float16_to_int8_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2514 int16_t float16_to_int16(float16 a
, float_status
*s
)
2516 return float16_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2519 int32_t float16_to_int32(float16 a
, float_status
*s
)
2521 return float16_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2524 int64_t float16_to_int64(float16 a
, float_status
*s
)
2526 return float16_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2529 int16_t float32_to_int16(float32 a
, float_status
*s
)
2531 return float32_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2534 int32_t float32_to_int32(float32 a
, float_status
*s
)
2536 return float32_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2539 int64_t float32_to_int64(float32 a
, float_status
*s
)
2541 return float32_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2544 int16_t float64_to_int16(float64 a
, float_status
*s
)
2546 return float64_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2549 int32_t float64_to_int32(float64 a
, float_status
*s
)
2551 return float64_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2554 int64_t float64_to_int64(float64 a
, float_status
*s
)
2556 return float64_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2559 int16_t float16_to_int16_round_to_zero(float16 a
, float_status
*s
)
2561 return float16_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2564 int32_t float16_to_int32_round_to_zero(float16 a
, float_status
*s
)
2566 return float16_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2569 int64_t float16_to_int64_round_to_zero(float16 a
, float_status
*s
)
2571 return float16_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2574 int16_t float32_to_int16_round_to_zero(float32 a
, float_status
*s
)
2576 return float32_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2579 int32_t float32_to_int32_round_to_zero(float32 a
, float_status
*s
)
2581 return float32_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2584 int64_t float32_to_int64_round_to_zero(float32 a
, float_status
*s
)
2586 return float32_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2589 int16_t float64_to_int16_round_to_zero(float64 a
, float_status
*s
)
2591 return float64_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2594 int32_t float64_to_int32_round_to_zero(float64 a
, float_status
*s
)
2596 return float64_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2599 int64_t float64_to_int64_round_to_zero(float64 a
, float_status
*s
)
2601 return float64_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2605 * Returns the result of converting the floating-point value `a' to
2606 * the two's complement integer format.
2609 int16_t bfloat16_to_int16_scalbn(bfloat16 a
, FloatRoundMode rmode
, int scale
,
2614 bfloat16_unpack_canonical(&p
, a
, s
);
2615 return round_to_int_and_pack(p
, rmode
, scale
, INT16_MIN
, INT16_MAX
, s
);
2618 int32_t bfloat16_to_int32_scalbn(bfloat16 a
, FloatRoundMode rmode
, int scale
,
2623 bfloat16_unpack_canonical(&p
, a
, s
);
2624 return round_to_int_and_pack(p
, rmode
, scale
, INT32_MIN
, INT32_MAX
, s
);
2627 int64_t bfloat16_to_int64_scalbn(bfloat16 a
, FloatRoundMode rmode
, int scale
,
2632 bfloat16_unpack_canonical(&p
, a
, s
);
2633 return round_to_int_and_pack(p
, rmode
, scale
, INT64_MIN
, INT64_MAX
, s
);
2636 int16_t bfloat16_to_int16(bfloat16 a
, float_status
*s
)
2638 return bfloat16_to_int16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2641 int32_t bfloat16_to_int32(bfloat16 a
, float_status
*s
)
2643 return bfloat16_to_int32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2646 int64_t bfloat16_to_int64(bfloat16 a
, float_status
*s
)
2648 return bfloat16_to_int64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2651 int16_t bfloat16_to_int16_round_to_zero(bfloat16 a
, float_status
*s
)
2653 return bfloat16_to_int16_scalbn(a
, float_round_to_zero
, 0, s
);
2656 int32_t bfloat16_to_int32_round_to_zero(bfloat16 a
, float_status
*s
)
2658 return bfloat16_to_int32_scalbn(a
, float_round_to_zero
, 0, s
);
2661 int64_t bfloat16_to_int64_round_to_zero(bfloat16 a
, float_status
*s
)
2663 return bfloat16_to_int64_scalbn(a
, float_round_to_zero
, 0, s
);
2667 * Returns the result of converting the floating-point value `a' to
2668 * the unsigned integer format. The conversion is performed according
2669 * to the IEC/IEEE Standard for Binary Floating-Point
2670 * Arithmetic---which means in particular that the conversion is
2671 * rounded according to the current rounding mode. If `a' is a NaN,
2672 * the largest unsigned integer is returned. Otherwise, if the
2673 * conversion overflows, the largest unsigned integer is returned. If
2674 * the 'a' is negative, the result is rounded and zero is returned;
2675 * values that do not round to zero will raise the inexact exception
2679 static uint64_t round_to_uint_and_pack(FloatParts64 p
, FloatRoundMode rmode
,
2680 int scale
, uint64_t max
,
2687 case float_class_snan
:
2688 case float_class_qnan
:
2689 flags
= float_flag_invalid
;
2693 case float_class_inf
:
2694 flags
= float_flag_invalid
;
2695 r
= p
.sign
? 0 : max
;
2698 case float_class_zero
:
2701 case float_class_normal
:
2702 /* TODO: 62 = N - 2, frac_size for rounding */
2703 if (parts_round_to_int_normal(&p
, rmode
, scale
, 62)) {
2704 flags
= float_flag_inexact
;
2705 if (p
.cls
== float_class_zero
) {
2712 flags
= float_flag_invalid
;
2714 } else if (p
.exp
> DECOMPOSED_BINARY_POINT
) {
2715 flags
= float_flag_invalid
;
2718 r
= p
.frac
>> (DECOMPOSED_BINARY_POINT
- p
.exp
);
2720 flags
= float_flag_invalid
;
2727 g_assert_not_reached();
2730 float_raise(flags
, s
);
2734 uint8_t float16_to_uint8_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2739 float16_unpack_canonical(&p
, a
, s
);
2740 return round_to_uint_and_pack(p
, rmode
, scale
, UINT8_MAX
, s
);
2743 uint16_t float16_to_uint16_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2748 float16_unpack_canonical(&p
, a
, s
);
2749 return round_to_uint_and_pack(p
, rmode
, scale
, UINT16_MAX
, s
);
2752 uint32_t float16_to_uint32_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2757 float16_unpack_canonical(&p
, a
, s
);
2758 return round_to_uint_and_pack(p
, rmode
, scale
, UINT32_MAX
, s
);
2761 uint64_t float16_to_uint64_scalbn(float16 a
, FloatRoundMode rmode
, int scale
,
2766 float16_unpack_canonical(&p
, a
, s
);
2767 return round_to_uint_and_pack(p
, rmode
, scale
, UINT64_MAX
, s
);
2770 uint16_t float32_to_uint16_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2775 float32_unpack_canonical(&p
, a
, s
);
2776 return round_to_uint_and_pack(p
, rmode
, scale
, UINT16_MAX
, s
);
2779 uint32_t float32_to_uint32_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2784 float32_unpack_canonical(&p
, a
, s
);
2785 return round_to_uint_and_pack(p
, rmode
, scale
, UINT32_MAX
, s
);
2788 uint64_t float32_to_uint64_scalbn(float32 a
, FloatRoundMode rmode
, int scale
,
2793 float32_unpack_canonical(&p
, a
, s
);
2794 return round_to_uint_and_pack(p
, rmode
, scale
, UINT64_MAX
, s
);
2797 uint16_t float64_to_uint16_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2802 float64_unpack_canonical(&p
, a
, s
);
2803 return round_to_uint_and_pack(p
, rmode
, scale
, UINT16_MAX
, s
);
2806 uint32_t float64_to_uint32_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2811 float64_unpack_canonical(&p
, a
, s
);
2812 return round_to_uint_and_pack(p
, rmode
, scale
, UINT32_MAX
, s
);
2815 uint64_t float64_to_uint64_scalbn(float64 a
, FloatRoundMode rmode
, int scale
,
2820 float64_unpack_canonical(&p
, a
, s
);
2821 return round_to_uint_and_pack(p
, rmode
, scale
, UINT64_MAX
, s
);
2824 uint8_t float16_to_uint8(float16 a
, float_status
*s
)
2826 return float16_to_uint8_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2829 uint16_t float16_to_uint16(float16 a
, float_status
*s
)
2831 return float16_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2834 uint32_t float16_to_uint32(float16 a
, float_status
*s
)
2836 return float16_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2839 uint64_t float16_to_uint64(float16 a
, float_status
*s
)
2841 return float16_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2844 uint16_t float32_to_uint16(float32 a
, float_status
*s
)
2846 return float32_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2849 uint32_t float32_to_uint32(float32 a
, float_status
*s
)
2851 return float32_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2854 uint64_t float32_to_uint64(float32 a
, float_status
*s
)
2856 return float32_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2859 uint16_t float64_to_uint16(float64 a
, float_status
*s
)
2861 return float64_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2864 uint32_t float64_to_uint32(float64 a
, float_status
*s
)
2866 return float64_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2869 uint64_t float64_to_uint64(float64 a
, float_status
*s
)
2871 return float64_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2874 uint16_t float16_to_uint16_round_to_zero(float16 a
, float_status
*s
)
2876 return float16_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
2879 uint32_t float16_to_uint32_round_to_zero(float16 a
, float_status
*s
)
2881 return float16_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
2884 uint64_t float16_to_uint64_round_to_zero(float16 a
, float_status
*s
)
2886 return float16_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
2889 uint16_t float32_to_uint16_round_to_zero(float32 a
, float_status
*s
)
2891 return float32_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
2894 uint32_t float32_to_uint32_round_to_zero(float32 a
, float_status
*s
)
2896 return float32_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
2899 uint64_t float32_to_uint64_round_to_zero(float32 a
, float_status
*s
)
2901 return float32_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
2904 uint16_t float64_to_uint16_round_to_zero(float64 a
, float_status
*s
)
2906 return float64_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
2909 uint32_t float64_to_uint32_round_to_zero(float64 a
, float_status
*s
)
2911 return float64_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
2914 uint64_t float64_to_uint64_round_to_zero(float64 a
, float_status
*s
)
2916 return float64_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
2920 * Returns the result of converting the bfloat16 value `a' to
2921 * the unsigned integer format.
2924 uint16_t bfloat16_to_uint16_scalbn(bfloat16 a
, FloatRoundMode rmode
,
2925 int scale
, float_status
*s
)
2929 bfloat16_unpack_canonical(&p
, a
, s
);
2930 return round_to_uint_and_pack(p
, rmode
, scale
, UINT16_MAX
, s
);
2933 uint32_t bfloat16_to_uint32_scalbn(bfloat16 a
, FloatRoundMode rmode
,
2934 int scale
, float_status
*s
)
2938 bfloat16_unpack_canonical(&p
, a
, s
);
2939 return round_to_uint_and_pack(p
, rmode
, scale
, UINT32_MAX
, s
);
2942 uint64_t bfloat16_to_uint64_scalbn(bfloat16 a
, FloatRoundMode rmode
,
2943 int scale
, float_status
*s
)
2947 bfloat16_unpack_canonical(&p
, a
, s
);
2948 return round_to_uint_and_pack(p
, rmode
, scale
, UINT64_MAX
, s
);
2951 uint16_t bfloat16_to_uint16(bfloat16 a
, float_status
*s
)
2953 return bfloat16_to_uint16_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2956 uint32_t bfloat16_to_uint32(bfloat16 a
, float_status
*s
)
2958 return bfloat16_to_uint32_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2961 uint64_t bfloat16_to_uint64(bfloat16 a
, float_status
*s
)
2963 return bfloat16_to_uint64_scalbn(a
, s
->float_rounding_mode
, 0, s
);
2966 uint16_t bfloat16_to_uint16_round_to_zero(bfloat16 a
, float_status
*s
)
2968 return bfloat16_to_uint16_scalbn(a
, float_round_to_zero
, 0, s
);
2971 uint32_t bfloat16_to_uint32_round_to_zero(bfloat16 a
, float_status
*s
)
2973 return bfloat16_to_uint32_scalbn(a
, float_round_to_zero
, 0, s
);
2976 uint64_t bfloat16_to_uint64_round_to_zero(bfloat16 a
, float_status
*s
)
2978 return bfloat16_to_uint64_scalbn(a
, float_round_to_zero
, 0, s
);
2982 * Integer to float conversions
2984 * Returns the result of converting the two's complement integer `a'
2985 * to the floating-point format. The conversion is performed according
2986 * to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2989 static FloatParts64
int_to_float(int64_t a
, int scale
, float_status
*status
)
2991 FloatParts64 r
= { .sign
= false };
2994 r
.cls
= float_class_zero
;
2999 r
.cls
= float_class_normal
;
3005 scale
= MIN(MAX(scale
, -0x10000), 0x10000);
3007 r
.exp
= DECOMPOSED_BINARY_POINT
- shift
+ scale
;
3008 r
.frac
= f
<< shift
;
3014 float16
int64_to_float16_scalbn(int64_t a
, int scale
, float_status
*status
)
3016 FloatParts64 pa
= int_to_float(a
, scale
, status
);
3017 return float16_round_pack_canonical(&pa
, status
);
3020 float16
int32_to_float16_scalbn(int32_t a
, int scale
, float_status
*status
)
3022 return int64_to_float16_scalbn(a
, scale
, status
);
3025 float16
int16_to_float16_scalbn(int16_t a
, int scale
, float_status
*status
)
3027 return int64_to_float16_scalbn(a
, scale
, status
);
3030 float16
int64_to_float16(int64_t a
, float_status
*status
)
3032 return int64_to_float16_scalbn(a
, 0, status
);
3035 float16
int32_to_float16(int32_t a
, float_status
*status
)
3037 return int64_to_float16_scalbn(a
, 0, status
);
3040 float16
int16_to_float16(int16_t a
, float_status
*status
)
3042 return int64_to_float16_scalbn(a
, 0, status
);
3045 float16
int8_to_float16(int8_t a
, float_status
*status
)
3047 return int64_to_float16_scalbn(a
, 0, status
);
3050 float32
int64_to_float32_scalbn(int64_t a
, int scale
, float_status
*status
)
3052 FloatParts64 pa
= int_to_float(a
, scale
, status
);
3053 return float32_round_pack_canonical(&pa
, status
);
3056 float32
int32_to_float32_scalbn(int32_t a
, int scale
, float_status
*status
)
3058 return int64_to_float32_scalbn(a
, scale
, status
);
3061 float32
int16_to_float32_scalbn(int16_t a
, int scale
, float_status
*status
)
3063 return int64_to_float32_scalbn(a
, scale
, status
);
3066 float32
int64_to_float32(int64_t a
, float_status
*status
)
3068 return int64_to_float32_scalbn(a
, 0, status
);
3071 float32
int32_to_float32(int32_t a
, float_status
*status
)
3073 return int64_to_float32_scalbn(a
, 0, status
);
3076 float32
int16_to_float32(int16_t a
, float_status
*status
)
3078 return int64_to_float32_scalbn(a
, 0, status
);
3081 float64
int64_to_float64_scalbn(int64_t a
, int scale
, float_status
*status
)
3083 FloatParts64 pa
= int_to_float(a
, scale
, status
);
3084 return float64_round_pack_canonical(&pa
, status
);
3087 float64
int32_to_float64_scalbn(int32_t a
, int scale
, float_status
*status
)
3089 return int64_to_float64_scalbn(a
, scale
, status
);
3092 float64
int16_to_float64_scalbn(int16_t a
, int scale
, float_status
*status
)
3094 return int64_to_float64_scalbn(a
, scale
, status
);
3097 float64
int64_to_float64(int64_t a
, float_status
*status
)
3099 return int64_to_float64_scalbn(a
, 0, status
);
3102 float64
int32_to_float64(int32_t a
, float_status
*status
)
3104 return int64_to_float64_scalbn(a
, 0, status
);
3107 float64
int16_to_float64(int16_t a
, float_status
*status
)
3109 return int64_to_float64_scalbn(a
, 0, status
);
3113 * Returns the result of converting the two's complement integer `a'
3114 * to the bfloat16 format.
3117 bfloat16
int64_to_bfloat16_scalbn(int64_t a
, int scale
, float_status
*status
)
3119 FloatParts64 pa
= int_to_float(a
, scale
, status
);
3120 return bfloat16_round_pack_canonical(&pa
, status
);
3123 bfloat16
int32_to_bfloat16_scalbn(int32_t a
, int scale
, float_status
*status
)
3125 return int64_to_bfloat16_scalbn(a
, scale
, status
);
3128 bfloat16
int16_to_bfloat16_scalbn(int16_t a
, int scale
, float_status
*status
)
3130 return int64_to_bfloat16_scalbn(a
, scale
, status
);
3133 bfloat16
int64_to_bfloat16(int64_t a
, float_status
*status
)
3135 return int64_to_bfloat16_scalbn(a
, 0, status
);
3138 bfloat16
int32_to_bfloat16(int32_t a
, float_status
*status
)
3140 return int64_to_bfloat16_scalbn(a
, 0, status
);
3143 bfloat16
int16_to_bfloat16(int16_t a
, float_status
*status
)
3145 return int64_to_bfloat16_scalbn(a
, 0, status
);
3149 * Unsigned Integer to float conversions
3151 * Returns the result of converting the unsigned integer `a' to the
3152 * floating-point format. The conversion is performed according to the
3153 * IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3156 static FloatParts64
uint_to_float(uint64_t a
, int scale
, float_status
*status
)
3158 FloatParts64 r
= { .sign
= false };
3162 r
.cls
= float_class_zero
;
3164 scale
= MIN(MAX(scale
, -0x10000), 0x10000);
3166 r
.cls
= float_class_normal
;
3167 r
.exp
= DECOMPOSED_BINARY_POINT
- shift
+ scale
;
3168 r
.frac
= a
<< shift
;
3174 float16
uint64_to_float16_scalbn(uint64_t a
, int scale
, float_status
*status
)
3176 FloatParts64 pa
= uint_to_float(a
, scale
, status
);
3177 return float16_round_pack_canonical(&pa
, status
);
3180 float16
uint32_to_float16_scalbn(uint32_t a
, int scale
, float_status
*status
)
3182 return uint64_to_float16_scalbn(a
, scale
, status
);
3185 float16
uint16_to_float16_scalbn(uint16_t a
, int scale
, float_status
*status
)
3187 return uint64_to_float16_scalbn(a
, scale
, status
);
3190 float16
uint64_to_float16(uint64_t a
, float_status
*status
)
3192 return uint64_to_float16_scalbn(a
, 0, status
);
3195 float16
uint32_to_float16(uint32_t a
, float_status
*status
)
3197 return uint64_to_float16_scalbn(a
, 0, status
);
3200 float16
uint16_to_float16(uint16_t a
, float_status
*status
)
3202 return uint64_to_float16_scalbn(a
, 0, status
);
3205 float16
uint8_to_float16(uint8_t a
, float_status
*status
)
3207 return uint64_to_float16_scalbn(a
, 0, status
);
3210 float32
uint64_to_float32_scalbn(uint64_t a
, int scale
, float_status
*status
)
3212 FloatParts64 pa
= uint_to_float(a
, scale
, status
);
3213 return float32_round_pack_canonical(&pa
, status
);
3216 float32
uint32_to_float32_scalbn(uint32_t a
, int scale
, float_status
*status
)
3218 return uint64_to_float32_scalbn(a
, scale
, status
);
3221 float32
uint16_to_float32_scalbn(uint16_t a
, int scale
, float_status
*status
)
3223 return uint64_to_float32_scalbn(a
, scale
, status
);
3226 float32
uint64_to_float32(uint64_t a
, float_status
*status
)
3228 return uint64_to_float32_scalbn(a
, 0, status
);
3231 float32
uint32_to_float32(uint32_t a
, float_status
*status
)
3233 return uint64_to_float32_scalbn(a
, 0, status
);
3236 float32
uint16_to_float32(uint16_t a
, float_status
*status
)
3238 return uint64_to_float32_scalbn(a
, 0, status
);
3241 float64
uint64_to_float64_scalbn(uint64_t a
, int scale
, float_status
*status
)
3243 FloatParts64 pa
= uint_to_float(a
, scale
, status
);
3244 return float64_round_pack_canonical(&pa
, status
);
3247 float64
uint32_to_float64_scalbn(uint32_t a
, int scale
, float_status
*status
)
3249 return uint64_to_float64_scalbn(a
, scale
, status
);
3252 float64
uint16_to_float64_scalbn(uint16_t a
, int scale
, float_status
*status
)
3254 return uint64_to_float64_scalbn(a
, scale
, status
);
3257 float64
uint64_to_float64(uint64_t a
, float_status
*status
)
3259 return uint64_to_float64_scalbn(a
, 0, status
);
3262 float64
uint32_to_float64(uint32_t a
, float_status
*status
)
3264 return uint64_to_float64_scalbn(a
, 0, status
);
3267 float64
uint16_to_float64(uint16_t a
, float_status
*status
)
3269 return uint64_to_float64_scalbn(a
, 0, status
);
3273 * Returns the result of converting the unsigned integer `a' to the
3277 bfloat16
uint64_to_bfloat16_scalbn(uint64_t a
, int scale
, float_status
*status
)
3279 FloatParts64 pa
= uint_to_float(a
, scale
, status
);
3280 return bfloat16_round_pack_canonical(&pa
, status
);
3283 bfloat16
uint32_to_bfloat16_scalbn(uint32_t a
, int scale
, float_status
*status
)
3285 return uint64_to_bfloat16_scalbn(a
, scale
, status
);
3288 bfloat16
uint16_to_bfloat16_scalbn(uint16_t a
, int scale
, float_status
*status
)
3290 return uint64_to_bfloat16_scalbn(a
, scale
, status
);
3293 bfloat16
uint64_to_bfloat16(uint64_t a
, float_status
*status
)
3295 return uint64_to_bfloat16_scalbn(a
, 0, status
);
3298 bfloat16
uint32_to_bfloat16(uint32_t a
, float_status
*status
)
3300 return uint64_to_bfloat16_scalbn(a
, 0, status
);
3303 bfloat16
uint16_to_bfloat16(uint16_t a
, float_status
*status
)
3305 return uint64_to_bfloat16_scalbn(a
, 0, status
);
3309 /* min() and max() functions. These can't be implemented as
3310 * 'compare and pick one input' because that would mishandle
3311 * NaNs and +0 vs -0.
3313 * minnum() and maxnum() functions. These are similar to the min()
3314 * and max() functions but if one of the arguments is a QNaN and
3315 * the other is numerical then the numerical argument is returned.
3316 * SNaNs will get quietened before being returned.
3317 * minnum() and maxnum correspond to the IEEE 754-2008 minNum()
3318 * and maxNum() operations. min() and max() are the typical min/max
3319 * semantics provided by many CPUs which predate that specification.
3321 * minnummag() and maxnummag() functions correspond to minNumMag()
3322 * and minNumMag() from the IEEE-754 2008.
3324 static FloatParts64
minmax_floats(FloatParts64 a
, FloatParts64 b
, bool ismin
,
3325 bool ieee
, bool ismag
, float_status
*s
)
3327 if (unlikely(is_nan(a
.cls
) || is_nan(b
.cls
))) {
3329 /* Takes two floating-point values `a' and `b', one of
3330 * which is a NaN, and returns the appropriate NaN
3331 * result. If either `a' or `b' is a signaling NaN,
3332 * the invalid exception is raised.
3334 if (is_snan(a
.cls
) || is_snan(b
.cls
)) {
3335 return *parts_pick_nan(&a
, &b
, s
);
3336 } else if (is_nan(a
.cls
) && !is_nan(b
.cls
)) {
3338 } else if (is_nan(b
.cls
) && !is_nan(a
.cls
)) {
3342 return *parts_pick_nan(&a
, &b
, s
);
3347 case float_class_normal
:
3350 case float_class_inf
:
3353 case float_class_zero
:
3357 g_assert_not_reached();
3361 case float_class_normal
:
3364 case float_class_inf
:
3367 case float_class_zero
:
3371 g_assert_not_reached();
3375 if (ismag
&& (a_exp
!= b_exp
|| a
.frac
!= b
.frac
)) {
3376 bool a_less
= a_exp
< b_exp
;
3377 if (a_exp
== b_exp
) {
3378 a_less
= a
.frac
< b
.frac
;
3380 return a_less
^ ismin
? b
: a
;
3383 if (a
.sign
== b
.sign
) {
3384 bool a_less
= a_exp
< b_exp
;
3385 if (a_exp
== b_exp
) {
3386 a_less
= a
.frac
< b
.frac
;
3388 return a
.sign
^ a_less
^ ismin
? b
: a
;
3390 return a
.sign
^ ismin
? b
: a
;
3395 #define MINMAX(sz, name, ismin, isiee, ismag) \
3396 float ## sz float ## sz ## _ ## name(float ## sz a, float ## sz b, \
3399 FloatParts64 pa, pb, pr; \
3400 float ## sz ## _unpack_canonical(&pa, a, s); \
3401 float ## sz ## _unpack_canonical(&pb, b, s); \
3402 pr = minmax_floats(pa, pb, ismin, isiee, ismag, s); \
3403 return float ## sz ## _round_pack_canonical(&pr, s); \
3406 MINMAX(16, min
, true, false, false)
3407 MINMAX(16, minnum
, true, true, false)
3408 MINMAX(16, minnummag
, true, true, true)
3409 MINMAX(16, max
, false, false, false)
3410 MINMAX(16, maxnum
, false, true, false)
3411 MINMAX(16, maxnummag
, false, true, true)
3413 MINMAX(32, min
, true, false, false)
3414 MINMAX(32, minnum
, true, true, false)
3415 MINMAX(32, minnummag
, true, true, true)
3416 MINMAX(32, max
, false, false, false)
3417 MINMAX(32, maxnum
, false, true, false)
3418 MINMAX(32, maxnummag
, false, true, true)
3420 MINMAX(64, min
, true, false, false)
3421 MINMAX(64, minnum
, true, true, false)
3422 MINMAX(64, minnummag
, true, true, true)
3423 MINMAX(64, max
, false, false, false)
3424 MINMAX(64, maxnum
, false, true, false)
3425 MINMAX(64, maxnummag
, false, true, true)
3429 #define BF16_MINMAX(name, ismin, isiee, ismag) \
3430 bfloat16 bfloat16_ ## name(bfloat16 a, bfloat16 b, float_status *s) \
3432 FloatParts64 pa, pb, pr; \
3433 bfloat16_unpack_canonical(&pa, a, s); \
3434 bfloat16_unpack_canonical(&pb, b, s); \
3435 pr = minmax_floats(pa, pb, ismin, isiee, ismag, s); \
3436 return bfloat16_round_pack_canonical(&pr, s); \
3439 BF16_MINMAX(min
, true, false, false)
3440 BF16_MINMAX(minnum
, true, true, false)
3441 BF16_MINMAX(minnummag
, true, true, true)
3442 BF16_MINMAX(max
, false, false, false)
3443 BF16_MINMAX(maxnum
, false, true, false)
3444 BF16_MINMAX(maxnummag
, false, true, true)
3448 /* Floating point compare */
3449 static FloatRelation
compare_floats(FloatParts64 a
, FloatParts64 b
, bool is_quiet
,
3452 if (is_nan(a
.cls
) || is_nan(b
.cls
)) {
3454 a
.cls
== float_class_snan
||
3455 b
.cls
== float_class_snan
) {
3456 float_raise(float_flag_invalid
, s
);
3458 return float_relation_unordered
;
3461 if (a
.cls
== float_class_zero
) {
3462 if (b
.cls
== float_class_zero
) {
3463 return float_relation_equal
;
3465 return b
.sign
? float_relation_greater
: float_relation_less
;
3466 } else if (b
.cls
== float_class_zero
) {
3467 return a
.sign
? float_relation_less
: float_relation_greater
;
3470 /* The only really important thing about infinity is its sign. If
3471 * both are infinities the sign marks the smallest of the two.
3473 if (a
.cls
== float_class_inf
) {
3474 if ((b
.cls
== float_class_inf
) && (a
.sign
== b
.sign
)) {
3475 return float_relation_equal
;
3477 return a
.sign
? float_relation_less
: float_relation_greater
;
3478 } else if (b
.cls
== float_class_inf
) {
3479 return b
.sign
? float_relation_greater
: float_relation_less
;
3482 if (a
.sign
!= b
.sign
) {
3483 return a
.sign
? float_relation_less
: float_relation_greater
;
3486 if (a
.exp
== b
.exp
) {
3487 if (a
.frac
== b
.frac
) {
3488 return float_relation_equal
;
3491 return a
.frac
> b
.frac
?
3492 float_relation_less
: float_relation_greater
;
3494 return a
.frac
> b
.frac
?
3495 float_relation_greater
: float_relation_less
;
3499 return a
.exp
> b
.exp
? float_relation_less
: float_relation_greater
;
3501 return a
.exp
> b
.exp
? float_relation_greater
: float_relation_less
;
3506 #define COMPARE(name, attr, sz) \
3508 name(float ## sz a, float ## sz b, bool is_quiet, float_status *s) \
3510 FloatParts64 pa, pb; \
3511 float ## sz ## _unpack_canonical(&pa, a, s); \
3512 float ## sz ## _unpack_canonical(&pb, b, s); \
3513 return compare_floats(pa, pb, is_quiet, s); \
3516 COMPARE(soft_f16_compare
, QEMU_FLATTEN
, 16)
3517 COMPARE(soft_f32_compare
, QEMU_SOFTFLOAT_ATTR
, 32)
3518 COMPARE(soft_f64_compare
, QEMU_SOFTFLOAT_ATTR
, 64)
3522 FloatRelation
float16_compare(float16 a
, float16 b
, float_status
*s
)
3524 return soft_f16_compare(a
, b
, false, s
);
3527 FloatRelation
float16_compare_quiet(float16 a
, float16 b
, float_status
*s
)
3529 return soft_f16_compare(a
, b
, true, s
);
3532 static FloatRelation QEMU_FLATTEN
3533 f32_compare(float32 xa
, float32 xb
, bool is_quiet
, float_status
*s
)
3535 union_float32 ua
, ub
;
3540 if (QEMU_NO_HARDFLOAT
) {
3544 float32_input_flush2(&ua
.s
, &ub
.s
, s
);
3545 if (isgreaterequal(ua
.h
, ub
.h
)) {
3546 if (isgreater(ua
.h
, ub
.h
)) {
3547 return float_relation_greater
;
3549 return float_relation_equal
;
3551 if (likely(isless(ua
.h
, ub
.h
))) {
3552 return float_relation_less
;
3554 /* The only condition remaining is unordered.
3555 * Fall through to set flags.
3558 return soft_f32_compare(ua
.s
, ub
.s
, is_quiet
, s
);
3561 FloatRelation
float32_compare(float32 a
, float32 b
, float_status
*s
)
3563 return f32_compare(a
, b
, false, s
);
3566 FloatRelation
float32_compare_quiet(float32 a
, float32 b
, float_status
*s
)
3568 return f32_compare(a
, b
, true, s
);
3571 static FloatRelation QEMU_FLATTEN
3572 f64_compare(float64 xa
, float64 xb
, bool is_quiet
, float_status
*s
)
3574 union_float64 ua
, ub
;
3579 if (QEMU_NO_HARDFLOAT
) {
3583 float64_input_flush2(&ua
.s
, &ub
.s
, s
);
3584 if (isgreaterequal(ua
.h
, ub
.h
)) {
3585 if (isgreater(ua
.h
, ub
.h
)) {
3586 return float_relation_greater
;
3588 return float_relation_equal
;
3590 if (likely(isless(ua
.h
, ub
.h
))) {
3591 return float_relation_less
;
3593 /* The only condition remaining is unordered.
3594 * Fall through to set flags.
3597 return soft_f64_compare(ua
.s
, ub
.s
, is_quiet
, s
);
3600 FloatRelation
float64_compare(float64 a
, float64 b
, float_status
*s
)
3602 return f64_compare(a
, b
, false, s
);
3605 FloatRelation
float64_compare_quiet(float64 a
, float64 b
, float_status
*s
)
3607 return f64_compare(a
, b
, true, s
);
3610 static FloatRelation QEMU_FLATTEN
3611 soft_bf16_compare(bfloat16 a
, bfloat16 b
, bool is_quiet
, float_status
*s
)
3613 FloatParts64 pa
, pb
;
3615 bfloat16_unpack_canonical(&pa
, a
, s
);
3616 bfloat16_unpack_canonical(&pb
, b
, s
);
3617 return compare_floats(pa
, pb
, is_quiet
, s
);
3620 FloatRelation
bfloat16_compare(bfloat16 a
, bfloat16 b
, float_status
*s
)
3622 return soft_bf16_compare(a
, b
, false, s
);
3625 FloatRelation
bfloat16_compare_quiet(bfloat16 a
, bfloat16 b
, float_status
*s
)
3627 return soft_bf16_compare(a
, b
, true, s
);
3630 /* Multiply A by 2 raised to the power N. */
3631 static FloatParts64
scalbn_decomposed(FloatParts64 a
, int n
, float_status
*s
)
3633 if (unlikely(is_nan(a
.cls
))) {
3634 parts_return_nan(&a
, s
);
3636 if (a
.cls
== float_class_normal
) {
3637 /* The largest float type (even though not supported by FloatParts64)
3638 * is float128, which has a 15 bit exponent. Bounding N to 16 bits
3639 * still allows rounding to infinity, without allowing overflow
3640 * within the int32_t that backs FloatParts64.exp.
3642 n
= MIN(MAX(n
, -0x10000), 0x10000);
3648 float16
float16_scalbn(float16 a
, int n
, float_status
*status
)
3650 FloatParts64 pa
, pr
;
3652 float16_unpack_canonical(&pa
, a
, status
);
3653 pr
= scalbn_decomposed(pa
, n
, status
);
3654 return float16_round_pack_canonical(&pr
, status
);
3657 float32
float32_scalbn(float32 a
, int n
, float_status
*status
)
3659 FloatParts64 pa
, pr
;
3661 float32_unpack_canonical(&pa
, a
, status
);
3662 pr
= scalbn_decomposed(pa
, n
, status
);
3663 return float32_round_pack_canonical(&pr
, status
);
3666 float64
float64_scalbn(float64 a
, int n
, float_status
*status
)
3668 FloatParts64 pa
, pr
;
3670 float64_unpack_canonical(&pa
, a
, status
);
3671 pr
= scalbn_decomposed(pa
, n
, status
);
3672 return float64_round_pack_canonical(&pr
, status
);
3675 bfloat16
bfloat16_scalbn(bfloat16 a
, int n
, float_status
*status
)
3677 FloatParts64 pa
, pr
;
3679 bfloat16_unpack_canonical(&pa
, a
, status
);
3680 pr
= scalbn_decomposed(pa
, n
, status
);
3681 return bfloat16_round_pack_canonical(&pr
, status
);
3687 * The old softfloat code did an approximation step before zeroing in
3688 * on the final result. However for simpleness we just compute the
3689 * square root by iterating down from the implicit bit to enough extra
3690 * bits to ensure we get a correctly rounded result.
3692 * This does mean however the calculation is slower than before,
3693 * especially for 64 bit floats.
3696 static FloatParts64
sqrt_float(FloatParts64 a
, float_status
*s
, const FloatFmt
*p
)
3698 uint64_t a_frac
, r_frac
, s_frac
;
3701 if (is_nan(a
.cls
)) {
3702 parts_return_nan(&a
, s
);
3705 if (a
.cls
== float_class_zero
) {
3706 return a
; /* sqrt(+-0) = +-0 */
3709 float_raise(float_flag_invalid
, s
);
3710 parts_default_nan(&a
, s
);
3713 if (a
.cls
== float_class_inf
) {
3714 return a
; /* sqrt(+inf) = +inf */
3717 assert(a
.cls
== float_class_normal
);
3719 /* We need two overflow bits at the top. Adding room for that is a
3720 * right shift. If the exponent is odd, we can discard the low bit
3721 * by multiplying the fraction by 2; that's a left shift. Combine
3722 * those and we shift right by 1 if the exponent is odd, otherwise 2.
3724 a_frac
= a
.frac
>> (2 - (a
.exp
& 1));
3727 /* Bit-by-bit computation of sqrt. */
3731 /* Iterate from implicit bit down to the 3 extra bits to compute a
3732 * properly rounded result. Remember we've inserted two more bits
3733 * at the top, so these positions are two less.
3735 bit
= DECOMPOSED_BINARY_POINT
- 2;
3736 last_bit
= MAX(p
->frac_shift
- 4, 0);
3738 uint64_t q
= 1ULL << bit
;
3739 uint64_t t_frac
= s_frac
+ q
;
3740 if (t_frac
<= a_frac
) {
3741 s_frac
= t_frac
+ q
;
3746 } while (--bit
>= last_bit
);
3748 /* Undo the right shift done above. If there is any remaining
3749 * fraction, the result is inexact. Set the sticky bit.
3751 a
.frac
= (r_frac
<< 2) + (a_frac
!= 0);
3756 float16 QEMU_FLATTEN
float16_sqrt(float16 a
, float_status
*status
)
3758 FloatParts64 pa
, pr
;
3760 float16_unpack_canonical(&pa
, a
, status
);
3761 pr
= sqrt_float(pa
, status
, &float16_params
);
3762 return float16_round_pack_canonical(&pr
, status
);
3765 static float32 QEMU_SOFTFLOAT_ATTR
3766 soft_f32_sqrt(float32 a
, float_status
*status
)
3768 FloatParts64 pa
, pr
;
3770 float32_unpack_canonical(&pa
, a
, status
);
3771 pr
= sqrt_float(pa
, status
, &float32_params
);
3772 return float32_round_pack_canonical(&pr
, status
);
3775 static float64 QEMU_SOFTFLOAT_ATTR
3776 soft_f64_sqrt(float64 a
, float_status
*status
)
3778 FloatParts64 pa
, pr
;
3780 float64_unpack_canonical(&pa
, a
, status
);
3781 pr
= sqrt_float(pa
, status
, &float64_params
);
3782 return float64_round_pack_canonical(&pr
, status
);
3785 float32 QEMU_FLATTEN
float32_sqrt(float32 xa
, float_status
*s
)
3787 union_float32 ua
, ur
;
3790 if (unlikely(!can_use_fpu(s
))) {
3794 float32_input_flush1(&ua
.s
, s
);
3795 if (QEMU_HARDFLOAT_1F32_USE_FP
) {
3796 if (unlikely(!(fpclassify(ua
.h
) == FP_NORMAL
||
3797 fpclassify(ua
.h
) == FP_ZERO
) ||
3801 } else if (unlikely(!float32_is_zero_or_normal(ua
.s
) ||
3802 float32_is_neg(ua
.s
))) {
3809 return soft_f32_sqrt(ua
.s
, s
);
3812 float64 QEMU_FLATTEN
float64_sqrt(float64 xa
, float_status
*s
)
3814 union_float64 ua
, ur
;
3817 if (unlikely(!can_use_fpu(s
))) {
3821 float64_input_flush1(&ua
.s
, s
);
3822 if (QEMU_HARDFLOAT_1F64_USE_FP
) {
3823 if (unlikely(!(fpclassify(ua
.h
) == FP_NORMAL
||
3824 fpclassify(ua
.h
) == FP_ZERO
) ||
3828 } else if (unlikely(!float64_is_zero_or_normal(ua
.s
) ||
3829 float64_is_neg(ua
.s
))) {
3836 return soft_f64_sqrt(ua
.s
, s
);
3839 bfloat16 QEMU_FLATTEN
bfloat16_sqrt(bfloat16 a
, float_status
*status
)
3841 FloatParts64 pa
, pr
;
3843 bfloat16_unpack_canonical(&pa
, a
, status
);
3844 pr
= sqrt_float(pa
, status
, &bfloat16_params
);
3845 return bfloat16_round_pack_canonical(&pr
, status
);
3848 /*----------------------------------------------------------------------------
3849 | The pattern for a default generated NaN.
3850 *----------------------------------------------------------------------------*/
3852 float16
float16_default_nan(float_status
*status
)
3856 parts_default_nan(&p
, status
);
3857 p
.frac
>>= float16_params
.frac_shift
;
3858 return float16_pack_raw(&p
);
3861 float32
float32_default_nan(float_status
*status
)
3865 parts_default_nan(&p
, status
);
3866 p
.frac
>>= float32_params
.frac_shift
;
3867 return float32_pack_raw(&p
);
3870 float64
float64_default_nan(float_status
*status
)
3874 parts_default_nan(&p
, status
);
3875 p
.frac
>>= float64_params
.frac_shift
;
3876 return float64_pack_raw(&p
);
3879 float128
float128_default_nan(float_status
*status
)
3883 parts_default_nan(&p
, status
);
3884 frac_shr(&p
, float128_params
.frac_shift
);
3885 return float128_pack_raw(&p
);
3888 bfloat16
bfloat16_default_nan(float_status
*status
)
3892 parts_default_nan(&p
, status
);
3893 p
.frac
>>= bfloat16_params
.frac_shift
;
3894 return bfloat16_pack_raw(&p
);
3897 /*----------------------------------------------------------------------------
3898 | Returns a quiet NaN from a signalling NaN for the floating point value `a'.
3899 *----------------------------------------------------------------------------*/
3901 float16
float16_silence_nan(float16 a
, float_status
*status
)
3905 float16_unpack_raw(&p
, a
);
3906 p
.frac
<<= float16_params
.frac_shift
;
3907 parts_silence_nan(&p
, status
);
3908 p
.frac
>>= float16_params
.frac_shift
;
3909 return float16_pack_raw(&p
);
3912 float32
float32_silence_nan(float32 a
, float_status
*status
)
3916 float32_unpack_raw(&p
, a
);
3917 p
.frac
<<= float32_params
.frac_shift
;
3918 parts_silence_nan(&p
, status
);
3919 p
.frac
>>= float32_params
.frac_shift
;
3920 return float32_pack_raw(&p
);
3923 float64
float64_silence_nan(float64 a
, float_status
*status
)
3927 float64_unpack_raw(&p
, a
);
3928 p
.frac
<<= float64_params
.frac_shift
;
3929 parts_silence_nan(&p
, status
);
3930 p
.frac
>>= float64_params
.frac_shift
;
3931 return float64_pack_raw(&p
);
3934 bfloat16
bfloat16_silence_nan(bfloat16 a
, float_status
*status
)
3938 bfloat16_unpack_raw(&p
, a
);
3939 p
.frac
<<= bfloat16_params
.frac_shift
;
3940 parts_silence_nan(&p
, status
);
3941 p
.frac
>>= bfloat16_params
.frac_shift
;
3942 return bfloat16_pack_raw(&p
);
3945 float128
float128_silence_nan(float128 a
, float_status
*status
)
3949 float128_unpack_raw(&p
, a
);
3950 frac_shl(&p
, float128_params
.frac_shift
);
3951 parts_silence_nan(&p
, status
);
3952 frac_shr(&p
, float128_params
.frac_shift
);
3953 return float128_pack_raw(&p
);
3956 /*----------------------------------------------------------------------------
3957 | If `a' is denormal and we are in flush-to-zero mode then set the
3958 | input-denormal exception and return zero. Otherwise just return the value.
3959 *----------------------------------------------------------------------------*/
3961 static bool parts_squash_denormal(FloatParts64 p
, float_status
*status
)
3963 if (p
.exp
== 0 && p
.frac
!= 0) {
3964 float_raise(float_flag_input_denormal
, status
);
3971 float16
float16_squash_input_denormal(float16 a
, float_status
*status
)
3973 if (status
->flush_inputs_to_zero
) {
3976 float16_unpack_raw(&p
, a
);
3977 if (parts_squash_denormal(p
, status
)) {
3978 return float16_set_sign(float16_zero
, p
.sign
);
3984 float32
float32_squash_input_denormal(float32 a
, float_status
*status
)
3986 if (status
->flush_inputs_to_zero
) {
3989 float32_unpack_raw(&p
, a
);
3990 if (parts_squash_denormal(p
, status
)) {
3991 return float32_set_sign(float32_zero
, p
.sign
);
3997 float64
float64_squash_input_denormal(float64 a
, float_status
*status
)
3999 if (status
->flush_inputs_to_zero
) {
4002 float64_unpack_raw(&p
, a
);
4003 if (parts_squash_denormal(p
, status
)) {
4004 return float64_set_sign(float64_zero
, p
.sign
);
4010 bfloat16
bfloat16_squash_input_denormal(bfloat16 a
, float_status
*status
)
4012 if (status
->flush_inputs_to_zero
) {
4015 bfloat16_unpack_raw(&p
, a
);
4016 if (parts_squash_denormal(p
, status
)) {
4017 return bfloat16_set_sign(bfloat16_zero
, p
.sign
);
4023 /*----------------------------------------------------------------------------
4024 | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
4025 | and 7, and returns the properly rounded 32-bit integer corresponding to the
4026 | input. If `zSign' is 1, the input is negated before being converted to an
4027 | integer. Bit 63 of `absZ' must be zero. Ordinarily, the fixed-point input
4028 | is simply rounded to an integer, with the inexact exception raised if the
4029 | input cannot be represented exactly as an integer. However, if the fixed-
4030 | point input is too large, the invalid exception is raised and the largest
4031 | positive or negative integer is returned.
4032 *----------------------------------------------------------------------------*/
4034 static int32_t roundAndPackInt32(bool zSign
, uint64_t absZ
,
4035 float_status
*status
)
4037 int8_t roundingMode
;
4038 bool roundNearestEven
;
4039 int8_t roundIncrement
, roundBits
;
4042 roundingMode
= status
->float_rounding_mode
;
4043 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4044 switch (roundingMode
) {
4045 case float_round_nearest_even
:
4046 case float_round_ties_away
:
4047 roundIncrement
= 0x40;
4049 case float_round_to_zero
:
4052 case float_round_up
:
4053 roundIncrement
= zSign
? 0 : 0x7f;
4055 case float_round_down
:
4056 roundIncrement
= zSign
? 0x7f : 0;
4058 case float_round_to_odd
:
4059 roundIncrement
= absZ
& 0x80 ? 0 : 0x7f;
4064 roundBits
= absZ
& 0x7F;
4065 absZ
= ( absZ
+ roundIncrement
)>>7;
4066 if (!(roundBits
^ 0x40) && roundNearestEven
) {
4070 if ( zSign
) z
= - z
;
4071 if ( ( absZ
>>32 ) || ( z
&& ( ( z
< 0 ) ^ zSign
) ) ) {
4072 float_raise(float_flag_invalid
, status
);
4073 return zSign
? INT32_MIN
: INT32_MAX
;
4076 float_raise(float_flag_inexact
, status
);
4082 /*----------------------------------------------------------------------------
4083 | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
4084 | `absZ1', with binary point between bits 63 and 64 (between the input words),
4085 | and returns the properly rounded 64-bit integer corresponding to the input.
4086 | If `zSign' is 1, the input is negated before being converted to an integer.
4087 | Ordinarily, the fixed-point input is simply rounded to an integer, with
4088 | the inexact exception raised if the input cannot be represented exactly as
4089 | an integer. However, if the fixed-point input is too large, the invalid
4090 | exception is raised and the largest positive or negative integer is
4092 *----------------------------------------------------------------------------*/
4094 static int64_t roundAndPackInt64(bool zSign
, uint64_t absZ0
, uint64_t absZ1
,
4095 float_status
*status
)
4097 int8_t roundingMode
;
4098 bool roundNearestEven
, increment
;
4101 roundingMode
= status
->float_rounding_mode
;
4102 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4103 switch (roundingMode
) {
4104 case float_round_nearest_even
:
4105 case float_round_ties_away
:
4106 increment
= ((int64_t) absZ1
< 0);
4108 case float_round_to_zero
:
4111 case float_round_up
:
4112 increment
= !zSign
&& absZ1
;
4114 case float_round_down
:
4115 increment
= zSign
&& absZ1
;
4117 case float_round_to_odd
:
4118 increment
= !(absZ0
& 1) && absZ1
;
4125 if ( absZ0
== 0 ) goto overflow
;
4126 if (!(absZ1
<< 1) && roundNearestEven
) {
4131 if ( zSign
) z
= - z
;
4132 if ( z
&& ( ( z
< 0 ) ^ zSign
) ) {
4134 float_raise(float_flag_invalid
, status
);
4135 return zSign
? INT64_MIN
: INT64_MAX
;
4138 float_raise(float_flag_inexact
, status
);
4144 /*----------------------------------------------------------------------------
4145 | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
4146 | `absZ1', with binary point between bits 63 and 64 (between the input words),
4147 | and returns the properly rounded 64-bit unsigned integer corresponding to the
4148 | input. Ordinarily, the fixed-point input is simply rounded to an integer,
4149 | with the inexact exception raised if the input cannot be represented exactly
4150 | as an integer. However, if the fixed-point input is too large, the invalid
4151 | exception is raised and the largest unsigned integer is returned.
4152 *----------------------------------------------------------------------------*/
4154 static int64_t roundAndPackUint64(bool zSign
, uint64_t absZ0
,
4155 uint64_t absZ1
, float_status
*status
)
4157 int8_t roundingMode
;
4158 bool roundNearestEven
, increment
;
4160 roundingMode
= status
->float_rounding_mode
;
4161 roundNearestEven
= (roundingMode
== float_round_nearest_even
);
4162 switch (roundingMode
) {
4163 case float_round_nearest_even
:
4164 case float_round_ties_away
:
4165 increment
= ((int64_t)absZ1
< 0);
4167 case float_round_to_zero
:
4170 case float_round_up
:
4171 increment
= !zSign
&& absZ1
;
4173 case float_round_down
:
4174 increment
= zSign
&& absZ1
;
4176 case float_round_to_odd
:
4177 increment
= !(absZ0
& 1) && absZ1
;
4185 float_raise(float_flag_invalid
, status
);
4188 if (!(absZ1
<< 1) && roundNearestEven
) {
4193 if (zSign
&& absZ0
) {
4194 float_raise(float_flag_invalid
, status
);
4199 float_raise(float_flag_inexact
, status
);
4204 /*----------------------------------------------------------------------------
4205 | Normalizes the subnormal single-precision floating-point value represented
4206 | by the denormalized significand `aSig'. The normalized exponent and
4207 | significand are stored at the locations pointed to by `zExpPtr' and
4208 | `zSigPtr', respectively.
4209 *----------------------------------------------------------------------------*/
4212 normalizeFloat32Subnormal(uint32_t aSig
, int *zExpPtr
, uint32_t *zSigPtr
)
4216 shiftCount
= clz32(aSig
) - 8;
4217 *zSigPtr
= aSig
<<shiftCount
;
4218 *zExpPtr
= 1 - shiftCount
;
4222 /*----------------------------------------------------------------------------
4223 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4224 | and significand `zSig', and returns the proper single-precision floating-
4225 | point value corresponding to the abstract input. Ordinarily, the abstract
4226 | value is simply rounded and packed into the single-precision format, with
4227 | the inexact exception raised if the abstract input cannot be represented
4228 | exactly. However, if the abstract value is too large, the overflow and
4229 | inexact exceptions are raised and an infinity or maximal finite value is
4230 | returned. If the abstract value is too small, the input value is rounded to
4231 | a subnormal number, and the underflow and inexact exceptions are raised if
4232 | the abstract input cannot be represented exactly as a subnormal single-
4233 | precision floating-point number.
4234 | The input significand `zSig' has its binary point between bits 30
4235 | and 29, which is 7 bits to the left of the usual location. This shifted
4236 | significand must be normalized or smaller. If `zSig' is not normalized,
4237 | `zExp' must be 0; in that case, the result returned is a subnormal number,
4238 | and it must not require rounding. In the usual case that `zSig' is
4239 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
4240 | The handling of underflow and overflow follows the IEC/IEEE Standard for
4241 | Binary Floating-Point Arithmetic.
4242 *----------------------------------------------------------------------------*/
4244 static float32
roundAndPackFloat32(bool zSign
, int zExp
, uint32_t zSig
,
4245 float_status
*status
)
4247 int8_t roundingMode
;
4248 bool roundNearestEven
;
4249 int8_t roundIncrement
, roundBits
;
4252 roundingMode
= status
->float_rounding_mode
;
4253 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4254 switch (roundingMode
) {
4255 case float_round_nearest_even
:
4256 case float_round_ties_away
:
4257 roundIncrement
= 0x40;
4259 case float_round_to_zero
:
4262 case float_round_up
:
4263 roundIncrement
= zSign
? 0 : 0x7f;
4265 case float_round_down
:
4266 roundIncrement
= zSign
? 0x7f : 0;
4268 case float_round_to_odd
:
4269 roundIncrement
= zSig
& 0x80 ? 0 : 0x7f;
4275 roundBits
= zSig
& 0x7F;
4276 if ( 0xFD <= (uint16_t) zExp
) {
4277 if ( ( 0xFD < zExp
)
4278 || ( ( zExp
== 0xFD )
4279 && ( (int32_t) ( zSig
+ roundIncrement
) < 0 ) )
4281 bool overflow_to_inf
= roundingMode
!= float_round_to_odd
&&
4282 roundIncrement
!= 0;
4283 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4284 return packFloat32(zSign
, 0xFF, -!overflow_to_inf
);
4287 if (status
->flush_to_zero
) {
4288 float_raise(float_flag_output_denormal
, status
);
4289 return packFloat32(zSign
, 0, 0);
4291 isTiny
= status
->tininess_before_rounding
4293 || (zSig
+ roundIncrement
< 0x80000000);
4294 shift32RightJamming( zSig
, - zExp
, &zSig
);
4296 roundBits
= zSig
& 0x7F;
4297 if (isTiny
&& roundBits
) {
4298 float_raise(float_flag_underflow
, status
);
4300 if (roundingMode
== float_round_to_odd
) {
4302 * For round-to-odd case, the roundIncrement depends on
4303 * zSig which just changed.
4305 roundIncrement
= zSig
& 0x80 ? 0 : 0x7f;
4310 float_raise(float_flag_inexact
, status
);
4312 zSig
= ( zSig
+ roundIncrement
)>>7;
4313 if (!(roundBits
^ 0x40) && roundNearestEven
) {
4316 if ( zSig
== 0 ) zExp
= 0;
4317 return packFloat32( zSign
, zExp
, zSig
);
4321 /*----------------------------------------------------------------------------
4322 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4323 | and significand `zSig', and returns the proper single-precision floating-
4324 | point value corresponding to the abstract input. This routine is just like
4325 | `roundAndPackFloat32' except that `zSig' does not have to be normalized.
4326 | Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
4327 | floating-point exponent.
4328 *----------------------------------------------------------------------------*/
4331 normalizeRoundAndPackFloat32(bool zSign
, int zExp
, uint32_t zSig
,
4332 float_status
*status
)
4336 shiftCount
= clz32(zSig
) - 1;
4337 return roundAndPackFloat32(zSign
, zExp
- shiftCount
, zSig
<<shiftCount
,
4342 /*----------------------------------------------------------------------------
4343 | Normalizes the subnormal double-precision floating-point value represented
4344 | by the denormalized significand `aSig'. The normalized exponent and
4345 | significand are stored at the locations pointed to by `zExpPtr' and
4346 | `zSigPtr', respectively.
4347 *----------------------------------------------------------------------------*/
4350 normalizeFloat64Subnormal(uint64_t aSig
, int *zExpPtr
, uint64_t *zSigPtr
)
4354 shiftCount
= clz64(aSig
) - 11;
4355 *zSigPtr
= aSig
<<shiftCount
;
4356 *zExpPtr
= 1 - shiftCount
;
4360 /*----------------------------------------------------------------------------
4361 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
4362 | double-precision floating-point value, returning the result. After being
4363 | shifted into the proper positions, the three fields are simply added
4364 | together to form the result. This means that any integer portion of `zSig'
4365 | will be added into the exponent. Since a properly normalized significand
4366 | will have an integer portion equal to 1, the `zExp' input should be 1 less
4367 | than the desired result exponent whenever `zSig' is a complete, normalized
4369 *----------------------------------------------------------------------------*/
4371 static inline float64
packFloat64(bool zSign
, int zExp
, uint64_t zSig
)
4374 return make_float64(
4375 ( ( (uint64_t) zSign
)<<63 ) + ( ( (uint64_t) zExp
)<<52 ) + zSig
);
4379 /*----------------------------------------------------------------------------
4380 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4381 | and significand `zSig', and returns the proper double-precision floating-
4382 | point value corresponding to the abstract input. Ordinarily, the abstract
4383 | value is simply rounded and packed into the double-precision format, with
4384 | the inexact exception raised if the abstract input cannot be represented
4385 | exactly. However, if the abstract value is too large, the overflow and
4386 | inexact exceptions are raised and an infinity or maximal finite value is
4387 | returned. If the abstract value is too small, the input value is rounded to
4388 | a subnormal number, and the underflow and inexact exceptions are raised if
4389 | the abstract input cannot be represented exactly as a subnormal double-
4390 | precision floating-point number.
4391 | The input significand `zSig' has its binary point between bits 62
4392 | and 61, which is 10 bits to the left of the usual location. This shifted
4393 | significand must be normalized or smaller. If `zSig' is not normalized,
4394 | `zExp' must be 0; in that case, the result returned is a subnormal number,
4395 | and it must not require rounding. In the usual case that `zSig' is
4396 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
4397 | The handling of underflow and overflow follows the IEC/IEEE Standard for
4398 | Binary Floating-Point Arithmetic.
4399 *----------------------------------------------------------------------------*/
4401 static float64
roundAndPackFloat64(bool zSign
, int zExp
, uint64_t zSig
,
4402 float_status
*status
)
4404 int8_t roundingMode
;
4405 bool roundNearestEven
;
4406 int roundIncrement
, roundBits
;
4409 roundingMode
= status
->float_rounding_mode
;
4410 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4411 switch (roundingMode
) {
4412 case float_round_nearest_even
:
4413 case float_round_ties_away
:
4414 roundIncrement
= 0x200;
4416 case float_round_to_zero
:
4419 case float_round_up
:
4420 roundIncrement
= zSign
? 0 : 0x3ff;
4422 case float_round_down
:
4423 roundIncrement
= zSign
? 0x3ff : 0;
4425 case float_round_to_odd
:
4426 roundIncrement
= (zSig
& 0x400) ? 0 : 0x3ff;
4431 roundBits
= zSig
& 0x3FF;
4432 if ( 0x7FD <= (uint16_t) zExp
) {
4433 if ( ( 0x7FD < zExp
)
4434 || ( ( zExp
== 0x7FD )
4435 && ( (int64_t) ( zSig
+ roundIncrement
) < 0 ) )
4437 bool overflow_to_inf
= roundingMode
!= float_round_to_odd
&&
4438 roundIncrement
!= 0;
4439 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4440 return packFloat64(zSign
, 0x7FF, -(!overflow_to_inf
));
4443 if (status
->flush_to_zero
) {
4444 float_raise(float_flag_output_denormal
, status
);
4445 return packFloat64(zSign
, 0, 0);
4447 isTiny
= status
->tininess_before_rounding
4449 || (zSig
+ roundIncrement
< UINT64_C(0x8000000000000000));
4450 shift64RightJamming( zSig
, - zExp
, &zSig
);
4452 roundBits
= zSig
& 0x3FF;
4453 if (isTiny
&& roundBits
) {
4454 float_raise(float_flag_underflow
, status
);
4456 if (roundingMode
== float_round_to_odd
) {
4458 * For round-to-odd case, the roundIncrement depends on
4459 * zSig which just changed.
4461 roundIncrement
= (zSig
& 0x400) ? 0 : 0x3ff;
4466 float_raise(float_flag_inexact
, status
);
4468 zSig
= ( zSig
+ roundIncrement
)>>10;
4469 if (!(roundBits
^ 0x200) && roundNearestEven
) {
4472 if ( zSig
== 0 ) zExp
= 0;
4473 return packFloat64( zSign
, zExp
, zSig
);
4477 /*----------------------------------------------------------------------------
4478 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4479 | and significand `zSig', and returns the proper double-precision floating-
4480 | point value corresponding to the abstract input. This routine is just like
4481 | `roundAndPackFloat64' except that `zSig' does not have to be normalized.
4482 | Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
4483 | floating-point exponent.
4484 *----------------------------------------------------------------------------*/
4487 normalizeRoundAndPackFloat64(bool zSign
, int zExp
, uint64_t zSig
,
4488 float_status
*status
)
4492 shiftCount
= clz64(zSig
) - 1;
4493 return roundAndPackFloat64(zSign
, zExp
- shiftCount
, zSig
<<shiftCount
,
4498 /*----------------------------------------------------------------------------
4499 | Normalizes the subnormal extended double-precision floating-point value
4500 | represented by the denormalized significand `aSig'. The normalized exponent
4501 | and significand are stored at the locations pointed to by `zExpPtr' and
4502 | `zSigPtr', respectively.
4503 *----------------------------------------------------------------------------*/
4505 void normalizeFloatx80Subnormal(uint64_t aSig
, int32_t *zExpPtr
,
4510 shiftCount
= clz64(aSig
);
4511 *zSigPtr
= aSig
<<shiftCount
;
4512 *zExpPtr
= 1 - shiftCount
;
4515 /*----------------------------------------------------------------------------
4516 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4517 | and extended significand formed by the concatenation of `zSig0' and `zSig1',
4518 | and returns the proper extended double-precision floating-point value
4519 | corresponding to the abstract input. Ordinarily, the abstract value is
4520 | rounded and packed into the extended double-precision format, with the
4521 | inexact exception raised if the abstract input cannot be represented
4522 | exactly. However, if the abstract value is too large, the overflow and
4523 | inexact exceptions are raised and an infinity or maximal finite value is
4524 | returned. If the abstract value is too small, the input value is rounded to
4525 | a subnormal number, and the underflow and inexact exceptions are raised if
4526 | the abstract input cannot be represented exactly as a subnormal extended
4527 | double-precision floating-point number.
4528 | If `roundingPrecision' is 32 or 64, the result is rounded to the same
4529 | number of bits as single or double precision, respectively. Otherwise, the
4530 | result is rounded to the full precision of the extended double-precision
4532 | The input significand must be normalized or smaller. If the input
4533 | significand is not normalized, `zExp' must be 0; in that case, the result
4534 | returned is a subnormal number, and it must not require rounding. The
4535 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
4536 | Floating-Point Arithmetic.
4537 *----------------------------------------------------------------------------*/
4539 floatx80
roundAndPackFloatx80(int8_t roundingPrecision
, bool zSign
,
4540 int32_t zExp
, uint64_t zSig0
, uint64_t zSig1
,
4541 float_status
*status
)
4543 int8_t roundingMode
;
4544 bool roundNearestEven
, increment
, isTiny
;
4545 int64_t roundIncrement
, roundMask
, roundBits
;
4547 roundingMode
= status
->float_rounding_mode
;
4548 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4549 if ( roundingPrecision
== 80 ) goto precision80
;
4550 if ( roundingPrecision
== 64 ) {
4551 roundIncrement
= UINT64_C(0x0000000000000400);
4552 roundMask
= UINT64_C(0x00000000000007FF);
4554 else if ( roundingPrecision
== 32 ) {
4555 roundIncrement
= UINT64_C(0x0000008000000000);
4556 roundMask
= UINT64_C(0x000000FFFFFFFFFF);
4561 zSig0
|= ( zSig1
!= 0 );
4562 switch (roundingMode
) {
4563 case float_round_nearest_even
:
4564 case float_round_ties_away
:
4566 case float_round_to_zero
:
4569 case float_round_up
:
4570 roundIncrement
= zSign
? 0 : roundMask
;
4572 case float_round_down
:
4573 roundIncrement
= zSign
? roundMask
: 0;
4578 roundBits
= zSig0
& roundMask
;
4579 if ( 0x7FFD <= (uint32_t) ( zExp
- 1 ) ) {
4580 if ( ( 0x7FFE < zExp
)
4581 || ( ( zExp
== 0x7FFE ) && ( zSig0
+ roundIncrement
< zSig0
) )
4586 if (status
->flush_to_zero
) {
4587 float_raise(float_flag_output_denormal
, status
);
4588 return packFloatx80(zSign
, 0, 0);
4590 isTiny
= status
->tininess_before_rounding
4592 || (zSig0
<= zSig0
+ roundIncrement
);
4593 shift64RightJamming( zSig0
, 1 - zExp
, &zSig0
);
4595 roundBits
= zSig0
& roundMask
;
4596 if (isTiny
&& roundBits
) {
4597 float_raise(float_flag_underflow
, status
);
4600 float_raise(float_flag_inexact
, status
);
4602 zSig0
+= roundIncrement
;
4603 if ( (int64_t) zSig0
< 0 ) zExp
= 1;
4604 roundIncrement
= roundMask
+ 1;
4605 if ( roundNearestEven
&& ( roundBits
<<1 == roundIncrement
) ) {
4606 roundMask
|= roundIncrement
;
4608 zSig0
&= ~ roundMask
;
4609 return packFloatx80( zSign
, zExp
, zSig0
);
4613 float_raise(float_flag_inexact
, status
);
4615 zSig0
+= roundIncrement
;
4616 if ( zSig0
< roundIncrement
) {
4618 zSig0
= UINT64_C(0x8000000000000000);
4620 roundIncrement
= roundMask
+ 1;
4621 if ( roundNearestEven
&& ( roundBits
<<1 == roundIncrement
) ) {
4622 roundMask
|= roundIncrement
;
4624 zSig0
&= ~ roundMask
;
4625 if ( zSig0
== 0 ) zExp
= 0;
4626 return packFloatx80( zSign
, zExp
, zSig0
);
4628 switch (roundingMode
) {
4629 case float_round_nearest_even
:
4630 case float_round_ties_away
:
4631 increment
= ((int64_t)zSig1
< 0);
4633 case float_round_to_zero
:
4636 case float_round_up
:
4637 increment
= !zSign
&& zSig1
;
4639 case float_round_down
:
4640 increment
= zSign
&& zSig1
;
4645 if ( 0x7FFD <= (uint32_t) ( zExp
- 1 ) ) {
4646 if ( ( 0x7FFE < zExp
)
4647 || ( ( zExp
== 0x7FFE )
4648 && ( zSig0
== UINT64_C(0xFFFFFFFFFFFFFFFF) )
4654 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4655 if ( ( roundingMode
== float_round_to_zero
)
4656 || ( zSign
&& ( roundingMode
== float_round_up
) )
4657 || ( ! zSign
&& ( roundingMode
== float_round_down
) )
4659 return packFloatx80( zSign
, 0x7FFE, ~ roundMask
);
4661 return packFloatx80(zSign
,
4662 floatx80_infinity_high
,
4663 floatx80_infinity_low
);
4666 isTiny
= status
->tininess_before_rounding
4669 || (zSig0
< UINT64_C(0xFFFFFFFFFFFFFFFF));
4670 shift64ExtraRightJamming( zSig0
, zSig1
, 1 - zExp
, &zSig0
, &zSig1
);
4672 if (isTiny
&& zSig1
) {
4673 float_raise(float_flag_underflow
, status
);
4676 float_raise(float_flag_inexact
, status
);
4678 switch (roundingMode
) {
4679 case float_round_nearest_even
:
4680 case float_round_ties_away
:
4681 increment
= ((int64_t)zSig1
< 0);
4683 case float_round_to_zero
:
4686 case float_round_up
:
4687 increment
= !zSign
&& zSig1
;
4689 case float_round_down
:
4690 increment
= zSign
&& zSig1
;
4697 if (!(zSig1
<< 1) && roundNearestEven
) {
4700 if ( (int64_t) zSig0
< 0 ) zExp
= 1;
4702 return packFloatx80( zSign
, zExp
, zSig0
);
4706 float_raise(float_flag_inexact
, status
);
4712 zSig0
= UINT64_C(0x8000000000000000);
4715 if (!(zSig1
<< 1) && roundNearestEven
) {
4721 if ( zSig0
== 0 ) zExp
= 0;
4723 return packFloatx80( zSign
, zExp
, zSig0
);
4727 /*----------------------------------------------------------------------------
4728 | Takes an abstract floating-point value having sign `zSign', exponent
4729 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
4730 | and returns the proper extended double-precision floating-point value
4731 | corresponding to the abstract input. This routine is just like
4732 | `roundAndPackFloatx80' except that the input significand does not have to be
4734 *----------------------------------------------------------------------------*/
4736 floatx80
normalizeRoundAndPackFloatx80(int8_t roundingPrecision
,
4737 bool zSign
, int32_t zExp
,
4738 uint64_t zSig0
, uint64_t zSig1
,
4739 float_status
*status
)
4748 shiftCount
= clz64(zSig0
);
4749 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
4751 return roundAndPackFloatx80(roundingPrecision
, zSign
, zExp
,
4752 zSig0
, zSig1
, status
);
4756 /*----------------------------------------------------------------------------
4757 | Returns the least-significant 64 fraction bits of the quadruple-precision
4758 | floating-point value `a'.
4759 *----------------------------------------------------------------------------*/
4761 static inline uint64_t extractFloat128Frac1( float128 a
)
4768 /*----------------------------------------------------------------------------
4769 | Returns the most-significant 48 fraction bits of the quadruple-precision
4770 | floating-point value `a'.
4771 *----------------------------------------------------------------------------*/
4773 static inline uint64_t extractFloat128Frac0( float128 a
)
4776 return a
.high
& UINT64_C(0x0000FFFFFFFFFFFF);
4780 /*----------------------------------------------------------------------------
4781 | Returns the exponent bits of the quadruple-precision floating-point value
4783 *----------------------------------------------------------------------------*/
4785 static inline int32_t extractFloat128Exp( float128 a
)
4788 return ( a
.high
>>48 ) & 0x7FFF;
4792 /*----------------------------------------------------------------------------
4793 | Returns the sign bit of the quadruple-precision floating-point value `a'.
4794 *----------------------------------------------------------------------------*/
4796 static inline bool extractFloat128Sign(float128 a
)
4798 return a
.high
>> 63;
4801 /*----------------------------------------------------------------------------
4802 | Normalizes the subnormal quadruple-precision floating-point value
4803 | represented by the denormalized significand formed by the concatenation of
4804 | `aSig0' and `aSig1'. The normalized exponent is stored at the location
4805 | pointed to by `zExpPtr'. The most significant 49 bits of the normalized
4806 | significand are stored at the location pointed to by `zSig0Ptr', and the
4807 | least significant 64 bits of the normalized significand are stored at the
4808 | location pointed to by `zSig1Ptr'.
4809 *----------------------------------------------------------------------------*/
4812 normalizeFloat128Subnormal(
4823 shiftCount
= clz64(aSig1
) - 15;
4824 if ( shiftCount
< 0 ) {
4825 *zSig0Ptr
= aSig1
>>( - shiftCount
);
4826 *zSig1Ptr
= aSig1
<<( shiftCount
& 63 );
4829 *zSig0Ptr
= aSig1
<<shiftCount
;
4832 *zExpPtr
= - shiftCount
- 63;
4835 shiftCount
= clz64(aSig0
) - 15;
4836 shortShift128Left( aSig0
, aSig1
, shiftCount
, zSig0Ptr
, zSig1Ptr
);
4837 *zExpPtr
= 1 - shiftCount
;
4842 /*----------------------------------------------------------------------------
4843 | Packs the sign `zSign', the exponent `zExp', and the significand formed
4844 | by the concatenation of `zSig0' and `zSig1' into a quadruple-precision
4845 | floating-point value, returning the result. After being shifted into the
4846 | proper positions, the three fields `zSign', `zExp', and `zSig0' are simply
4847 | added together to form the most significant 32 bits of the result. This
4848 | means that any integer portion of `zSig0' will be added into the exponent.
4849 | Since a properly normalized significand will have an integer portion equal
4850 | to 1, the `zExp' input should be 1 less than the desired result exponent
4851 | whenever `zSig0' and `zSig1' concatenated form a complete, normalized
4853 *----------------------------------------------------------------------------*/
4855 static inline float128
4856 packFloat128(bool zSign
, int32_t zExp
, uint64_t zSig0
, uint64_t zSig1
)
4861 z
.high
= ((uint64_t)zSign
<< 63) + ((uint64_t)zExp
<< 48) + zSig0
;
4865 /*----------------------------------------------------------------------------
4866 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
4867 | and extended significand formed by the concatenation of `zSig0', `zSig1',
4868 | and `zSig2', and returns the proper quadruple-precision floating-point value
4869 | corresponding to the abstract input. Ordinarily, the abstract value is
4870 | simply rounded and packed into the quadruple-precision format, with the
4871 | inexact exception raised if the abstract input cannot be represented
4872 | exactly. However, if the abstract value is too large, the overflow and
4873 | inexact exceptions are raised and an infinity or maximal finite value is
4874 | returned. If the abstract value is too small, the input value is rounded to
4875 | a subnormal number, and the underflow and inexact exceptions are raised if
4876 | the abstract input cannot be represented exactly as a subnormal quadruple-
4877 | precision floating-point number.
4878 | The input significand must be normalized or smaller. If the input
4879 | significand is not normalized, `zExp' must be 0; in that case, the result
4880 | returned is a subnormal number, and it must not require rounding. In the
4881 | usual case that the input significand is normalized, `zExp' must be 1 less
4882 | than the ``true'' floating-point exponent. The handling of underflow and
4883 | overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4884 *----------------------------------------------------------------------------*/
4886 static float128
roundAndPackFloat128(bool zSign
, int32_t zExp
,
4887 uint64_t zSig0
, uint64_t zSig1
,
4888 uint64_t zSig2
, float_status
*status
)
4890 int8_t roundingMode
;
4891 bool roundNearestEven
, increment
, isTiny
;
4893 roundingMode
= status
->float_rounding_mode
;
4894 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
4895 switch (roundingMode
) {
4896 case float_round_nearest_even
:
4897 case float_round_ties_away
:
4898 increment
= ((int64_t)zSig2
< 0);
4900 case float_round_to_zero
:
4903 case float_round_up
:
4904 increment
= !zSign
&& zSig2
;
4906 case float_round_down
:
4907 increment
= zSign
&& zSig2
;
4909 case float_round_to_odd
:
4910 increment
= !(zSig1
& 0x1) && zSig2
;
4915 if ( 0x7FFD <= (uint32_t) zExp
) {
4916 if ( ( 0x7FFD < zExp
)
4917 || ( ( zExp
== 0x7FFD )
4919 UINT64_C(0x0001FFFFFFFFFFFF),
4920 UINT64_C(0xFFFFFFFFFFFFFFFF),
4927 float_raise(float_flag_overflow
| float_flag_inexact
, status
);
4928 if ( ( roundingMode
== float_round_to_zero
)
4929 || ( zSign
&& ( roundingMode
== float_round_up
) )
4930 || ( ! zSign
&& ( roundingMode
== float_round_down
) )
4931 || (roundingMode
== float_round_to_odd
)
4937 UINT64_C(0x0000FFFFFFFFFFFF),
4938 UINT64_C(0xFFFFFFFFFFFFFFFF)
4941 return packFloat128( zSign
, 0x7FFF, 0, 0 );
4944 if (status
->flush_to_zero
) {
4945 float_raise(float_flag_output_denormal
, status
);
4946 return packFloat128(zSign
, 0, 0, 0);
4948 isTiny
= status
->tininess_before_rounding
4951 || lt128(zSig0
, zSig1
,
4952 UINT64_C(0x0001FFFFFFFFFFFF),
4953 UINT64_C(0xFFFFFFFFFFFFFFFF));
4954 shift128ExtraRightJamming(
4955 zSig0
, zSig1
, zSig2
, - zExp
, &zSig0
, &zSig1
, &zSig2
);
4957 if (isTiny
&& zSig2
) {
4958 float_raise(float_flag_underflow
, status
);
4960 switch (roundingMode
) {
4961 case float_round_nearest_even
:
4962 case float_round_ties_away
:
4963 increment
= ((int64_t)zSig2
< 0);
4965 case float_round_to_zero
:
4968 case float_round_up
:
4969 increment
= !zSign
&& zSig2
;
4971 case float_round_down
:
4972 increment
= zSign
&& zSig2
;
4974 case float_round_to_odd
:
4975 increment
= !(zSig1
& 0x1) && zSig2
;
4983 float_raise(float_flag_inexact
, status
);
4986 add128( zSig0
, zSig1
, 0, 1, &zSig0
, &zSig1
);
4987 if ((zSig2
+ zSig2
== 0) && roundNearestEven
) {
4992 if ( ( zSig0
| zSig1
) == 0 ) zExp
= 0;
4994 return packFloat128( zSign
, zExp
, zSig0
, zSig1
);
4998 /*----------------------------------------------------------------------------
4999 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
5000 | and significand formed by the concatenation of `zSig0' and `zSig1', and
5001 | returns the proper quadruple-precision floating-point value corresponding
5002 | to the abstract input. This routine is just like `roundAndPackFloat128'
5003 | except that the input significand has fewer bits and does not have to be
5004 | normalized. In all cases, `zExp' must be 1 less than the ``true'' floating-
5006 *----------------------------------------------------------------------------*/
5008 static float128
normalizeRoundAndPackFloat128(bool zSign
, int32_t zExp
,
5009 uint64_t zSig0
, uint64_t zSig1
,
5010 float_status
*status
)
5020 shiftCount
= clz64(zSig0
) - 15;
5021 if ( 0 <= shiftCount
) {
5023 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
5026 shift128ExtraRightJamming(
5027 zSig0
, zSig1
, 0, - shiftCount
, &zSig0
, &zSig1
, &zSig2
);
5030 return roundAndPackFloat128(zSign
, zExp
, zSig0
, zSig1
, zSig2
, status
);
5035 /*----------------------------------------------------------------------------
5036 | Returns the result of converting the 32-bit two's complement integer `a'
5037 | to the extended double-precision floating-point format. The conversion
5038 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5040 *----------------------------------------------------------------------------*/
5042 floatx80
int32_to_floatx80(int32_t a
, float_status
*status
)
5049 if ( a
== 0 ) return packFloatx80( 0, 0, 0 );
5051 absA
= zSign
? - a
: a
;
5052 shiftCount
= clz32(absA
) + 32;
5054 return packFloatx80( zSign
, 0x403E - shiftCount
, zSig
<<shiftCount
);
5058 /*----------------------------------------------------------------------------
5059 | Returns the result of converting the 32-bit two's complement integer `a' to
5060 | the quadruple-precision floating-point format. The conversion is performed
5061 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5062 *----------------------------------------------------------------------------*/
5064 float128
int32_to_float128(int32_t a
, float_status
*status
)
5071 if ( a
== 0 ) return packFloat128( 0, 0, 0, 0 );
5073 absA
= zSign
? - a
: a
;
5074 shiftCount
= clz32(absA
) + 17;
5076 return packFloat128( zSign
, 0x402E - shiftCount
, zSig0
<<shiftCount
, 0 );
5080 /*----------------------------------------------------------------------------
5081 | Returns the result of converting the 64-bit two's complement integer `a'
5082 | to the extended double-precision floating-point format. The conversion
5083 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5085 *----------------------------------------------------------------------------*/
5087 floatx80
int64_to_floatx80(int64_t a
, float_status
*status
)
5093 if ( a
== 0 ) return packFloatx80( 0, 0, 0 );
5095 absA
= zSign
? - a
: a
;
5096 shiftCount
= clz64(absA
);
5097 return packFloatx80( zSign
, 0x403E - shiftCount
, absA
<<shiftCount
);
5101 /*----------------------------------------------------------------------------
5102 | Returns the result of converting the 64-bit two's complement integer `a' to
5103 | the quadruple-precision floating-point format. The conversion is performed
5104 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5105 *----------------------------------------------------------------------------*/
5107 float128
int64_to_float128(int64_t a
, float_status
*status
)
5113 uint64_t zSig0
, zSig1
;
5115 if ( a
== 0 ) return packFloat128( 0, 0, 0, 0 );
5117 absA
= zSign
? - a
: a
;
5118 shiftCount
= clz64(absA
) + 49;
5119 zExp
= 0x406E - shiftCount
;
5120 if ( 64 <= shiftCount
) {
5129 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
5130 return packFloat128( zSign
, zExp
, zSig0
, zSig1
);
5134 /*----------------------------------------------------------------------------
5135 | Returns the result of converting the 64-bit unsigned integer `a'
5136 | to the quadruple-precision floating-point format. The conversion is performed
5137 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5138 *----------------------------------------------------------------------------*/
5140 float128
uint64_to_float128(uint64_t a
, float_status
*status
)
5143 return float128_zero
;
5145 return normalizeRoundAndPackFloat128(0, 0x406E, 0, a
, status
);
5148 /*----------------------------------------------------------------------------
5149 | Returns the result of converting the single-precision floating-point value
5150 | `a' to the extended double-precision floating-point format. The conversion
5151 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5153 *----------------------------------------------------------------------------*/
5155 floatx80
float32_to_floatx80(float32 a
, float_status
*status
)
5161 a
= float32_squash_input_denormal(a
, status
);
5162 aSig
= extractFloat32Frac( a
);
5163 aExp
= extractFloat32Exp( a
);
5164 aSign
= extractFloat32Sign( a
);
5165 if ( aExp
== 0xFF ) {
5167 floatx80 res
= commonNaNToFloatx80(float32ToCommonNaN(a
, status
),
5169 return floatx80_silence_nan(res
, status
);
5171 return packFloatx80(aSign
,
5172 floatx80_infinity_high
,
5173 floatx80_infinity_low
);
5176 if ( aSig
== 0 ) return packFloatx80( aSign
, 0, 0 );
5177 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
5180 return packFloatx80( aSign
, aExp
+ 0x3F80, ( (uint64_t) aSig
)<<40 );
5184 /*----------------------------------------------------------------------------
5185 | Returns the remainder of the single-precision floating-point value `a'
5186 | with respect to the corresponding value `b'. The operation is performed
5187 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5188 *----------------------------------------------------------------------------*/
5190 float32
float32_rem(float32 a
, float32 b
, float_status
*status
)
5193 int aExp
, bExp
, expDiff
;
5194 uint32_t aSig
, bSig
;
5196 uint64_t aSig64
, bSig64
, q64
;
5197 uint32_t alternateASig
;
5199 a
= float32_squash_input_denormal(a
, status
);
5200 b
= float32_squash_input_denormal(b
, status
);
5202 aSig
= extractFloat32Frac( a
);
5203 aExp
= extractFloat32Exp( a
);
5204 aSign
= extractFloat32Sign( a
);
5205 bSig
= extractFloat32Frac( b
);
5206 bExp
= extractFloat32Exp( b
);
5207 if ( aExp
== 0xFF ) {
5208 if ( aSig
|| ( ( bExp
== 0xFF ) && bSig
) ) {
5209 return propagateFloat32NaN(a
, b
, status
);
5211 float_raise(float_flag_invalid
, status
);
5212 return float32_default_nan(status
);
5214 if ( bExp
== 0xFF ) {
5216 return propagateFloat32NaN(a
, b
, status
);
5222 float_raise(float_flag_invalid
, status
);
5223 return float32_default_nan(status
);
5225 normalizeFloat32Subnormal( bSig
, &bExp
, &bSig
);
5228 if ( aSig
== 0 ) return a
;
5229 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
5231 expDiff
= aExp
- bExp
;
5234 if ( expDiff
< 32 ) {
5237 if ( expDiff
< 0 ) {
5238 if ( expDiff
< -1 ) return a
;
5241 q
= ( bSig
<= aSig
);
5242 if ( q
) aSig
-= bSig
;
5243 if ( 0 < expDiff
) {
5244 q
= ( ( (uint64_t) aSig
)<<32 ) / bSig
;
5247 aSig
= ( ( aSig
>>1 )<<( expDiff
- 1 ) ) - bSig
* q
;
5255 if ( bSig
<= aSig
) aSig
-= bSig
;
5256 aSig64
= ( (uint64_t) aSig
)<<40;
5257 bSig64
= ( (uint64_t) bSig
)<<40;
5259 while ( 0 < expDiff
) {
5260 q64
= estimateDiv128To64( aSig64
, 0, bSig64
);
5261 q64
= ( 2 < q64
) ? q64
- 2 : 0;
5262 aSig64
= - ( ( bSig
* q64
)<<38 );
5266 q64
= estimateDiv128To64( aSig64
, 0, bSig64
);
5267 q64
= ( 2 < q64
) ? q64
- 2 : 0;
5268 q
= q64
>>( 64 - expDiff
);
5270 aSig
= ( ( aSig64
>>33 )<<( expDiff
- 1 ) ) - bSig
* q
;
5273 alternateASig
= aSig
;
5276 } while ( 0 <= (int32_t) aSig
);
5277 sigMean
= aSig
+ alternateASig
;
5278 if ( ( sigMean
< 0 ) || ( ( sigMean
== 0 ) && ( q
& 1 ) ) ) {
5279 aSig
= alternateASig
;
5281 zSign
= ( (int32_t) aSig
< 0 );
5282 if ( zSign
) aSig
= - aSig
;
5283 return normalizeRoundAndPackFloat32(aSign
^ zSign
, bExp
, aSig
, status
);
5288 /*----------------------------------------------------------------------------
5289 | Returns the binary exponential of the single-precision floating-point value
5290 | `a'. The operation is performed according to the IEC/IEEE Standard for
5291 | Binary Floating-Point Arithmetic.
5293 | Uses the following identities:
5295 | 1. -------------------------------------------------------------------------
5299 | 2. -------------------------------------------------------------------------
5302 | e = 1 + --- + --- + --- + --- + --- + ... + --- + ...
5304 *----------------------------------------------------------------------------*/
5306 static const float64 float32_exp2_coefficients
[15] =
5308 const_float64( 0x3ff0000000000000ll
), /* 1 */
5309 const_float64( 0x3fe0000000000000ll
), /* 2 */
5310 const_float64( 0x3fc5555555555555ll
), /* 3 */
5311 const_float64( 0x3fa5555555555555ll
), /* 4 */
5312 const_float64( 0x3f81111111111111ll
), /* 5 */
5313 const_float64( 0x3f56c16c16c16c17ll
), /* 6 */
5314 const_float64( 0x3f2a01a01a01a01all
), /* 7 */
5315 const_float64( 0x3efa01a01a01a01all
), /* 8 */
5316 const_float64( 0x3ec71de3a556c734ll
), /* 9 */
5317 const_float64( 0x3e927e4fb7789f5cll
), /* 10 */
5318 const_float64( 0x3e5ae64567f544e4ll
), /* 11 */
5319 const_float64( 0x3e21eed8eff8d898ll
), /* 12 */
5320 const_float64( 0x3de6124613a86d09ll
), /* 13 */
5321 const_float64( 0x3da93974a8c07c9dll
), /* 14 */
5322 const_float64( 0x3d6ae7f3e733b81fll
), /* 15 */
5325 float32
float32_exp2(float32 a
, float_status
*status
)
5332 a
= float32_squash_input_denormal(a
, status
);
5334 aSig
= extractFloat32Frac( a
);
5335 aExp
= extractFloat32Exp( a
);
5336 aSign
= extractFloat32Sign( a
);
5338 if ( aExp
== 0xFF) {
5340 return propagateFloat32NaN(a
, float32_zero
, status
);
5342 return (aSign
) ? float32_zero
: a
;
5345 if (aSig
== 0) return float32_one
;
5348 float_raise(float_flag_inexact
, status
);
5350 /* ******************************* */
5351 /* using float64 for approximation */
5352 /* ******************************* */
5353 x
= float32_to_float64(a
, status
);
5354 x
= float64_mul(x
, float64_ln2
, status
);
5358 for (i
= 0 ; i
< 15 ; i
++) {
5361 f
= float64_mul(xn
, float32_exp2_coefficients
[i
], status
);
5362 r
= float64_add(r
, f
, status
);
5364 xn
= float64_mul(xn
, x
, status
);
5367 return float64_to_float32(r
, status
);
5370 /*----------------------------------------------------------------------------
5371 | Returns the binary log of the single-precision floating-point value `a'.
5372 | The operation is performed according to the IEC/IEEE Standard for Binary
5373 | Floating-Point Arithmetic.
5374 *----------------------------------------------------------------------------*/
5375 float32
float32_log2(float32 a
, float_status
*status
)
5379 uint32_t aSig
, zSig
, i
;
5381 a
= float32_squash_input_denormal(a
, status
);
5382 aSig
= extractFloat32Frac( a
);
5383 aExp
= extractFloat32Exp( a
);
5384 aSign
= extractFloat32Sign( a
);
5387 if ( aSig
== 0 ) return packFloat32( 1, 0xFF, 0 );
5388 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
5391 float_raise(float_flag_invalid
, status
);
5392 return float32_default_nan(status
);
5394 if ( aExp
== 0xFF ) {
5396 return propagateFloat32NaN(a
, float32_zero
, status
);
5406 for (i
= 1 << 22; i
> 0; i
>>= 1) {
5407 aSig
= ( (uint64_t)aSig
* aSig
) >> 23;
5408 if ( aSig
& 0x01000000 ) {
5417 return normalizeRoundAndPackFloat32(zSign
, 0x85, zSig
, status
);
5420 /*----------------------------------------------------------------------------
5421 | Returns the result of converting the double-precision floating-point value
5422 | `a' to the extended double-precision floating-point format. The conversion
5423 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5425 *----------------------------------------------------------------------------*/
5427 floatx80
float64_to_floatx80(float64 a
, float_status
*status
)
5433 a
= float64_squash_input_denormal(a
, status
);
5434 aSig
= extractFloat64Frac( a
);
5435 aExp
= extractFloat64Exp( a
);
5436 aSign
= extractFloat64Sign( a
);
5437 if ( aExp
== 0x7FF ) {
5439 floatx80 res
= commonNaNToFloatx80(float64ToCommonNaN(a
, status
),
5441 return floatx80_silence_nan(res
, status
);
5443 return packFloatx80(aSign
,
5444 floatx80_infinity_high
,
5445 floatx80_infinity_low
);
5448 if ( aSig
== 0 ) return packFloatx80( aSign
, 0, 0 );
5449 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
5453 aSign
, aExp
+ 0x3C00, (aSig
| UINT64_C(0x0010000000000000)) << 11);
5457 /*----------------------------------------------------------------------------
5458 | Returns the remainder of the double-precision floating-point value `a'
5459 | with respect to the corresponding value `b'. The operation is performed
5460 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5461 *----------------------------------------------------------------------------*/
5463 float64
float64_rem(float64 a
, float64 b
, float_status
*status
)
5466 int aExp
, bExp
, expDiff
;
5467 uint64_t aSig
, bSig
;
5468 uint64_t q
, alternateASig
;
5471 a
= float64_squash_input_denormal(a
, status
);
5472 b
= float64_squash_input_denormal(b
, status
);
5473 aSig
= extractFloat64Frac( a
);
5474 aExp
= extractFloat64Exp( a
);
5475 aSign
= extractFloat64Sign( a
);
5476 bSig
= extractFloat64Frac( b
);
5477 bExp
= extractFloat64Exp( b
);
5478 if ( aExp
== 0x7FF ) {
5479 if ( aSig
|| ( ( bExp
== 0x7FF ) && bSig
) ) {
5480 return propagateFloat64NaN(a
, b
, status
);
5482 float_raise(float_flag_invalid
, status
);
5483 return float64_default_nan(status
);
5485 if ( bExp
== 0x7FF ) {
5487 return propagateFloat64NaN(a
, b
, status
);
5493 float_raise(float_flag_invalid
, status
);
5494 return float64_default_nan(status
);
5496 normalizeFloat64Subnormal( bSig
, &bExp
, &bSig
);
5499 if ( aSig
== 0 ) return a
;
5500 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
5502 expDiff
= aExp
- bExp
;
5503 aSig
= (aSig
| UINT64_C(0x0010000000000000)) << 11;
5504 bSig
= (bSig
| UINT64_C(0x0010000000000000)) << 11;
5505 if ( expDiff
< 0 ) {
5506 if ( expDiff
< -1 ) return a
;
5509 q
= ( bSig
<= aSig
);
5510 if ( q
) aSig
-= bSig
;
5512 while ( 0 < expDiff
) {
5513 q
= estimateDiv128To64( aSig
, 0, bSig
);
5514 q
= ( 2 < q
) ? q
- 2 : 0;
5515 aSig
= - ( ( bSig
>>2 ) * q
);
5519 if ( 0 < expDiff
) {
5520 q
= estimateDiv128To64( aSig
, 0, bSig
);
5521 q
= ( 2 < q
) ? q
- 2 : 0;
5524 aSig
= ( ( aSig
>>1 )<<( expDiff
- 1 ) ) - bSig
* q
;
5531 alternateASig
= aSig
;
5534 } while ( 0 <= (int64_t) aSig
);
5535 sigMean
= aSig
+ alternateASig
;
5536 if ( ( sigMean
< 0 ) || ( ( sigMean
== 0 ) && ( q
& 1 ) ) ) {
5537 aSig
= alternateASig
;
5539 zSign
= ( (int64_t) aSig
< 0 );
5540 if ( zSign
) aSig
= - aSig
;
5541 return normalizeRoundAndPackFloat64(aSign
^ zSign
, bExp
, aSig
, status
);
5545 /*----------------------------------------------------------------------------
5546 | Returns the binary log of the double-precision floating-point value `a'.
5547 | The operation is performed according to the IEC/IEEE Standard for Binary
5548 | Floating-Point Arithmetic.
5549 *----------------------------------------------------------------------------*/
5550 float64
float64_log2(float64 a
, float_status
*status
)
5554 uint64_t aSig
, aSig0
, aSig1
, zSig
, i
;
5555 a
= float64_squash_input_denormal(a
, status
);
5557 aSig
= extractFloat64Frac( a
);
5558 aExp
= extractFloat64Exp( a
);
5559 aSign
= extractFloat64Sign( a
);
5562 if ( aSig
== 0 ) return packFloat64( 1, 0x7FF, 0 );
5563 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
5566 float_raise(float_flag_invalid
, status
);
5567 return float64_default_nan(status
);
5569 if ( aExp
== 0x7FF ) {
5571 return propagateFloat64NaN(a
, float64_zero
, status
);
5577 aSig
|= UINT64_C(0x0010000000000000);
5579 zSig
= (uint64_t)aExp
<< 52;
5580 for (i
= 1LL << 51; i
> 0; i
>>= 1) {
5581 mul64To128( aSig
, aSig
, &aSig0
, &aSig1
);
5582 aSig
= ( aSig0
<< 12 ) | ( aSig1
>> 52 );
5583 if ( aSig
& UINT64_C(0x0020000000000000) ) {
5591 return normalizeRoundAndPackFloat64(zSign
, 0x408, zSig
, status
);
5594 /*----------------------------------------------------------------------------
5595 | Returns the result of converting the extended double-precision floating-
5596 | point value `a' to the 32-bit two's complement integer format. The
5597 | conversion is performed according to the IEC/IEEE Standard for Binary
5598 | Floating-Point Arithmetic---which means in particular that the conversion
5599 | is rounded according to the current rounding mode. If `a' is a NaN, the
5600 | largest positive integer is returned. Otherwise, if the conversion
5601 | overflows, the largest integer with the same sign as `a' is returned.
5602 *----------------------------------------------------------------------------*/
5604 int32_t floatx80_to_int32(floatx80 a
, float_status
*status
)
5607 int32_t aExp
, shiftCount
;
5610 if (floatx80_invalid_encoding(a
)) {
5611 float_raise(float_flag_invalid
, status
);
5614 aSig
= extractFloatx80Frac( a
);
5615 aExp
= extractFloatx80Exp( a
);
5616 aSign
= extractFloatx80Sign( a
);
5617 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) aSign
= 0;
5618 shiftCount
= 0x4037 - aExp
;
5619 if ( shiftCount
<= 0 ) shiftCount
= 1;
5620 shift64RightJamming( aSig
, shiftCount
, &aSig
);
5621 return roundAndPackInt32(aSign
, aSig
, status
);
5625 /*----------------------------------------------------------------------------
5626 | Returns the result of converting the extended double-precision floating-
5627 | point value `a' to the 32-bit two's complement integer format. The
5628 | conversion is performed according to the IEC/IEEE Standard for Binary
5629 | Floating-Point Arithmetic, except that the conversion is always rounded
5630 | toward zero. If `a' is a NaN, the largest positive integer is returned.
5631 | Otherwise, if the conversion overflows, the largest integer with the same
5632 | sign as `a' is returned.
5633 *----------------------------------------------------------------------------*/
5635 int32_t floatx80_to_int32_round_to_zero(floatx80 a
, float_status
*status
)
5638 int32_t aExp
, shiftCount
;
5639 uint64_t aSig
, savedASig
;
5642 if (floatx80_invalid_encoding(a
)) {
5643 float_raise(float_flag_invalid
, status
);
5646 aSig
= extractFloatx80Frac( a
);
5647 aExp
= extractFloatx80Exp( a
);
5648 aSign
= extractFloatx80Sign( a
);
5649 if ( 0x401E < aExp
) {
5650 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) aSign
= 0;
5653 else if ( aExp
< 0x3FFF ) {
5655 float_raise(float_flag_inexact
, status
);
5659 shiftCount
= 0x403E - aExp
;
5661 aSig
>>= shiftCount
;
5663 if ( aSign
) z
= - z
;
5664 if ( ( z
< 0 ) ^ aSign
) {
5666 float_raise(float_flag_invalid
, status
);
5667 return aSign
? (int32_t) 0x80000000 : 0x7FFFFFFF;
5669 if ( ( aSig
<<shiftCount
) != savedASig
) {
5670 float_raise(float_flag_inexact
, status
);
5676 /*----------------------------------------------------------------------------
5677 | Returns the result of converting the extended double-precision floating-
5678 | point value `a' to the 64-bit two's complement integer format. The
5679 | conversion is performed according to the IEC/IEEE Standard for Binary
5680 | Floating-Point Arithmetic---which means in particular that the conversion
5681 | is rounded according to the current rounding mode. If `a' is a NaN,
5682 | the largest positive integer is returned. Otherwise, if the conversion
5683 | overflows, the largest integer with the same sign as `a' is returned.
5684 *----------------------------------------------------------------------------*/
5686 int64_t floatx80_to_int64(floatx80 a
, float_status
*status
)
5689 int32_t aExp
, shiftCount
;
5690 uint64_t aSig
, aSigExtra
;
5692 if (floatx80_invalid_encoding(a
)) {
5693 float_raise(float_flag_invalid
, status
);
5696 aSig
= extractFloatx80Frac( a
);
5697 aExp
= extractFloatx80Exp( a
);
5698 aSign
= extractFloatx80Sign( a
);
5699 shiftCount
= 0x403E - aExp
;
5700 if ( shiftCount
<= 0 ) {
5702 float_raise(float_flag_invalid
, status
);
5703 if (!aSign
|| floatx80_is_any_nan(a
)) {
5711 shift64ExtraRightJamming( aSig
, 0, shiftCount
, &aSig
, &aSigExtra
);
5713 return roundAndPackInt64(aSign
, aSig
, aSigExtra
, status
);
5717 /*----------------------------------------------------------------------------
5718 | Returns the result of converting the extended double-precision floating-
5719 | point value `a' to the 64-bit two's complement integer format. The
5720 | conversion is performed according to the IEC/IEEE Standard for Binary
5721 | Floating-Point Arithmetic, except that the conversion is always rounded
5722 | toward zero. If `a' is a NaN, the largest positive integer is returned.
5723 | Otherwise, if the conversion overflows, the largest integer with the same
5724 | sign as `a' is returned.
5725 *----------------------------------------------------------------------------*/
5727 int64_t floatx80_to_int64_round_to_zero(floatx80 a
, float_status
*status
)
5730 int32_t aExp
, shiftCount
;
5734 if (floatx80_invalid_encoding(a
)) {
5735 float_raise(float_flag_invalid
, status
);
5738 aSig
= extractFloatx80Frac( a
);
5739 aExp
= extractFloatx80Exp( a
);
5740 aSign
= extractFloatx80Sign( a
);
5741 shiftCount
= aExp
- 0x403E;
5742 if ( 0 <= shiftCount
) {
5743 aSig
&= UINT64_C(0x7FFFFFFFFFFFFFFF);
5744 if ( ( a
.high
!= 0xC03E ) || aSig
) {
5745 float_raise(float_flag_invalid
, status
);
5746 if ( ! aSign
|| ( ( aExp
== 0x7FFF ) && aSig
) ) {
5752 else if ( aExp
< 0x3FFF ) {
5754 float_raise(float_flag_inexact
, status
);
5758 z
= aSig
>>( - shiftCount
);
5759 if ( (uint64_t) ( aSig
<<( shiftCount
& 63 ) ) ) {
5760 float_raise(float_flag_inexact
, status
);
5762 if ( aSign
) z
= - z
;
5767 /*----------------------------------------------------------------------------
5768 | Returns the result of converting the extended double-precision floating-
5769 | point value `a' to the single-precision floating-point format. The
5770 | conversion is performed according to the IEC/IEEE Standard for Binary
5771 | Floating-Point Arithmetic.
5772 *----------------------------------------------------------------------------*/
5774 float32
floatx80_to_float32(floatx80 a
, float_status
*status
)
5780 if (floatx80_invalid_encoding(a
)) {
5781 float_raise(float_flag_invalid
, status
);
5782 return float32_default_nan(status
);
5784 aSig
= extractFloatx80Frac( a
);
5785 aExp
= extractFloatx80Exp( a
);
5786 aSign
= extractFloatx80Sign( a
);
5787 if ( aExp
== 0x7FFF ) {
5788 if ( (uint64_t) ( aSig
<<1 ) ) {
5789 float32 res
= commonNaNToFloat32(floatx80ToCommonNaN(a
, status
),
5791 return float32_silence_nan(res
, status
);
5793 return packFloat32( aSign
, 0xFF, 0 );
5795 shift64RightJamming( aSig
, 33, &aSig
);
5796 if ( aExp
|| aSig
) aExp
-= 0x3F81;
5797 return roundAndPackFloat32(aSign
, aExp
, aSig
, status
);
5801 /*----------------------------------------------------------------------------
5802 | Returns the result of converting the extended double-precision floating-
5803 | point value `a' to the double-precision floating-point format. The
5804 | conversion is performed according to the IEC/IEEE Standard for Binary
5805 | Floating-Point Arithmetic.
5806 *----------------------------------------------------------------------------*/
5808 float64
floatx80_to_float64(floatx80 a
, float_status
*status
)
5812 uint64_t aSig
, zSig
;
5814 if (floatx80_invalid_encoding(a
)) {
5815 float_raise(float_flag_invalid
, status
);
5816 return float64_default_nan(status
);
5818 aSig
= extractFloatx80Frac( a
);
5819 aExp
= extractFloatx80Exp( a
);
5820 aSign
= extractFloatx80Sign( a
);
5821 if ( aExp
== 0x7FFF ) {
5822 if ( (uint64_t) ( aSig
<<1 ) ) {
5823 float64 res
= commonNaNToFloat64(floatx80ToCommonNaN(a
, status
),
5825 return float64_silence_nan(res
, status
);
5827 return packFloat64( aSign
, 0x7FF, 0 );
5829 shift64RightJamming( aSig
, 1, &zSig
);
5830 if ( aExp
|| aSig
) aExp
-= 0x3C01;
5831 return roundAndPackFloat64(aSign
, aExp
, zSig
, status
);
5835 /*----------------------------------------------------------------------------
5836 | Returns the result of converting the extended double-precision floating-
5837 | point value `a' to the quadruple-precision floating-point format. The
5838 | conversion is performed according to the IEC/IEEE Standard for Binary
5839 | Floating-Point Arithmetic.
5840 *----------------------------------------------------------------------------*/
5842 float128
floatx80_to_float128(floatx80 a
, float_status
*status
)
5846 uint64_t aSig
, zSig0
, zSig1
;
5848 if (floatx80_invalid_encoding(a
)) {
5849 float_raise(float_flag_invalid
, status
);
5850 return float128_default_nan(status
);
5852 aSig
= extractFloatx80Frac( a
);
5853 aExp
= extractFloatx80Exp( a
);
5854 aSign
= extractFloatx80Sign( a
);
5855 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) {
5856 float128 res
= commonNaNToFloat128(floatx80ToCommonNaN(a
, status
),
5858 return float128_silence_nan(res
, status
);
5860 shift128Right( aSig
<<1, 0, 16, &zSig0
, &zSig1
);
5861 return packFloat128( aSign
, aExp
, zSig0
, zSig1
);
5865 /*----------------------------------------------------------------------------
5866 | Rounds the extended double-precision floating-point value `a'
5867 | to the precision provided by floatx80_rounding_precision and returns the
5868 | result as an extended double-precision floating-point value.
5869 | The operation is performed according to the IEC/IEEE Standard for Binary
5870 | Floating-Point Arithmetic.
5871 *----------------------------------------------------------------------------*/
5873 floatx80
floatx80_round(floatx80 a
, float_status
*status
)
5875 return roundAndPackFloatx80(status
->floatx80_rounding_precision
,
5876 extractFloatx80Sign(a
),
5877 extractFloatx80Exp(a
),
5878 extractFloatx80Frac(a
), 0, status
);
5881 /*----------------------------------------------------------------------------
5882 | Rounds the extended double-precision floating-point value `a' to an integer,
5883 | and returns the result as an extended quadruple-precision floating-point
5884 | value. The operation is performed according to the IEC/IEEE Standard for
5885 | Binary Floating-Point Arithmetic.
5886 *----------------------------------------------------------------------------*/
5888 floatx80
floatx80_round_to_int(floatx80 a
, float_status
*status
)
5892 uint64_t lastBitMask
, roundBitsMask
;
5895 if (floatx80_invalid_encoding(a
)) {
5896 float_raise(float_flag_invalid
, status
);
5897 return floatx80_default_nan(status
);
5899 aExp
= extractFloatx80Exp( a
);
5900 if ( 0x403E <= aExp
) {
5901 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) ) {
5902 return propagateFloatx80NaN(a
, a
, status
);
5906 if ( aExp
< 0x3FFF ) {
5908 && ( (uint64_t) ( extractFloatx80Frac( a
) ) == 0 ) ) {
5911 float_raise(float_flag_inexact
, status
);
5912 aSign
= extractFloatx80Sign( a
);
5913 switch (status
->float_rounding_mode
) {
5914 case float_round_nearest_even
:
5915 if ( ( aExp
== 0x3FFE ) && (uint64_t) ( extractFloatx80Frac( a
)<<1 )
5918 packFloatx80( aSign
, 0x3FFF, UINT64_C(0x8000000000000000));
5921 case float_round_ties_away
:
5922 if (aExp
== 0x3FFE) {
5923 return packFloatx80(aSign
, 0x3FFF, UINT64_C(0x8000000000000000));
5926 case float_round_down
:
5929 packFloatx80( 1, 0x3FFF, UINT64_C(0x8000000000000000))
5930 : packFloatx80( 0, 0, 0 );
5931 case float_round_up
:
5933 aSign
? packFloatx80( 1, 0, 0 )
5934 : packFloatx80( 0, 0x3FFF, UINT64_C(0x8000000000000000));
5936 case float_round_to_zero
:
5939 g_assert_not_reached();
5941 return packFloatx80( aSign
, 0, 0 );
5944 lastBitMask
<<= 0x403E - aExp
;
5945 roundBitsMask
= lastBitMask
- 1;
5947 switch (status
->float_rounding_mode
) {
5948 case float_round_nearest_even
:
5949 z
.low
+= lastBitMask
>>1;
5950 if ((z
.low
& roundBitsMask
) == 0) {
5951 z
.low
&= ~lastBitMask
;
5954 case float_round_ties_away
:
5955 z
.low
+= lastBitMask
>> 1;
5957 case float_round_to_zero
:
5959 case float_round_up
:
5960 if (!extractFloatx80Sign(z
)) {
5961 z
.low
+= roundBitsMask
;
5964 case float_round_down
:
5965 if (extractFloatx80Sign(z
)) {
5966 z
.low
+= roundBitsMask
;
5972 z
.low
&= ~ roundBitsMask
;
5975 z
.low
= UINT64_C(0x8000000000000000);
5977 if (z
.low
!= a
.low
) {
5978 float_raise(float_flag_inexact
, status
);
5984 /*----------------------------------------------------------------------------
5985 | Returns the result of adding the absolute values of the extended double-
5986 | precision floating-point values `a' and `b'. If `zSign' is 1, the sum is
5987 | negated before being returned. `zSign' is ignored if the result is a NaN.
5988 | The addition is performed according to the IEC/IEEE Standard for Binary
5989 | Floating-Point Arithmetic.
5990 *----------------------------------------------------------------------------*/
5992 static floatx80
addFloatx80Sigs(floatx80 a
, floatx80 b
, bool zSign
,
5993 float_status
*status
)
5995 int32_t aExp
, bExp
, zExp
;
5996 uint64_t aSig
, bSig
, zSig0
, zSig1
;
5999 aSig
= extractFloatx80Frac( a
);
6000 aExp
= extractFloatx80Exp( a
);
6001 bSig
= extractFloatx80Frac( b
);
6002 bExp
= extractFloatx80Exp( b
);
6003 expDiff
= aExp
- bExp
;
6004 if ( 0 < expDiff
) {
6005 if ( aExp
== 0x7FFF ) {
6006 if ((uint64_t)(aSig
<< 1)) {
6007 return propagateFloatx80NaN(a
, b
, status
);
6011 if ( bExp
== 0 ) --expDiff
;
6012 shift64ExtraRightJamming( bSig
, 0, expDiff
, &bSig
, &zSig1
);
6015 else if ( expDiff
< 0 ) {
6016 if ( bExp
== 0x7FFF ) {
6017 if ((uint64_t)(bSig
<< 1)) {
6018 return propagateFloatx80NaN(a
, b
, status
);
6020 return packFloatx80(zSign
,
6021 floatx80_infinity_high
,
6022 floatx80_infinity_low
);
6024 if ( aExp
== 0 ) ++expDiff
;
6025 shift64ExtraRightJamming( aSig
, 0, - expDiff
, &aSig
, &zSig1
);
6029 if ( aExp
== 0x7FFF ) {
6030 if ( (uint64_t) ( ( aSig
| bSig
)<<1 ) ) {
6031 return propagateFloatx80NaN(a
, b
, status
);
6036 zSig0
= aSig
+ bSig
;
6038 if ((aSig
| bSig
) & UINT64_C(0x8000000000000000) && zSig0
< aSig
) {
6039 /* At least one of the values is a pseudo-denormal,
6040 * and there is a carry out of the result. */
6045 return packFloatx80(zSign
, 0, 0);
6047 normalizeFloatx80Subnormal( zSig0
, &zExp
, &zSig0
);
6053 zSig0
= aSig
+ bSig
;
6054 if ( (int64_t) zSig0
< 0 ) goto roundAndPack
;
6056 shift64ExtraRightJamming( zSig0
, zSig1
, 1, &zSig0
, &zSig1
);
6057 zSig0
|= UINT64_C(0x8000000000000000);
6060 return roundAndPackFloatx80(status
->floatx80_rounding_precision
,
6061 zSign
, zExp
, zSig0
, zSig1
, status
);
6064 /*----------------------------------------------------------------------------
6065 | Returns the result of subtracting the absolute values of the extended
6066 | double-precision floating-point values `a' and `b'. If `zSign' is 1, the
6067 | difference is negated before being returned. `zSign' is ignored if the
6068 | result is a NaN. The subtraction is performed according to the IEC/IEEE
6069 | Standard for Binary Floating-Point Arithmetic.
6070 *----------------------------------------------------------------------------*/
6072 static floatx80
subFloatx80Sigs(floatx80 a
, floatx80 b
, bool zSign
,
6073 float_status
*status
)
6075 int32_t aExp
, bExp
, zExp
;
6076 uint64_t aSig
, bSig
, zSig0
, zSig1
;
6079 aSig
= extractFloatx80Frac( a
);
6080 aExp
= extractFloatx80Exp( a
);
6081 bSig
= extractFloatx80Frac( b
);
6082 bExp
= extractFloatx80Exp( b
);
6083 expDiff
= aExp
- bExp
;
6084 if ( 0 < expDiff
) goto aExpBigger
;
6085 if ( expDiff
< 0 ) goto bExpBigger
;
6086 if ( aExp
== 0x7FFF ) {
6087 if ( (uint64_t) ( ( aSig
| bSig
)<<1 ) ) {
6088 return propagateFloatx80NaN(a
, b
, status
);
6090 float_raise(float_flag_invalid
, status
);
6091 return floatx80_default_nan(status
);
6098 if ( bSig
< aSig
) goto aBigger
;
6099 if ( aSig
< bSig
) goto bBigger
;
6100 return packFloatx80(status
->float_rounding_mode
== float_round_down
, 0, 0);
6102 if ( bExp
== 0x7FFF ) {
6103 if ((uint64_t)(bSig
<< 1)) {
6104 return propagateFloatx80NaN(a
, b
, status
);
6106 return packFloatx80(zSign
^ 1, floatx80_infinity_high
,
6107 floatx80_infinity_low
);
6109 if ( aExp
== 0 ) ++expDiff
;
6110 shift128RightJamming( aSig
, 0, - expDiff
, &aSig
, &zSig1
);
6112 sub128( bSig
, 0, aSig
, zSig1
, &zSig0
, &zSig1
);
6115 goto normalizeRoundAndPack
;
6117 if ( aExp
== 0x7FFF ) {
6118 if ((uint64_t)(aSig
<< 1)) {
6119 return propagateFloatx80NaN(a
, b
, status
);
6123 if ( bExp
== 0 ) --expDiff
;
6124 shift128RightJamming( bSig
, 0, expDiff
, &bSig
, &zSig1
);
6126 sub128( aSig
, 0, bSig
, zSig1
, &zSig0
, &zSig1
);
6128 normalizeRoundAndPack
:
6129 return normalizeRoundAndPackFloatx80(status
->floatx80_rounding_precision
,
6130 zSign
, zExp
, zSig0
, zSig1
, status
);
6133 /*----------------------------------------------------------------------------
6134 | Returns the result of adding the extended double-precision floating-point
6135 | values `a' and `b'. The operation is performed according to the IEC/IEEE
6136 | Standard for Binary Floating-Point Arithmetic.
6137 *----------------------------------------------------------------------------*/
6139 floatx80
floatx80_add(floatx80 a
, floatx80 b
, float_status
*status
)
6143 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
6144 float_raise(float_flag_invalid
, status
);
6145 return floatx80_default_nan(status
);
6147 aSign
= extractFloatx80Sign( a
);
6148 bSign
= extractFloatx80Sign( b
);
6149 if ( aSign
== bSign
) {
6150 return addFloatx80Sigs(a
, b
, aSign
, status
);
6153 return subFloatx80Sigs(a
, b
, aSign
, status
);
6158 /*----------------------------------------------------------------------------
6159 | Returns the result of subtracting the extended double-precision floating-
6160 | point values `a' and `b'. The operation is performed according to the
6161 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6162 *----------------------------------------------------------------------------*/
6164 floatx80
floatx80_sub(floatx80 a
, floatx80 b
, float_status
*status
)
6168 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
6169 float_raise(float_flag_invalid
, status
);
6170 return floatx80_default_nan(status
);
6172 aSign
= extractFloatx80Sign( a
);
6173 bSign
= extractFloatx80Sign( b
);
6174 if ( aSign
== bSign
) {
6175 return subFloatx80Sigs(a
, b
, aSign
, status
);
6178 return addFloatx80Sigs(a
, b
, aSign
, status
);
6183 /*----------------------------------------------------------------------------
6184 | Returns the result of multiplying the extended double-precision floating-
6185 | point values `a' and `b'. The operation is performed according to the
6186 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6187 *----------------------------------------------------------------------------*/
6189 floatx80
floatx80_mul(floatx80 a
, floatx80 b
, float_status
*status
)
6191 bool aSign
, bSign
, zSign
;
6192 int32_t aExp
, bExp
, zExp
;
6193 uint64_t aSig
, bSig
, zSig0
, zSig1
;
6195 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
6196 float_raise(float_flag_invalid
, status
);
6197 return floatx80_default_nan(status
);
6199 aSig
= extractFloatx80Frac( a
);
6200 aExp
= extractFloatx80Exp( a
);
6201 aSign
= extractFloatx80Sign( a
);
6202 bSig
= extractFloatx80Frac( b
);
6203 bExp
= extractFloatx80Exp( b
);
6204 bSign
= extractFloatx80Sign( b
);
6205 zSign
= aSign
^ bSign
;
6206 if ( aExp
== 0x7FFF ) {
6207 if ( (uint64_t) ( aSig
<<1 )
6208 || ( ( bExp
== 0x7FFF ) && (uint64_t) ( bSig
<<1 ) ) ) {
6209 return propagateFloatx80NaN(a
, b
, status
);
6211 if ( ( bExp
| bSig
) == 0 ) goto invalid
;
6212 return packFloatx80(zSign
, floatx80_infinity_high
,
6213 floatx80_infinity_low
);
6215 if ( bExp
== 0x7FFF ) {
6216 if ((uint64_t)(bSig
<< 1)) {
6217 return propagateFloatx80NaN(a
, b
, status
);
6219 if ( ( aExp
| aSig
) == 0 ) {
6221 float_raise(float_flag_invalid
, status
);
6222 return floatx80_default_nan(status
);
6224 return packFloatx80(zSign
, floatx80_infinity_high
,
6225 floatx80_infinity_low
);
6228 if ( aSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
6229 normalizeFloatx80Subnormal( aSig
, &aExp
, &aSig
);
6232 if ( bSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
6233 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
6235 zExp
= aExp
+ bExp
- 0x3FFE;
6236 mul64To128( aSig
, bSig
, &zSig0
, &zSig1
);
6237 if ( 0 < (int64_t) zSig0
) {
6238 shortShift128Left( zSig0
, zSig1
, 1, &zSig0
, &zSig1
);
6241 return roundAndPackFloatx80(status
->floatx80_rounding_precision
,
6242 zSign
, zExp
, zSig0
, zSig1
, status
);
6245 /*----------------------------------------------------------------------------
6246 | Returns the result of dividing the extended double-precision floating-point
6247 | value `a' by the corresponding value `b'. The operation is performed
6248 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6249 *----------------------------------------------------------------------------*/
6251 floatx80
floatx80_div(floatx80 a
, floatx80 b
, float_status
*status
)
6253 bool aSign
, bSign
, zSign
;
6254 int32_t aExp
, bExp
, zExp
;
6255 uint64_t aSig
, bSig
, zSig0
, zSig1
;
6256 uint64_t rem0
, rem1
, rem2
, term0
, term1
, term2
;
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 return propagateFloatx80NaN(a
, b
, status
);
6273 if ( bExp
== 0x7FFF ) {
6274 if ((uint64_t)(bSig
<< 1)) {
6275 return propagateFloatx80NaN(a
, b
, status
);
6279 return packFloatx80(zSign
, floatx80_infinity_high
,
6280 floatx80_infinity_low
);
6282 if ( bExp
== 0x7FFF ) {
6283 if ((uint64_t)(bSig
<< 1)) {
6284 return propagateFloatx80NaN(a
, b
, status
);
6286 return packFloatx80( zSign
, 0, 0 );
6290 if ( ( aExp
| aSig
) == 0 ) {
6292 float_raise(float_flag_invalid
, status
);
6293 return floatx80_default_nan(status
);
6295 float_raise(float_flag_divbyzero
, status
);
6296 return packFloatx80(zSign
, floatx80_infinity_high
,
6297 floatx80_infinity_low
);
6299 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
6302 if ( aSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
6303 normalizeFloatx80Subnormal( aSig
, &aExp
, &aSig
);
6305 zExp
= aExp
- bExp
+ 0x3FFE;
6307 if ( bSig
<= aSig
) {
6308 shift128Right( aSig
, 0, 1, &aSig
, &rem1
);
6311 zSig0
= estimateDiv128To64( aSig
, rem1
, bSig
);
6312 mul64To128( bSig
, zSig0
, &term0
, &term1
);
6313 sub128( aSig
, rem1
, term0
, term1
, &rem0
, &rem1
);
6314 while ( (int64_t) rem0
< 0 ) {
6316 add128( rem0
, rem1
, 0, bSig
, &rem0
, &rem1
);
6318 zSig1
= estimateDiv128To64( rem1
, 0, bSig
);
6319 if ( (uint64_t) ( zSig1
<<1 ) <= 8 ) {
6320 mul64To128( bSig
, zSig1
, &term1
, &term2
);
6321 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
6322 while ( (int64_t) rem1
< 0 ) {
6324 add128( rem1
, rem2
, 0, bSig
, &rem1
, &rem2
);
6326 zSig1
|= ( ( rem1
| rem2
) != 0 );
6328 return roundAndPackFloatx80(status
->floatx80_rounding_precision
,
6329 zSign
, zExp
, zSig0
, zSig1
, status
);
6332 /*----------------------------------------------------------------------------
6333 | Returns the remainder of the extended double-precision floating-point value
6334 | `a' with respect to the corresponding value `b'. The operation is performed
6335 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic,
6336 | if 'mod' is false; if 'mod' is true, return the remainder based on truncating
6337 | the quotient toward zero instead. '*quotient' is set to the low 64 bits of
6338 | the absolute value of the integer quotient.
6339 *----------------------------------------------------------------------------*/
6341 floatx80
floatx80_modrem(floatx80 a
, floatx80 b
, bool mod
, uint64_t *quotient
,
6342 float_status
*status
)
6345 int32_t aExp
, bExp
, expDiff
, aExpOrig
;
6346 uint64_t aSig0
, aSig1
, bSig
;
6347 uint64_t q
, term0
, term1
, alternateASig0
, alternateASig1
;
6350 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
6351 float_raise(float_flag_invalid
, status
);
6352 return floatx80_default_nan(status
);
6354 aSig0
= extractFloatx80Frac( a
);
6355 aExpOrig
= aExp
= extractFloatx80Exp( a
);
6356 aSign
= extractFloatx80Sign( a
);
6357 bSig
= extractFloatx80Frac( b
);
6358 bExp
= extractFloatx80Exp( b
);
6359 if ( aExp
== 0x7FFF ) {
6360 if ( (uint64_t) ( aSig0
<<1 )
6361 || ( ( bExp
== 0x7FFF ) && (uint64_t) ( bSig
<<1 ) ) ) {
6362 return propagateFloatx80NaN(a
, b
, status
);
6366 if ( bExp
== 0x7FFF ) {
6367 if ((uint64_t)(bSig
<< 1)) {
6368 return propagateFloatx80NaN(a
, b
, status
);
6370 if (aExp
== 0 && aSig0
>> 63) {
6372 * Pseudo-denormal argument must be returned in normalized
6375 return packFloatx80(aSign
, 1, aSig0
);
6382 float_raise(float_flag_invalid
, status
);
6383 return floatx80_default_nan(status
);
6385 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
6388 if ( aSig0
== 0 ) return a
;
6389 normalizeFloatx80Subnormal( aSig0
, &aExp
, &aSig0
);
6392 expDiff
= aExp
- bExp
;
6394 if ( expDiff
< 0 ) {
6395 if ( mod
|| expDiff
< -1 ) {
6396 if (aExp
== 1 && aExpOrig
== 0) {
6398 * Pseudo-denormal argument must be returned in
6401 return packFloatx80(aSign
, aExp
, aSig0
);
6405 shift128Right( aSig0
, 0, 1, &aSig0
, &aSig1
);
6408 *quotient
= q
= ( bSig
<= aSig0
);
6409 if ( q
) aSig0
-= bSig
;
6411 while ( 0 < expDiff
) {
6412 q
= estimateDiv128To64( aSig0
, aSig1
, bSig
);
6413 q
= ( 2 < q
) ? q
- 2 : 0;
6414 mul64To128( bSig
, q
, &term0
, &term1
);
6415 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
6416 shortShift128Left( aSig0
, aSig1
, 62, &aSig0
, &aSig1
);
6422 if ( 0 < expDiff
) {
6423 q
= estimateDiv128To64( aSig0
, aSig1
, bSig
);
6424 q
= ( 2 < q
) ? q
- 2 : 0;
6426 mul64To128( bSig
, q
<<( 64 - expDiff
), &term0
, &term1
);
6427 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
6428 shortShift128Left( 0, bSig
, 64 - expDiff
, &term0
, &term1
);
6429 while ( le128( term0
, term1
, aSig0
, aSig1
) ) {
6431 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
6434 *quotient
<<= expDiff
;
6445 sub128( term0
, term1
, aSig0
, aSig1
, &alternateASig0
, &alternateASig1
);
6446 if ( lt128( alternateASig0
, alternateASig1
, aSig0
, aSig1
)
6447 || ( eq128( alternateASig0
, alternateASig1
, aSig0
, aSig1
)
6450 aSig0
= alternateASig0
;
6451 aSig1
= alternateASig1
;
6457 normalizeRoundAndPackFloatx80(
6458 80, zSign
, bExp
+ expDiff
, aSig0
, aSig1
, status
);
6462 /*----------------------------------------------------------------------------
6463 | Returns the remainder of the extended double-precision floating-point value
6464 | `a' with respect to the corresponding value `b'. The operation is performed
6465 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6466 *----------------------------------------------------------------------------*/
6468 floatx80
floatx80_rem(floatx80 a
, floatx80 b
, float_status
*status
)
6471 return floatx80_modrem(a
, b
, false, "ient
, status
);
6474 /*----------------------------------------------------------------------------
6475 | Returns the remainder of the extended double-precision floating-point value
6476 | `a' with respect to the corresponding value `b', with the quotient truncated
6478 *----------------------------------------------------------------------------*/
6480 floatx80
floatx80_mod(floatx80 a
, floatx80 b
, float_status
*status
)
6483 return floatx80_modrem(a
, b
, true, "ient
, status
);
6486 /*----------------------------------------------------------------------------
6487 | Returns the square root of the extended double-precision floating-point
6488 | value `a'. The operation is performed according to the IEC/IEEE Standard
6489 | for Binary Floating-Point Arithmetic.
6490 *----------------------------------------------------------------------------*/
6492 floatx80
floatx80_sqrt(floatx80 a
, float_status
*status
)
6496 uint64_t aSig0
, aSig1
, zSig0
, zSig1
, doubleZSig0
;
6497 uint64_t rem0
, rem1
, rem2
, rem3
, term0
, term1
, term2
, term3
;
6499 if (floatx80_invalid_encoding(a
)) {
6500 float_raise(float_flag_invalid
, status
);
6501 return floatx80_default_nan(status
);
6503 aSig0
= extractFloatx80Frac( a
);
6504 aExp
= extractFloatx80Exp( a
);
6505 aSign
= extractFloatx80Sign( a
);
6506 if ( aExp
== 0x7FFF ) {
6507 if ((uint64_t)(aSig0
<< 1)) {
6508 return propagateFloatx80NaN(a
, a
, status
);
6510 if ( ! aSign
) return a
;
6514 if ( ( aExp
| aSig0
) == 0 ) return a
;
6516 float_raise(float_flag_invalid
, status
);
6517 return floatx80_default_nan(status
);
6520 if ( aSig0
== 0 ) return packFloatx80( 0, 0, 0 );
6521 normalizeFloatx80Subnormal( aSig0
, &aExp
, &aSig0
);
6523 zExp
= ( ( aExp
- 0x3FFF )>>1 ) + 0x3FFF;
6524 zSig0
= estimateSqrt32( aExp
, aSig0
>>32 );
6525 shift128Right( aSig0
, 0, 2 + ( aExp
& 1 ), &aSig0
, &aSig1
);
6526 zSig0
= estimateDiv128To64( aSig0
, aSig1
, zSig0
<<32 ) + ( zSig0
<<30 );
6527 doubleZSig0
= zSig0
<<1;
6528 mul64To128( zSig0
, zSig0
, &term0
, &term1
);
6529 sub128( aSig0
, aSig1
, term0
, term1
, &rem0
, &rem1
);
6530 while ( (int64_t) rem0
< 0 ) {
6533 add128( rem0
, rem1
, zSig0
>>63, doubleZSig0
| 1, &rem0
, &rem1
);
6535 zSig1
= estimateDiv128To64( rem1
, 0, doubleZSig0
);
6536 if ( ( zSig1
& UINT64_C(0x3FFFFFFFFFFFFFFF) ) <= 5 ) {
6537 if ( zSig1
== 0 ) zSig1
= 1;
6538 mul64To128( doubleZSig0
, zSig1
, &term1
, &term2
);
6539 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
6540 mul64To128( zSig1
, zSig1
, &term2
, &term3
);
6541 sub192( rem1
, rem2
, 0, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
6542 while ( (int64_t) rem1
< 0 ) {
6544 shortShift128Left( 0, zSig1
, 1, &term2
, &term3
);
6546 term2
|= doubleZSig0
;
6547 add192( rem1
, rem2
, rem3
, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
6549 zSig1
|= ( ( rem1
| rem2
| rem3
) != 0 );
6551 shortShift128Left( 0, zSig1
, 1, &zSig0
, &zSig1
);
6552 zSig0
|= doubleZSig0
;
6553 return roundAndPackFloatx80(status
->floatx80_rounding_precision
,
6554 0, zExp
, zSig0
, zSig1
, status
);
6557 /*----------------------------------------------------------------------------
6558 | Returns the result of converting the quadruple-precision floating-point
6559 | value `a' to the 32-bit two's complement integer format. The conversion
6560 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
6561 | Arithmetic---which means in particular that the conversion is rounded
6562 | according to the current rounding mode. If `a' is a NaN, the largest
6563 | positive integer is returned. Otherwise, if the conversion overflows, the
6564 | largest integer with the same sign as `a' is returned.
6565 *----------------------------------------------------------------------------*/
6567 int32_t float128_to_int32(float128 a
, float_status
*status
)
6570 int32_t aExp
, shiftCount
;
6571 uint64_t aSig0
, aSig1
;
6573 aSig1
= extractFloat128Frac1( a
);
6574 aSig0
= extractFloat128Frac0( a
);
6575 aExp
= extractFloat128Exp( a
);
6576 aSign
= extractFloat128Sign( a
);
6577 if ( ( aExp
== 0x7FFF ) && ( aSig0
| aSig1
) ) aSign
= 0;
6578 if ( aExp
) aSig0
|= UINT64_C(0x0001000000000000);
6579 aSig0
|= ( aSig1
!= 0 );
6580 shiftCount
= 0x4028 - aExp
;
6581 if ( 0 < shiftCount
) shift64RightJamming( aSig0
, shiftCount
, &aSig0
);
6582 return roundAndPackInt32(aSign
, aSig0
, status
);
6586 /*----------------------------------------------------------------------------
6587 | Returns the result of converting the quadruple-precision floating-point
6588 | value `a' to the 32-bit two's complement integer format. The conversion
6589 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
6590 | Arithmetic, except that the conversion is always rounded toward zero. If
6591 | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
6592 | conversion overflows, the largest integer with the same sign as `a' is
6594 *----------------------------------------------------------------------------*/
6596 int32_t float128_to_int32_round_to_zero(float128 a
, float_status
*status
)
6599 int32_t aExp
, shiftCount
;
6600 uint64_t aSig0
, aSig1
, savedASig
;
6603 aSig1
= extractFloat128Frac1( a
);
6604 aSig0
= extractFloat128Frac0( a
);
6605 aExp
= extractFloat128Exp( a
);
6606 aSign
= extractFloat128Sign( a
);
6607 aSig0
|= ( aSig1
!= 0 );
6608 if ( 0x401E < aExp
) {
6609 if ( ( aExp
== 0x7FFF ) && aSig0
) aSign
= 0;
6612 else if ( aExp
< 0x3FFF ) {
6613 if (aExp
|| aSig0
) {
6614 float_raise(float_flag_inexact
, status
);
6618 aSig0
|= UINT64_C(0x0001000000000000);
6619 shiftCount
= 0x402F - aExp
;
6621 aSig0
>>= shiftCount
;
6623 if ( aSign
) z
= - z
;
6624 if ( ( z
< 0 ) ^ aSign
) {
6626 float_raise(float_flag_invalid
, status
);
6627 return aSign
? INT32_MIN
: INT32_MAX
;
6629 if ( ( aSig0
<<shiftCount
) != savedASig
) {
6630 float_raise(float_flag_inexact
, status
);
6636 /*----------------------------------------------------------------------------
6637 | Returns the result of converting the quadruple-precision floating-point
6638 | value `a' to the 64-bit two's complement integer format. The conversion
6639 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
6640 | Arithmetic---which means in particular that the conversion is rounded
6641 | according to the current rounding mode. If `a' is a NaN, the largest
6642 | positive integer is returned. Otherwise, if the conversion overflows, the
6643 | largest integer with the same sign as `a' is returned.
6644 *----------------------------------------------------------------------------*/
6646 int64_t float128_to_int64(float128 a
, float_status
*status
)
6649 int32_t aExp
, shiftCount
;
6650 uint64_t aSig0
, aSig1
;
6652 aSig1
= extractFloat128Frac1( a
);
6653 aSig0
= extractFloat128Frac0( a
);
6654 aExp
= extractFloat128Exp( a
);
6655 aSign
= extractFloat128Sign( a
);
6656 if ( aExp
) aSig0
|= UINT64_C(0x0001000000000000);
6657 shiftCount
= 0x402F - aExp
;
6658 if ( shiftCount
<= 0 ) {
6659 if ( 0x403E < aExp
) {
6660 float_raise(float_flag_invalid
, status
);
6662 || ( ( aExp
== 0x7FFF )
6663 && ( aSig1
|| ( aSig0
!= UINT64_C(0x0001000000000000) ) )
6670 shortShift128Left( aSig0
, aSig1
, - shiftCount
, &aSig0
, &aSig1
);
6673 shift64ExtraRightJamming( aSig0
, aSig1
, shiftCount
, &aSig0
, &aSig1
);
6675 return roundAndPackInt64(aSign
, aSig0
, aSig1
, status
);
6679 /*----------------------------------------------------------------------------
6680 | Returns the result of converting the quadruple-precision floating-point
6681 | value `a' to the 64-bit two's complement integer format. The conversion
6682 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
6683 | Arithmetic, except that the conversion is always rounded toward zero.
6684 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
6685 | the conversion overflows, the largest integer with the same sign as `a' is
6687 *----------------------------------------------------------------------------*/
6689 int64_t float128_to_int64_round_to_zero(float128 a
, float_status
*status
)
6692 int32_t aExp
, shiftCount
;
6693 uint64_t aSig0
, aSig1
;
6696 aSig1
= extractFloat128Frac1( a
);
6697 aSig0
= extractFloat128Frac0( a
);
6698 aExp
= extractFloat128Exp( a
);
6699 aSign
= extractFloat128Sign( a
);
6700 if ( aExp
) aSig0
|= UINT64_C(0x0001000000000000);
6701 shiftCount
= aExp
- 0x402F;
6702 if ( 0 < shiftCount
) {
6703 if ( 0x403E <= aExp
) {
6704 aSig0
&= UINT64_C(0x0000FFFFFFFFFFFF);
6705 if ( ( a
.high
== UINT64_C(0xC03E000000000000) )
6706 && ( aSig1
< UINT64_C(0x0002000000000000) ) ) {
6708 float_raise(float_flag_inexact
, status
);
6712 float_raise(float_flag_invalid
, status
);
6713 if ( ! aSign
|| ( ( aExp
== 0x7FFF ) && ( aSig0
| aSig1
) ) ) {
6719 z
= ( aSig0
<<shiftCount
) | ( aSig1
>>( ( - shiftCount
) & 63 ) );
6720 if ( (uint64_t) ( aSig1
<<shiftCount
) ) {
6721 float_raise(float_flag_inexact
, status
);
6725 if ( aExp
< 0x3FFF ) {
6726 if ( aExp
| aSig0
| aSig1
) {
6727 float_raise(float_flag_inexact
, status
);
6731 z
= aSig0
>>( - shiftCount
);
6733 || ( shiftCount
&& (uint64_t) ( aSig0
<<( shiftCount
& 63 ) ) ) ) {
6734 float_raise(float_flag_inexact
, status
);
6737 if ( aSign
) z
= - z
;
6742 /*----------------------------------------------------------------------------
6743 | Returns the result of converting the quadruple-precision floating-point value
6744 | `a' to the 64-bit unsigned integer format. The conversion is
6745 | performed according to the IEC/IEEE Standard for Binary Floating-Point
6746 | Arithmetic---which means in particular that the conversion is rounded
6747 | according to the current rounding mode. If `a' is a NaN, the largest
6748 | positive integer is returned. If the conversion overflows, the
6749 | largest unsigned integer is returned. If 'a' is negative, the value is
6750 | rounded and zero is returned; negative values that do not round to zero
6751 | will raise the inexact exception.
6752 *----------------------------------------------------------------------------*/
6754 uint64_t float128_to_uint64(float128 a
, float_status
*status
)
6759 uint64_t aSig0
, aSig1
;
6761 aSig0
= extractFloat128Frac0(a
);
6762 aSig1
= extractFloat128Frac1(a
);
6763 aExp
= extractFloat128Exp(a
);
6764 aSign
= extractFloat128Sign(a
);
6765 if (aSign
&& (aExp
> 0x3FFE)) {
6766 float_raise(float_flag_invalid
, status
);
6767 if (float128_is_any_nan(a
)) {
6774 aSig0
|= UINT64_C(0x0001000000000000);
6776 shiftCount
= 0x402F - aExp
;
6777 if (shiftCount
<= 0) {
6778 if (0x403E < aExp
) {
6779 float_raise(float_flag_invalid
, status
);
6782 shortShift128Left(aSig0
, aSig1
, -shiftCount
, &aSig0
, &aSig1
);
6784 shift64ExtraRightJamming(aSig0
, aSig1
, shiftCount
, &aSig0
, &aSig1
);
6786 return roundAndPackUint64(aSign
, aSig0
, aSig1
, status
);
6789 uint64_t float128_to_uint64_round_to_zero(float128 a
, float_status
*status
)
6792 signed char current_rounding_mode
= status
->float_rounding_mode
;
6794 set_float_rounding_mode(float_round_to_zero
, status
);
6795 v
= float128_to_uint64(a
, status
);
6796 set_float_rounding_mode(current_rounding_mode
, status
);
6801 /*----------------------------------------------------------------------------
6802 | Returns the result of converting the quadruple-precision floating-point
6803 | value `a' to the 32-bit unsigned integer format. The conversion
6804 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
6805 | Arithmetic except that the conversion is always rounded toward zero.
6806 | If `a' is a NaN, the largest positive integer is returned. Otherwise,
6807 | if the conversion overflows, the largest unsigned integer is returned.
6808 | If 'a' is negative, the value is rounded and zero is returned; negative
6809 | values that do not round to zero will raise the inexact exception.
6810 *----------------------------------------------------------------------------*/
6812 uint32_t float128_to_uint32_round_to_zero(float128 a
, float_status
*status
)
6816 int old_exc_flags
= get_float_exception_flags(status
);
6818 v
= float128_to_uint64_round_to_zero(a
, status
);
6819 if (v
> 0xffffffff) {
6824 set_float_exception_flags(old_exc_flags
, status
);
6825 float_raise(float_flag_invalid
, status
);
6829 /*----------------------------------------------------------------------------
6830 | Returns the result of converting the quadruple-precision floating-point value
6831 | `a' to the 32-bit unsigned integer format. The conversion is
6832 | performed according to the IEC/IEEE Standard for Binary Floating-Point
6833 | Arithmetic---which means in particular that the conversion is rounded
6834 | according to the current rounding mode. If `a' is a NaN, the largest
6835 | positive integer is returned. If the conversion overflows, the
6836 | largest unsigned integer is returned. If 'a' is negative, the value is
6837 | rounded and zero is returned; negative values that do not round to zero
6838 | will raise the inexact exception.
6839 *----------------------------------------------------------------------------*/
6841 uint32_t float128_to_uint32(float128 a
, float_status
*status
)
6845 int old_exc_flags
= get_float_exception_flags(status
);
6847 v
= float128_to_uint64(a
, status
);
6848 if (v
> 0xffffffff) {
6853 set_float_exception_flags(old_exc_flags
, status
);
6854 float_raise(float_flag_invalid
, status
);
6858 /*----------------------------------------------------------------------------
6859 | Returns the result of converting the quadruple-precision floating-point
6860 | value `a' to the extended double-precision floating-point format. The
6861 | conversion is performed according to the IEC/IEEE Standard for Binary
6862 | Floating-Point Arithmetic.
6863 *----------------------------------------------------------------------------*/
6865 floatx80
float128_to_floatx80(float128 a
, float_status
*status
)
6869 uint64_t aSig0
, aSig1
;
6871 aSig1
= extractFloat128Frac1( a
);
6872 aSig0
= extractFloat128Frac0( a
);
6873 aExp
= extractFloat128Exp( a
);
6874 aSign
= extractFloat128Sign( a
);
6875 if ( aExp
== 0x7FFF ) {
6876 if ( aSig0
| aSig1
) {
6877 floatx80 res
= commonNaNToFloatx80(float128ToCommonNaN(a
, status
),
6879 return floatx80_silence_nan(res
, status
);
6881 return packFloatx80(aSign
, floatx80_infinity_high
,
6882 floatx80_infinity_low
);
6885 if ( ( aSig0
| aSig1
) == 0 ) return packFloatx80( aSign
, 0, 0 );
6886 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
6889 aSig0
|= UINT64_C(0x0001000000000000);
6891 shortShift128Left( aSig0
, aSig1
, 15, &aSig0
, &aSig1
);
6892 return roundAndPackFloatx80(80, aSign
, aExp
, aSig0
, aSig1
, status
);
6896 /*----------------------------------------------------------------------------
6897 | Returns the remainder of the quadruple-precision floating-point value `a'
6898 | with respect to the corresponding value `b'. The operation is performed
6899 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6900 *----------------------------------------------------------------------------*/
6902 float128
float128_rem(float128 a
, float128 b
, float_status
*status
)
6905 int32_t aExp
, bExp
, expDiff
;
6906 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, q
, term0
, term1
, term2
;
6907 uint64_t allZero
, alternateASig0
, alternateASig1
, sigMean1
;
6910 aSig1
= extractFloat128Frac1( a
);
6911 aSig0
= extractFloat128Frac0( a
);
6912 aExp
= extractFloat128Exp( a
);
6913 aSign
= extractFloat128Sign( a
);
6914 bSig1
= extractFloat128Frac1( b
);
6915 bSig0
= extractFloat128Frac0( b
);
6916 bExp
= extractFloat128Exp( b
);
6917 if ( aExp
== 0x7FFF ) {
6918 if ( ( aSig0
| aSig1
)
6919 || ( ( bExp
== 0x7FFF ) && ( bSig0
| bSig1
) ) ) {
6920 return propagateFloat128NaN(a
, b
, status
);
6924 if ( bExp
== 0x7FFF ) {
6925 if (bSig0
| bSig1
) {
6926 return propagateFloat128NaN(a
, b
, status
);
6931 if ( ( bSig0
| bSig1
) == 0 ) {
6933 float_raise(float_flag_invalid
, status
);
6934 return float128_default_nan(status
);
6936 normalizeFloat128Subnormal( bSig0
, bSig1
, &bExp
, &bSig0
, &bSig1
);
6939 if ( ( aSig0
| aSig1
) == 0 ) return a
;
6940 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
6942 expDiff
= aExp
- bExp
;
6943 if ( expDiff
< -1 ) return a
;
6945 aSig0
| UINT64_C(0x0001000000000000),
6947 15 - ( expDiff
< 0 ),
6952 bSig0
| UINT64_C(0x0001000000000000), bSig1
, 15, &bSig0
, &bSig1
);
6953 q
= le128( bSig0
, bSig1
, aSig0
, aSig1
);
6954 if ( q
) sub128( aSig0
, aSig1
, bSig0
, bSig1
, &aSig0
, &aSig1
);
6956 while ( 0 < expDiff
) {
6957 q
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
6958 q
= ( 4 < q
) ? q
- 4 : 0;
6959 mul128By64To192( bSig0
, bSig1
, q
, &term0
, &term1
, &term2
);
6960 shortShift192Left( term0
, term1
, term2
, 61, &term1
, &term2
, &allZero
);
6961 shortShift128Left( aSig0
, aSig1
, 61, &aSig0
, &allZero
);
6962 sub128( aSig0
, 0, term1
, term2
, &aSig0
, &aSig1
);
6965 if ( -64 < expDiff
) {
6966 q
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
6967 q
= ( 4 < q
) ? q
- 4 : 0;
6969 shift128Right( bSig0
, bSig1
, 12, &bSig0
, &bSig1
);
6971 if ( expDiff
< 0 ) {
6972 shift128Right( aSig0
, aSig1
, - expDiff
, &aSig0
, &aSig1
);
6975 shortShift128Left( aSig0
, aSig1
, expDiff
, &aSig0
, &aSig1
);
6977 mul128By64To192( bSig0
, bSig1
, q
, &term0
, &term1
, &term2
);
6978 sub128( aSig0
, aSig1
, term1
, term2
, &aSig0
, &aSig1
);
6981 shift128Right( aSig0
, aSig1
, 12, &aSig0
, &aSig1
);
6982 shift128Right( bSig0
, bSig1
, 12, &bSig0
, &bSig1
);
6985 alternateASig0
= aSig0
;
6986 alternateASig1
= aSig1
;
6988 sub128( aSig0
, aSig1
, bSig0
, bSig1
, &aSig0
, &aSig1
);
6989 } while ( 0 <= (int64_t) aSig0
);
6991 aSig0
, aSig1
, alternateASig0
, alternateASig1
, (uint64_t *)&sigMean0
, &sigMean1
);
6992 if ( ( sigMean0
< 0 )
6993 || ( ( ( sigMean0
| sigMean1
) == 0 ) && ( q
& 1 ) ) ) {
6994 aSig0
= alternateASig0
;
6995 aSig1
= alternateASig1
;
6997 zSign
= ( (int64_t) aSig0
< 0 );
6998 if ( zSign
) sub128( 0, 0, aSig0
, aSig1
, &aSig0
, &aSig1
);
6999 return normalizeRoundAndPackFloat128(aSign
^ zSign
, bExp
- 4, aSig0
, aSig1
,
7003 /*----------------------------------------------------------------------------
7004 | Returns the square root of the quadruple-precision floating-point value `a'.
7005 | The operation is performed according to the IEC/IEEE Standard for Binary
7006 | Floating-Point Arithmetic.
7007 *----------------------------------------------------------------------------*/
7009 float128
float128_sqrt(float128 a
, float_status
*status
)
7013 uint64_t aSig0
, aSig1
, zSig0
, zSig1
, zSig2
, doubleZSig0
;
7014 uint64_t rem0
, rem1
, rem2
, rem3
, term0
, term1
, term2
, term3
;
7016 aSig1
= extractFloat128Frac1( a
);
7017 aSig0
= extractFloat128Frac0( a
);
7018 aExp
= extractFloat128Exp( a
);
7019 aSign
= extractFloat128Sign( a
);
7020 if ( aExp
== 0x7FFF ) {
7021 if (aSig0
| aSig1
) {
7022 return propagateFloat128NaN(a
, a
, status
);
7024 if ( ! aSign
) return a
;
7028 if ( ( aExp
| aSig0
| aSig1
) == 0 ) return a
;
7030 float_raise(float_flag_invalid
, status
);
7031 return float128_default_nan(status
);
7034 if ( ( aSig0
| aSig1
) == 0 ) return packFloat128( 0, 0, 0, 0 );
7035 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
7037 zExp
= ( ( aExp
- 0x3FFF )>>1 ) + 0x3FFE;
7038 aSig0
|= UINT64_C(0x0001000000000000);
7039 zSig0
= estimateSqrt32( aExp
, aSig0
>>17 );
7040 shortShift128Left( aSig0
, aSig1
, 13 - ( aExp
& 1 ), &aSig0
, &aSig1
);
7041 zSig0
= estimateDiv128To64( aSig0
, aSig1
, zSig0
<<32 ) + ( zSig0
<<30 );
7042 doubleZSig0
= zSig0
<<1;
7043 mul64To128( zSig0
, zSig0
, &term0
, &term1
);
7044 sub128( aSig0
, aSig1
, term0
, term1
, &rem0
, &rem1
);
7045 while ( (int64_t) rem0
< 0 ) {
7048 add128( rem0
, rem1
, zSig0
>>63, doubleZSig0
| 1, &rem0
, &rem1
);
7050 zSig1
= estimateDiv128To64( rem1
, 0, doubleZSig0
);
7051 if ( ( zSig1
& 0x1FFF ) <= 5 ) {
7052 if ( zSig1
== 0 ) zSig1
= 1;
7053 mul64To128( doubleZSig0
, zSig1
, &term1
, &term2
);
7054 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
7055 mul64To128( zSig1
, zSig1
, &term2
, &term3
);
7056 sub192( rem1
, rem2
, 0, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
7057 while ( (int64_t) rem1
< 0 ) {
7059 shortShift128Left( 0, zSig1
, 1, &term2
, &term3
);
7061 term2
|= doubleZSig0
;
7062 add192( rem1
, rem2
, rem3
, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
7064 zSig1
|= ( ( rem1
| rem2
| rem3
) != 0 );
7066 shift128ExtraRightJamming( zSig0
, zSig1
, 0, 14, &zSig0
, &zSig1
, &zSig2
);
7067 return roundAndPackFloat128(0, zExp
, zSig0
, zSig1
, zSig2
, status
);
7071 static inline FloatRelation
7072 floatx80_compare_internal(floatx80 a
, floatx80 b
, bool is_quiet
,
7073 float_status
*status
)
7077 if (floatx80_invalid_encoding(a
) || floatx80_invalid_encoding(b
)) {
7078 float_raise(float_flag_invalid
, status
);
7079 return float_relation_unordered
;
7081 if (( ( extractFloatx80Exp( a
) == 0x7fff ) &&
7082 ( extractFloatx80Frac( a
)<<1 ) ) ||
7083 ( ( extractFloatx80Exp( b
) == 0x7fff ) &&
7084 ( extractFloatx80Frac( b
)<<1 ) )) {
7086 floatx80_is_signaling_nan(a
, status
) ||
7087 floatx80_is_signaling_nan(b
, status
)) {
7088 float_raise(float_flag_invalid
, status
);
7090 return float_relation_unordered
;
7092 aSign
= extractFloatx80Sign( a
);
7093 bSign
= extractFloatx80Sign( b
);
7094 if ( aSign
!= bSign
) {
7096 if ( ( ( (uint16_t) ( ( a
.high
| b
.high
) << 1 ) ) == 0) &&
7097 ( ( a
.low
| b
.low
) == 0 ) ) {
7099 return float_relation_equal
;
7101 return 1 - (2 * aSign
);
7104 /* Normalize pseudo-denormals before comparison. */
7105 if ((a
.high
& 0x7fff) == 0 && a
.low
& UINT64_C(0x8000000000000000)) {
7108 if ((b
.high
& 0x7fff) == 0 && b
.low
& UINT64_C(0x8000000000000000)) {
7111 if (a
.low
== b
.low
&& a
.high
== b
.high
) {
7112 return float_relation_equal
;
7114 return 1 - 2 * (aSign
^ ( lt128( a
.high
, a
.low
, b
.high
, b
.low
) ));
7119 FloatRelation
floatx80_compare(floatx80 a
, floatx80 b
, float_status
*status
)
7121 return floatx80_compare_internal(a
, b
, 0, status
);
7124 FloatRelation
floatx80_compare_quiet(floatx80 a
, floatx80 b
,
7125 float_status
*status
)
7127 return floatx80_compare_internal(a
, b
, 1, status
);
7130 static inline FloatRelation
7131 float128_compare_internal(float128 a
, float128 b
, bool is_quiet
,
7132 float_status
*status
)
7136 if (( ( extractFloat128Exp( a
) == 0x7fff ) &&
7137 ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) ) ||
7138 ( ( extractFloat128Exp( b
) == 0x7fff ) &&
7139 ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )) {
7141 float128_is_signaling_nan(a
, status
) ||
7142 float128_is_signaling_nan(b
, status
)) {
7143 float_raise(float_flag_invalid
, status
);
7145 return float_relation_unordered
;
7147 aSign
= extractFloat128Sign( a
);
7148 bSign
= extractFloat128Sign( b
);
7149 if ( aSign
!= bSign
) {
7150 if ( ( ( ( a
.high
| b
.high
)<<1 ) | a
.low
| b
.low
) == 0 ) {
7152 return float_relation_equal
;
7154 return 1 - (2 * aSign
);
7157 if (a
.low
== b
.low
&& a
.high
== b
.high
) {
7158 return float_relation_equal
;
7160 return 1 - 2 * (aSign
^ ( lt128( a
.high
, a
.low
, b
.high
, b
.low
) ));
7165 FloatRelation
float128_compare(float128 a
, float128 b
, float_status
*status
)
7167 return float128_compare_internal(a
, b
, 0, status
);
7170 FloatRelation
float128_compare_quiet(float128 a
, float128 b
,
7171 float_status
*status
)
7173 return float128_compare_internal(a
, b
, 1, status
);
7176 floatx80
floatx80_scalbn(floatx80 a
, int n
, float_status
*status
)
7182 if (floatx80_invalid_encoding(a
)) {
7183 float_raise(float_flag_invalid
, status
);
7184 return floatx80_default_nan(status
);
7186 aSig
= extractFloatx80Frac( a
);
7187 aExp
= extractFloatx80Exp( a
);
7188 aSign
= extractFloatx80Sign( a
);
7190 if ( aExp
== 0x7FFF ) {
7192 return propagateFloatx80NaN(a
, a
, status
);
7206 } else if (n
< -0x10000) {
7211 return normalizeRoundAndPackFloatx80(status
->floatx80_rounding_precision
,
7212 aSign
, aExp
, aSig
, 0, status
);
7215 float128
float128_scalbn(float128 a
, int n
, float_status
*status
)
7219 uint64_t aSig0
, aSig1
;
7221 aSig1
= extractFloat128Frac1( a
);
7222 aSig0
= extractFloat128Frac0( a
);
7223 aExp
= extractFloat128Exp( a
);
7224 aSign
= extractFloat128Sign( a
);
7225 if ( aExp
== 0x7FFF ) {
7226 if ( aSig0
| aSig1
) {
7227 return propagateFloat128NaN(a
, a
, status
);
7232 aSig0
|= UINT64_C(0x0001000000000000);
7233 } else if (aSig0
== 0 && aSig1
== 0) {
7241 } else if (n
< -0x10000) {
7246 return normalizeRoundAndPackFloat128( aSign
, aExp
, aSig0
, aSig1
7251 static void __attribute__((constructor
)) softfloat_init(void)
7253 union_float64 ua
, ub
, uc
, ur
;
7255 if (QEMU_NO_HARDFLOAT
) {
7259 * Test that the host's FMA is not obviously broken. For example,
7260 * glibc < 2.23 can perform an incorrect FMA on certain hosts; see
7261 * https://sourceware.org/bugzilla/show_bug.cgi?id=13304
7263 ua
.s
= 0x0020000000000001ULL
;
7264 ub
.s
= 0x3ca0000000000000ULL
;
7265 uc
.s
= 0x0020000000000000ULL
;
7266 ur
.h
= fma(ua
.h
, ub
.h
, uc
.h
);
7267 if (ur
.s
!= 0x0020000000000001ULL
) {
7268 force_soft_fma
= true;