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