if (unlikely(soft_t ## _is_denormal(*a))) { \
*a = soft_t ## _set_sign(soft_t ## _zero, \
soft_t ## _is_neg(*a)); \
- s->float_exception_flags |= float_flag_input_denormal; \
+ float_raise(float_flag_input_denormal, s); \
} \
}
ur.h = hard(ua.h, ub.h);
if (unlikely(f32_is_inf(ur))) {
- s->float_exception_flags |= float_flag_overflow;
+ float_raise(float_flag_overflow, s);
} else if (unlikely(fabsf(ur.h) <= FLT_MIN) && post(ua, ub)) {
goto soft;
}
ur.h = hard(ua.h, ub.h);
if (unlikely(f64_is_inf(ur))) {
- s->float_exception_flags |= float_flag_overflow;
+ float_raise(float_flag_overflow, s);
} else if (unlikely(fabs(ur.h) <= DBL_MIN) && post(ua, ub)) {
goto soft;
}
float_class_snan,
} FloatClass;
+#define float_cmask(bit) (1u << (bit))
+
+enum {
+ float_cmask_zero = float_cmask(float_class_zero),
+ float_cmask_normal = float_cmask(float_class_normal),
+ float_cmask_inf = float_cmask(float_class_inf),
+ float_cmask_qnan = float_cmask(float_class_qnan),
+ float_cmask_snan = float_cmask(float_class_snan),
+
+ float_cmask_infzero = float_cmask_zero | float_cmask_inf,
+ float_cmask_anynan = float_cmask_qnan | float_cmask_snan,
+};
+
+
/* Simple helpers for checking if, or what kind of, NaN we have */
static inline __attribute__((unused)) bool is_nan(FloatClass c)
{
}
/*
- * Structure holding all of the decomposed parts of a float. The
- * exponent is unbiased and the fraction is normalized. All
- * calculations are done with a 64 bit fraction and then rounded as
- * appropriate for the final format.
+ * Structure holding all of the decomposed parts of a float.
+ * The exponent is unbiased and the fraction is normalized.
*
- * Thanks to the packed FloatClass a decent compiler should be able to
- * fit the whole structure into registers and avoid using the stack
- * for parameter passing.
+ * The fraction words are stored in big-endian word ordering,
+ * so that truncation from a larger format to a smaller format
+ * can be done simply by ignoring subsequent elements.
*/
typedef struct {
- uint64_t frac;
- int32_t exp;
FloatClass cls;
bool sign;
-} FloatParts;
+ int32_t exp;
+ union {
+ /* Routines that know the structure may reference the singular name. */
+ uint64_t frac;
+ /*
+ * Routines expanded with multiple structures reference "hi" and "lo"
+ * depending on the operation. In FloatParts64, "hi" and "lo" are
+ * both the same word and aliased here.
+ */
+ uint64_t frac_hi;
+ uint64_t frac_lo;
+ };
+} FloatParts64;
+
+typedef struct {
+ FloatClass cls;
+ bool sign;
+ int32_t exp;
+ uint64_t frac_hi;
+ uint64_t frac_lo;
+} FloatParts128;
-#define DECOMPOSED_BINARY_POINT (64 - 2)
+typedef struct {
+ FloatClass cls;
+ bool sign;
+ int32_t exp;
+ uint64_t frac_hi;
+ uint64_t frac_hm; /* high-middle */
+ uint64_t frac_lm; /* low-middle */
+ uint64_t frac_lo;
+} FloatParts256;
+
+/* These apply to the most significant word of each FloatPartsN. */
+#define DECOMPOSED_BINARY_POINT 63
#define DECOMPOSED_IMPLICIT_BIT (1ull << DECOMPOSED_BINARY_POINT)
-#define DECOMPOSED_OVERFLOW_BIT (DECOMPOSED_IMPLICIT_BIT << 1)
/* Structure holding all of the relevant parameters for a format.
* exp_size: the size of the exponent field
.exp_bias = ((1 << E) - 1) >> 1, \
.exp_max = (1 << E) - 1, \
.frac_size = F, \
- .frac_shift = DECOMPOSED_BINARY_POINT - F, \
- .frac_lsb = 1ull << (DECOMPOSED_BINARY_POINT - F), \
- .frac_lsbm1 = 1ull << ((DECOMPOSED_BINARY_POINT - F) - 1), \
- .round_mask = (1ull << (DECOMPOSED_BINARY_POINT - F)) - 1, \
- .roundeven_mask = (2ull << (DECOMPOSED_BINARY_POINT - F)) - 1
+ .frac_shift = (-F - 1) & 63, \
+ .frac_lsb = 1ull << ((-F - 1) & 63), \
+ .frac_lsbm1 = 1ull << ((-F - 2) & 63), \
+ .round_mask = (1ull << ((-F - 1) & 63)) - 1, \
+ .roundeven_mask = (2ull << ((-F - 1) & 63)) - 1
static const FloatFmt float16_params = {
FLOAT_PARAMS(5, 10)
FLOAT_PARAMS(11, 52)
};
+static const FloatFmt float128_params = {
+ FLOAT_PARAMS(15, 112)
+};
+
/* Unpack a float to parts, but do not canonicalize. */
-static inline FloatParts unpack_raw(FloatFmt fmt, uint64_t raw)
+static void unpack_raw64(FloatParts64 *r, const FloatFmt *fmt, uint64_t raw)
{
- const int sign_pos = fmt.frac_size + fmt.exp_size;
+ const int f_size = fmt->frac_size;
+ const int e_size = fmt->exp_size;
- return (FloatParts) {
+ *r = (FloatParts64) {
.cls = float_class_unclassified,
- .sign = extract64(raw, sign_pos, 1),
- .exp = extract64(raw, fmt.frac_size, fmt.exp_size),
- .frac = extract64(raw, 0, fmt.frac_size),
+ .sign = extract64(raw, f_size + e_size, 1),
+ .exp = extract64(raw, f_size, e_size),
+ .frac = extract64(raw, 0, f_size)
};
}
-static inline FloatParts float16_unpack_raw(float16 f)
+static inline void float16_unpack_raw(FloatParts64 *p, float16 f)
+{
+ unpack_raw64(p, &float16_params, f);
+}
+
+static inline void bfloat16_unpack_raw(FloatParts64 *p, bfloat16 f)
{
- return unpack_raw(float16_params, f);
+ unpack_raw64(p, &bfloat16_params, f);
}
-static inline FloatParts bfloat16_unpack_raw(bfloat16 f)
+static inline void float32_unpack_raw(FloatParts64 *p, float32 f)
{
- return unpack_raw(bfloat16_params, f);
+ unpack_raw64(p, &float32_params, f);
}
-static inline FloatParts float32_unpack_raw(float32 f)
+static inline void float64_unpack_raw(FloatParts64 *p, float64 f)
{
- return unpack_raw(float32_params, f);
+ unpack_raw64(p, &float64_params, f);
}
-static inline FloatParts float64_unpack_raw(float64 f)
+static void float128_unpack_raw(FloatParts128 *p, float128 f)
{
- return unpack_raw(float64_params, f);
+ const int f_size = float128_params.frac_size - 64;
+ const int e_size = float128_params.exp_size;
+
+ *p = (FloatParts128) {
+ .cls = float_class_unclassified,
+ .sign = extract64(f.high, f_size + e_size, 1),
+ .exp = extract64(f.high, f_size, e_size),
+ .frac_hi = extract64(f.high, 0, f_size),
+ .frac_lo = f.low,
+ };
}
/* Pack a float from parts, but do not canonicalize. */
-static inline uint64_t pack_raw(FloatFmt fmt, FloatParts p)
+static uint64_t pack_raw64(const FloatParts64 *p, const FloatFmt *fmt)
+{
+ const int f_size = fmt->frac_size;
+ const int e_size = fmt->exp_size;
+ uint64_t ret;
+
+ ret = (uint64_t)p->sign << (f_size + e_size);
+ ret = deposit64(ret, f_size, e_size, p->exp);
+ ret = deposit64(ret, 0, f_size, p->frac);
+ return ret;
+}
+
+static inline float16 float16_pack_raw(const FloatParts64 *p)
{
- const int sign_pos = fmt.frac_size + fmt.exp_size;
- uint64_t ret = deposit64(p.frac, fmt.frac_size, fmt.exp_size, p.exp);
- return deposit64(ret, sign_pos, 1, p.sign);
+ return make_float16(pack_raw64(p, &float16_params));
}
-static inline float16 float16_pack_raw(FloatParts p)
+static inline bfloat16 bfloat16_pack_raw(const FloatParts64 *p)
{
- return make_float16(pack_raw(float16_params, p));
+ return pack_raw64(p, &bfloat16_params);
}
-static inline bfloat16 bfloat16_pack_raw(FloatParts p)
+static inline float32 float32_pack_raw(const FloatParts64 *p)
{
- return pack_raw(bfloat16_params, p);
+ return make_float32(pack_raw64(p, &float32_params));
}
-static inline float32 float32_pack_raw(FloatParts p)
+static inline float64 float64_pack_raw(const FloatParts64 *p)
{
- return make_float32(pack_raw(float32_params, p));
+ return make_float64(pack_raw64(p, &float64_params));
}
-static inline float64 float64_pack_raw(FloatParts p)
+static float128 float128_pack_raw(const FloatParts128 *p)
{
- return make_float64(pack_raw(float64_params, p));
+ const int f_size = float128_params.frac_size - 64;
+ const int e_size = float128_params.exp_size;
+ uint64_t hi;
+
+ hi = (uint64_t)p->sign << (f_size + e_size);
+ hi = deposit64(hi, f_size, e_size, p->exp);
+ hi = deposit64(hi, 0, f_size, p->frac_hi);
+ return make_float128(hi, p->frac_lo);
}
/*----------------------------------------------------------------------------
*----------------------------------------------------------------------------*/
#include "softfloat-specialize.c.inc"
-/* Canonicalize EXP and FRAC, setting CLS. */
-static FloatParts sf_canonicalize(FloatParts part, const FloatFmt *parm,
- float_status *status)
-{
- if (part.exp == parm->exp_max && !parm->arm_althp) {
- if (part.frac == 0) {
- part.cls = float_class_inf;
- } else {
- part.frac <<= parm->frac_shift;
- part.cls = (parts_is_snan_frac(part.frac, status)
- ? float_class_snan : float_class_qnan);
- }
- } else if (part.exp == 0) {
- if (likely(part.frac == 0)) {
- part.cls = float_class_zero;
- } else if (status->flush_inputs_to_zero) {
- float_raise(float_flag_input_denormal, status);
- part.cls = float_class_zero;
- part.frac = 0;
- } else {
- int shift = clz64(part.frac) - 1;
- part.cls = float_class_normal;
- part.exp = parm->frac_shift - parm->exp_bias - shift + 1;
- part.frac <<= shift;
- }
- } else {
- part.cls = float_class_normal;
- part.exp -= parm->exp_bias;
- part.frac = DECOMPOSED_IMPLICIT_BIT + (part.frac << parm->frac_shift);
- }
- return part;
-}
+#define PARTS_GENERIC_64_128(NAME, P) \
+ QEMU_GENERIC(P, (FloatParts128 *, parts128_##NAME), parts64_##NAME)
-/* Round and uncanonicalize a floating-point number by parts. There
- * are FRAC_SHIFT bits that may require rounding at the bottom of the
- * fraction; these bits will be removed. The exponent will be biased
- * by EXP_BIAS and must be bounded by [EXP_MAX-1, 0].
- */
+#define PARTS_GENERIC_64_128_256(NAME, P) \
+ QEMU_GENERIC(P, (FloatParts256 *, parts256_##NAME), \
+ (FloatParts128 *, parts128_##NAME), parts64_##NAME)
-static FloatParts round_canonical(FloatParts p, float_status *s,
- const FloatFmt *parm)
-{
- const uint64_t frac_lsb = parm->frac_lsb;
- const uint64_t frac_lsbm1 = parm->frac_lsbm1;
- const uint64_t round_mask = parm->round_mask;
- const uint64_t roundeven_mask = parm->roundeven_mask;
- const int exp_max = parm->exp_max;
- const int frac_shift = parm->frac_shift;
- uint64_t frac, inc;
- int exp, flags = 0;
- bool overflow_norm;
+#define parts_default_nan(P, S) PARTS_GENERIC_64_128(default_nan, P)(P, S)
+#define parts_silence_nan(P, S) PARTS_GENERIC_64_128(silence_nan, P)(P, S)
- frac = p.frac;
- exp = p.exp;
+static void parts64_return_nan(FloatParts64 *a, float_status *s);
+static void parts128_return_nan(FloatParts128 *a, float_status *s);
- switch (p.cls) {
- case float_class_normal:
- switch (s->float_rounding_mode) {
- case float_round_nearest_even:
- overflow_norm = false;
- inc = ((frac & roundeven_mask) != frac_lsbm1 ? frac_lsbm1 : 0);
- break;
- case float_round_ties_away:
- overflow_norm = false;
- inc = frac_lsbm1;
- break;
- case float_round_to_zero:
- overflow_norm = true;
- inc = 0;
- break;
- case float_round_up:
- inc = p.sign ? 0 : round_mask;
- overflow_norm = p.sign;
- break;
- case float_round_down:
- inc = p.sign ? round_mask : 0;
- overflow_norm = !p.sign;
- break;
- case float_round_to_odd:
- overflow_norm = true;
- inc = frac & frac_lsb ? 0 : round_mask;
- break;
- default:
- g_assert_not_reached();
- }
+#define parts_return_nan(P, S) PARTS_GENERIC_64_128(return_nan, P)(P, S)
- exp += parm->exp_bias;
- if (likely(exp > 0)) {
- if (frac & round_mask) {
- flags |= float_flag_inexact;
- frac += inc;
- if (frac & DECOMPOSED_OVERFLOW_BIT) {
- frac >>= 1;
- exp++;
- }
- }
- frac >>= frac_shift;
-
- if (parm->arm_althp) {
- /* ARM Alt HP eschews Inf and NaN for a wider exponent. */
- if (unlikely(exp > exp_max)) {
- /* Overflow. Return the maximum normal. */
- flags = float_flag_invalid;
- exp = exp_max;
- frac = -1;
- }
- } else if (unlikely(exp >= exp_max)) {
- flags |= float_flag_overflow | float_flag_inexact;
- if (overflow_norm) {
- exp = exp_max - 1;
- frac = -1;
- } else {
- p.cls = float_class_inf;
- goto do_inf;
- }
- }
- } else if (s->flush_to_zero) {
- flags |= float_flag_output_denormal;
- p.cls = float_class_zero;
- goto do_zero;
- } else {
- bool is_tiny = s->tininess_before_rounding
- || (exp < 0)
- || !((frac + inc) & DECOMPOSED_OVERFLOW_BIT);
-
- shift64RightJamming(frac, 1 - exp, &frac);
- if (frac & round_mask) {
- /* Need to recompute round-to-even. */
- switch (s->float_rounding_mode) {
- case float_round_nearest_even:
- inc = ((frac & roundeven_mask) != frac_lsbm1
- ? frac_lsbm1 : 0);
- break;
- case float_round_to_odd:
- inc = frac & frac_lsb ? 0 : round_mask;
- break;
- default:
- break;
- }
- flags |= float_flag_inexact;
- frac += inc;
- }
+static FloatParts64 *parts64_pick_nan(FloatParts64 *a, FloatParts64 *b,
+ float_status *s);
+static FloatParts128 *parts128_pick_nan(FloatParts128 *a, FloatParts128 *b,
+ float_status *s);
- exp = (frac & DECOMPOSED_IMPLICIT_BIT ? 1 : 0);
- frac >>= frac_shift;
+#define parts_pick_nan(A, B, S) PARTS_GENERIC_64_128(pick_nan, A)(A, B, S)
- if (is_tiny && (flags & float_flag_inexact)) {
- flags |= float_flag_underflow;
- }
- if (exp == 0 && frac == 0) {
- p.cls = float_class_zero;
- }
- }
- break;
+static FloatParts64 *parts64_pick_nan_muladd(FloatParts64 *a, FloatParts64 *b,
+ FloatParts64 *c, float_status *s,
+ int ab_mask, int abc_mask);
+static FloatParts128 *parts128_pick_nan_muladd(FloatParts128 *a,
+ FloatParts128 *b,
+ FloatParts128 *c,
+ float_status *s,
+ int ab_mask, int abc_mask);
- case float_class_zero:
- do_zero:
- exp = 0;
- frac = 0;
- break;
+#define parts_pick_nan_muladd(A, B, C, S, ABM, ABCM) \
+ PARTS_GENERIC_64_128(pick_nan_muladd, A)(A, B, C, S, ABM, ABCM)
- case float_class_inf:
- do_inf:
- assert(!parm->arm_althp);
- exp = exp_max;
- frac = 0;
- break;
+static void parts64_canonicalize(FloatParts64 *p, float_status *status,
+ const FloatFmt *fmt);
+static void parts128_canonicalize(FloatParts128 *p, float_status *status,
+ const FloatFmt *fmt);
- case float_class_qnan:
- case float_class_snan:
- assert(!parm->arm_althp);
- exp = exp_max;
- frac >>= parm->frac_shift;
- break;
+#define parts_canonicalize(A, S, F) \
+ PARTS_GENERIC_64_128(canonicalize, A)(A, S, F)
- default:
- g_assert_not_reached();
- }
+static void parts64_uncanon(FloatParts64 *p, float_status *status,
+ const FloatFmt *fmt);
+static void parts128_uncanon(FloatParts128 *p, float_status *status,
+ const FloatFmt *fmt);
- float_raise(flags, s);
- p.exp = exp;
- p.frac = frac;
- return p;
-}
+#define parts_uncanon(A, S, F) \
+ PARTS_GENERIC_64_128(uncanon, A)(A, S, F)
+
+static void parts64_add_normal(FloatParts64 *a, FloatParts64 *b);
+static void parts128_add_normal(FloatParts128 *a, FloatParts128 *b);
+static void parts256_add_normal(FloatParts256 *a, FloatParts256 *b);
+
+#define parts_add_normal(A, B) \
+ PARTS_GENERIC_64_128_256(add_normal, A)(A, B)
+
+static bool parts64_sub_normal(FloatParts64 *a, FloatParts64 *b);
+static bool parts128_sub_normal(FloatParts128 *a, FloatParts128 *b);
+static bool parts256_sub_normal(FloatParts256 *a, FloatParts256 *b);
+
+#define parts_sub_normal(A, B) \
+ PARTS_GENERIC_64_128_256(sub_normal, A)(A, B)
+
+static FloatParts64 *parts64_addsub(FloatParts64 *a, FloatParts64 *b,
+ float_status *s, bool subtract);
+static FloatParts128 *parts128_addsub(FloatParts128 *a, FloatParts128 *b,
+ float_status *s, bool subtract);
+
+#define parts_addsub(A, B, S, Z) \
+ PARTS_GENERIC_64_128(addsub, A)(A, B, S, Z)
+
+static FloatParts64 *parts64_mul(FloatParts64 *a, FloatParts64 *b,
+ float_status *s);
+static FloatParts128 *parts128_mul(FloatParts128 *a, FloatParts128 *b,
+ float_status *s);
+
+#define parts_mul(A, B, S) \
+ PARTS_GENERIC_64_128(mul, A)(A, B, S)
+
+static FloatParts64 *parts64_muladd(FloatParts64 *a, FloatParts64 *b,
+ FloatParts64 *c, int flags,
+ float_status *s);
+static FloatParts128 *parts128_muladd(FloatParts128 *a, FloatParts128 *b,
+ FloatParts128 *c, int flags,
+ float_status *s);
+
+#define parts_muladd(A, B, C, Z, S) \
+ PARTS_GENERIC_64_128(muladd, A)(A, B, C, Z, S)
+
+static FloatParts64 *parts64_div(FloatParts64 *a, FloatParts64 *b,
+ float_status *s);
+static FloatParts128 *parts128_div(FloatParts128 *a, FloatParts128 *b,
+ float_status *s);
+
+#define parts_div(A, B, S) \
+ PARTS_GENERIC_64_128(div, A)(A, B, S)
+
+static bool parts64_round_to_int_normal(FloatParts64 *a, FloatRoundMode rm,
+ int scale, int frac_size);
+static bool parts128_round_to_int_normal(FloatParts128 *a, FloatRoundMode r,
+ int scale, int frac_size);
+
+#define parts_round_to_int_normal(A, R, C, F) \
+ PARTS_GENERIC_64_128(round_to_int_normal, A)(A, R, C, F)
+
+static void parts64_round_to_int(FloatParts64 *a, FloatRoundMode rm,
+ int scale, float_status *s,
+ const FloatFmt *fmt);
+static void parts128_round_to_int(FloatParts128 *a, FloatRoundMode r,
+ int scale, float_status *s,
+ const FloatFmt *fmt);
+
+#define parts_round_to_int(A, R, C, S, F) \
+ PARTS_GENERIC_64_128(round_to_int, A)(A, R, C, S, F)
+
+static int64_t parts64_float_to_sint(FloatParts64 *p, FloatRoundMode rmode,
+ int scale, int64_t min, int64_t max,
+ float_status *s);
+static int64_t parts128_float_to_sint(FloatParts128 *p, FloatRoundMode rmode,
+ int scale, int64_t min, int64_t max,
+ float_status *s);
+
+#define parts_float_to_sint(P, R, Z, MN, MX, S) \
+ PARTS_GENERIC_64_128(float_to_sint, P)(P, R, Z, MN, MX, S)
+
+/*
+ * Helper functions for softfloat-parts.c.inc, per-size operations.
+ */
-/* Explicit FloatFmt version */
-static FloatParts float16a_unpack_canonical(float16 f, float_status *s,
- const FloatFmt *params)
+#define FRAC_GENERIC_64_128(NAME, P) \
+ QEMU_GENERIC(P, (FloatParts128 *, frac128_##NAME), frac64_##NAME)
+
+#define FRAC_GENERIC_64_128_256(NAME, P) \
+ QEMU_GENERIC(P, (FloatParts256 *, frac256_##NAME), \
+ (FloatParts128 *, frac128_##NAME), frac64_##NAME)
+
+static bool frac64_add(FloatParts64 *r, FloatParts64 *a, FloatParts64 *b)
{
- return sf_canonicalize(float16_unpack_raw(f), params, s);
+ return uadd64_overflow(a->frac, b->frac, &r->frac);
}
-static FloatParts float16_unpack_canonical(float16 f, float_status *s)
+static bool frac128_add(FloatParts128 *r, FloatParts128 *a, FloatParts128 *b)
{
- return float16a_unpack_canonical(f, s, &float16_params);
+ bool c = 0;
+ r->frac_lo = uadd64_carry(a->frac_lo, b->frac_lo, &c);
+ r->frac_hi = uadd64_carry(a->frac_hi, b->frac_hi, &c);
+ return c;
}
-static FloatParts bfloat16_unpack_canonical(bfloat16 f, float_status *s)
+static bool frac256_add(FloatParts256 *r, FloatParts256 *a, FloatParts256 *b)
{
- return sf_canonicalize(bfloat16_unpack_raw(f), &bfloat16_params, s);
+ bool c = 0;
+ r->frac_lo = uadd64_carry(a->frac_lo, b->frac_lo, &c);
+ r->frac_lm = uadd64_carry(a->frac_lm, b->frac_lm, &c);
+ r->frac_hm = uadd64_carry(a->frac_hm, b->frac_hm, &c);
+ r->frac_hi = uadd64_carry(a->frac_hi, b->frac_hi, &c);
+ return c;
}
-static float16 float16a_round_pack_canonical(FloatParts p, float_status *s,
- const FloatFmt *params)
+#define frac_add(R, A, B) FRAC_GENERIC_64_128_256(add, R)(R, A, B)
+
+static bool frac64_addi(FloatParts64 *r, FloatParts64 *a, uint64_t c)
{
- return float16_pack_raw(round_canonical(p, s, params));
+ return uadd64_overflow(a->frac, c, &r->frac);
}
-static float16 float16_round_pack_canonical(FloatParts p, float_status *s)
+static bool frac128_addi(FloatParts128 *r, FloatParts128 *a, uint64_t c)
{
- return float16a_round_pack_canonical(p, s, &float16_params);
+ c = uadd64_overflow(a->frac_lo, c, &r->frac_lo);
+ return uadd64_overflow(a->frac_hi, c, &r->frac_hi);
}
-static bfloat16 bfloat16_round_pack_canonical(FloatParts p, float_status *s)
+#define frac_addi(R, A, C) FRAC_GENERIC_64_128(addi, R)(R, A, C)
+
+static void frac64_allones(FloatParts64 *a)
{
- return bfloat16_pack_raw(round_canonical(p, s, &bfloat16_params));
+ a->frac = -1;
}
-static FloatParts float32_unpack_canonical(float32 f, float_status *s)
+static void frac128_allones(FloatParts128 *a)
{
- return sf_canonicalize(float32_unpack_raw(f), &float32_params, s);
+ a->frac_hi = a->frac_lo = -1;
}
-static float32 float32_round_pack_canonical(FloatParts p, float_status *s)
+#define frac_allones(A) FRAC_GENERIC_64_128(allones, A)(A)
+
+static int frac64_cmp(FloatParts64 *a, FloatParts64 *b)
{
- return float32_pack_raw(round_canonical(p, s, &float32_params));
+ return a->frac == b->frac ? 0 : a->frac < b->frac ? -1 : 1;
}
-static FloatParts float64_unpack_canonical(float64 f, float_status *s)
+static int frac128_cmp(FloatParts128 *a, FloatParts128 *b)
{
- return sf_canonicalize(float64_unpack_raw(f), &float64_params, s);
+ uint64_t ta = a->frac_hi, tb = b->frac_hi;
+ if (ta == tb) {
+ ta = a->frac_lo, tb = b->frac_lo;
+ if (ta == tb) {
+ return 0;
+ }
+ }
+ return ta < tb ? -1 : 1;
}
-static float64 float64_round_pack_canonical(FloatParts p, float_status *s)
+#define frac_cmp(A, B) FRAC_GENERIC_64_128(cmp, A)(A, B)
+
+static void frac64_clear(FloatParts64 *a)
{
- return float64_pack_raw(round_canonical(p, s, &float64_params));
+ a->frac = 0;
}
-static FloatParts return_nan(FloatParts a, float_status *s)
+static void frac128_clear(FloatParts128 *a)
{
- switch (a.cls) {
- case float_class_snan:
- s->float_exception_flags |= float_flag_invalid;
- a = parts_silence_nan(a, s);
- /* fall through */
- case float_class_qnan:
- if (s->default_nan_mode) {
- return parts_default_nan(s);
- }
- break;
-
- default:
- g_assert_not_reached();
- }
- return a;
+ a->frac_hi = a->frac_lo = 0;
}
-static FloatParts pick_nan(FloatParts a, FloatParts b, float_status *s)
+#define frac_clear(A) FRAC_GENERIC_64_128(clear, A)(A)
+
+static bool frac64_div(FloatParts64 *a, FloatParts64 *b)
{
- if (is_snan(a.cls) || is_snan(b.cls)) {
- s->float_exception_flags |= float_flag_invalid;
- }
+ uint64_t n1, n0, r, q;
+ bool ret;
- if (s->default_nan_mode) {
- return parts_default_nan(s);
+ /*
+ * We want a 2*N / N-bit division to produce exactly an N-bit
+ * result, so that we do not lose any precision and so that we
+ * do not have to renormalize afterward. If A.frac < B.frac,
+ * then division would produce an (N-1)-bit result; shift A left
+ * by one to produce the an N-bit result, and return true to
+ * decrement the exponent to match.
+ *
+ * The udiv_qrnnd algorithm that we're using requires normalization,
+ * i.e. the msb of the denominator must be set, which is already true.
+ */
+ ret = a->frac < b->frac;
+ if (ret) {
+ n0 = a->frac;
+ n1 = 0;
} else {
- if (pickNaN(a.cls, b.cls,
- a.frac > b.frac ||
- (a.frac == b.frac && a.sign < b.sign), s)) {
- a = b;
- }
- if (is_snan(a.cls)) {
- return parts_silence_nan(a, s);
- }
+ n0 = a->frac >> 1;
+ n1 = a->frac << 63;
}
- return a;
+ q = udiv_qrnnd(&r, n0, n1, b->frac);
+
+ /* Set lsb if there is a remainder, to set inexact. */
+ a->frac = q | (r != 0);
+
+ return ret;
}
-static FloatParts pick_nan_muladd(FloatParts a, FloatParts b, FloatParts c,
- bool inf_zero, float_status *s)
+static bool frac128_div(FloatParts128 *a, FloatParts128 *b)
{
- int which;
+ uint64_t q0, q1, a0, a1, b0, b1;
+ uint64_t r0, r1, r2, r3, t0, t1, t2, t3;
+ bool ret = false;
+
+ a0 = a->frac_hi, a1 = a->frac_lo;
+ b0 = b->frac_hi, b1 = b->frac_lo;
- if (is_snan(a.cls) || is_snan(b.cls) || is_snan(c.cls)) {
- s->float_exception_flags |= float_flag_invalid;
+ ret = lt128(a0, a1, b0, b1);
+ if (!ret) {
+ a1 = shr_double(a0, a1, 1);
+ a0 = a0 >> 1;
}
- which = pickNaNMulAdd(a.cls, b.cls, c.cls, inf_zero, s);
+ /* Use 128/64 -> 64 division as estimate for 192/128 -> 128 division. */
+ q0 = estimateDiv128To64(a0, a1, b0);
- if (s->default_nan_mode) {
- /* Note that this check is after pickNaNMulAdd so that function
- * has an opportunity to set the Invalid flag.
- */
- which = 3;
+ /*
+ * Estimate is high because B1 was not included (unless B1 == 0).
+ * Reduce quotient and increase remainder until remainder is non-negative.
+ * This loop will execute 0 to 2 times.
+ */
+ mul128By64To192(b0, b1, q0, &t0, &t1, &t2);
+ sub192(a0, a1, 0, t0, t1, t2, &r0, &r1, &r2);
+ while (r0 != 0) {
+ q0--;
+ add192(r0, r1, r2, 0, b0, b1, &r0, &r1, &r2);
}
- switch (which) {
- case 0:
- break;
- case 1:
- a = b;
- break;
- case 2:
- a = c;
- break;
- case 3:
- return parts_default_nan(s);
- default:
- g_assert_not_reached();
+ /* Repeat using the remainder, producing a second word of quotient. */
+ q1 = estimateDiv128To64(r1, r2, b0);
+ mul128By64To192(b0, b1, q1, &t1, &t2, &t3);
+ sub192(r1, r2, 0, t1, t2, t3, &r1, &r2, &r3);
+ while (r1 != 0) {
+ q1--;
+ add192(r1, r2, r3, 0, b0, b1, &r1, &r2, &r3);
}
- if (is_snan(a.cls)) {
- return parts_silence_nan(a, s);
- }
- return a;
+ /* Any remainder indicates inexact; set sticky bit. */
+ q1 |= (r2 | r3) != 0;
+
+ a->frac_hi = q0;
+ a->frac_lo = q1;
+ return ret;
}
-/*
- * Returns the result of adding or subtracting the values of the
- * floating-point values `a' and `b'. The operation is performed
- * according to the IEC/IEEE Standard for Binary Floating-Point
- * Arithmetic.
- */
+#define frac_div(A, B) FRAC_GENERIC_64_128(div, A)(A, B)
-static FloatParts addsub_floats(FloatParts a, FloatParts b, bool subtract,
- float_status *s)
+static bool frac64_eqz(FloatParts64 *a)
{
- bool a_sign = a.sign;
- bool b_sign = b.sign ^ subtract;
-
- if (a_sign != b_sign) {
- /* Subtraction */
-
- if (a.cls == float_class_normal && b.cls == float_class_normal) {
- if (a.exp > b.exp || (a.exp == b.exp && a.frac >= b.frac)) {
- shift64RightJamming(b.frac, a.exp - b.exp, &b.frac);
- a.frac = a.frac - b.frac;
- } else {
- shift64RightJamming(a.frac, b.exp - a.exp, &a.frac);
- a.frac = b.frac - a.frac;
- a.exp = b.exp;
- a_sign ^= 1;
- }
+ return a->frac == 0;
+}
- if (a.frac == 0) {
- a.cls = float_class_zero;
- a.sign = s->float_rounding_mode == float_round_down;
- } else {
- int shift = clz64(a.frac) - 1;
- a.frac = a.frac << shift;
- a.exp = a.exp - shift;
- a.sign = a_sign;
- }
- return a;
- }
- if (is_nan(a.cls) || is_nan(b.cls)) {
- return pick_nan(a, b, s);
- }
- if (a.cls == float_class_inf) {
- if (b.cls == float_class_inf) {
- float_raise(float_flag_invalid, s);
- return parts_default_nan(s);
- }
- return a;
- }
- if (a.cls == float_class_zero && b.cls == float_class_zero) {
- a.sign = s->float_rounding_mode == float_round_down;
- return a;
- }
- if (a.cls == float_class_zero || b.cls == float_class_inf) {
- b.sign = a_sign ^ 1;
- return b;
- }
- if (b.cls == float_class_zero) {
- return a;
- }
- } else {
- /* Addition */
- if (a.cls == float_class_normal && b.cls == float_class_normal) {
- if (a.exp > b.exp) {
- shift64RightJamming(b.frac, a.exp - b.exp, &b.frac);
- } else if (a.exp < b.exp) {
- shift64RightJamming(a.frac, b.exp - a.exp, &a.frac);
- a.exp = b.exp;
- }
- a.frac += b.frac;
- if (a.frac & DECOMPOSED_OVERFLOW_BIT) {
- shift64RightJamming(a.frac, 1, &a.frac);
- a.exp += 1;
- }
- return a;
- }
- if (is_nan(a.cls) || is_nan(b.cls)) {
- return pick_nan(a, b, s);
- }
- if (a.cls == float_class_inf || b.cls == float_class_zero) {
- return a;
- }
- if (b.cls == float_class_inf || a.cls == float_class_zero) {
- b.sign = b_sign;
- return b;
- }
- }
- g_assert_not_reached();
+static bool frac128_eqz(FloatParts128 *a)
+{
+ return (a->frac_hi | a->frac_lo) == 0;
}
-/*
- * Returns the result of adding or subtracting the floating-point
- * values `a' and `b'. The operation is performed according to the
- * IEC/IEEE Standard for Binary Floating-Point Arithmetic.
- */
+#define frac_eqz(A) FRAC_GENERIC_64_128(eqz, A)(A)
-float16 QEMU_FLATTEN float16_add(float16 a, float16 b, float_status *status)
+static void frac64_mulw(FloatParts128 *r, FloatParts64 *a, FloatParts64 *b)
{
- FloatParts pa = float16_unpack_canonical(a, status);
- FloatParts pb = float16_unpack_canonical(b, status);
- FloatParts pr = addsub_floats(pa, pb, false, status);
-
- return float16_round_pack_canonical(pr, status);
+ mulu64(&r->frac_lo, &r->frac_hi, a->frac, b->frac);
}
-float16 QEMU_FLATTEN float16_sub(float16 a, float16 b, float_status *status)
+static void frac128_mulw(FloatParts256 *r, FloatParts128 *a, FloatParts128 *b)
{
- FloatParts pa = float16_unpack_canonical(a, status);
- FloatParts pb = float16_unpack_canonical(b, status);
- FloatParts pr = addsub_floats(pa, pb, true, status);
-
- return float16_round_pack_canonical(pr, status);
+ mul128To256(a->frac_hi, a->frac_lo, b->frac_hi, b->frac_lo,
+ &r->frac_hi, &r->frac_hm, &r->frac_lm, &r->frac_lo);
}
-static float32 QEMU_SOFTFLOAT_ATTR
-soft_f32_addsub(float32 a, float32 b, bool subtract, float_status *status)
-{
- FloatParts pa = float32_unpack_canonical(a, status);
- FloatParts pb = float32_unpack_canonical(b, status);
- FloatParts pr = addsub_floats(pa, pb, subtract, status);
+#define frac_mulw(R, A, B) FRAC_GENERIC_64_128(mulw, A)(R, A, B)
- return float32_round_pack_canonical(pr, status);
+static void frac64_neg(FloatParts64 *a)
+{
+ a->frac = -a->frac;
}
-static inline float32 soft_f32_add(float32 a, float32 b, float_status *status)
+static void frac128_neg(FloatParts128 *a)
{
- return soft_f32_addsub(a, b, false, status);
+ bool c = 0;
+ a->frac_lo = usub64_borrow(0, a->frac_lo, &c);
+ a->frac_hi = usub64_borrow(0, a->frac_hi, &c);
}
-static inline float32 soft_f32_sub(float32 a, float32 b, float_status *status)
+static void frac256_neg(FloatParts256 *a)
{
- return soft_f32_addsub(a, b, true, status);
+ bool c = 0;
+ a->frac_lo = usub64_borrow(0, a->frac_lo, &c);
+ a->frac_lm = usub64_borrow(0, a->frac_lm, &c);
+ a->frac_hm = usub64_borrow(0, a->frac_hm, &c);
+ a->frac_hi = usub64_borrow(0, a->frac_hi, &c);
}
-static float64 QEMU_SOFTFLOAT_ATTR
-soft_f64_addsub(float64 a, float64 b, bool subtract, float_status *status)
+#define frac_neg(A) FRAC_GENERIC_64_128_256(neg, A)(A)
+
+static int frac64_normalize(FloatParts64 *a)
{
- FloatParts pa = float64_unpack_canonical(a, status);
- FloatParts pb = float64_unpack_canonical(b, status);
- FloatParts pr = addsub_floats(pa, pb, subtract, status);
+ if (a->frac) {
+ int shift = clz64(a->frac);
+ a->frac <<= shift;
+ return shift;
+ }
+ return 64;
+}
- return float64_round_pack_canonical(pr, status);
+static int frac128_normalize(FloatParts128 *a)
+{
+ if (a->frac_hi) {
+ int shl = clz64(a->frac_hi);
+ a->frac_hi = shl_double(a->frac_hi, a->frac_lo, shl);
+ a->frac_lo <<= shl;
+ return shl;
+ } else if (a->frac_lo) {
+ int shl = clz64(a->frac_lo);
+ a->frac_hi = a->frac_lo << shl;
+ a->frac_lo = 0;
+ return shl + 64;
+ }
+ return 128;
}
-static inline float64 soft_f64_add(float64 a, float64 b, float_status *status)
+static int frac256_normalize(FloatParts256 *a)
{
- return soft_f64_addsub(a, b, false, status);
+ uint64_t a0 = a->frac_hi, a1 = a->frac_hm;
+ uint64_t a2 = a->frac_lm, a3 = a->frac_lo;
+ int ret, shl;
+
+ if (likely(a0)) {
+ shl = clz64(a0);
+ if (shl == 0) {
+ return 0;
+ }
+ ret = shl;
+ } else {
+ if (a1) {
+ ret = 64;
+ a0 = a1, a1 = a2, a2 = a3, a3 = 0;
+ } else if (a2) {
+ ret = 128;
+ a0 = a2, a1 = a3, a2 = 0, a3 = 0;
+ } else if (a3) {
+ ret = 192;
+ a0 = a3, a1 = 0, a2 = 0, a3 = 0;
+ } else {
+ ret = 256;
+ a0 = 0, a1 = 0, a2 = 0, a3 = 0;
+ goto done;
+ }
+ shl = clz64(a0);
+ if (shl == 0) {
+ goto done;
+ }
+ ret += shl;
+ }
+
+ a0 = shl_double(a0, a1, shl);
+ a1 = shl_double(a1, a2, shl);
+ a2 = shl_double(a2, a3, shl);
+ a3 <<= shl;
+
+ done:
+ a->frac_hi = a0;
+ a->frac_hm = a1;
+ a->frac_lm = a2;
+ a->frac_lo = a3;
+ return ret;
}
-static inline float64 soft_f64_sub(float64 a, float64 b, float_status *status)
+#define frac_normalize(A) FRAC_GENERIC_64_128_256(normalize, A)(A)
+
+static void frac64_shl(FloatParts64 *a, int c)
{
- return soft_f64_addsub(a, b, true, status);
+ a->frac <<= c;
}
-static float hard_f32_add(float a, float b)
+static void frac128_shl(FloatParts128 *a, int c)
{
- return a + b;
+ uint64_t a0 = a->frac_hi, a1 = a->frac_lo;
+
+ if (c & 64) {
+ a0 = a1, a1 = 0;
+ }
+
+ c &= 63;
+ if (c) {
+ a0 = shl_double(a0, a1, c);
+ a1 = a1 << c;
+ }
+
+ a->frac_hi = a0;
+ a->frac_lo = a1;
}
-static float hard_f32_sub(float a, float b)
+#define frac_shl(A, C) FRAC_GENERIC_64_128(shl, A)(A, C)
+
+static void frac64_shr(FloatParts64 *a, int c)
{
- return a - b;
+ a->frac >>= c;
}
-static double hard_f64_add(double a, double b)
+static void frac128_shr(FloatParts128 *a, int c)
{
- return a + b;
+ uint64_t a0 = a->frac_hi, a1 = a->frac_lo;
+
+ if (c & 64) {
+ a1 = a0, a0 = 0;
+ }
+
+ c &= 63;
+ if (c) {
+ a1 = shr_double(a0, a1, c);
+ a0 = a0 >> c;
+ }
+
+ a->frac_hi = a0;
+ a->frac_lo = a1;
}
-static double hard_f64_sub(double a, double b)
+#define frac_shr(A, C) FRAC_GENERIC_64_128(shr, A)(A, C)
+
+static void frac64_shrjam(FloatParts64 *a, int c)
{
- return a - b;
+ uint64_t a0 = a->frac;
+
+ if (likely(c != 0)) {
+ if (likely(c < 64)) {
+ a0 = (a0 >> c) | (shr_double(a0, 0, c) != 0);
+ } else {
+ a0 = a0 != 0;
+ }
+ a->frac = a0;
+ }
}
-static bool f32_addsubmul_post(union_float32 a, union_float32 b)
+static void frac128_shrjam(FloatParts128 *a, int c)
{
- if (QEMU_HARDFLOAT_2F32_USE_FP) {
- return !(fpclassify(a.h) == FP_ZERO && fpclassify(b.h) == FP_ZERO);
+ uint64_t a0 = a->frac_hi, a1 = a->frac_lo;
+ uint64_t sticky = 0;
+
+ if (unlikely(c == 0)) {
+ return;
+ } else if (likely(c < 64)) {
+ /* nothing */
+ } else if (likely(c < 128)) {
+ sticky = a1;
+ a1 = a0;
+ a0 = 0;
+ c &= 63;
+ if (c == 0) {
+ goto done;
+ }
+ } else {
+ sticky = a0 | a1;
+ a0 = a1 = 0;
+ goto done;
}
- return !(float32_is_zero(a.s) && float32_is_zero(b.s));
+
+ sticky |= shr_double(a1, 0, c);
+ a1 = shr_double(a0, a1, c);
+ a0 = a0 >> c;
+
+ done:
+ a->frac_lo = a1 | (sticky != 0);
+ a->frac_hi = a0;
}
-static bool f64_addsubmul_post(union_float64 a, union_float64 b)
+static void frac256_shrjam(FloatParts256 *a, int c)
{
- if (QEMU_HARDFLOAT_2F64_USE_FP) {
- return !(fpclassify(a.h) == FP_ZERO && fpclassify(b.h) == FP_ZERO);
+ uint64_t a0 = a->frac_hi, a1 = a->frac_hm;
+ uint64_t a2 = a->frac_lm, a3 = a->frac_lo;
+ uint64_t sticky = 0;
+
+ if (unlikely(c == 0)) {
+ return;
+ } else if (likely(c < 64)) {
+ /* nothing */
+ } else if (likely(c < 256)) {
+ if (unlikely(c & 128)) {
+ sticky |= a2 | a3;
+ a3 = a1, a2 = a0, a1 = 0, a0 = 0;
+ }
+ if (unlikely(c & 64)) {
+ sticky |= a3;
+ a3 = a2, a2 = a1, a1 = a0, a0 = 0;
+ }
+ c &= 63;
+ if (c == 0) {
+ goto done;
+ }
} else {
- return !(float64_is_zero(a.s) && float64_is_zero(b.s));
+ sticky = a0 | a1 | a2 | a3;
+ a0 = a1 = a2 = a3 = 0;
+ goto done;
}
+
+ sticky |= shr_double(a3, 0, c);
+ a3 = shr_double(a2, a3, c);
+ a2 = shr_double(a1, a2, c);
+ a1 = shr_double(a0, a1, c);
+ a0 = a0 >> c;
+
+ done:
+ a->frac_lo = a3 | (sticky != 0);
+ a->frac_lm = a2;
+ a->frac_hm = a1;
+ a->frac_hi = a0;
}
-static float32 float32_addsub(float32 a, float32 b, float_status *s,
- hard_f32_op2_fn hard, soft_f32_op2_fn soft)
+#define frac_shrjam(A, C) FRAC_GENERIC_64_128_256(shrjam, A)(A, C)
+
+static bool frac64_sub(FloatParts64 *r, FloatParts64 *a, FloatParts64 *b)
{
- return float32_gen2(a, b, s, hard, soft,
- f32_is_zon2, f32_addsubmul_post);
+ return usub64_overflow(a->frac, b->frac, &r->frac);
}
-static float64 float64_addsub(float64 a, float64 b, float_status *s,
- hard_f64_op2_fn hard, soft_f64_op2_fn soft)
+static bool frac128_sub(FloatParts128 *r, FloatParts128 *a, FloatParts128 *b)
{
- return float64_gen2(a, b, s, hard, soft,
- f64_is_zon2, f64_addsubmul_post);
+ bool c = 0;
+ r->frac_lo = usub64_borrow(a->frac_lo, b->frac_lo, &c);
+ r->frac_hi = usub64_borrow(a->frac_hi, b->frac_hi, &c);
+ return c;
}
-float32 QEMU_FLATTEN
-float32_add(float32 a, float32 b, float_status *s)
+static bool frac256_sub(FloatParts256 *r, FloatParts256 *a, FloatParts256 *b)
{
- return float32_addsub(a, b, s, hard_f32_add, soft_f32_add);
+ bool c = 0;
+ r->frac_lo = usub64_borrow(a->frac_lo, b->frac_lo, &c);
+ r->frac_lm = usub64_borrow(a->frac_lm, b->frac_lm, &c);
+ r->frac_hm = usub64_borrow(a->frac_hm, b->frac_hm, &c);
+ r->frac_hi = usub64_borrow(a->frac_hi, b->frac_hi, &c);
+ return c;
}
-float32 QEMU_FLATTEN
-float32_sub(float32 a, float32 b, float_status *s)
+#define frac_sub(R, A, B) FRAC_GENERIC_64_128_256(sub, R)(R, A, B)
+
+static void frac64_truncjam(FloatParts64 *r, FloatParts128 *a)
{
- return float32_addsub(a, b, s, hard_f32_sub, soft_f32_sub);
+ r->frac = a->frac_hi | (a->frac_lo != 0);
}
-float64 QEMU_FLATTEN
-float64_add(float64 a, float64 b, float_status *s)
+static void frac128_truncjam(FloatParts128 *r, FloatParts256 *a)
{
- return float64_addsub(a, b, s, hard_f64_add, soft_f64_add);
+ r->frac_hi = a->frac_hi;
+ r->frac_lo = a->frac_hm | ((a->frac_lm | a->frac_lo) != 0);
}
-float64 QEMU_FLATTEN
-float64_sub(float64 a, float64 b, float_status *s)
+#define frac_truncjam(R, A) FRAC_GENERIC_64_128(truncjam, R)(R, A)
+
+static void frac64_widen(FloatParts128 *r, FloatParts64 *a)
{
- return float64_addsub(a, b, s, hard_f64_sub, soft_f64_sub);
+ r->frac_hi = a->frac;
+ r->frac_lo = 0;
+}
+
+static void frac128_widen(FloatParts256 *r, FloatParts128 *a)
+{
+ r->frac_hi = a->frac_hi;
+ r->frac_hm = a->frac_lo;
+ r->frac_lm = 0;
+ r->frac_lo = 0;
}
+#define frac_widen(A, B) FRAC_GENERIC_64_128(widen, B)(A, B)
+
+#define partsN(NAME) glue(glue(glue(parts,N),_),NAME)
+#define FloatPartsN glue(FloatParts,N)
+#define FloatPartsW glue(FloatParts,W)
+
+#define N 64
+#define W 128
+
+#include "softfloat-parts-addsub.c.inc"
+#include "softfloat-parts.c.inc"
+
+#undef N
+#undef W
+#define N 128
+#define W 256
+
+#include "softfloat-parts-addsub.c.inc"
+#include "softfloat-parts.c.inc"
+
+#undef N
+#undef W
+#define N 256
+
+#include "softfloat-parts-addsub.c.inc"
+
+#undef N
+#undef W
+#undef partsN
+#undef FloatPartsN
+#undef FloatPartsW
+
/*
- * Returns the result of adding or subtracting the bfloat16
- * values `a' and `b'.
+ * Pack/unpack routines with a specific FloatFmt.
*/
-bfloat16 QEMU_FLATTEN bfloat16_add(bfloat16 a, bfloat16 b, float_status *status)
+
+static void float16a_unpack_canonical(FloatParts64 *p, float16 f,
+ float_status *s, const FloatFmt *params)
{
- FloatParts pa = bfloat16_unpack_canonical(a, status);
- FloatParts pb = bfloat16_unpack_canonical(b, status);
- FloatParts pr = addsub_floats(pa, pb, false, status);
+ float16_unpack_raw(p, f);
+ parts_canonicalize(p, s, params);
+}
- return bfloat16_round_pack_canonical(pr, status);
+static void float16_unpack_canonical(FloatParts64 *p, float16 f,
+ float_status *s)
+{
+ float16a_unpack_canonical(p, f, s, &float16_params);
}
-bfloat16 QEMU_FLATTEN bfloat16_sub(bfloat16 a, bfloat16 b, float_status *status)
+static void bfloat16_unpack_canonical(FloatParts64 *p, bfloat16 f,
+ float_status *s)
{
- FloatParts pa = bfloat16_unpack_canonical(a, status);
- FloatParts pb = bfloat16_unpack_canonical(b, status);
- FloatParts pr = addsub_floats(pa, pb, true, status);
+ bfloat16_unpack_raw(p, f);
+ parts_canonicalize(p, s, &bfloat16_params);
+}
- return bfloat16_round_pack_canonical(pr, status);
+static float16 float16a_round_pack_canonical(FloatParts64 *p,
+ float_status *s,
+ const FloatFmt *params)
+{
+ parts_uncanon(p, s, params);
+ return float16_pack_raw(p);
}
-/*
- * Returns the result of multiplying the floating-point values `a' and
- * `b'. The operation is performed according to the IEC/IEEE Standard
- * for Binary Floating-Point Arithmetic.
- */
+static float16 float16_round_pack_canonical(FloatParts64 *p,
+ float_status *s)
+{
+ return float16a_round_pack_canonical(p, s, &float16_params);
+}
-static FloatParts mul_floats(FloatParts a, FloatParts b, float_status *s)
+static bfloat16 bfloat16_round_pack_canonical(FloatParts64 *p,
+ float_status *s)
{
- bool sign = a.sign ^ b.sign;
+ parts_uncanon(p, s, &bfloat16_params);
+ return bfloat16_pack_raw(p);
+}
- if (a.cls == float_class_normal && b.cls == float_class_normal) {
- uint64_t hi, lo;
- int exp = a.exp + b.exp;
+static void float32_unpack_canonical(FloatParts64 *p, float32 f,
+ float_status *s)
+{
+ float32_unpack_raw(p, f);
+ parts_canonicalize(p, s, &float32_params);
+}
- mul64To128(a.frac, b.frac, &hi, &lo);
- shift128RightJamming(hi, lo, DECOMPOSED_BINARY_POINT, &hi, &lo);
- if (lo & DECOMPOSED_OVERFLOW_BIT) {
- shift64RightJamming(lo, 1, &lo);
- exp += 1;
- }
+static float32 float32_round_pack_canonical(FloatParts64 *p,
+ float_status *s)
+{
+ parts_uncanon(p, s, &float32_params);
+ return float32_pack_raw(p);
+}
- /* Re-use a */
- a.exp = exp;
- a.sign = sign;
- a.frac = lo;
- return a;
- }
- /* handle all the NaN cases */
- if (is_nan(a.cls) || is_nan(b.cls)) {
- return pick_nan(a, b, s);
- }
- /* Inf * Zero == NaN */
- if ((a.cls == float_class_inf && b.cls == float_class_zero) ||
- (a.cls == float_class_zero && b.cls == float_class_inf)) {
- s->float_exception_flags |= float_flag_invalid;
- return parts_default_nan(s);
- }
- /* Multiply by 0 or Inf */
- if (a.cls == float_class_inf || a.cls == float_class_zero) {
- a.sign = sign;
- return a;
- }
- if (b.cls == float_class_inf || b.cls == float_class_zero) {
- b.sign = sign;
- return b;
- }
- g_assert_not_reached();
+static void float64_unpack_canonical(FloatParts64 *p, float64 f,
+ float_status *s)
+{
+ float64_unpack_raw(p, f);
+ parts_canonicalize(p, s, &float64_params);
}
-float16 QEMU_FLATTEN float16_mul(float16 a, float16 b, float_status *status)
+static float64 float64_round_pack_canonical(FloatParts64 *p,
+ float_status *s)
+{
+ parts_uncanon(p, s, &float64_params);
+ return float64_pack_raw(p);
+}
+
+static void float128_unpack_canonical(FloatParts128 *p, float128 f,
+ float_status *s)
+{
+ float128_unpack_raw(p, f);
+ parts_canonicalize(p, s, &float128_params);
+}
+
+static float128 float128_round_pack_canonical(FloatParts128 *p,
+ float_status *s)
+{
+ parts_uncanon(p, s, &float128_params);
+ return float128_pack_raw(p);
+}
+
+/*
+ * Addition and subtraction
+ */
+
+static float16 QEMU_FLATTEN
+float16_addsub(float16 a, float16 b, float_status *status, bool subtract)
{
- FloatParts pa = float16_unpack_canonical(a, status);
- FloatParts pb = float16_unpack_canonical(b, status);
- FloatParts pr = mul_floats(pa, pb, status);
+ FloatParts64 pa, pb, *pr;
+
+ float16_unpack_canonical(&pa, a, status);
+ float16_unpack_canonical(&pb, b, status);
+ pr = parts_addsub(&pa, &pb, status, subtract);
return float16_round_pack_canonical(pr, status);
}
+float16 float16_add(float16 a, float16 b, float_status *status)
+{
+ return float16_addsub(a, b, status, false);
+}
+
+float16 float16_sub(float16 a, float16 b, float_status *status)
+{
+ return float16_addsub(a, b, status, true);
+}
+
static float32 QEMU_SOFTFLOAT_ATTR
-soft_f32_mul(float32 a, float32 b, float_status *status)
+soft_f32_addsub(float32 a, float32 b, float_status *status, bool subtract)
{
- FloatParts pa = float32_unpack_canonical(a, status);
- FloatParts pb = float32_unpack_canonical(b, status);
- FloatParts pr = mul_floats(pa, pb, status);
+ FloatParts64 pa, pb, *pr;
+
+ float32_unpack_canonical(&pa, a, status);
+ float32_unpack_canonical(&pb, b, status);
+ pr = parts_addsub(&pa, &pb, status, subtract);
return float32_round_pack_canonical(pr, status);
}
+static float32 soft_f32_add(float32 a, float32 b, float_status *status)
+{
+ return soft_f32_addsub(a, b, status, false);
+}
+
+static float32 soft_f32_sub(float32 a, float32 b, float_status *status)
+{
+ return soft_f32_addsub(a, b, status, true);
+}
+
static float64 QEMU_SOFTFLOAT_ATTR
-soft_f64_mul(float64 a, float64 b, float_status *status)
+soft_f64_addsub(float64 a, float64 b, float_status *status, bool subtract)
{
- FloatParts pa = float64_unpack_canonical(a, status);
- FloatParts pb = float64_unpack_canonical(b, status);
- FloatParts pr = mul_floats(pa, pb, status);
+ FloatParts64 pa, pb, *pr;
+
+ float64_unpack_canonical(&pa, a, status);
+ float64_unpack_canonical(&pb, b, status);
+ pr = parts_addsub(&pa, &pb, status, subtract);
return float64_round_pack_canonical(pr, status);
}
-static float hard_f32_mul(float a, float b)
+static float64 soft_f64_add(float64 a, float64 b, float_status *status)
{
- return a * b;
+ return soft_f64_addsub(a, b, status, false);
}
-static double hard_f64_mul(double a, double b)
+static float64 soft_f64_sub(float64 a, float64 b, float_status *status)
{
- return a * b;
+ return soft_f64_addsub(a, b, status, true);
}
-float32 QEMU_FLATTEN
-float32_mul(float32 a, float32 b, float_status *s)
+static float hard_f32_add(float a, float b)
{
- return float32_gen2(a, b, s, hard_f32_mul, soft_f32_mul,
- f32_is_zon2, f32_addsubmul_post);
+ return a + b;
}
-float64 QEMU_FLATTEN
-float64_mul(float64 a, float64 b, float_status *s)
+static float hard_f32_sub(float a, float b)
{
- return float64_gen2(a, b, s, hard_f64_mul, soft_f64_mul,
- f64_is_zon2, f64_addsubmul_post);
+ return a - b;
}
-/*
- * Returns the result of multiplying the bfloat16
- * values `a' and `b'.
- */
-
-bfloat16 QEMU_FLATTEN bfloat16_mul(bfloat16 a, bfloat16 b, float_status *status)
+static double hard_f64_add(double a, double b)
{
- FloatParts pa = bfloat16_unpack_canonical(a, status);
- FloatParts pb = bfloat16_unpack_canonical(b, status);
- FloatParts pr = mul_floats(pa, pb, status);
-
- return bfloat16_round_pack_canonical(pr, status);
+ return a + b;
}
-/*
- * Returns the result of multiplying the floating-point values `a' and
- * `b' then adding 'c', with no intermediate rounding step after the
- * multiplication. The operation is performed according to the
- * IEC/IEEE Standard for Binary Floating-Point Arithmetic 754-2008.
- * The flags argument allows the caller to select negation of the
- * addend, the intermediate product, or the final result. (The
- * difference between this and having the caller do a separate
- * negation is that negating externally will flip the sign bit on
- * NaNs.)
- */
-
-static FloatParts muladd_floats(FloatParts a, FloatParts b, FloatParts c,
- int flags, float_status *s)
-{
- bool inf_zero = ((1 << a.cls) | (1 << b.cls)) ==
- ((1 << float_class_inf) | (1 << float_class_zero));
- bool p_sign;
- bool sign_flip = flags & float_muladd_negate_result;
- FloatClass p_class;
- uint64_t hi, lo;
- int p_exp;
-
- /* It is implementation-defined whether the cases of (0,inf,qnan)
- * and (inf,0,qnan) raise InvalidOperation or not (and what QNaN
- * they return if they do), so we have to hand this information
- * off to the target-specific pick-a-NaN routine.
- */
- if (is_nan(a.cls) || is_nan(b.cls) || is_nan(c.cls)) {
- return pick_nan_muladd(a, b, c, inf_zero, s);
- }
+static double hard_f64_sub(double a, double b)
+{
+ return a - b;
+}
- if (inf_zero) {
- s->float_exception_flags |= float_flag_invalid;
- return parts_default_nan(s);
+static bool f32_addsubmul_post(union_float32 a, union_float32 b)
+{
+ if (QEMU_HARDFLOAT_2F32_USE_FP) {
+ return !(fpclassify(a.h) == FP_ZERO && fpclassify(b.h) == FP_ZERO);
}
+ return !(float32_is_zero(a.s) && float32_is_zero(b.s));
+}
- if (flags & float_muladd_negate_c) {
- c.sign ^= 1;
+static bool f64_addsubmul_post(union_float64 a, union_float64 b)
+{
+ if (QEMU_HARDFLOAT_2F64_USE_FP) {
+ return !(fpclassify(a.h) == FP_ZERO && fpclassify(b.h) == FP_ZERO);
+ } else {
+ return !(float64_is_zero(a.s) && float64_is_zero(b.s));
}
+}
- p_sign = a.sign ^ b.sign;
-
- if (flags & float_muladd_negate_product) {
- p_sign ^= 1;
- }
+static float32 float32_addsub(float32 a, float32 b, float_status *s,
+ hard_f32_op2_fn hard, soft_f32_op2_fn soft)
+{
+ return float32_gen2(a, b, s, hard, soft,
+ f32_is_zon2, f32_addsubmul_post);
+}
- if (a.cls == float_class_inf || b.cls == float_class_inf) {
- p_class = float_class_inf;
- } else if (a.cls == float_class_zero || b.cls == float_class_zero) {
- p_class = float_class_zero;
- } else {
- p_class = float_class_normal;
- }
+static float64 float64_addsub(float64 a, float64 b, float_status *s,
+ hard_f64_op2_fn hard, soft_f64_op2_fn soft)
+{
+ return float64_gen2(a, b, s, hard, soft,
+ f64_is_zon2, f64_addsubmul_post);
+}
- if (c.cls == float_class_inf) {
- if (p_class == float_class_inf && p_sign != c.sign) {
- s->float_exception_flags |= float_flag_invalid;
- return parts_default_nan(s);
- } else {
- a.cls = float_class_inf;
- a.sign = c.sign ^ sign_flip;
- return a;
- }
- }
+float32 QEMU_FLATTEN
+float32_add(float32 a, float32 b, float_status *s)
+{
+ return float32_addsub(a, b, s, hard_f32_add, soft_f32_add);
+}
- if (p_class == float_class_inf) {
- a.cls = float_class_inf;
- a.sign = p_sign ^ sign_flip;
- return a;
- }
+float32 QEMU_FLATTEN
+float32_sub(float32 a, float32 b, float_status *s)
+{
+ return float32_addsub(a, b, s, hard_f32_sub, soft_f32_sub);
+}
- if (p_class == float_class_zero) {
- if (c.cls == float_class_zero) {
- if (p_sign != c.sign) {
- p_sign = s->float_rounding_mode == float_round_down;
- }
- c.sign = p_sign;
- } else if (flags & float_muladd_halve_result) {
- c.exp -= 1;
- }
- c.sign ^= sign_flip;
- return c;
- }
+float64 QEMU_FLATTEN
+float64_add(float64 a, float64 b, float_status *s)
+{
+ return float64_addsub(a, b, s, hard_f64_add, soft_f64_add);
+}
- /* a & b should be normals now... */
- assert(a.cls == float_class_normal &&
- b.cls == float_class_normal);
+float64 QEMU_FLATTEN
+float64_sub(float64 a, float64 b, float_status *s)
+{
+ return float64_addsub(a, b, s, hard_f64_sub, soft_f64_sub);
+}
- p_exp = a.exp + b.exp;
+static bfloat16 QEMU_FLATTEN
+bfloat16_addsub(bfloat16 a, bfloat16 b, float_status *status, bool subtract)
+{
+ FloatParts64 pa, pb, *pr;
- /* Multiply of 2 62-bit numbers produces a (2*62) == 124-bit
- * result.
- */
- mul64To128(a.frac, b.frac, &hi, &lo);
- /* binary point now at bit 124 */
+ bfloat16_unpack_canonical(&pa, a, status);
+ bfloat16_unpack_canonical(&pb, b, status);
+ pr = parts_addsub(&pa, &pb, status, subtract);
- /* check for overflow */
- if (hi & (1ULL << (DECOMPOSED_BINARY_POINT * 2 + 1 - 64))) {
- shift128RightJamming(hi, lo, 1, &hi, &lo);
- p_exp += 1;
- }
+ return bfloat16_round_pack_canonical(pr, status);
+}
- /* + add/sub */
- if (c.cls == float_class_zero) {
- /* move binary point back to 62 */
- shift128RightJamming(hi, lo, DECOMPOSED_BINARY_POINT, &hi, &lo);
- } else {
- int exp_diff = p_exp - c.exp;
- if (p_sign == c.sign) {
- /* Addition */
- if (exp_diff <= 0) {
- shift128RightJamming(hi, lo,
- DECOMPOSED_BINARY_POINT - exp_diff,
- &hi, &lo);
- lo += c.frac;
- p_exp = c.exp;
- } else {
- uint64_t c_hi, c_lo;
- /* shift c to the same binary point as the product (124) */
- c_hi = c.frac >> 2;
- c_lo = 0;
- shift128RightJamming(c_hi, c_lo,
- exp_diff,
- &c_hi, &c_lo);
- add128(hi, lo, c_hi, c_lo, &hi, &lo);
- /* move binary point back to 62 */
- shift128RightJamming(hi, lo, DECOMPOSED_BINARY_POINT, &hi, &lo);
- }
+bfloat16 bfloat16_add(bfloat16 a, bfloat16 b, float_status *status)
+{
+ return bfloat16_addsub(a, b, status, false);
+}
- if (lo & DECOMPOSED_OVERFLOW_BIT) {
- shift64RightJamming(lo, 1, &lo);
- p_exp += 1;
- }
+bfloat16 bfloat16_sub(bfloat16 a, bfloat16 b, float_status *status)
+{
+ return bfloat16_addsub(a, b, status, true);
+}
- } else {
- /* Subtraction */
- uint64_t c_hi, c_lo;
- /* make C binary point match product at bit 124 */
- c_hi = c.frac >> 2;
- c_lo = 0;
-
- if (exp_diff <= 0) {
- shift128RightJamming(hi, lo, -exp_diff, &hi, &lo);
- if (exp_diff == 0
- &&
- (hi > c_hi || (hi == c_hi && lo >= c_lo))) {
- sub128(hi, lo, c_hi, c_lo, &hi, &lo);
- } else {
- sub128(c_hi, c_lo, hi, lo, &hi, &lo);
- p_sign ^= 1;
- p_exp = c.exp;
- }
- } else {
- shift128RightJamming(c_hi, c_lo,
- exp_diff,
- &c_hi, &c_lo);
- sub128(hi, lo, c_hi, c_lo, &hi, &lo);
- }
+static float128 QEMU_FLATTEN
+float128_addsub(float128 a, float128 b, float_status *status, bool subtract)
+{
+ FloatParts128 pa, pb, *pr;
- if (hi == 0 && lo == 0) {
- a.cls = float_class_zero;
- a.sign = s->float_rounding_mode == float_round_down;
- a.sign ^= sign_flip;
- return a;
- } else {
- int shift;
- if (hi != 0) {
- shift = clz64(hi);
- } else {
- shift = clz64(lo) + 64;
- }
- /* Normalizing to a binary point of 124 is the
- correct adjust for the exponent. However since we're
- shifting, we might as well put the binary point back
- at 62 where we really want it. Therefore shift as
- if we're leaving 1 bit at the top of the word, but
- adjust the exponent as if we're leaving 3 bits. */
- shift -= 1;
- if (shift >= 64) {
- lo = lo << (shift - 64);
- } else {
- hi = (hi << shift) | (lo >> (64 - shift));
- lo = hi | ((lo << shift) != 0);
- }
- p_exp -= shift - 2;
- }
- }
- }
+ float128_unpack_canonical(&pa, a, status);
+ float128_unpack_canonical(&pb, b, status);
+ pr = parts_addsub(&pa, &pb, status, subtract);
- if (flags & float_muladd_halve_result) {
- p_exp -= 1;
- }
+ return float128_round_pack_canonical(pr, status);
+}
- /* finally prepare our result */
- a.cls = float_class_normal;
- a.sign = p_sign ^ sign_flip;
- a.exp = p_exp;
- a.frac = lo;
+float128 float128_add(float128 a, float128 b, float_status *status)
+{
+ return float128_addsub(a, b, status, false);
+}
- return a;
+float128 float128_sub(float128 a, float128 b, float_status *status)
+{
+ return float128_addsub(a, b, status, true);
}
-float16 QEMU_FLATTEN float16_muladd(float16 a, float16 b, float16 c,
- int flags, float_status *status)
+/*
+ * Multiplication
+ */
+
+float16 QEMU_FLATTEN float16_mul(float16 a, float16 b, float_status *status)
{
- FloatParts pa = float16_unpack_canonical(a, status);
- FloatParts pb = float16_unpack_canonical(b, status);
- FloatParts pc = float16_unpack_canonical(c, status);
- FloatParts pr = muladd_floats(pa, pb, pc, flags, status);
+ FloatParts64 pa, pb, *pr;
+
+ float16_unpack_canonical(&pa, a, status);
+ float16_unpack_canonical(&pb, b, status);
+ pr = parts_mul(&pa, &pb, status);
return float16_round_pack_canonical(pr, status);
}
static float32 QEMU_SOFTFLOAT_ATTR
-soft_f32_muladd(float32 a, float32 b, float32 c, int flags,
- float_status *status)
+soft_f32_mul(float32 a, float32 b, float_status *status)
{
- FloatParts pa = float32_unpack_canonical(a, status);
- FloatParts pb = float32_unpack_canonical(b, status);
- FloatParts pc = float32_unpack_canonical(c, status);
- FloatParts pr = muladd_floats(pa, pb, pc, flags, status);
+ FloatParts64 pa, pb, *pr;
+
+ float32_unpack_canonical(&pa, a, status);
+ float32_unpack_canonical(&pb, b, status);
+ pr = parts_mul(&pa, &pb, status);
return float32_round_pack_canonical(pr, status);
}
static float64 QEMU_SOFTFLOAT_ATTR
-soft_f64_muladd(float64 a, float64 b, float64 c, int flags,
- float_status *status)
+soft_f64_mul(float64 a, float64 b, float_status *status)
{
- FloatParts pa = float64_unpack_canonical(a, status);
- FloatParts pb = float64_unpack_canonical(b, status);
- FloatParts pc = float64_unpack_canonical(c, status);
- FloatParts pr = muladd_floats(pa, pb, pc, flags, status);
+ FloatParts64 pa, pb, *pr;
+
+ float64_unpack_canonical(&pa, a, status);
+ float64_unpack_canonical(&pb, b, status);
+ pr = parts_mul(&pa, &pb, status);
return float64_round_pack_canonical(pr, status);
}
-static bool force_soft_fma;
+static float hard_f32_mul(float a, float b)
+{
+ return a * b;
+}
+
+static double hard_f64_mul(double a, double b)
+{
+ return a * b;
+}
float32 QEMU_FLATTEN
-float32_muladd(float32 xa, float32 xb, float32 xc, int flags, float_status *s)
+float32_mul(float32 a, float32 b, float_status *s)
{
- union_float32 ua, ub, uc, ur;
+ return float32_gen2(a, b, s, hard_f32_mul, soft_f32_mul,
+ f32_is_zon2, f32_addsubmul_post);
+}
- ua.s = xa;
- ub.s = xb;
- uc.s = xc;
+float64 QEMU_FLATTEN
+float64_mul(float64 a, float64 b, float_status *s)
+{
+ return float64_gen2(a, b, s, hard_f64_mul, soft_f64_mul,
+ f64_is_zon2, f64_addsubmul_post);
+}
- if (unlikely(!can_use_fpu(s))) {
- goto soft;
- }
- if (unlikely(flags & float_muladd_halve_result)) {
- goto soft;
- }
+bfloat16 QEMU_FLATTEN
+bfloat16_mul(bfloat16 a, bfloat16 b, float_status *status)
+{
+ FloatParts64 pa, pb, *pr;
- float32_input_flush3(&ua.s, &ub.s, &uc.s, s);
- if (unlikely(!f32_is_zon3(ua, ub, uc))) {
- goto soft;
- }
+ bfloat16_unpack_canonical(&pa, a, status);
+ bfloat16_unpack_canonical(&pb, b, status);
+ pr = parts_mul(&pa, &pb, status);
- if (unlikely(force_soft_fma)) {
- goto soft;
+ return bfloat16_round_pack_canonical(pr, status);
+}
+
+float128 QEMU_FLATTEN
+float128_mul(float128 a, float128 b, float_status *status)
+{
+ FloatParts128 pa, pb, *pr;
+
+ float128_unpack_canonical(&pa, a, status);
+ float128_unpack_canonical(&pb, b, status);
+ pr = parts_mul(&pa, &pb, status);
+
+ return float128_round_pack_canonical(pr, status);
+}
+
+/*
+ * Fused multiply-add
+ */
+
+float16 QEMU_FLATTEN float16_muladd(float16 a, float16 b, float16 c,
+ int flags, float_status *status)
+{
+ FloatParts64 pa, pb, pc, *pr;
+
+ float16_unpack_canonical(&pa, a, status);
+ float16_unpack_canonical(&pb, b, status);
+ float16_unpack_canonical(&pc, c, status);
+ pr = parts_muladd(&pa, &pb, &pc, flags, status);
+
+ return float16_round_pack_canonical(pr, status);
+}
+
+static float32 QEMU_SOFTFLOAT_ATTR
+soft_f32_muladd(float32 a, float32 b, float32 c, int flags,
+ float_status *status)
+{
+ FloatParts64 pa, pb, pc, *pr;
+
+ float32_unpack_canonical(&pa, a, status);
+ float32_unpack_canonical(&pb, b, status);
+ float32_unpack_canonical(&pc, c, status);
+ pr = parts_muladd(&pa, &pb, &pc, flags, status);
+
+ return float32_round_pack_canonical(pr, status);
+}
+
+static float64 QEMU_SOFTFLOAT_ATTR
+soft_f64_muladd(float64 a, float64 b, float64 c, int flags,
+ float_status *status)
+{
+ FloatParts64 pa, pb, pc, *pr;
+
+ float64_unpack_canonical(&pa, a, status);
+ float64_unpack_canonical(&pb, b, status);
+ float64_unpack_canonical(&pc, c, status);
+ pr = parts_muladd(&pa, &pb, &pc, flags, status);
+
+ return float64_round_pack_canonical(pr, status);
+}
+
+static bool force_soft_fma;
+
+float32 QEMU_FLATTEN
+float32_muladd(float32 xa, float32 xb, float32 xc, int flags, float_status *s)
+{
+ union_float32 ua, ub, uc, ur;
+
+ ua.s = xa;
+ ub.s = xb;
+ uc.s = xc;
+
+ if (unlikely(!can_use_fpu(s))) {
+ goto soft;
+ }
+ if (unlikely(flags & float_muladd_halve_result)) {
+ goto soft;
+ }
+
+ float32_input_flush3(&ua.s, &ub.s, &uc.s, s);
+ if (unlikely(!f32_is_zon3(ua, ub, uc))) {
+ goto soft;
+ }
+
+ if (unlikely(force_soft_fma)) {
+ goto soft;
}
/*
ur.h = fmaf(ua.h, ub.h, uc.h);
if (unlikely(f32_is_inf(ur))) {
- s->float_exception_flags |= float_flag_overflow;
+ float_raise(float_flag_overflow, s);
} else if (unlikely(fabsf(ur.h) <= FLT_MIN)) {
ua = ua_orig;
uc = uc_orig;
ur.h = fma(ua.h, ub.h, uc.h);
if (unlikely(f64_is_inf(ur))) {
- s->float_exception_flags |= float_flag_overflow;
+ float_raise(float_flag_overflow, s);
} else if (unlikely(fabs(ur.h) <= FLT_MIN)) {
ua = ua_orig;
uc = uc_orig;
return soft_f64_muladd(ua.s, ub.s, uc.s, flags, s);
}
-/*
- * Returns the result of multiplying the bfloat16 values `a'
- * and `b' then adding 'c', with no intermediate rounding step after the
- * multiplication.
- */
-
bfloat16 QEMU_FLATTEN bfloat16_muladd(bfloat16 a, bfloat16 b, bfloat16 c,
int flags, float_status *status)
{
- FloatParts pa = bfloat16_unpack_canonical(a, status);
- FloatParts pb = bfloat16_unpack_canonical(b, status);
- FloatParts pc = bfloat16_unpack_canonical(c, status);
- FloatParts pr = muladd_floats(pa, pb, pc, flags, status);
+ FloatParts64 pa, pb, pc, *pr;
+
+ bfloat16_unpack_canonical(&pa, a, status);
+ bfloat16_unpack_canonical(&pb, b, status);
+ bfloat16_unpack_canonical(&pc, c, status);
+ pr = parts_muladd(&pa, &pb, &pc, flags, status);
return bfloat16_round_pack_canonical(pr, status);
}
-/*
- * Returns the result of dividing the floating-point value `a' by the
- * corresponding value `b'. The operation is performed according to
- * the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
- */
-
-static FloatParts div_floats(FloatParts a, FloatParts b, float_status *s)
+float128 QEMU_FLATTEN float128_muladd(float128 a, float128 b, float128 c,
+ int flags, float_status *status)
{
- bool sign = a.sign ^ b.sign;
-
- if (a.cls == float_class_normal && b.cls == float_class_normal) {
- uint64_t n0, n1, q, r;
- int exp = a.exp - b.exp;
+ FloatParts128 pa, pb, pc, *pr;
- /*
- * We want a 2*N / N-bit division to produce exactly an N-bit
- * result, so that we do not lose any precision and so that we
- * do not have to renormalize afterward. If A.frac < B.frac,
- * then division would produce an (N-1)-bit result; shift A left
- * by one to produce the an N-bit result, and decrement the
- * exponent to match.
- *
- * The udiv_qrnnd algorithm that we're using requires normalization,
- * i.e. the msb of the denominator must be set. Since we know that
- * DECOMPOSED_BINARY_POINT is msb-1, the inputs must be shifted left
- * by one (more), and the remainder must be shifted right by one.
- */
- if (a.frac < b.frac) {
- exp -= 1;
- shift128Left(0, a.frac, DECOMPOSED_BINARY_POINT + 2, &n1, &n0);
- } else {
- shift128Left(0, a.frac, DECOMPOSED_BINARY_POINT + 1, &n1, &n0);
- }
- q = udiv_qrnnd(&r, n1, n0, b.frac << 1);
+ float128_unpack_canonical(&pa, a, status);
+ float128_unpack_canonical(&pb, b, status);
+ float128_unpack_canonical(&pc, c, status);
+ pr = parts_muladd(&pa, &pb, &pc, flags, status);
- /*
- * Set lsb if there is a remainder, to set inexact.
- * As mentioned above, to find the actual value of the remainder we
- * would need to shift right, but (1) we are only concerned about
- * non-zero-ness, and (2) the remainder will always be even because
- * both inputs to the division primitive are even.
- */
- a.frac = q | (r != 0);
- a.sign = sign;
- a.exp = exp;
- return a;
- }
- /* handle all the NaN cases */
- if (is_nan(a.cls) || is_nan(b.cls)) {
- return pick_nan(a, b, s);
- }
- /* 0/0 or Inf/Inf */
- if (a.cls == b.cls
- &&
- (a.cls == float_class_inf || a.cls == float_class_zero)) {
- s->float_exception_flags |= float_flag_invalid;
- return parts_default_nan(s);
- }
- /* Inf / x or 0 / x */
- if (a.cls == float_class_inf || a.cls == float_class_zero) {
- a.sign = sign;
- return a;
- }
- /* Div 0 => Inf */
- if (b.cls == float_class_zero) {
- s->float_exception_flags |= float_flag_divbyzero;
- a.cls = float_class_inf;
- a.sign = sign;
- return a;
- }
- /* Div by Inf */
- if (b.cls == float_class_inf) {
- a.cls = float_class_zero;
- a.sign = sign;
- return a;
- }
- g_assert_not_reached();
+ return float128_round_pack_canonical(pr, status);
}
+/*
+ * Division
+ */
+
float16 float16_div(float16 a, float16 b, float_status *status)
{
- FloatParts pa = float16_unpack_canonical(a, status);
- FloatParts pb = float16_unpack_canonical(b, status);
- FloatParts pr = div_floats(pa, pb, status);
+ FloatParts64 pa, pb, *pr;
+
+ float16_unpack_canonical(&pa, a, status);
+ float16_unpack_canonical(&pb, b, status);
+ pr = parts_div(&pa, &pb, status);
return float16_round_pack_canonical(pr, status);
}
static float32 QEMU_SOFTFLOAT_ATTR
soft_f32_div(float32 a, float32 b, float_status *status)
{
- FloatParts pa = float32_unpack_canonical(a, status);
- FloatParts pb = float32_unpack_canonical(b, status);
- FloatParts pr = div_floats(pa, pb, status);
+ FloatParts64 pa, pb, *pr;
+
+ float32_unpack_canonical(&pa, a, status);
+ float32_unpack_canonical(&pb, b, status);
+ pr = parts_div(&pa, &pb, status);
return float32_round_pack_canonical(pr, status);
}
static float64 QEMU_SOFTFLOAT_ATTR
soft_f64_div(float64 a, float64 b, float_status *status)
{
- FloatParts pa = float64_unpack_canonical(a, status);
- FloatParts pb = float64_unpack_canonical(b, status);
- FloatParts pr = div_floats(pa, pb, status);
+ FloatParts64 pa, pb, *pr;
+
+ float64_unpack_canonical(&pa, a, status);
+ float64_unpack_canonical(&pb, b, status);
+ pr = parts_div(&pa, &pb, status);
return float64_round_pack_canonical(pr, status);
}
f64_div_pre, f64_div_post);
}
-/*
- * Returns the result of dividing the bfloat16
- * value `a' by the corresponding value `b'.
- */
-
-bfloat16 bfloat16_div(bfloat16 a, bfloat16 b, float_status *status)
+bfloat16 QEMU_FLATTEN
+bfloat16_div(bfloat16 a, bfloat16 b, float_status *status)
{
- FloatParts pa = bfloat16_unpack_canonical(a, status);
- FloatParts pb = bfloat16_unpack_canonical(b, status);
- FloatParts pr = div_floats(pa, pb, status);
+ FloatParts64 pa, pb, *pr;
+
+ bfloat16_unpack_canonical(&pa, a, status);
+ bfloat16_unpack_canonical(&pb, b, status);
+ pr = parts_div(&pa, &pb, status);
return bfloat16_round_pack_canonical(pr, status);
}
+float128 QEMU_FLATTEN
+float128_div(float128 a, float128 b, float_status *status)
+{
+ FloatParts128 pa, pb, *pr;
+
+ float128_unpack_canonical(&pa, a, status);
+ float128_unpack_canonical(&pb, b, status);
+ pr = parts_div(&pa, &pb, status);
+
+ return float128_round_pack_canonical(pr, status);
+}
+
/*
* Float to Float conversions
*
* conversion is performed according to the IEC/IEEE Standard for
* Binary Floating-Point Arithmetic.
*
- * The float_to_float helper only needs to take care of raising
- * invalid exceptions and handling the conversion on NaNs.
+ * Usually this only needs to take care of raising invalid exceptions
+ * and handling the conversion on NaNs.
*/
-static FloatParts float_to_float(FloatParts a, const FloatFmt *dstf,
- float_status *s)
+static void parts_float_to_ahp(FloatParts64 *a, float_status *s)
{
- if (dstf->arm_althp) {
- switch (a.cls) {
- case float_class_qnan:
- case float_class_snan:
- /* There is no NaN in the destination format. Raise Invalid
- * and return a zero with the sign of the input NaN.
- */
- s->float_exception_flags |= float_flag_invalid;
- a.cls = float_class_zero;
- a.frac = 0;
- a.exp = 0;
- break;
+ switch (a->cls) {
+ case float_class_qnan:
+ case float_class_snan:
+ /*
+ * There is no NaN in the destination format. Raise Invalid
+ * and return a zero with the sign of the input NaN.
+ */
+ float_raise(float_flag_invalid, s);
+ a->cls = float_class_zero;
+ break;
- case float_class_inf:
- /* There is no Inf in the destination format. Raise Invalid
- * and return the maximum normal with the correct sign.
- */
- s->float_exception_flags |= float_flag_invalid;
- a.cls = float_class_normal;
- a.exp = dstf->exp_max;
- a.frac = ((1ull << dstf->frac_size) - 1) << dstf->frac_shift;
- break;
+ case float_class_inf:
+ /*
+ * There is no Inf in the destination format. Raise Invalid
+ * and return the maximum normal with the correct sign.
+ */
+ float_raise(float_flag_invalid, s);
+ a->cls = float_class_normal;
+ a->exp = float16_params_ahp.exp_max;
+ a->frac = MAKE_64BIT_MASK(float16_params_ahp.frac_shift,
+ float16_params_ahp.frac_size + 1);
+ break;
- default:
- break;
- }
- } else if (is_nan(a.cls)) {
- if (is_snan(a.cls)) {
- s->float_exception_flags |= float_flag_invalid;
- a = parts_silence_nan(a, s);
- }
- if (s->default_nan_mode) {
- return parts_default_nan(s);
- }
+ case float_class_normal:
+ case float_class_zero:
+ break;
+
+ default:
+ g_assert_not_reached();
+ }
+}
+
+static void parts64_float_to_float(FloatParts64 *a, float_status *s)
+{
+ if (is_nan(a->cls)) {
+ parts_return_nan(a, s);
+ }
+}
+
+static void parts128_float_to_float(FloatParts128 *a, float_status *s)
+{
+ if (is_nan(a->cls)) {
+ parts_return_nan(a, s);
+ }
+}
+
+#define parts_float_to_float(P, S) \
+ PARTS_GENERIC_64_128(float_to_float, P)(P, S)
+
+static void parts_float_to_float_narrow(FloatParts64 *a, FloatParts128 *b,
+ float_status *s)
+{
+ a->cls = b->cls;
+ a->sign = b->sign;
+ a->exp = b->exp;
+
+ if (a->cls == float_class_normal) {
+ frac_truncjam(a, b);
+ } else if (is_nan(a->cls)) {
+ /* Discard the low bits of the NaN. */
+ a->frac = b->frac_hi;
+ parts_return_nan(a, s);
+ }
+}
+
+static void parts_float_to_float_widen(FloatParts128 *a, FloatParts64 *b,
+ float_status *s)
+{
+ a->cls = b->cls;
+ a->sign = b->sign;
+ a->exp = b->exp;
+ frac_widen(a, b);
+
+ if (is_nan(a->cls)) {
+ parts_return_nan(a, s);
}
- return a;
}
float32 float16_to_float32(float16 a, bool ieee, float_status *s)
{
const FloatFmt *fmt16 = ieee ? &float16_params : &float16_params_ahp;
- FloatParts p = float16a_unpack_canonical(a, s, fmt16);
- FloatParts pr = float_to_float(p, &float32_params, s);
- return float32_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ float16a_unpack_canonical(&p, a, s, fmt16);
+ parts_float_to_float(&p, s);
+ return float32_round_pack_canonical(&p, s);
}
float64 float16_to_float64(float16 a, bool ieee, float_status *s)
{
const FloatFmt *fmt16 = ieee ? &float16_params : &float16_params_ahp;
- FloatParts p = float16a_unpack_canonical(a, s, fmt16);
- FloatParts pr = float_to_float(p, &float64_params, s);
- return float64_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ float16a_unpack_canonical(&p, a, s, fmt16);
+ parts_float_to_float(&p, s);
+ return float64_round_pack_canonical(&p, s);
}
float16 float32_to_float16(float32 a, bool ieee, float_status *s)
{
- const FloatFmt *fmt16 = ieee ? &float16_params : &float16_params_ahp;
- FloatParts p = float32_unpack_canonical(a, s);
- FloatParts pr = float_to_float(p, fmt16, s);
- return float16a_round_pack_canonical(pr, s, fmt16);
+ FloatParts64 p;
+ const FloatFmt *fmt;
+
+ float32_unpack_canonical(&p, a, s);
+ if (ieee) {
+ parts_float_to_float(&p, s);
+ fmt = &float16_params;
+ } else {
+ parts_float_to_ahp(&p, s);
+ fmt = &float16_params_ahp;
+ }
+ return float16a_round_pack_canonical(&p, s, fmt);
}
static float64 QEMU_SOFTFLOAT_ATTR
soft_float32_to_float64(float32 a, float_status *s)
{
- FloatParts p = float32_unpack_canonical(a, s);
- FloatParts pr = float_to_float(p, &float64_params, s);
- return float64_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ float32_unpack_canonical(&p, a, s);
+ parts_float_to_float(&p, s);
+ return float64_round_pack_canonical(&p, s);
}
float64 float32_to_float64(float32 a, float_status *s)
float16 float64_to_float16(float64 a, bool ieee, float_status *s)
{
- const FloatFmt *fmt16 = ieee ? &float16_params : &float16_params_ahp;
- FloatParts p = float64_unpack_canonical(a, s);
- FloatParts pr = float_to_float(p, fmt16, s);
- return float16a_round_pack_canonical(pr, s, fmt16);
+ FloatParts64 p;
+ const FloatFmt *fmt;
+
+ float64_unpack_canonical(&p, a, s);
+ if (ieee) {
+ parts_float_to_float(&p, s);
+ fmt = &float16_params;
+ } else {
+ parts_float_to_ahp(&p, s);
+ fmt = &float16_params_ahp;
+ }
+ return float16a_round_pack_canonical(&p, s, fmt);
}
float32 float64_to_float32(float64 a, float_status *s)
{
- FloatParts p = float64_unpack_canonical(a, s);
- FloatParts pr = float_to_float(p, &float32_params, s);
- return float32_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ float64_unpack_canonical(&p, a, s);
+ parts_float_to_float(&p, s);
+ return float32_round_pack_canonical(&p, s);
}
float32 bfloat16_to_float32(bfloat16 a, float_status *s)
{
- FloatParts p = bfloat16_unpack_canonical(a, s);
- FloatParts pr = float_to_float(p, &float32_params, s);
- return float32_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ bfloat16_unpack_canonical(&p, a, s);
+ parts_float_to_float(&p, s);
+ return float32_round_pack_canonical(&p, s);
}
float64 bfloat16_to_float64(bfloat16 a, float_status *s)
{
- FloatParts p = bfloat16_unpack_canonical(a, s);
- FloatParts pr = float_to_float(p, &float64_params, s);
- return float64_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ bfloat16_unpack_canonical(&p, a, s);
+ parts_float_to_float(&p, s);
+ return float64_round_pack_canonical(&p, s);
}
bfloat16 float32_to_bfloat16(float32 a, float_status *s)
{
- FloatParts p = float32_unpack_canonical(a, s);
- FloatParts pr = float_to_float(p, &bfloat16_params, s);
- return bfloat16_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ float32_unpack_canonical(&p, a, s);
+ parts_float_to_float(&p, s);
+ return bfloat16_round_pack_canonical(&p, s);
}
bfloat16 float64_to_bfloat16(float64 a, float_status *s)
{
- FloatParts p = float64_unpack_canonical(a, s);
- FloatParts pr = float_to_float(p, &bfloat16_params, s);
- return bfloat16_round_pack_canonical(pr, s);
-}
+ FloatParts64 p;
-/*
- * Rounds the floating-point value `a' to an integer, and returns the
- * result as a floating-point value. The operation is performed
- * according to the IEC/IEEE Standard for Binary Floating-Point
- * Arithmetic.
- */
+ float64_unpack_canonical(&p, a, s);
+ parts_float_to_float(&p, s);
+ return bfloat16_round_pack_canonical(&p, s);
+}
-static FloatParts round_to_int(FloatParts a, FloatRoundMode rmode,
- int scale, float_status *s)
+float32 float128_to_float32(float128 a, float_status *s)
{
- switch (a.cls) {
- case float_class_qnan:
- case float_class_snan:
- return return_nan(a, s);
+ FloatParts64 p64;
+ FloatParts128 p128;
- case float_class_zero:
- case float_class_inf:
- /* already "integral" */
- break;
+ float128_unpack_canonical(&p128, a, s);
+ parts_float_to_float_narrow(&p64, &p128, s);
+ return float32_round_pack_canonical(&p64, s);
+}
- case float_class_normal:
- scale = MIN(MAX(scale, -0x10000), 0x10000);
- a.exp += scale;
+float64 float128_to_float64(float128 a, float_status *s)
+{
+ FloatParts64 p64;
+ FloatParts128 p128;
- if (a.exp >= DECOMPOSED_BINARY_POINT) {
- /* already integral */
- break;
- }
- if (a.exp < 0) {
- bool one;
- /* all fractional */
- s->float_exception_flags |= float_flag_inexact;
- switch (rmode) {
- case float_round_nearest_even:
- one = a.exp == -1 && a.frac > DECOMPOSED_IMPLICIT_BIT;
- break;
- case float_round_ties_away:
- one = a.exp == -1 && a.frac >= DECOMPOSED_IMPLICIT_BIT;
- break;
- case float_round_to_zero:
- one = false;
- break;
- case float_round_up:
- one = !a.sign;
- break;
- case float_round_down:
- one = a.sign;
- break;
- case float_round_to_odd:
- one = true;
- break;
- default:
- g_assert_not_reached();
- }
+ float128_unpack_canonical(&p128, a, s);
+ parts_float_to_float_narrow(&p64, &p128, s);
+ return float64_round_pack_canonical(&p64, s);
+}
- if (one) {
- a.frac = DECOMPOSED_IMPLICIT_BIT;
- a.exp = 0;
- } else {
- a.cls = float_class_zero;
- }
- } else {
- uint64_t frac_lsb = DECOMPOSED_IMPLICIT_BIT >> a.exp;
- uint64_t frac_lsbm1 = frac_lsb >> 1;
- uint64_t rnd_even_mask = (frac_lsb - 1) | frac_lsb;
- uint64_t rnd_mask = rnd_even_mask >> 1;
- uint64_t inc;
+float128 float32_to_float128(float32 a, float_status *s)
+{
+ FloatParts64 p64;
+ FloatParts128 p128;
- switch (rmode) {
- case float_round_nearest_even:
- inc = ((a.frac & rnd_even_mask) != frac_lsbm1 ? frac_lsbm1 : 0);
- break;
- case float_round_ties_away:
- inc = frac_lsbm1;
- break;
- case float_round_to_zero:
- inc = 0;
- break;
- case float_round_up:
- inc = a.sign ? 0 : rnd_mask;
- break;
- case float_round_down:
- inc = a.sign ? rnd_mask : 0;
- break;
- case float_round_to_odd:
- inc = a.frac & frac_lsb ? 0 : rnd_mask;
- break;
- default:
- g_assert_not_reached();
- }
+ float32_unpack_canonical(&p64, a, s);
+ parts_float_to_float_widen(&p128, &p64, s);
+ return float128_round_pack_canonical(&p128, s);
+}
- if (a.frac & rnd_mask) {
- s->float_exception_flags |= float_flag_inexact;
- a.frac += inc;
- a.frac &= ~rnd_mask;
- if (a.frac & DECOMPOSED_OVERFLOW_BIT) {
- a.frac >>= 1;
- a.exp++;
- }
- }
- }
- break;
- default:
- g_assert_not_reached();
- }
- return a;
+float128 float64_to_float128(float64 a, float_status *s)
+{
+ FloatParts64 p64;
+ FloatParts128 p128;
+
+ float64_unpack_canonical(&p64, a, s);
+ parts_float_to_float_widen(&p128, &p64, s);
+ return float128_round_pack_canonical(&p128, s);
}
+/*
+ * Round to integral value
+ */
+
float16 float16_round_to_int(float16 a, float_status *s)
{
- FloatParts pa = float16_unpack_canonical(a, s);
- FloatParts pr = round_to_int(pa, s->float_rounding_mode, 0, s);
- return float16_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ float16_unpack_canonical(&p, a, s);
+ parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float16_params);
+ return float16_round_pack_canonical(&p, s);
}
float32 float32_round_to_int(float32 a, float_status *s)
{
- FloatParts pa = float32_unpack_canonical(a, s);
- FloatParts pr = round_to_int(pa, s->float_rounding_mode, 0, s);
- return float32_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ float32_unpack_canonical(&p, a, s);
+ parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float32_params);
+ return float32_round_pack_canonical(&p, s);
}
float64 float64_round_to_int(float64 a, float_status *s)
{
- FloatParts pa = float64_unpack_canonical(a, s);
- FloatParts pr = round_to_int(pa, s->float_rounding_mode, 0, s);
- return float64_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ float64_unpack_canonical(&p, a, s);
+ parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float64_params);
+ return float64_round_pack_canonical(&p, s);
+}
+
+bfloat16 bfloat16_round_to_int(bfloat16 a, float_status *s)
+{
+ FloatParts64 p;
+
+ bfloat16_unpack_canonical(&p, a, s);
+ parts_round_to_int(&p, s->float_rounding_mode, 0, s, &bfloat16_params);
+ return bfloat16_round_pack_canonical(&p, s);
+}
+
+float128 float128_round_to_int(float128 a, float_status *s)
+{
+ FloatParts128 p;
+
+ float128_unpack_canonical(&p, a, s);
+ parts_round_to_int(&p, s->float_rounding_mode, 0, s, &float128_params);
+ return float128_round_pack_canonical(&p, s);
}
/*
- * Rounds the bfloat16 value `a' to an integer, and returns the
- * result as a bfloat16 value.
+ * Floating-point to signed integer conversions
*/
-bfloat16 bfloat16_round_to_int(bfloat16 a, float_status *s)
+int8_t float16_to_int8_scalbn(float16 a, FloatRoundMode rmode, int scale,
+ float_status *s)
{
- FloatParts pa = bfloat16_unpack_canonical(a, s);
- FloatParts pr = round_to_int(pa, s->float_rounding_mode, 0, s);
- return bfloat16_round_pack_canonical(pr, s);
+ FloatParts64 p;
+
+ float16_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT8_MIN, INT8_MAX, s);
}
-/*
- * Returns the result of converting the floating-point value `a' to
- * the two's complement integer format. The conversion is performed
- * according to the IEC/IEEE Standard for Binary Floating-Point
- * Arithmetic---which means in particular that the conversion is
- * rounded according to the current rounding mode. If `a' is a NaN,
- * the largest positive integer is returned. Otherwise, if the
- * conversion overflows, the largest integer with the same sign as `a'
- * is returned.
-*/
+int16_t float16_to_int16_scalbn(float16 a, FloatRoundMode rmode, int scale,
+ float_status *s)
+{
+ FloatParts64 p;
-static int64_t round_to_int_and_pack(FloatParts in, FloatRoundMode rmode,
- int scale, int64_t min, int64_t max,
- float_status *s)
+ float16_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
+}
+
+int32_t float16_to_int32_scalbn(float16 a, FloatRoundMode rmode, int scale,
+ float_status *s)
{
- uint64_t r;
- int orig_flags = get_float_exception_flags(s);
- FloatParts p = round_to_int(in, rmode, scale, s);
+ FloatParts64 p;
- switch (p.cls) {
- case float_class_snan:
- case float_class_qnan:
- s->float_exception_flags = orig_flags | float_flag_invalid;
- return max;
- case float_class_inf:
- s->float_exception_flags = orig_flags | float_flag_invalid;
- return p.sign ? min : max;
- case float_class_zero:
- return 0;
- case float_class_normal:
- if (p.exp < DECOMPOSED_BINARY_POINT) {
- r = p.frac >> (DECOMPOSED_BINARY_POINT - p.exp);
- } else if (p.exp - DECOMPOSED_BINARY_POINT < 2) {
- r = p.frac << (p.exp - DECOMPOSED_BINARY_POINT);
- } else {
- r = UINT64_MAX;
- }
- if (p.sign) {
- if (r <= -(uint64_t) min) {
- return -r;
- } else {
- s->float_exception_flags = orig_flags | float_flag_invalid;
- return min;
- }
- } else {
- if (r <= max) {
- return r;
- } else {
- s->float_exception_flags = orig_flags | float_flag_invalid;
- return max;
- }
- }
- default:
- g_assert_not_reached();
- }
+ float16_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
}
-int8_t float16_to_int8_scalbn(float16 a, FloatRoundMode rmode, int scale,
- float_status *s)
+int64_t float16_to_int64_scalbn(float16 a, FloatRoundMode rmode, int scale,
+ float_status *s)
{
- return round_to_int_and_pack(float16_unpack_canonical(a, s),
- rmode, scale, INT8_MIN, INT8_MAX, s);
+ FloatParts64 p;
+
+ float16_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
}
-int16_t float16_to_int16_scalbn(float16 a, FloatRoundMode rmode, int scale,
+int16_t float32_to_int16_scalbn(float32 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_int_and_pack(float16_unpack_canonical(a, s),
- rmode, scale, INT16_MIN, INT16_MAX, s);
+ FloatParts64 p;
+
+ float32_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
}
-int32_t float16_to_int32_scalbn(float16 a, FloatRoundMode rmode, int scale,
+int32_t float32_to_int32_scalbn(float32 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_int_and_pack(float16_unpack_canonical(a, s),
- rmode, scale, INT32_MIN, INT32_MAX, s);
+ FloatParts64 p;
+
+ float32_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
}
-int64_t float16_to_int64_scalbn(float16 a, FloatRoundMode rmode, int scale,
+int64_t float32_to_int64_scalbn(float32 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_int_and_pack(float16_unpack_canonical(a, s),
- rmode, scale, INT64_MIN, INT64_MAX, s);
+ FloatParts64 p;
+
+ float32_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
}
-int16_t float32_to_int16_scalbn(float32 a, FloatRoundMode rmode, int scale,
+int16_t float64_to_int16_scalbn(float64 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_int_and_pack(float32_unpack_canonical(a, s),
- rmode, scale, INT16_MIN, INT16_MAX, s);
+ FloatParts64 p;
+
+ float64_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
}
-int32_t float32_to_int32_scalbn(float32 a, FloatRoundMode rmode, int scale,
+int32_t float64_to_int32_scalbn(float64 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_int_and_pack(float32_unpack_canonical(a, s),
- rmode, scale, INT32_MIN, INT32_MAX, s);
+ FloatParts64 p;
+
+ float64_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
}
-int64_t float32_to_int64_scalbn(float32 a, FloatRoundMode rmode, int scale,
+int64_t float64_to_int64_scalbn(float64 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_int_and_pack(float32_unpack_canonical(a, s),
- rmode, scale, INT64_MIN, INT64_MAX, s);
+ FloatParts64 p;
+
+ float64_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
}
-int16_t float64_to_int16_scalbn(float64 a, FloatRoundMode rmode, int scale,
- float_status *s)
+int16_t bfloat16_to_int16_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
+ float_status *s)
{
- return round_to_int_and_pack(float64_unpack_canonical(a, s),
- rmode, scale, INT16_MIN, INT16_MAX, s);
+ FloatParts64 p;
+
+ bfloat16_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT16_MIN, INT16_MAX, s);
}
-int32_t float64_to_int32_scalbn(float64 a, FloatRoundMode rmode, int scale,
- float_status *s)
+int32_t bfloat16_to_int32_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
+ float_status *s)
{
- return round_to_int_and_pack(float64_unpack_canonical(a, s),
- rmode, scale, INT32_MIN, INT32_MAX, s);
+ FloatParts64 p;
+
+ bfloat16_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
}
-int64_t float64_to_int64_scalbn(float64 a, FloatRoundMode rmode, int scale,
- float_status *s)
+int64_t bfloat16_to_int64_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
+ float_status *s)
+{
+ FloatParts64 p;
+
+ bfloat16_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
+}
+
+static int32_t float128_to_int32_scalbn(float128 a, FloatRoundMode rmode,
+ int scale, float_status *s)
+{
+ FloatParts128 p;
+
+ float128_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT32_MIN, INT32_MAX, s);
+}
+
+static int64_t float128_to_int64_scalbn(float128 a, FloatRoundMode rmode,
+ int scale, float_status *s)
{
- return round_to_int_and_pack(float64_unpack_canonical(a, s),
- rmode, scale, INT64_MIN, INT64_MAX, s);
+ FloatParts128 p;
+
+ float128_unpack_canonical(&p, a, s);
+ return parts_float_to_sint(&p, rmode, scale, INT64_MIN, INT64_MAX, s);
}
int8_t float16_to_int8(float16 a, float_status *s)
return float64_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
}
+int32_t float128_to_int32(float128 a, float_status *s)
+{
+ return float128_to_int32_scalbn(a, s->float_rounding_mode, 0, s);
+}
+
+int64_t float128_to_int64(float128 a, float_status *s)
+{
+ return float128_to_int64_scalbn(a, s->float_rounding_mode, 0, s);
+}
+
int16_t float16_to_int16_round_to_zero(float16 a, float_status *s)
{
return float16_to_int16_scalbn(a, float_round_to_zero, 0, s);
return float64_to_int64_scalbn(a, float_round_to_zero, 0, s);
}
-/*
- * Returns the result of converting the floating-point value `a' to
- * the two's complement integer format.
- */
-
-int16_t bfloat16_to_int16_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
- float_status *s)
-{
- return round_to_int_and_pack(bfloat16_unpack_canonical(a, s),
- rmode, scale, INT16_MIN, INT16_MAX, s);
-}
-
-int32_t bfloat16_to_int32_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
- float_status *s)
+int32_t float128_to_int32_round_to_zero(float128 a, float_status *s)
{
- return round_to_int_and_pack(bfloat16_unpack_canonical(a, s),
- rmode, scale, INT32_MIN, INT32_MAX, s);
+ return float128_to_int32_scalbn(a, float_round_to_zero, 0, s);
}
-int64_t bfloat16_to_int64_scalbn(bfloat16 a, FloatRoundMode rmode, int scale,
- float_status *s)
+int64_t float128_to_int64_round_to_zero(float128 a, float_status *s)
{
- return round_to_int_and_pack(bfloat16_unpack_canonical(a, s),
- rmode, scale, INT64_MIN, INT64_MAX, s);
+ return float128_to_int64_scalbn(a, float_round_to_zero, 0, s);
}
int16_t bfloat16_to_int16(bfloat16 a, float_status *s)
* flag.
*/
-static uint64_t round_to_uint_and_pack(FloatParts in, FloatRoundMode rmode,
+static uint64_t round_to_uint_and_pack(FloatParts64 p, FloatRoundMode rmode,
int scale, uint64_t max,
float_status *s)
{
- int orig_flags = get_float_exception_flags(s);
- FloatParts p = round_to_int(in, rmode, scale, s);
+ int flags = 0;
uint64_t r;
switch (p.cls) {
case float_class_snan:
case float_class_qnan:
- s->float_exception_flags = orig_flags | float_flag_invalid;
- return max;
+ flags = float_flag_invalid;
+ r = max;
+ break;
+
case float_class_inf:
- s->float_exception_flags = orig_flags | float_flag_invalid;
- return p.sign ? 0 : max;
+ flags = float_flag_invalid;
+ r = p.sign ? 0 : max;
+ break;
+
case float_class_zero:
return 0;
+
case float_class_normal:
- if (p.sign) {
- s->float_exception_flags = orig_flags | float_flag_invalid;
- return 0;
+ /* TODO: 62 = N - 2, frac_size for rounding */
+ if (parts_round_to_int_normal(&p, rmode, scale, 62)) {
+ flags = float_flag_inexact;
+ if (p.cls == float_class_zero) {
+ r = 0;
+ break;
+ }
}
- if (p.exp < DECOMPOSED_BINARY_POINT) {
- r = p.frac >> (DECOMPOSED_BINARY_POINT - p.exp);
- } else if (p.exp - DECOMPOSED_BINARY_POINT < 2) {
- r = p.frac << (p.exp - DECOMPOSED_BINARY_POINT);
+ if (p.sign) {
+ flags = float_flag_invalid;
+ r = 0;
+ } else if (p.exp > DECOMPOSED_BINARY_POINT) {
+ flags = float_flag_invalid;
+ r = max;
} else {
- s->float_exception_flags = orig_flags | float_flag_invalid;
- return max;
+ r = p.frac >> (DECOMPOSED_BINARY_POINT - p.exp);
+ if (r > max) {
+ flags = float_flag_invalid;
+ r = max;
+ }
}
+ break;
- /* For uint64 this will never trip, but if p.exp is too large
- * to shift a decomposed fraction we shall have exited via the
- * 3rd leg above.
- */
- if (r > max) {
- s->float_exception_flags = orig_flags | float_flag_invalid;
- return max;
- }
- return r;
default:
g_assert_not_reached();
}
+
+ float_raise(flags, s);
+ return r;
}
uint8_t float16_to_uint8_scalbn(float16 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_uint_and_pack(float16_unpack_canonical(a, s),
- rmode, scale, UINT8_MAX, s);
+ FloatParts64 p;
+
+ float16_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT8_MAX, s);
}
uint16_t float16_to_uint16_scalbn(float16 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_uint_and_pack(float16_unpack_canonical(a, s),
- rmode, scale, UINT16_MAX, s);
+ FloatParts64 p;
+
+ float16_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT16_MAX, s);
}
uint32_t float16_to_uint32_scalbn(float16 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_uint_and_pack(float16_unpack_canonical(a, s),
- rmode, scale, UINT32_MAX, s);
+ FloatParts64 p;
+
+ float16_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT32_MAX, s);
}
uint64_t float16_to_uint64_scalbn(float16 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_uint_and_pack(float16_unpack_canonical(a, s),
- rmode, scale, UINT64_MAX, s);
+ FloatParts64 p;
+
+ float16_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT64_MAX, s);
}
uint16_t float32_to_uint16_scalbn(float32 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_uint_and_pack(float32_unpack_canonical(a, s),
- rmode, scale, UINT16_MAX, s);
+ FloatParts64 p;
+
+ float32_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT16_MAX, s);
}
uint32_t float32_to_uint32_scalbn(float32 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_uint_and_pack(float32_unpack_canonical(a, s),
- rmode, scale, UINT32_MAX, s);
+ FloatParts64 p;
+
+ float32_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT32_MAX, s);
}
uint64_t float32_to_uint64_scalbn(float32 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_uint_and_pack(float32_unpack_canonical(a, s),
- rmode, scale, UINT64_MAX, s);
+ FloatParts64 p;
+
+ float32_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT64_MAX, s);
}
uint16_t float64_to_uint16_scalbn(float64 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_uint_and_pack(float64_unpack_canonical(a, s),
- rmode, scale, UINT16_MAX, s);
+ FloatParts64 p;
+
+ float64_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT16_MAX, s);
}
uint32_t float64_to_uint32_scalbn(float64 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_uint_and_pack(float64_unpack_canonical(a, s),
- rmode, scale, UINT32_MAX, s);
+ FloatParts64 p;
+
+ float64_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT32_MAX, s);
}
uint64_t float64_to_uint64_scalbn(float64 a, FloatRoundMode rmode, int scale,
float_status *s)
{
- return round_to_uint_and_pack(float64_unpack_canonical(a, s),
- rmode, scale, UINT64_MAX, s);
+ FloatParts64 p;
+
+ float64_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT64_MAX, s);
}
uint8_t float16_to_uint8(float16 a, float_status *s)
uint16_t bfloat16_to_uint16_scalbn(bfloat16 a, FloatRoundMode rmode,
int scale, float_status *s)
{
- return round_to_uint_and_pack(bfloat16_unpack_canonical(a, s),
- rmode, scale, UINT16_MAX, s);
+ FloatParts64 p;
+
+ bfloat16_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT16_MAX, s);
}
uint32_t bfloat16_to_uint32_scalbn(bfloat16 a, FloatRoundMode rmode,
int scale, float_status *s)
{
- return round_to_uint_and_pack(bfloat16_unpack_canonical(a, s),
- rmode, scale, UINT32_MAX, s);
+ FloatParts64 p;
+
+ bfloat16_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT32_MAX, s);
}
uint64_t bfloat16_to_uint64_scalbn(bfloat16 a, FloatRoundMode rmode,
int scale, float_status *s)
{
- return round_to_uint_and_pack(bfloat16_unpack_canonical(a, s),
- rmode, scale, UINT64_MAX, s);
+ FloatParts64 p;
+
+ bfloat16_unpack_canonical(&p, a, s);
+ return round_to_uint_and_pack(p, rmode, scale, UINT64_MAX, s);
}
uint16_t bfloat16_to_uint16(bfloat16 a, float_status *s)
* to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
*/
-static FloatParts int_to_float(int64_t a, int scale, float_status *status)
+static FloatParts64 int_to_float(int64_t a, int scale, float_status *status)
{
- FloatParts r = { .sign = false };
+ FloatParts64 r = { .sign = false };
if (a == 0) {
r.cls = float_class_zero;
f = -f;
r.sign = true;
}
- shift = clz64(f) - 1;
+ shift = clz64(f);
scale = MIN(MAX(scale, -0x10000), 0x10000);
r.exp = DECOMPOSED_BINARY_POINT - shift + scale;
- r.frac = (shift < 0 ? DECOMPOSED_IMPLICIT_BIT : f << shift);
+ r.frac = f << shift;
}
return r;
float16 int64_to_float16_scalbn(int64_t a, int scale, float_status *status)
{
- FloatParts pa = int_to_float(a, scale, status);
- return float16_round_pack_canonical(pa, status);
+ FloatParts64 pa = int_to_float(a, scale, status);
+ return float16_round_pack_canonical(&pa, status);
}
float16 int32_to_float16_scalbn(int32_t a, int scale, float_status *status)
float32 int64_to_float32_scalbn(int64_t a, int scale, float_status *status)
{
- FloatParts pa = int_to_float(a, scale, status);
- return float32_round_pack_canonical(pa, status);
+ FloatParts64 pa = int_to_float(a, scale, status);
+ return float32_round_pack_canonical(&pa, status);
}
float32 int32_to_float32_scalbn(int32_t a, int scale, float_status *status)
float64 int64_to_float64_scalbn(int64_t a, int scale, float_status *status)
{
- FloatParts pa = int_to_float(a, scale, status);
- return float64_round_pack_canonical(pa, status);
+ FloatParts64 pa = int_to_float(a, scale, status);
+ return float64_round_pack_canonical(&pa, status);
}
float64 int32_to_float64_scalbn(int32_t a, int scale, float_status *status)
bfloat16 int64_to_bfloat16_scalbn(int64_t a, int scale, float_status *status)
{
- FloatParts pa = int_to_float(a, scale, status);
- return bfloat16_round_pack_canonical(pa, status);
+ FloatParts64 pa = int_to_float(a, scale, status);
+ return bfloat16_round_pack_canonical(&pa, status);
}
bfloat16 int32_to_bfloat16_scalbn(int32_t a, int scale, float_status *status)
* IEC/IEEE Standard for Binary Floating-Point Arithmetic.
*/
-static FloatParts uint_to_float(uint64_t a, int scale, float_status *status)
+static FloatParts64 uint_to_float(uint64_t a, int scale, float_status *status)
{
- FloatParts r = { .sign = false };
+ FloatParts64 r = { .sign = false };
+ int shift;
if (a == 0) {
r.cls = float_class_zero;
} else {
scale = MIN(MAX(scale, -0x10000), 0x10000);
+ shift = clz64(a);
r.cls = float_class_normal;
- if ((int64_t)a < 0) {
- r.exp = DECOMPOSED_BINARY_POINT + 1 + scale;
- shift64RightJamming(a, 1, &a);
- r.frac = a;
- } else {
- int shift = clz64(a) - 1;
- r.exp = DECOMPOSED_BINARY_POINT - shift + scale;
- r.frac = a << shift;
- }
+ r.exp = DECOMPOSED_BINARY_POINT - shift + scale;
+ r.frac = a << shift;
}
return r;
float16 uint64_to_float16_scalbn(uint64_t a, int scale, float_status *status)
{
- FloatParts pa = uint_to_float(a, scale, status);
- return float16_round_pack_canonical(pa, status);
+ FloatParts64 pa = uint_to_float(a, scale, status);
+ return float16_round_pack_canonical(&pa, status);
}
float16 uint32_to_float16_scalbn(uint32_t a, int scale, float_status *status)
float32 uint64_to_float32_scalbn(uint64_t a, int scale, float_status *status)
{
- FloatParts pa = uint_to_float(a, scale, status);
- return float32_round_pack_canonical(pa, status);
+ FloatParts64 pa = uint_to_float(a, scale, status);
+ return float32_round_pack_canonical(&pa, status);
}
float32 uint32_to_float32_scalbn(uint32_t a, int scale, float_status *status)
float64 uint64_to_float64_scalbn(uint64_t a, int scale, float_status *status)
{
- FloatParts pa = uint_to_float(a, scale, status);
- return float64_round_pack_canonical(pa, status);
+ FloatParts64 pa = uint_to_float(a, scale, status);
+ return float64_round_pack_canonical(&pa, status);
}
float64 uint32_to_float64_scalbn(uint32_t a, int scale, float_status *status)
bfloat16 uint64_to_bfloat16_scalbn(uint64_t a, int scale, float_status *status)
{
- FloatParts pa = uint_to_float(a, scale, status);
- return bfloat16_round_pack_canonical(pa, status);
+ FloatParts64 pa = uint_to_float(a, scale, status);
+ return bfloat16_round_pack_canonical(&pa, status);
}
bfloat16 uint32_to_bfloat16_scalbn(uint32_t a, int scale, float_status *status)
* minnummag() and maxnummag() functions correspond to minNumMag()
* and minNumMag() from the IEEE-754 2008.
*/
-static FloatParts minmax_floats(FloatParts a, FloatParts b, bool ismin,
+static FloatParts64 minmax_floats(FloatParts64 a, FloatParts64 b, bool ismin,
bool ieee, bool ismag, float_status *s)
{
if (unlikely(is_nan(a.cls) || is_nan(b.cls))) {
* the invalid exception is raised.
*/
if (is_snan(a.cls) || is_snan(b.cls)) {
- return pick_nan(a, b, s);
+ return *parts_pick_nan(&a, &b, s);
} else if (is_nan(a.cls) && !is_nan(b.cls)) {
return b;
} else if (is_nan(b.cls) && !is_nan(a.cls)) {
return a;
}
}
- return pick_nan(a, b, s);
+ return *parts_pick_nan(&a, &b, s);
} else {
int a_exp, b_exp;
float ## sz float ## sz ## _ ## name(float ## sz a, float ## sz b, \
float_status *s) \
{ \
- FloatParts pa = float ## sz ## _unpack_canonical(a, s); \
- FloatParts pb = float ## sz ## _unpack_canonical(b, s); \
- FloatParts pr = minmax_floats(pa, pb, ismin, isiee, ismag, s); \
- \
- return float ## sz ## _round_pack_canonical(pr, s); \
+ FloatParts64 pa, pb, pr; \
+ float ## sz ## _unpack_canonical(&pa, a, s); \
+ float ## sz ## _unpack_canonical(&pb, b, s); \
+ pr = minmax_floats(pa, pb, ismin, isiee, ismag, s); \
+ return float ## sz ## _round_pack_canonical(&pr, s); \
}
MINMAX(16, min, true, false, false)
#define BF16_MINMAX(name, ismin, isiee, ismag) \
bfloat16 bfloat16_ ## name(bfloat16 a, bfloat16 b, float_status *s) \
{ \
- FloatParts pa = bfloat16_unpack_canonical(a, s); \
- FloatParts pb = bfloat16_unpack_canonical(b, s); \
- FloatParts pr = minmax_floats(pa, pb, ismin, isiee, ismag, s); \
- \
- return bfloat16_round_pack_canonical(pr, s); \
+ FloatParts64 pa, pb, pr; \
+ bfloat16_unpack_canonical(&pa, a, s); \
+ bfloat16_unpack_canonical(&pb, b, s); \
+ pr = minmax_floats(pa, pb, ismin, isiee, ismag, s); \
+ return bfloat16_round_pack_canonical(&pr, s); \
}
BF16_MINMAX(min, true, false, false)
#undef BF16_MINMAX
/* Floating point compare */
-static FloatRelation compare_floats(FloatParts a, FloatParts b, bool is_quiet,
+static FloatRelation compare_floats(FloatParts64 a, FloatParts64 b, bool is_quiet,
float_status *s)
{
if (is_nan(a.cls) || is_nan(b.cls)) {
if (!is_quiet ||
a.cls == float_class_snan ||
b.cls == float_class_snan) {
- s->float_exception_flags |= float_flag_invalid;
+ float_raise(float_flag_invalid, s);
}
return float_relation_unordered;
}
static int attr \
name(float ## sz a, float ## sz b, bool is_quiet, float_status *s) \
{ \
- FloatParts pa = float ## sz ## _unpack_canonical(a, s); \
- FloatParts pb = float ## sz ## _unpack_canonical(b, s); \
+ FloatParts64 pa, pb; \
+ float ## sz ## _unpack_canonical(&pa, a, s); \
+ float ## sz ## _unpack_canonical(&pb, b, s); \
return compare_floats(pa, pb, is_quiet, s); \
}
static FloatRelation QEMU_FLATTEN
soft_bf16_compare(bfloat16 a, bfloat16 b, bool is_quiet, float_status *s)
{
- FloatParts pa = bfloat16_unpack_canonical(a, s);
- FloatParts pb = bfloat16_unpack_canonical(b, s);
+ FloatParts64 pa, pb;
+
+ bfloat16_unpack_canonical(&pa, a, s);
+ bfloat16_unpack_canonical(&pb, b, s);
return compare_floats(pa, pb, is_quiet, s);
}
}
/* Multiply A by 2 raised to the power N. */
-static FloatParts scalbn_decomposed(FloatParts a, int n, float_status *s)
+static FloatParts64 scalbn_decomposed(FloatParts64 a, int n, float_status *s)
{
if (unlikely(is_nan(a.cls))) {
- return return_nan(a, s);
+ parts_return_nan(&a, s);
}
if (a.cls == float_class_normal) {
- /* The largest float type (even though not supported by FloatParts)
+ /* The largest float type (even though not supported by FloatParts64)
* is float128, which has a 15 bit exponent. Bounding N to 16 bits
* still allows rounding to infinity, without allowing overflow
- * within the int32_t that backs FloatParts.exp.
+ * within the int32_t that backs FloatParts64.exp.
*/
n = MIN(MAX(n, -0x10000), 0x10000);
a.exp += n;
float16 float16_scalbn(float16 a, int n, float_status *status)
{
- FloatParts pa = float16_unpack_canonical(a, status);
- FloatParts pr = scalbn_decomposed(pa, n, status);
- return float16_round_pack_canonical(pr, status);
+ FloatParts64 pa, pr;
+
+ float16_unpack_canonical(&pa, a, status);
+ pr = scalbn_decomposed(pa, n, status);
+ return float16_round_pack_canonical(&pr, status);
}
float32 float32_scalbn(float32 a, int n, float_status *status)
{
- FloatParts pa = float32_unpack_canonical(a, status);
- FloatParts pr = scalbn_decomposed(pa, n, status);
- return float32_round_pack_canonical(pr, status);
+ FloatParts64 pa, pr;
+
+ float32_unpack_canonical(&pa, a, status);
+ pr = scalbn_decomposed(pa, n, status);
+ return float32_round_pack_canonical(&pr, status);
}
float64 float64_scalbn(float64 a, int n, float_status *status)
{
- FloatParts pa = float64_unpack_canonical(a, status);
- FloatParts pr = scalbn_decomposed(pa, n, status);
- return float64_round_pack_canonical(pr, status);
+ FloatParts64 pa, pr;
+
+ float64_unpack_canonical(&pa, a, status);
+ pr = scalbn_decomposed(pa, n, status);
+ return float64_round_pack_canonical(&pr, status);
}
bfloat16 bfloat16_scalbn(bfloat16 a, int n, float_status *status)
{
- FloatParts pa = bfloat16_unpack_canonical(a, status);
- FloatParts pr = scalbn_decomposed(pa, n, status);
- return bfloat16_round_pack_canonical(pr, status);
+ FloatParts64 pa, pr;
+
+ bfloat16_unpack_canonical(&pa, a, status);
+ pr = scalbn_decomposed(pa, n, status);
+ return bfloat16_round_pack_canonical(&pr, status);
}
/*
* especially for 64 bit floats.
*/
-static FloatParts sqrt_float(FloatParts a, float_status *s, const FloatFmt *p)
+static FloatParts64 sqrt_float(FloatParts64 a, float_status *s, const FloatFmt *p)
{
uint64_t a_frac, r_frac, s_frac;
int bit, last_bit;
if (is_nan(a.cls)) {
- return return_nan(a, s);
+ parts_return_nan(&a, s);
+ return a;
}
if (a.cls == float_class_zero) {
return a; /* sqrt(+-0) = +-0 */
}
if (a.sign) {
- s->float_exception_flags |= float_flag_invalid;
- return parts_default_nan(s);
+ float_raise(float_flag_invalid, s);
+ parts_default_nan(&a, s);
+ return a;
}
if (a.cls == float_class_inf) {
return a; /* sqrt(+inf) = +inf */
/* We need two overflow bits at the top. Adding room for that is a
* right shift. If the exponent is odd, we can discard the low bit
* by multiplying the fraction by 2; that's a left shift. Combine
- * those and we shift right if the exponent is even.
+ * those and we shift right by 1 if the exponent is odd, otherwise 2.
*/
- a_frac = a.frac;
- if (!(a.exp & 1)) {
- a_frac >>= 1;
- }
+ a_frac = a.frac >> (2 - (a.exp & 1));
a.exp >>= 1;
/* Bit-by-bit computation of sqrt. */
s_frac = 0;
/* Iterate from implicit bit down to the 3 extra bits to compute a
- * properly rounded result. Remember we've inserted one more bit
- * at the top, so these positions are one less.
+ * properly rounded result. Remember we've inserted two more bits
+ * at the top, so these positions are two less.
*/
- bit = DECOMPOSED_BINARY_POINT - 1;
+ bit = DECOMPOSED_BINARY_POINT - 2;
last_bit = MAX(p->frac_shift - 4, 0);
do {
uint64_t q = 1ULL << bit;
/* Undo the right shift done above. If there is any remaining
* fraction, the result is inexact. Set the sticky bit.
*/
- a.frac = (r_frac << 1) + (a_frac != 0);
+ a.frac = (r_frac << 2) + (a_frac != 0);
return a;
}
float16 QEMU_FLATTEN float16_sqrt(float16 a, float_status *status)
{
- FloatParts pa = float16_unpack_canonical(a, status);
- FloatParts pr = sqrt_float(pa, status, &float16_params);
- return float16_round_pack_canonical(pr, status);
+ FloatParts64 pa, pr;
+
+ float16_unpack_canonical(&pa, a, status);
+ pr = sqrt_float(pa, status, &float16_params);
+ return float16_round_pack_canonical(&pr, status);
}
static float32 QEMU_SOFTFLOAT_ATTR
soft_f32_sqrt(float32 a, float_status *status)
{
- FloatParts pa = float32_unpack_canonical(a, status);
- FloatParts pr = sqrt_float(pa, status, &float32_params);
- return float32_round_pack_canonical(pr, status);
+ FloatParts64 pa, pr;
+
+ float32_unpack_canonical(&pa, a, status);
+ pr = sqrt_float(pa, status, &float32_params);
+ return float32_round_pack_canonical(&pr, status);
}
static float64 QEMU_SOFTFLOAT_ATTR
soft_f64_sqrt(float64 a, float_status *status)
{
- FloatParts pa = float64_unpack_canonical(a, status);
- FloatParts pr = sqrt_float(pa, status, &float64_params);
- return float64_round_pack_canonical(pr, status);
+ FloatParts64 pa, pr;
+
+ float64_unpack_canonical(&pa, a, status);
+ pr = sqrt_float(pa, status, &float64_params);
+ return float64_round_pack_canonical(&pr, status);
}
float32 QEMU_FLATTEN float32_sqrt(float32 xa, float_status *s)
bfloat16 QEMU_FLATTEN bfloat16_sqrt(bfloat16 a, float_status *status)
{
- FloatParts pa = bfloat16_unpack_canonical(a, status);
- FloatParts pr = sqrt_float(pa, status, &bfloat16_params);
- return bfloat16_round_pack_canonical(pr, status);
+ FloatParts64 pa, pr;
+
+ bfloat16_unpack_canonical(&pa, a, status);
+ pr = sqrt_float(pa, status, &bfloat16_params);
+ return bfloat16_round_pack_canonical(&pr, status);
}
/*----------------------------------------------------------------------------
float16 float16_default_nan(float_status *status)
{
- FloatParts p = parts_default_nan(status);
+ FloatParts64 p;
+
+ parts_default_nan(&p, status);
p.frac >>= float16_params.frac_shift;
- return float16_pack_raw(p);
+ return float16_pack_raw(&p);
}
float32 float32_default_nan(float_status *status)
{
- FloatParts p = parts_default_nan(status);
+ FloatParts64 p;
+
+ parts_default_nan(&p, status);
p.frac >>= float32_params.frac_shift;
- return float32_pack_raw(p);
+ return float32_pack_raw(&p);
}
float64 float64_default_nan(float_status *status)
{
- FloatParts p = parts_default_nan(status);
+ FloatParts64 p;
+
+ parts_default_nan(&p, status);
p.frac >>= float64_params.frac_shift;
- return float64_pack_raw(p);
+ return float64_pack_raw(&p);
}
float128 float128_default_nan(float_status *status)
{
- FloatParts p = parts_default_nan(status);
- float128 r;
-
- /* Extrapolate from the choices made by parts_default_nan to fill
- * in the quad-floating format. If the low bit is set, assume we
- * want to set all non-snan bits.
- */
- r.low = -(p.frac & 1);
- r.high = p.frac >> (DECOMPOSED_BINARY_POINT - 48);
- r.high |= UINT64_C(0x7FFF000000000000);
- r.high |= (uint64_t)p.sign << 63;
+ FloatParts128 p;
- return r;
+ parts_default_nan(&p, status);
+ frac_shr(&p, float128_params.frac_shift);
+ return float128_pack_raw(&p);
}
bfloat16 bfloat16_default_nan(float_status *status)
{
- FloatParts p = parts_default_nan(status);
+ FloatParts64 p;
+
+ parts_default_nan(&p, status);
p.frac >>= bfloat16_params.frac_shift;
- return bfloat16_pack_raw(p);
+ return bfloat16_pack_raw(&p);
}
/*----------------------------------------------------------------------------
float16 float16_silence_nan(float16 a, float_status *status)
{
- FloatParts p = float16_unpack_raw(a);
+ FloatParts64 p;
+
+ float16_unpack_raw(&p, a);
p.frac <<= float16_params.frac_shift;
- p = parts_silence_nan(p, status);
+ parts_silence_nan(&p, status);
p.frac >>= float16_params.frac_shift;
- return float16_pack_raw(p);
+ return float16_pack_raw(&p);
}
float32 float32_silence_nan(float32 a, float_status *status)
{
- FloatParts p = float32_unpack_raw(a);
+ FloatParts64 p;
+
+ float32_unpack_raw(&p, a);
p.frac <<= float32_params.frac_shift;
- p = parts_silence_nan(p, status);
+ parts_silence_nan(&p, status);
p.frac >>= float32_params.frac_shift;
- return float32_pack_raw(p);
+ return float32_pack_raw(&p);
}
float64 float64_silence_nan(float64 a, float_status *status)
{
- FloatParts p = float64_unpack_raw(a);
+ FloatParts64 p;
+
+ float64_unpack_raw(&p, a);
p.frac <<= float64_params.frac_shift;
- p = parts_silence_nan(p, status);
+ parts_silence_nan(&p, status);
p.frac >>= float64_params.frac_shift;
- return float64_pack_raw(p);
+ return float64_pack_raw(&p);
}
bfloat16 bfloat16_silence_nan(bfloat16 a, float_status *status)
{
- FloatParts p = bfloat16_unpack_raw(a);
+ FloatParts64 p;
+
+ bfloat16_unpack_raw(&p, a);
p.frac <<= bfloat16_params.frac_shift;
- p = parts_silence_nan(p, status);
+ parts_silence_nan(&p, status);
p.frac >>= bfloat16_params.frac_shift;
- return bfloat16_pack_raw(p);
+ return bfloat16_pack_raw(&p);
+}
+
+float128 float128_silence_nan(float128 a, float_status *status)
+{
+ FloatParts128 p;
+
+ float128_unpack_raw(&p, a);
+ frac_shl(&p, float128_params.frac_shift);
+ parts_silence_nan(&p, status);
+ frac_shr(&p, float128_params.frac_shift);
+ return float128_pack_raw(&p);
}
/*----------------------------------------------------------------------------
| input-denormal exception and return zero. Otherwise just return the value.
*----------------------------------------------------------------------------*/
-static bool parts_squash_denormal(FloatParts p, float_status *status)
+static bool parts_squash_denormal(FloatParts64 p, float_status *status)
{
if (p.exp == 0 && p.frac != 0) {
float_raise(float_flag_input_denormal, status);
float16 float16_squash_input_denormal(float16 a, float_status *status)
{
if (status->flush_inputs_to_zero) {
- FloatParts p = float16_unpack_raw(a);
+ FloatParts64 p;
+
+ float16_unpack_raw(&p, a);
if (parts_squash_denormal(p, status)) {
return float16_set_sign(float16_zero, p.sign);
}
float32 float32_squash_input_denormal(float32 a, float_status *status)
{
if (status->flush_inputs_to_zero) {
- FloatParts p = float32_unpack_raw(a);
+ FloatParts64 p;
+
+ float32_unpack_raw(&p, a);
if (parts_squash_denormal(p, status)) {
return float32_set_sign(float32_zero, p.sign);
}
float64 float64_squash_input_denormal(float64 a, float_status *status)
{
if (status->flush_inputs_to_zero) {
- FloatParts p = float64_unpack_raw(a);
+ FloatParts64 p;
+
+ float64_unpack_raw(&p, a);
if (parts_squash_denormal(p, status)) {
return float64_set_sign(float64_zero, p.sign);
}
bfloat16 bfloat16_squash_input_denormal(bfloat16 a, float_status *status)
{
if (status->flush_inputs_to_zero) {
- FloatParts p = bfloat16_unpack_raw(a);
+ FloatParts64 p;
+
+ bfloat16_unpack_raw(&p, a);
if (parts_squash_denormal(p, status)) {
return bfloat16_set_sign(bfloat16_zero, p.sign);
}
return zSign ? INT32_MIN : INT32_MAX;
}
if (roundBits) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
return z;
return zSign ? INT64_MIN : INT64_MAX;
}
if (absZ1) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
return z;
}
if (absZ1) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
return absZ0;
}
}
}
if (roundBits) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
zSig = ( zSig + roundIncrement )>>7;
if (!(roundBits ^ 0x40) && roundNearestEven) {
}
}
if (roundBits) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
zSig = ( zSig + roundIncrement )>>10;
if (!(roundBits ^ 0x200) && roundNearestEven) {
float_raise(float_flag_underflow, status);
}
if (roundBits) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
zSig0 += roundIncrement;
if ( (int64_t) zSig0 < 0 ) zExp = 1;
}
}
if (roundBits) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
zSig0 += roundIncrement;
if ( zSig0 < roundIncrement ) {
float_raise(float_flag_underflow, status);
}
if (zSig1) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
switch (roundingMode) {
case float_round_nearest_even:
}
}
if (zSig1) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
if ( increment ) {
++zSig0;
}
}
if (zSig2) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
if ( increment ) {
add128( zSig0, zSig1, 0, 1, &zSig0, &zSig1 );
}
-/*----------------------------------------------------------------------------
-| Returns the result of converting the single-precision floating-point value
-| `a' to the double-precision floating-point format. The conversion is
-| performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float128 float32_to_float128(float32 a, float_status *status)
-{
- bool aSign;
- int aExp;
- uint32_t aSig;
-
- a = float32_squash_input_denormal(a, status);
- aSig = extractFloat32Frac( a );
- aExp = extractFloat32Exp( a );
- aSign = extractFloat32Sign( a );
- if ( aExp == 0xFF ) {
- if (aSig) {
- return commonNaNToFloat128(float32ToCommonNaN(a, status), status);
- }
- return packFloat128( aSign, 0x7FFF, 0, 0 );
- }
- if ( aExp == 0 ) {
- if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
- normalizeFloat32Subnormal( aSig, &aExp, &aSig );
- --aExp;
- }
- return packFloat128( aSign, aExp + 0x3F80, ( (uint64_t) aSig )<<25, 0 );
-
-}
-
/*----------------------------------------------------------------------------
| Returns the remainder of the single-precision floating-point value `a'
| with respect to the corresponding value `b'. The operation is performed
}
-/*----------------------------------------------------------------------------
-| Returns the result of converting the double-precision floating-point value
-| `a' to the quadruple-precision floating-point format. The conversion is
-| performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float128 float64_to_float128(float64 a, float_status *status)
-{
- bool aSign;
- int aExp;
- uint64_t aSig, zSig0, zSig1;
-
- a = float64_squash_input_denormal(a, status);
- aSig = extractFloat64Frac( a );
- aExp = extractFloat64Exp( a );
- aSign = extractFloat64Sign( a );
- if ( aExp == 0x7FF ) {
- if (aSig) {
- return commonNaNToFloat128(float64ToCommonNaN(a, status), status);
- }
- return packFloat128( aSign, 0x7FFF, 0, 0 );
- }
- if ( aExp == 0 ) {
- if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
- normalizeFloat64Subnormal( aSig, &aExp, &aSig );
- --aExp;
- }
- shift128Right( aSig, 0, 4, &zSig0, &zSig1 );
- return packFloat128( aSign, aExp + 0x3C00, zSig0, zSig1 );
-
-}
-
-
/*----------------------------------------------------------------------------
| Returns the remainder of the double-precision floating-point value `a'
| with respect to the corresponding value `b'. The operation is performed
}
else if ( aExp < 0x3FFF ) {
if (aExp || aSig) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
return 0;
}
return aSign ? (int32_t) 0x80000000 : 0x7FFFFFFF;
}
if ( ( aSig<<shiftCount ) != savedASig ) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
return z;
}
else if ( aExp < 0x3FFF ) {
if (aExp | aSig) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
return 0;
}
z = aSig>>( - shiftCount );
if ( (uint64_t) ( aSig<<( shiftCount & 63 ) ) ) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
if ( aSign ) z = - z;
return z;
&& ( (uint64_t) ( extractFloatx80Frac( a ) ) == 0 ) ) {
return a;
}
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
aSign = extractFloatx80Sign( a );
switch (status->float_rounding_mode) {
case float_round_nearest_even:
z.low = UINT64_C(0x8000000000000000);
}
if (z.low != a.low) {
- status->float_exception_flags |= float_flag_inexact;
+ float_raise(float_flag_inexact, status);
}
return z;
}
/*----------------------------------------------------------------------------
-| Returns the result of converting the quadruple-precision floating-point
-| value `a' to the 32-bit two's complement integer format. The conversion
-| is performed according to the IEC/IEEE Standard for Binary Floating-Point
+| Returns the result of converting the quadruple-precision floating-point value
+| `a' to the 64-bit unsigned integer format. The conversion is
+| performed according to the IEC/IEEE Standard for Binary Floating-Point
| Arithmetic---which means in particular that the conversion is rounded
| according to the current rounding mode. If `a' is a NaN, the largest
-| positive integer is returned. Otherwise, if the conversion overflows, the
-| largest integer with the same sign as `a' is returned.
+| positive integer is returned. If the conversion overflows, the
+| largest unsigned integer is returned. If 'a' is negative, the value is
+| rounded and zero is returned; negative values that do not round to zero
+| will raise the inexact exception.
*----------------------------------------------------------------------------*/
-int32_t float128_to_int32(float128 a, float_status *status)
+uint64_t float128_to_uint64(float128 a, float_status *status)
{
bool aSign;
- int32_t aExp, shiftCount;
+ int aExp;
+ int shiftCount;
uint64_t aSig0, aSig1;
- aSig1 = extractFloat128Frac1( a );
- aSig0 = extractFloat128Frac0( a );
- aExp = extractFloat128Exp( a );
- aSign = extractFloat128Sign( a );
- if ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) aSign = 0;
- if ( aExp ) aSig0 |= UINT64_C(0x0001000000000000);
- aSig0 |= ( aSig1 != 0 );
- shiftCount = 0x4028 - aExp;
- if ( 0 < shiftCount ) shift64RightJamming( aSig0, shiftCount, &aSig0 );
- return roundAndPackInt32(aSign, aSig0, status);
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of converting the quadruple-precision floating-point
-| value `a' to the 32-bit two's complement integer format. The conversion
-| is performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic, except that the conversion is always rounded toward zero. If
-| `a' is a NaN, the largest positive integer is returned. Otherwise, if the
-| conversion overflows, the largest integer with the same sign as `a' is
-| returned.
-*----------------------------------------------------------------------------*/
-
-int32_t float128_to_int32_round_to_zero(float128 a, float_status *status)
-{
- bool aSign;
- int32_t aExp, shiftCount;
- uint64_t aSig0, aSig1, savedASig;
- int32_t z;
-
- aSig1 = extractFloat128Frac1( a );
- aSig0 = extractFloat128Frac0( a );
- aExp = extractFloat128Exp( a );
- aSign = extractFloat128Sign( a );
- aSig0 |= ( aSig1 != 0 );
- if ( 0x401E < aExp ) {
- if ( ( aExp == 0x7FFF ) && aSig0 ) aSign = 0;
- goto invalid;
- }
- else if ( aExp < 0x3FFF ) {
- if (aExp || aSig0) {
- status->float_exception_flags |= float_flag_inexact;
- }
- return 0;
- }
- aSig0 |= UINT64_C(0x0001000000000000);
- shiftCount = 0x402F - aExp;
- savedASig = aSig0;
- aSig0 >>= shiftCount;
- z = aSig0;
- if ( aSign ) z = - z;
- if ( ( z < 0 ) ^ aSign ) {
- invalid:
- float_raise(float_flag_invalid, status);
- return aSign ? INT32_MIN : INT32_MAX;
- }
- if ( ( aSig0<<shiftCount ) != savedASig ) {
- status->float_exception_flags |= float_flag_inexact;
- }
- return z;
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of converting the quadruple-precision floating-point
-| value `a' to the 64-bit two's complement integer format. The conversion
-| is performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic---which means in particular that the conversion is rounded
-| according to the current rounding mode. If `a' is a NaN, the largest
-| positive integer is returned. Otherwise, if the conversion overflows, the
-| largest integer with the same sign as `a' is returned.
-*----------------------------------------------------------------------------*/
-
-int64_t float128_to_int64(float128 a, float_status *status)
-{
- bool aSign;
- int32_t aExp, shiftCount;
- uint64_t aSig0, aSig1;
-
- aSig1 = extractFloat128Frac1( a );
- aSig0 = extractFloat128Frac0( a );
- aExp = extractFloat128Exp( a );
- aSign = extractFloat128Sign( a );
- if ( aExp ) aSig0 |= UINT64_C(0x0001000000000000);
- shiftCount = 0x402F - aExp;
- if ( shiftCount <= 0 ) {
- if ( 0x403E < aExp ) {
- float_raise(float_flag_invalid, status);
- if ( ! aSign
- || ( ( aExp == 0x7FFF )
- && ( aSig1 || ( aSig0 != UINT64_C(0x0001000000000000) ) )
- )
- ) {
- return INT64_MAX;
- }
- return INT64_MIN;
- }
- shortShift128Left( aSig0, aSig1, - shiftCount, &aSig0, &aSig1 );
- }
- else {
- shift64ExtraRightJamming( aSig0, aSig1, shiftCount, &aSig0, &aSig1 );
- }
- return roundAndPackInt64(aSign, aSig0, aSig1, status);
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of converting the quadruple-precision floating-point
-| value `a' to the 64-bit two's complement integer format. The conversion
-| is performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic, except that the conversion is always rounded toward zero.
-| If `a' is a NaN, the largest positive integer is returned. Otherwise, if
-| the conversion overflows, the largest integer with the same sign as `a' is
-| returned.
-*----------------------------------------------------------------------------*/
-
-int64_t float128_to_int64_round_to_zero(float128 a, float_status *status)
-{
- bool aSign;
- int32_t aExp, shiftCount;
- uint64_t aSig0, aSig1;
- int64_t z;
-
- aSig1 = extractFloat128Frac1( a );
- aSig0 = extractFloat128Frac0( a );
- aExp = extractFloat128Exp( a );
- aSign = extractFloat128Sign( a );
- if ( aExp ) aSig0 |= UINT64_C(0x0001000000000000);
- shiftCount = aExp - 0x402F;
- if ( 0 < shiftCount ) {
- if ( 0x403E <= aExp ) {
- aSig0 &= UINT64_C(0x0000FFFFFFFFFFFF);
- if ( ( a.high == UINT64_C(0xC03E000000000000) )
- && ( aSig1 < UINT64_C(0x0002000000000000) ) ) {
- if (aSig1) {
- status->float_exception_flags |= float_flag_inexact;
- }
- }
- else {
- float_raise(float_flag_invalid, status);
- if ( ! aSign || ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) ) {
- return INT64_MAX;
- }
- }
- return INT64_MIN;
- }
- z = ( aSig0<<shiftCount ) | ( aSig1>>( ( - shiftCount ) & 63 ) );
- if ( (uint64_t) ( aSig1<<shiftCount ) ) {
- status->float_exception_flags |= float_flag_inexact;
- }
- }
- else {
- if ( aExp < 0x3FFF ) {
- if ( aExp | aSig0 | aSig1 ) {
- status->float_exception_flags |= float_flag_inexact;
- }
- return 0;
- }
- z = aSig0>>( - shiftCount );
- if ( aSig1
- || ( shiftCount && (uint64_t) ( aSig0<<( shiftCount & 63 ) ) ) ) {
- status->float_exception_flags |= float_flag_inexact;
- }
- }
- if ( aSign ) z = - z;
- return z;
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of converting the quadruple-precision floating-point value
-| `a' to the 64-bit unsigned integer format. The conversion is
-| performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic---which means in particular that the conversion is rounded
-| according to the current rounding mode. If `a' is a NaN, the largest
-| positive integer is returned. If the conversion overflows, the
-| largest unsigned integer is returned. If 'a' is negative, the value is
-| rounded and zero is returned; negative values that do not round to zero
-| will raise the inexact exception.
-*----------------------------------------------------------------------------*/
-
-uint64_t float128_to_uint64(float128 a, float_status *status)
-{
- bool aSign;
- int aExp;
- int shiftCount;
- uint64_t aSig0, aSig1;
-
- aSig0 = extractFloat128Frac0(a);
- aSig1 = extractFloat128Frac1(a);
- aExp = extractFloat128Exp(a);
- aSign = extractFloat128Sign(a);
- if (aSign && (aExp > 0x3FFE)) {
- float_raise(float_flag_invalid, status);
- if (float128_is_any_nan(a)) {
- return UINT64_MAX;
- } else {
- return 0;
- }
- }
- if (aExp) {
- aSig0 |= UINT64_C(0x0001000000000000);
- }
- shiftCount = 0x402F - aExp;
- if (shiftCount <= 0) {
- if (0x403E < aExp) {
- float_raise(float_flag_invalid, status);
- return UINT64_MAX;
- }
- shortShift128Left(aSig0, aSig1, -shiftCount, &aSig0, &aSig1);
- } else {
- shift64ExtraRightJamming(aSig0, aSig1, shiftCount, &aSig0, &aSig1);
- }
- return roundAndPackUint64(aSign, aSig0, aSig1, status);
+ aSig0 = extractFloat128Frac0(a);
+ aSig1 = extractFloat128Frac1(a);
+ aExp = extractFloat128Exp(a);
+ aSign = extractFloat128Sign(a);
+ if (aSign && (aExp > 0x3FFE)) {
+ float_raise(float_flag_invalid, status);
+ if (float128_is_any_nan(a)) {
+ return UINT64_MAX;
+ } else {
+ return 0;
+ }
+ }
+ if (aExp) {
+ aSig0 |= UINT64_C(0x0001000000000000);
+ }
+ shiftCount = 0x402F - aExp;
+ if (shiftCount <= 0) {
+ if (0x403E < aExp) {
+ float_raise(float_flag_invalid, status);
+ return UINT64_MAX;
+ }
+ shortShift128Left(aSig0, aSig1, -shiftCount, &aSig0, &aSig1);
+ } else {
+ shift64ExtraRightJamming(aSig0, aSig1, shiftCount, &aSig0, &aSig1);
+ }
+ return roundAndPackUint64(aSign, aSig0, aSig1, status);
}
uint64_t float128_to_uint64_round_to_zero(float128 a, float_status *status)
return res;
}
-/*----------------------------------------------------------------------------
-| Returns the result of converting the quadruple-precision floating-point
-| value `a' to the single-precision floating-point format. The conversion
-| is performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float32 float128_to_float32(float128 a, float_status *status)
-{
- bool aSign;
- int32_t aExp;
- uint64_t aSig0, aSig1;
- uint32_t zSig;
-
- aSig1 = extractFloat128Frac1( a );
- aSig0 = extractFloat128Frac0( a );
- aExp = extractFloat128Exp( a );
- aSign = extractFloat128Sign( a );
- if ( aExp == 0x7FFF ) {
- if ( aSig0 | aSig1 ) {
- return commonNaNToFloat32(float128ToCommonNaN(a, status), status);
- }
- return packFloat32( aSign, 0xFF, 0 );
- }
- aSig0 |= ( aSig1 != 0 );
- shift64RightJamming( aSig0, 18, &aSig0 );
- zSig = aSig0;
- if ( aExp || zSig ) {
- zSig |= 0x40000000;
- aExp -= 0x3F81;
- }
- return roundAndPackFloat32(aSign, aExp, zSig, status);
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of converting the quadruple-precision floating-point
-| value `a' to the double-precision floating-point format. The conversion
-| is performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float64 float128_to_float64(float128 a, float_status *status)
-{
- bool aSign;
- int32_t aExp;
- uint64_t aSig0, aSig1;
-
- aSig1 = extractFloat128Frac1( a );
- aSig0 = extractFloat128Frac0( a );
- aExp = extractFloat128Exp( a );
- aSign = extractFloat128Sign( a );
- if ( aExp == 0x7FFF ) {
- if ( aSig0 | aSig1 ) {
- return commonNaNToFloat64(float128ToCommonNaN(a, status), status);
- }
- return packFloat64( aSign, 0x7FF, 0 );
- }
- shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
- aSig0 |= ( aSig1 != 0 );
- if ( aExp || aSig0 ) {
- aSig0 |= UINT64_C(0x4000000000000000);
- aExp -= 0x3C01;
- }
- return roundAndPackFloat64(aSign, aExp, aSig0, status);
-
-}
-
/*----------------------------------------------------------------------------
| Returns the result of converting the quadruple-precision floating-point
| value `a' to the extended double-precision floating-point format. The
}
-/*----------------------------------------------------------------------------
-| Rounds the quadruple-precision floating-point value `a' to an integer, and
-| returns the result as a quadruple-precision floating-point value. The
-| operation is performed according to the IEC/IEEE Standard for Binary
-| Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float128 float128_round_to_int(float128 a, float_status *status)
-{
- bool aSign;
- int32_t aExp;
- uint64_t lastBitMask, roundBitsMask;
- float128 z;
-
- aExp = extractFloat128Exp( a );
- if ( 0x402F <= aExp ) {
- if ( 0x406F <= aExp ) {
- if ( ( aExp == 0x7FFF )
- && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) )
- ) {
- return propagateFloat128NaN(a, a, status);
- }
- return a;
- }
- lastBitMask = 1;
- lastBitMask = ( lastBitMask<<( 0x406E - aExp ) )<<1;
- roundBitsMask = lastBitMask - 1;
- z = a;
- switch (status->float_rounding_mode) {
- case float_round_nearest_even:
- if ( lastBitMask ) {
- add128( z.high, z.low, 0, lastBitMask>>1, &z.high, &z.low );
- if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
- }
- else {
- if ( (int64_t) z.low < 0 ) {
- ++z.high;
- if ( (uint64_t) ( z.low<<1 ) == 0 ) z.high &= ~1;
- }
- }
- break;
- case float_round_ties_away:
- if (lastBitMask) {
- add128(z.high, z.low, 0, lastBitMask >> 1, &z.high, &z.low);
- } else {
- if ((int64_t) z.low < 0) {
- ++z.high;
- }
- }
- break;
- case float_round_to_zero:
- break;
- case float_round_up:
- if (!extractFloat128Sign(z)) {
- add128(z.high, z.low, 0, roundBitsMask, &z.high, &z.low);
- }
- break;
- case float_round_down:
- if (extractFloat128Sign(z)) {
- add128(z.high, z.low, 0, roundBitsMask, &z.high, &z.low);
- }
- break;
- case float_round_to_odd:
- /*
- * Note that if lastBitMask == 0, the last bit is the lsb
- * of high, and roundBitsMask == -1.
- */
- if ((lastBitMask ? z.low & lastBitMask : z.high & 1) == 0) {
- add128(z.high, z.low, 0, roundBitsMask, &z.high, &z.low);
- }
- break;
- default:
- abort();
- }
- z.low &= ~ roundBitsMask;
- }
- else {
- if ( aExp < 0x3FFF ) {
- if ( ( ( (uint64_t) ( a.high<<1 ) ) | a.low ) == 0 ) return a;
- status->float_exception_flags |= float_flag_inexact;
- aSign = extractFloat128Sign( a );
- switch (status->float_rounding_mode) {
- case float_round_nearest_even:
- if ( ( aExp == 0x3FFE )
- && ( extractFloat128Frac0( a )
- | extractFloat128Frac1( a ) )
- ) {
- return packFloat128( aSign, 0x3FFF, 0, 0 );
- }
- break;
- case float_round_ties_away:
- if (aExp == 0x3FFE) {
- return packFloat128(aSign, 0x3FFF, 0, 0);
- }
- break;
- case float_round_down:
- return
- aSign ? packFloat128( 1, 0x3FFF, 0, 0 )
- : packFloat128( 0, 0, 0, 0 );
- case float_round_up:
- return
- aSign ? packFloat128( 1, 0, 0, 0 )
- : packFloat128( 0, 0x3FFF, 0, 0 );
-
- case float_round_to_odd:
- return packFloat128(aSign, 0x3FFF, 0, 0);
-
- case float_round_to_zero:
- break;
- }
- return packFloat128( aSign, 0, 0, 0 );
- }
- lastBitMask = 1;
- lastBitMask <<= 0x402F - aExp;
- roundBitsMask = lastBitMask - 1;
- z.low = 0;
- z.high = a.high;
- switch (status->float_rounding_mode) {
- case float_round_nearest_even:
- z.high += lastBitMask>>1;
- if ( ( ( z.high & roundBitsMask ) | a.low ) == 0 ) {
- z.high &= ~ lastBitMask;
- }
- break;
- case float_round_ties_away:
- z.high += lastBitMask>>1;
- break;
- case float_round_to_zero:
- break;
- case float_round_up:
- if (!extractFloat128Sign(z)) {
- z.high |= ( a.low != 0 );
- z.high += roundBitsMask;
- }
- break;
- case float_round_down:
- if (extractFloat128Sign(z)) {
- z.high |= (a.low != 0);
- z.high += roundBitsMask;
- }
- break;
- case float_round_to_odd:
- if ((z.high & lastBitMask) == 0) {
- z.high |= (a.low != 0);
- z.high += roundBitsMask;
- }
- break;
- default:
- abort();
- }
- z.high &= ~ roundBitsMask;
- }
- if ( ( z.low != a.low ) || ( z.high != a.high ) ) {
- status->float_exception_flags |= float_flag_inexact;
- }
- return z;
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of adding the absolute values of the quadruple-precision
-| floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
-| before being returned. `zSign' is ignored if the result is a NaN.
-| The addition is performed according to the IEC/IEEE Standard for Binary
-| Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-static float128 addFloat128Sigs(float128 a, float128 b, bool zSign,
- float_status *status)
-{
- int32_t aExp, bExp, zExp;
- uint64_t aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
- int32_t expDiff;
-
- aSig1 = extractFloat128Frac1( a );
- aSig0 = extractFloat128Frac0( a );
- aExp = extractFloat128Exp( a );
- bSig1 = extractFloat128Frac1( b );
- bSig0 = extractFloat128Frac0( b );
- bExp = extractFloat128Exp( b );
- expDiff = aExp - bExp;
- if ( 0 < expDiff ) {
- if ( aExp == 0x7FFF ) {
- if (aSig0 | aSig1) {
- return propagateFloat128NaN(a, b, status);
- }
- return a;
- }
- if ( bExp == 0 ) {
- --expDiff;
- }
- else {
- bSig0 |= UINT64_C(0x0001000000000000);
- }
- shift128ExtraRightJamming(
- bSig0, bSig1, 0, expDiff, &bSig0, &bSig1, &zSig2 );
- zExp = aExp;
- }
- else if ( expDiff < 0 ) {
- if ( bExp == 0x7FFF ) {
- if (bSig0 | bSig1) {
- return propagateFloat128NaN(a, b, status);
- }
- return packFloat128( zSign, 0x7FFF, 0, 0 );
- }
- if ( aExp == 0 ) {
- ++expDiff;
- }
- else {
- aSig0 |= UINT64_C(0x0001000000000000);
- }
- shift128ExtraRightJamming(
- aSig0, aSig1, 0, - expDiff, &aSig0, &aSig1, &zSig2 );
- zExp = bExp;
- }
- else {
- if ( aExp == 0x7FFF ) {
- if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
- return propagateFloat128NaN(a, b, status);
- }
- return a;
- }
- add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
- if ( aExp == 0 ) {
- if (status->flush_to_zero) {
- if (zSig0 | zSig1) {
- float_raise(float_flag_output_denormal, status);
- }
- return packFloat128(zSign, 0, 0, 0);
- }
- return packFloat128( zSign, 0, zSig0, zSig1 );
- }
- zSig2 = 0;
- zSig0 |= UINT64_C(0x0002000000000000);
- zExp = aExp;
- goto shiftRight1;
- }
- aSig0 |= UINT64_C(0x0001000000000000);
- add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
- --zExp;
- if ( zSig0 < UINT64_C(0x0002000000000000) ) goto roundAndPack;
- ++zExp;
- shiftRight1:
- shift128ExtraRightJamming(
- zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
- roundAndPack:
- return roundAndPackFloat128(zSign, zExp, zSig0, zSig1, zSig2, status);
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of subtracting the absolute values of the quadruple-
-| precision floating-point values `a' and `b'. If `zSign' is 1, the
-| difference is negated before being returned. `zSign' is ignored if the
-| result is a NaN. The subtraction is performed according to the IEC/IEEE
-| Standard for Binary Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-static float128 subFloat128Sigs(float128 a, float128 b, bool zSign,
- float_status *status)
-{
- int32_t aExp, bExp, zExp;
- uint64_t aSig0, aSig1, bSig0, bSig1, zSig0, zSig1;
- int32_t expDiff;
-
- aSig1 = extractFloat128Frac1( a );
- aSig0 = extractFloat128Frac0( a );
- aExp = extractFloat128Exp( a );
- bSig1 = extractFloat128Frac1( b );
- bSig0 = extractFloat128Frac0( b );
- bExp = extractFloat128Exp( b );
- expDiff = aExp - bExp;
- shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
- shortShift128Left( bSig0, bSig1, 14, &bSig0, &bSig1 );
- if ( 0 < expDiff ) goto aExpBigger;
- if ( expDiff < 0 ) goto bExpBigger;
- if ( aExp == 0x7FFF ) {
- if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
- return propagateFloat128NaN(a, b, status);
- }
- float_raise(float_flag_invalid, status);
- return float128_default_nan(status);
- }
- if ( aExp == 0 ) {
- aExp = 1;
- bExp = 1;
- }
- if ( bSig0 < aSig0 ) goto aBigger;
- if ( aSig0 < bSig0 ) goto bBigger;
- if ( bSig1 < aSig1 ) goto aBigger;
- if ( aSig1 < bSig1 ) goto bBigger;
- return packFloat128(status->float_rounding_mode == float_round_down,
- 0, 0, 0);
- bExpBigger:
- if ( bExp == 0x7FFF ) {
- if (bSig0 | bSig1) {
- return propagateFloat128NaN(a, b, status);
- }
- return packFloat128( zSign ^ 1, 0x7FFF, 0, 0 );
- }
- if ( aExp == 0 ) {
- ++expDiff;
- }
- else {
- aSig0 |= UINT64_C(0x4000000000000000);
- }
- shift128RightJamming( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
- bSig0 |= UINT64_C(0x4000000000000000);
- bBigger:
- sub128( bSig0, bSig1, aSig0, aSig1, &zSig0, &zSig1 );
- zExp = bExp;
- zSign ^= 1;
- goto normalizeRoundAndPack;
- aExpBigger:
- if ( aExp == 0x7FFF ) {
- if (aSig0 | aSig1) {
- return propagateFloat128NaN(a, b, status);
- }
- return a;
- }
- if ( bExp == 0 ) {
- --expDiff;
- }
- else {
- bSig0 |= UINT64_C(0x4000000000000000);
- }
- shift128RightJamming( bSig0, bSig1, expDiff, &bSig0, &bSig1 );
- aSig0 |= UINT64_C(0x4000000000000000);
- aBigger:
- sub128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
- zExp = aExp;
- normalizeRoundAndPack:
- --zExp;
- return normalizeRoundAndPackFloat128(zSign, zExp - 14, zSig0, zSig1,
- status);
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of adding the quadruple-precision floating-point values
-| `a' and `b'. The operation is performed according to the IEC/IEEE Standard
-| for Binary Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float128 float128_add(float128 a, float128 b, float_status *status)
-{
- bool aSign, bSign;
-
- aSign = extractFloat128Sign( a );
- bSign = extractFloat128Sign( b );
- if ( aSign == bSign ) {
- return addFloat128Sigs(a, b, aSign, status);
- }
- else {
- return subFloat128Sigs(a, b, aSign, status);
- }
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of subtracting the quadruple-precision floating-point
-| values `a' and `b'. The operation is performed according to the IEC/IEEE
-| Standard for Binary Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float128 float128_sub(float128 a, float128 b, float_status *status)
-{
- bool aSign, bSign;
-
- aSign = extractFloat128Sign( a );
- bSign = extractFloat128Sign( b );
- if ( aSign == bSign ) {
- return subFloat128Sigs(a, b, aSign, status);
- }
- else {
- return addFloat128Sigs(a, b, aSign, status);
- }
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of multiplying the quadruple-precision floating-point
-| values `a' and `b'. The operation is performed according to the IEC/IEEE
-| Standard for Binary Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float128 float128_mul(float128 a, float128 b, float_status *status)
-{
- bool aSign, bSign, zSign;
- int32_t aExp, bExp, zExp;
- uint64_t aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2, zSig3;
-
- aSig1 = extractFloat128Frac1( a );
- aSig0 = extractFloat128Frac0( a );
- aExp = extractFloat128Exp( a );
- aSign = extractFloat128Sign( a );
- bSig1 = extractFloat128Frac1( b );
- bSig0 = extractFloat128Frac0( b );
- bExp = extractFloat128Exp( b );
- bSign = extractFloat128Sign( b );
- zSign = aSign ^ bSign;
- if ( aExp == 0x7FFF ) {
- if ( ( aSig0 | aSig1 )
- || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
- return propagateFloat128NaN(a, b, status);
- }
- if ( ( bExp | bSig0 | bSig1 ) == 0 ) goto invalid;
- return packFloat128( zSign, 0x7FFF, 0, 0 );
- }
- if ( bExp == 0x7FFF ) {
- if (bSig0 | bSig1) {
- return propagateFloat128NaN(a, b, status);
- }
- if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
- invalid:
- float_raise(float_flag_invalid, status);
- return float128_default_nan(status);
- }
- return packFloat128( zSign, 0x7FFF, 0, 0 );
- }
- if ( aExp == 0 ) {
- if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
- normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
- }
- if ( bExp == 0 ) {
- if ( ( bSig0 | bSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
- normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
- }
- zExp = aExp + bExp - 0x4000;
- aSig0 |= UINT64_C(0x0001000000000000);
- shortShift128Left( bSig0, bSig1, 16, &bSig0, &bSig1 );
- mul128To256( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1, &zSig2, &zSig3 );
- add128( zSig0, zSig1, aSig0, aSig1, &zSig0, &zSig1 );
- zSig2 |= ( zSig3 != 0 );
- if (UINT64_C( 0x0002000000000000) <= zSig0 ) {
- shift128ExtraRightJamming(
- zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
- ++zExp;
- }
- return roundAndPackFloat128(zSign, zExp, zSig0, zSig1, zSig2, status);
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of dividing the quadruple-precision floating-point value
-| `a' by the corresponding value `b'. The operation is performed according to
-| the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float128 float128_div(float128 a, float128 b, float_status *status)
-{
- bool aSign, bSign, zSign;
- int32_t aExp, bExp, zExp;
- uint64_t aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
- uint64_t rem0, rem1, rem2, rem3, term0, term1, term2, term3;
-
- aSig1 = extractFloat128Frac1( a );
- aSig0 = extractFloat128Frac0( a );
- aExp = extractFloat128Exp( a );
- aSign = extractFloat128Sign( a );
- bSig1 = extractFloat128Frac1( b );
- bSig0 = extractFloat128Frac0( b );
- bExp = extractFloat128Exp( b );
- bSign = extractFloat128Sign( b );
- zSign = aSign ^ bSign;
- if ( aExp == 0x7FFF ) {
- if (aSig0 | aSig1) {
- return propagateFloat128NaN(a, b, status);
- }
- if ( bExp == 0x7FFF ) {
- if (bSig0 | bSig1) {
- return propagateFloat128NaN(a, b, status);
- }
- goto invalid;
- }
- return packFloat128( zSign, 0x7FFF, 0, 0 );
- }
- if ( bExp == 0x7FFF ) {
- if (bSig0 | bSig1) {
- return propagateFloat128NaN(a, b, status);
- }
- return packFloat128( zSign, 0, 0, 0 );
- }
- if ( bExp == 0 ) {
- if ( ( bSig0 | bSig1 ) == 0 ) {
- if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
- invalid:
- float_raise(float_flag_invalid, status);
- return float128_default_nan(status);
- }
- float_raise(float_flag_divbyzero, status);
- return packFloat128( zSign, 0x7FFF, 0, 0 );
- }
- normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
- }
- if ( aExp == 0 ) {
- if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
- normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
- }
- zExp = aExp - bExp + 0x3FFD;
- shortShift128Left(
- aSig0 | UINT64_C(0x0001000000000000), aSig1, 15, &aSig0, &aSig1 );
- shortShift128Left(
- bSig0 | UINT64_C(0x0001000000000000), bSig1, 15, &bSig0, &bSig1 );
- if ( le128( bSig0, bSig1, aSig0, aSig1 ) ) {
- shift128Right( aSig0, aSig1, 1, &aSig0, &aSig1 );
- ++zExp;
- }
- zSig0 = estimateDiv128To64( aSig0, aSig1, bSig0 );
- mul128By64To192( bSig0, bSig1, zSig0, &term0, &term1, &term2 );
- sub192( aSig0, aSig1, 0, term0, term1, term2, &rem0, &rem1, &rem2 );
- while ( (int64_t) rem0 < 0 ) {
- --zSig0;
- add192( rem0, rem1, rem2, 0, bSig0, bSig1, &rem0, &rem1, &rem2 );
- }
- zSig1 = estimateDiv128To64( rem1, rem2, bSig0 );
- if ( ( zSig1 & 0x3FFF ) <= 4 ) {
- mul128By64To192( bSig0, bSig1, zSig1, &term1, &term2, &term3 );
- sub192( rem1, rem2, 0, term1, term2, term3, &rem1, &rem2, &rem3 );
- while ( (int64_t) rem1 < 0 ) {
- --zSig1;
- add192( rem1, rem2, rem3, 0, bSig0, bSig1, &rem1, &rem2, &rem3 );
- }
- zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
- }
- shift128ExtraRightJamming( zSig0, zSig1, 0, 15, &zSig0, &zSig1, &zSig2 );
- return roundAndPackFloat128(zSign, zExp, zSig0, zSig1, zSig2, status);
-
-}
-
/*----------------------------------------------------------------------------
| Returns the remainder of the quadruple-precision floating-point value `a'
| with respect to the corresponding value `b'. The operation is performed
projects / mirror_qemu.git / commitdiff