Based-on: 20260507160518.444797-1-richard.henderson@linaro.org ("[PULL 00/29] fpu patch queue") No changes for v3 other than merging 29 patches. r~ Richard Henderson (11): fpu: Reorganize partsN(muladd) fpu: Return struct from parts{64,128}_muladd fpu: Hoist nan check in partsN_addsub fpu: Simplify 0 +/- N case in parts_addsub fpu: Use parts64_round_to_int in parts_s390_divide_to_integer target/s390x: Move float{32,64}_s390_divide_to_integer target/arm: Use FloatParts64 in bfdotadd_ebf target/arm: Drop oddstatus from is_ebf and bfdotadd_ebf target/arm: Use FloatParts64 in f16_dotadd fpu: Saturate the exponent in uncanon_normal fpu: Introduce exp_scalbn include/fpu/softfloat-parts.h | 9 ++ include/fpu/softfloat.h | 11 -- target/arm/tcg/vec_internal.h | 12 +- fpu/softfloat.c | 224 ++++++++-------------------------- target/arm/tcg/sme_helper.c | 102 +++++++--------- target/arm/tcg/vec_helper.c | 127 ++++++++++--------- target/s390x/tcg/fpu_helper.c | 135 ++++++++++++++++++++ fpu/softfloat-parts.c.inc | 223 ++++++++++++++++----------------- 8 files changed, 407 insertions(+), 436 deletions(-) -- 2.43.0 Check the likely case of normal product and normal or zero addend first; shift NaN and infinity detection down; end with zero product + addend. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat-parts.c.inc | 155 +++++++++++++++++--------------------- 1 file changed, 70 insertions(+), 85 deletions(-) diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, FloatPartsN *c, int flags, float_status *s) { - int ab_mask, abc_mask; - FloatPartsW p_widen, c_widen; + int ab_mask = float_cmask(a->cls) | float_cmask(b->cls); + int c_mask = float_cmask(c->cls); + int abc_mask = ab_mask | c_mask; + bool c_sign = c->sign ^ !!(flags & float_muladd_negate_c); + bool p_sign = a->sign ^ b->sign ^ !!(flags & float_muladd_negate_product); - ab_mask = float_cmask(a->cls) | float_cmask(b->cls); - abc_mask = float_cmask(c->cls) | ab_mask; + /* + * The "likely" case is A and B normal, so that the product is normal, + * and C normal or zero so that the result is normal. + */ + int likely_mask = ab_mask | (c_mask & ~float_cmask_zero); + if (likely(cmask_is_only_normals(likely_mask))) { + record_denormals_used(abc_mask, s); + + /* Perform the multiplication step. */ + FloatPartsW p_widen = { .sign = p_sign, .exp = a->exp + b->exp + 1 }; + fracN(mulw)(&p_widen, a, b); + if (!(p_widen.frac_hi & DECOMPOSED_IMPLICIT_BIT)) { + fracW(add)(&p_widen, &p_widen, &p_widen); + p_widen.exp -= 1; + } + + /* Perform the addition step. */ + if (!(c_mask & float_cmask_zero)) { + /* Zero-extend C to less significant bits. */ + FloatPartsW c_widen = { .sign = c_sign, .exp = c->exp }; + fracN(widen)(&c_widen, c); + + if (p_sign == c_sign) { + partsW(add_normal)(&p_widen, &c_widen); + } else if (!partsW(sub_normal)(&p_widen, &c_widen)) { + goto return_sub_zero; + } + } + + /* Narrow with sticky bit, for proper rounding later. */ + fracN(truncjam)(a, &p_widen); + a->sign = p_widen.sign; + a->exp = p_widen.exp; + return a; + } /* * It is implementation-defined whether the cases of (0,inf,qnan) @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, return a; } - if (flags & float_muladd_negate_c) { - c->sign ^= 1; + if (unlikely(ab_mask == float_cmask_infzero)) { + /* Inf * Zero == NaN */ + float_raise(float_flag_invalid | float_flag_invalid_imz, s); + goto d_nan; } - /* Compute the sign of the product into A. */ - a->sign ^= b->sign; - if (flags & float_muladd_negate_product) { - a->sign ^= 1; - } - - if (unlikely(!cmask_is_only_normals(ab_mask))) { - if (unlikely(ab_mask == float_cmask_infzero)) { - float_raise(float_flag_invalid | float_flag_invalid_imz, s); + if (unlikely(ab_mask & float_cmask_inf)) { + if ((c_mask & float_cmask_inf) && p_sign != c_sign) { + /* Inf - Inf == NaN */ + float_raise(float_flag_invalid | float_flag_invalid_isi, s); goto d_nan; } - - if (ab_mask & float_cmask_inf) { - if (c->cls == float_class_inf && a->sign != c->sign) { - float_raise(float_flag_invalid | float_flag_invalid_isi, s); - goto d_nan; - } - goto return_inf; - } - - g_assert(ab_mask & float_cmask_zero); - if (is_anynorm(c->cls)) { - *a = *c; - goto finish_sign; - } - if (c->cls == float_class_zero) { - if (flags & float_muladd_suppress_add_product_zero) { - a->sign = c->sign; - } else if (a->sign != c->sign) { - goto return_sub_zero; - } - goto return_zero; - } - g_assert(c->cls == float_class_inf); + /* Inf + C == Inf */ + record_denormals_used(abc_mask, s); + a->sign = p_sign; + a->cls = float_class_inf; + return a; } - - if (unlikely(c->cls == float_class_inf)) { - a->sign = c->sign; - goto return_inf; - } - - /* Perform the multiplication step. */ - p_widen.sign = a->sign; - p_widen.exp = a->exp + b->exp + 1; - fracN(mulw)(&p_widen, a, b); - if (!(p_widen.frac_hi & DECOMPOSED_IMPLICIT_BIT)) { - fracW(add)(&p_widen, &p_widen, &p_widen); - p_widen.exp -= 1; - } - - /* Perform the addition step. */ - if (c->cls != float_class_zero) { - /* Zero-extend C to less significant bits. */ - fracN(widen)(&c_widen, c); - c_widen.exp = c->exp; - - if (a->sign == c->sign) { - partsW(add_normal)(&p_widen, &c_widen); - } else if (!partsW(sub_normal)(&p_widen, &c_widen)) { - goto return_sub_zero; - } - } - - /* Narrow with sticky bit, for proper rounding later. */ - fracN(truncjam)(a, &p_widen); - a->sign = p_widen.sign; - a->exp = p_widen.exp; - - finish_sign: - /* - * All result types except for "return the default NaN - * because this is an Invalid Operation" go through here; - * this matches the set of cases where we consumed a - * denormal input. - */ record_denormals_used(abc_mask, s); - return a; + + /* Only remaining case is zero product. */ + assert(ab_mask & float_cmask_zero); + + /* + * 0 + C == C, + * except for 0 - 0, which needs special rounding, + * except for when we want to suppress this addition step. + */ + if (!(c_mask & float_cmask_zero) + || p_sign == c_sign + || (flags & float_muladd_suppress_add_product_zero)) { + c->sign = c_sign; + return c; + } return_sub_zero: + /* 0 - 0 == -0 for round_down, +0 otherwise. */ a->sign = s->float_rounding_mode == float_round_down; - return_zero: a->cls = float_class_zero; - goto finish_sign; - - return_inf: - a->cls = float_class_inf; - goto finish_sign; + return a; d_nan: *a = partsN(default_nan)(s); -- 2.43.0 At the same time, export. Reviewed-by: Philippe Mathieu-Daudé <philmd@linaro.org> Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- include/fpu/softfloat-parts.h | 9 ++++ fpu/softfloat.c | 83 +++++++++++++++++------------------ fpu/softfloat-parts.c.inc | 38 ++++++++-------- 3 files changed, 69 insertions(+), 61 deletions(-) diff --git a/include/fpu/softfloat-parts.h b/include/fpu/softfloat-parts.h index XXXXXXX..XXXXXXX 100644 --- a/include/fpu/softfloat-parts.h +++ b/include/fpu/softfloat-parts.h @@ -XXX,XX +XXX,XX @@ FloatParts64 parts64_mul(const FloatParts64 *a, const FloatParts64 *b, FloatParts128 parts128_mul(const FloatParts128 *a, const FloatParts128 *b, float_status *s); +FloatParts64 parts64_muladd(const FloatParts64 *a, + const FloatParts64 *b, + const FloatParts64 *c, + int flags, float_status *s); +FloatParts128 parts128_muladd(const FloatParts128 *a, + const FloatParts128 *b, + const FloatParts128 *c, + int flags, float_status *s); + FloatParts64 parts64_round_to_int(const FloatParts64 *a, FloatRoundMode rmode, int scale, float_status *s, diff --git a/fpu/softfloat.c b/fpu/softfloat.c index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -XXX,XX +XXX,XX @@ float16 float16_muladd_scalbn(float16 a, float16 b, float16 c, FloatParts64 pa = float16_unpack_canonical(a, status); FloatParts64 pb = float16_unpack_canonical(b, status); FloatParts64 pc = float16_unpack_canonical(c, status); - FloatParts64 *pr = parts64_muladd(&pa, &pb, &pc, flags, status); + FloatParts64 pr = parts64_muladd(&pa, &pb, &pc, flags, status); /* Before rounding, scale. */ if (scale) { - parts64_scalbn(pr, scale, status); + parts64_scalbn(&pr, scale, status); } - parts64_uncanon(pr, status, &float16_params, false); + parts64_uncanon(&pr, status, &float16_params, false); /* After rounding, apply negate result, especially for -0.0. */ - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return pack_raw64(pr, &float16_params); + return pack_raw64(&pr, &float16_params); } float16 float16_muladd(float16 a, float16 b, float16 c, @@ -XXX,XX +XXX,XX @@ float32_muladd_scalbn(float32 a, float32 b, float32 c, FloatParts64 pa = float32_unpack_canonical(a, status); FloatParts64 pb = float32_unpack_canonical(b, status); FloatParts64 pc = float32_unpack_canonical(c, status); - FloatParts64 *pr = parts64_muladd(&pa, &pb, &pc, flags, status); + FloatParts64 pr = parts64_muladd(&pa, &pb, &pc, flags, status); /* Before rounding, scale. */ if (scale) { - parts64_scalbn(pr, scale, status); + parts64_scalbn(&pr, scale, status); } - parts64_uncanon(pr, status, &float32_params, false); + parts64_uncanon(&pr, status, &float32_params, false); /* After rounding, apply negate result, especially for -0.0. */ - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return pack_raw64(pr, &float32_params); + return pack_raw64(&pr, &float32_params); } float64 QEMU_SOFTFLOAT_ATTR @@ -XXX,XX +XXX,XX @@ float64_muladd_scalbn(float64 a, float64 b, float64 c, FloatParts64 pa = float64_unpack_canonical(a, status); FloatParts64 pb = float64_unpack_canonical(b, status); FloatParts64 pc = float64_unpack_canonical(c, status); - FloatParts64 *pr = parts64_muladd(&pa, &pb, &pc, flags, status); + FloatParts64 pr = parts64_muladd(&pa, &pb, &pc, flags, status); /* Before rounding, scale. */ if (scale) { - parts64_scalbn(pr, scale, status); + parts64_scalbn(&pr, scale, status); } - parts64_uncanon(pr, status, &float64_params, false); + parts64_uncanon(&pr, status, &float64_params, false); /* After rounding, apply negate result, especially for -0.0. */ - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return pack_raw64(pr, &float64_params); + return pack_raw64(&pr, &float64_params); } static bool force_soft_fma; @@ -XXX,XX +XXX,XX @@ float64 float64r32_muladd(float64 a, float64 b, float64 c, FloatParts64 pa = float64_unpack_canonical(a, status); FloatParts64 pb = float64_unpack_canonical(b, status); FloatParts64 pc = float64_unpack_canonical(c, status); - FloatParts64 *pr = parts64_muladd(&pa, &pb, &pc, flags, status); + FloatParts64 pr = parts64_muladd(&pa, &pb, &pc, flags, status); /* Round before applying negate result. */ - parts64_uncanon(pr, status, &float32_params, false); - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + parts64_uncanon(&pr, status, &float32_params, false); + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return float64r32_pack_raw(pr); + return float64r32_pack_raw(&pr); } bfloat16 bfloat16_muladd(bfloat16 a, bfloat16 b, bfloat16 c, @@ -XXX,XX +XXX,XX @@ bfloat16 bfloat16_muladd(bfloat16 a, bfloat16 b, bfloat16 c, FloatParts64 pa = bfloat16_unpack_canonical(a, status); FloatParts64 pb = bfloat16_unpack_canonical(b, status); FloatParts64 pc = bfloat16_unpack_canonical(c, status); - FloatParts64 *pr = parts64_muladd(&pa, &pb, &pc, flags, status); + FloatParts64 pr = parts64_muladd(&pa, &pb, &pc, flags, status); /* Round before applying negate result. */ - parts64_uncanon(pr, status, &bfloat16_params, false); - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + parts64_uncanon(&pr, status, &bfloat16_params, false); + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return pack_raw64(pr, &bfloat16_params); + return pack_raw64(&pr, &bfloat16_params); } float128 float128_muladd(float128 a, float128 b, float128 c, @@ -XXX,XX +XXX,XX @@ float128 float128_muladd(float128 a, float128 b, float128 c, FloatParts128 pa = float128_unpack_canonical(a, status); FloatParts128 pb = float128_unpack_canonical(b, status); FloatParts128 pc = float128_unpack_canonical(c, status); - FloatParts128 *pr = parts128_muladd(&pa, &pb, &pc, flags, status); + FloatParts128 pr = parts128_muladd(&pa, &pb, &pc, flags, status); /* Round before applying negate result. */ - parts128_uncanon(pr, status, &float128_params, false); - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + parts128_uncanon(&pr, status, &float128_params, false); + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return float128_pack_raw(pr); + return float128_pack_raw(&pr); } /* @@ -XXX,XX +XXX,XX @@ float32 float32_exp2(float32 a, float_status *status) rp = float64_unpack_canonical(float64_one, status); for (int i = 0; i < 15; i++) { tp = float64_unpack_canonical(float32_exp2_coefficients[i], status); - rp = *parts64_muladd(&tp, &xnp, &rp, 0, status); + rp = parts64_muladd(&tp, &xnp, &rp, 0, status); xnp = parts64_mul(&xnp, &xp, status); } @@ -XXX,XX +XXX,XX @@ static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, n->sign = a->sign ^ b->sign; *cc = 0; } else { - FloatParts64 *q, q_buf, *r_precise, r_precise_buf; + FloatParts64 *q, q_buf, r_precise; int float_exception_flags = 0; bool is_q_smallish; uint32_t r_flags; @@ -XXX,XX +XXX,XX @@ static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, 0, fmt->frac_size); /* Compute precise remainder */ - r_precise_buf = *b; - r_precise = parts64_muladd(&r_precise_buf, n, a, + r_precise = parts64_muladd(b, n, a, float_muladd_negate_product, status); /* Round remainder to the target format */ - *r = *r_precise; + *r = r_precise; status->float_exception_flags = 0; *r = parts64_round_to_fmt(r, status, fmt); r_flags = status->float_exception_flags; @@ -XXX,XX +XXX,XX @@ static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, * toward zero) or incremented. */ saved_r_sign = r->sign; - saved_r_precise_sign = r_precise->sign; + saved_r_precise_sign = r_precise.sign; r->sign = false; - r_precise->sign = false; - if (parts64_compare(r, r_precise, status, true) < + r_precise.sign = false; + if (parts64_compare(r, &r_precise, status, true) < float_relation_equal) { *dxc = 0x8; } else { *dxc = 0xc; } r->sign = saved_r_sign; - r_precise->sign = saved_r_precise_sign; + r_precise.sign = saved_r_precise_sign; } } } diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(mul)(const FloatPartsN *a, const FloatPartsN *b, * Requires A and C extracted into a double-sized structure to provide the * extra space for the widening multiply. */ -static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, - FloatPartsN *c, - int flags, float_status *s) +FloatPartsN partsN(muladd)(const FloatPartsN *a, const FloatPartsN *b, + const FloatPartsN *c, int flags, float_status *s) { int ab_mask = float_cmask(a->cls) | float_cmask(b->cls); int c_mask = float_cmask(c->cls); @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, } /* Narrow with sticky bit, for proper rounding later. */ - fracN(truncjam)(a, &p_widen); - a->sign = p_widen.sign; - a->exp = p_widen.exp; - return a; + FloatPartsN r = { + .sign = p_widen.sign, + .exp = p_widen.exp, + .cls = float_class_normal, + }; + fracN(truncjam)(&r, &p_widen); + return r; } /* @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, * off to the target-specific pick-a-NaN routine. */ if (unlikely(abc_mask & float_cmask_anynan)) { - *a = partsN(pick_nan_muladd)(a, b, c, s, ab_mask, abc_mask); - return a; + return partsN(pick_nan_muladd)(a, b, c, s, ab_mask, abc_mask); } if (unlikely(ab_mask == float_cmask_infzero)) { @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, } /* Inf + C == Inf */ record_denormals_used(abc_mask, s); - a->sign = p_sign; - a->cls = float_class_inf; - return a; + return (FloatPartsN){ .sign = p_sign, .cls = float_class_inf }; } record_denormals_used(abc_mask, s); @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, if (!(c_mask & float_cmask_zero) || p_sign == c_sign || (flags & float_muladd_suppress_add_product_zero)) { - c->sign = c_sign; - return c; + FloatPartsN r = *c; + r.sign = c_sign; + return r; } return_sub_zero: /* 0 - 0 == -0 for round_down, +0 otherwise. */ - a->sign = s->float_rounding_mode == float_round_down; - a->cls = float_class_zero; - return a; + return (FloatPartsN){ + .sign = s->float_rounding_mode == float_round_down, + .cls = float_class_zero + }; d_nan: - *a = partsN(default_nan)(s); - return a; + return partsN(default_nan)(s); } /* -- 2.43.0 The nan test had been down below because it was unlikely. But if we have to have one anyway because of denormals, we might as well take care of them right away. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat-parts.c.inc | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, { int ab_mask = float_cmask(a_orig->cls) | float_cmask(b_orig->cls); + if (unlikely(ab_mask & float_cmask_anynan)) { + return partsN(pick_nan)(a_orig, b_orig, s); + } + /* * For addition and subtraction, we will consume an * input denormal unless the other input is a NaN. */ - if (!(ab_mask & float_cmask_anynan)) { - record_denormals_used(ab_mask, s); - } + record_denormals_used(ab_mask, s); FloatPartsN a = *a_orig; FloatPartsN b = *b_orig; @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, return a; } - if (unlikely(ab_mask & float_cmask_anynan)) { - goto p_nan; - } - if (ab_mask & float_cmask_inf) { if (a.cls != float_class_inf) { /* N - Inf */ @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, return a; } - if (unlikely(ab_mask & float_cmask_anynan)) { - goto p_nan; - } - if (ab_mask & float_cmask_inf) { a.cls = float_class_inf; return a; @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, g_assert(a.cls == float_class_zero); g_assert(is_anynorm(b.cls)); return b; - - p_nan: - return partsN(pick_nan)(a_orig, b_orig, s); } /* -- 2.43.0 Consolidate the tests for zero and anynorm. Add comments for a few cases. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat-parts.c.inc | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, } if (ab_mask == float_cmask_zero) { + /* 0 - 0 */ a.sign = s->float_rounding_mode == float_round_down; return a; } @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, } if (ab_mask == float_cmask_zero) { + /* 0 + 0 */ return a; } if (ab_mask & float_cmask_inf) { + /* N + Inf or Inf + N */ a.cls = float_class_inf; return a; } } - if (b.cls == float_class_zero) { - g_assert(is_anynorm(a.cls)); - return a; - } - - g_assert(a.cls == float_class_zero); - g_assert(is_anynorm(b.cls)); - return b; + /* 0 +/- N or N +/- 0 */ + assert((ab_mask & float_cmask_zero) && (ab_mask & float_cmask_anynorm)); + return b.cls == float_class_zero ? a : b; } /* -- 2.43.0 We will not expose parts_round_to_int_normal. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat.c | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/fpu/softfloat.c b/fpu/softfloat.c index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -XXX,XX +XXX,XX @@ static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, * Rounding of partial quotient may be inexact. This is the whole point * of distinguishing partial quotients, so ignore the exception. */ - *n = *q; - parts64_round_to_int_normal(n, - is_q_smallish - ? final_quotient_rounding_mode - : float_round_to_zero, - 0, fmt->frac_size); + *n = parts64_round_to_int(q, + is_q_smallish + ? final_quotient_rounding_mode + : float_round_to_zero, + 0, status, fmt); /* Compute precise remainder */ r_precise = parts64_muladd(b, n, a, -- 2.43.0 Now that we've exposed enough infrastructure, this can be implemented in the backend that needs it. Reviewed-by: Philippe Mathieu-Daudé <philmd@linaro.org> Reviewed-by: Ilya Leoshkevich <iii@linux.ibm.com> Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- include/fpu/softfloat.h | 11 --- fpu/softfloat.c | 137 ---------------------------------- target/s390x/tcg/fpu_helper.c | 135 +++++++++++++++++++++++++++++++++ 3 files changed, 135 insertions(+), 148 deletions(-) diff --git a/include/fpu/softfloat.h b/include/fpu/softfloat.h index XXXXXXX..XXXXXXX 100644 --- a/include/fpu/softfloat.h +++ b/include/fpu/softfloat.h @@ -XXX,XX +XXX,XX @@ static inline bool float128_unordered_quiet(float128 a, float128 b, *----------------------------------------------------------------------------*/ float128 float128_default_nan(float_status *status); -#define DECLARE_S390_DIVIDE_TO_INTEGER(floatN) \ -void floatN ## _s390_divide_to_integer(floatN a, floatN b, \ - int final_quotient_rounding_mode, \ - bool mask_underflow, bool mask_inexact, \ - floatN *r, floatN *n, \ - uint32_t *cc, int *dxc, \ - float_status *status) -DECLARE_S390_DIVIDE_TO_INTEGER(float32); -DECLARE_S390_DIVIDE_TO_INTEGER(float64); - - #endif /* SOFTFLOAT_H */ diff --git a/fpu/softfloat.c b/fpu/softfloat.c index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -XXX,XX +XXX,XX @@ floatx80 floatx80_round(floatx80 a, float_status *status) return floatx80_round_pack_canonical(&p, status); } -static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, - int final_quotient_rounding_mode, - bool mask_underflow, bool mask_inexact, - const FloatFmt *fmt, - FloatParts64 *r, FloatParts64 *n, - uint32_t *cc, int *dxc, - float_status *status) -{ - /* POp table "Results: DIVIDE TO INTEGER (Part 1 of 2)" */ - if ((float_cmask(a->cls) | float_cmask(b->cls)) & float_cmask_anynan) { - *r = parts64_pick_nan(a, b, status); - *n = *r; - *cc = 1; - } else if (a->cls == float_class_inf || b->cls == float_class_zero) { - *r = parts64_default_nan(status); - *n = *r; - *cc = 1; - status->float_exception_flags |= float_flag_invalid; - } else if (b->cls == float_class_inf) { - *r = *a; - n->cls = float_class_zero; - n->sign = a->sign ^ b->sign; - *cc = 0; - } else { - FloatParts64 *q, q_buf, r_precise; - int float_exception_flags = 0; - bool is_q_smallish; - uint32_t r_flags; - - /* Compute precise quotient */ - q_buf = parts64_div(a, b, status); - q = &q_buf; - - /* - * Check whether two closest integers can be precisely represented, - * i.e., all their bits fit into the fractional part. - */ - is_q_smallish = q->exp < (fmt->frac_size + 1); - - /* - * Final quotient is rounded using final-quotient-rounding method, and - * partial quotient is rounded toward zero. - * - * Rounding of partial quotient may be inexact. This is the whole point - * of distinguishing partial quotients, so ignore the exception. - */ - *n = parts64_round_to_int(q, - is_q_smallish - ? final_quotient_rounding_mode - : float_round_to_zero, - 0, status, fmt); - - /* Compute precise remainder */ - r_precise = parts64_muladd(b, n, a, - float_muladd_negate_product, status); - - /* Round remainder to the target format */ - *r = r_precise; - status->float_exception_flags = 0; - *r = parts64_round_to_fmt(r, status, fmt); - r_flags = status->float_exception_flags; - - /* POp table "Results: DIVIDE TO INTEGER (Part 2 of 2)" */ - if (is_q_smallish) { - if (r->cls != float_class_zero) { - if (r->exp < 2 - (1 << (fmt->exp_size - 1))) { - if (mask_underflow) { - float_exception_flags |= float_flag_underflow; - *dxc = 0x10; - r->exp += fmt->exp_re_bias; - } - } else if (r_flags & float_flag_inexact) { - float_exception_flags |= float_flag_inexact; - if (mask_inexact) { - bool saved_r_sign, saved_r_precise_sign; - - /* - * Check whether remainder was truncated (rounded - * toward zero) or incremented. - */ - saved_r_sign = r->sign; - saved_r_precise_sign = r_precise.sign; - r->sign = false; - r_precise.sign = false; - if (parts64_compare(r, &r_precise, status, true) < - float_relation_equal) { - *dxc = 0x8; - } else { - *dxc = 0xc; - } - r->sign = saved_r_sign; - r_precise.sign = saved_r_precise_sign; - } - } - } - *cc = 0; - } else if (n->exp > (1 << (fmt->exp_size - 1)) - 1) { - n->exp -= fmt->exp_re_bias; - *cc = r->cls == float_class_zero ? 1 : 3; - } else { - *cc = r->cls == float_class_zero ? 0 : 2; - } - - /* Adjust signs of zero results */ - if (r->cls == float_class_zero) { - r->sign = a->sign; - } - if (n->cls == float_class_zero) { - n->sign = a->sign ^ b->sign; - } - - status->float_exception_flags = float_exception_flags; - } -} - -#define DEFINE_S390_DIVIDE_TO_INTEGER(floatN) \ -void floatN ## _s390_divide_to_integer(floatN a, floatN b, \ - int final_quotient_rounding_mode, \ - bool mask_underflow, bool mask_inexact, \ - floatN *r, floatN *n, \ - uint32_t *cc, int *dxc, \ - float_status *status) \ -{ \ - FloatParts64 pa = floatN ## _unpack_canonical(a, status); \ - FloatParts64 pb = floatN ## _unpack_canonical(b, status); \ - FloatParts64 pr, pn; \ - parts_s390_divide_to_integer(&pa, &pb, final_quotient_rounding_mode, \ - mask_underflow, mask_inexact, \ - &floatN ## _params, \ - &pr, &pn, cc, dxc, status); \ - *r = floatN ## _round_pack_canonical(&pr, status); \ - *n = floatN ## _round_pack_canonical(&pn, status); \ -} - -DEFINE_S390_DIVIDE_TO_INTEGER(float32) -DEFINE_S390_DIVIDE_TO_INTEGER(float64) - static void __attribute__((constructor)) softfloat_init(void) { union_float64 ua, ub, uc, ur; diff --git a/target/s390x/tcg/fpu_helper.c b/target/s390x/tcg/fpu_helper.c index XXXXXXX..XXXXXXX 100644 --- a/target/s390x/tcg/fpu_helper.c +++ b/target/s390x/tcg/fpu_helper.c @@ -XXX,XX +XXX,XX @@ #include "tcg_s390x.h" #include "exec/helper-proto.h" #include "fpu/softfloat.h" +#include "fpu/softfloat-parts.h" /* #define DEBUG_HELPER */ #ifdef DEBUG_HELPER @@ -XXX,XX +XXX,XX @@ Int128 HELPER(dxb)(CPUS390XState *env, Int128 a, Int128 b) return RET128(ret); } +static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, + int final_quotient_rounding_mode, + bool mask_underflow, bool mask_inexact, + const FloatFmt *fmt, + FloatParts64 *r, FloatParts64 *n, + uint32_t *cc, int *dxc, + float_status *status) +{ + /* POp table "Results: DIVIDE TO INTEGER (Part 1 of 2)" */ + if ((float_cmask(a->cls) | float_cmask(b->cls)) & float_cmask_anynan) { + *r = parts64_pick_nan(a, b, status); + *n = *r; + *cc = 1; + } else if (a->cls == float_class_inf || b->cls == float_class_zero) { + *r = parts64_default_nan(status); + *n = *r; + *cc = 1; + status->float_exception_flags |= float_flag_invalid; + } else if (b->cls == float_class_inf) { + *r = *a; + n->cls = float_class_zero; + n->sign = a->sign ^ b->sign; + *cc = 0; + } else { + FloatParts64 *q, q_buf, r_precise; + int float_exception_flags = 0; + bool is_q_smallish; + uint32_t r_flags; + + /* Compute precise quotient */ + q_buf = parts64_div(a, b, status); + q = &q_buf; + + /* + * Check whether two closest integers can be precisely represented, + * i.e., all their bits fit into the fractional part. + */ + is_q_smallish = q->exp < (fmt->frac_size + 1); + + /* + * Final quotient is rounded using final-quotient-rounding method, and + * partial quotient is rounded toward zero. + * + * Rounding of partial quotient may be inexact. This is the whole point + * of distinguishing partial quotients, so ignore the exception. + */ + *n = parts64_round_to_int(q, + is_q_smallish + ? final_quotient_rounding_mode + : float_round_to_zero, + 0, status, fmt); + + /* Compute precise remainder */ + r_precise = parts64_muladd(b, n, a, + float_muladd_negate_product, status); + + /* Round remainder to the target format */ + *r = r_precise; + status->float_exception_flags = 0; + *r = parts64_round_to_fmt(r, status, fmt); + r_flags = status->float_exception_flags; + + /* POp table "Results: DIVIDE TO INTEGER (Part 2 of 2)" */ + if (is_q_smallish) { + if (r->cls != float_class_zero) { + if (r->exp < 2 - (1 << (fmt->exp_size - 1))) { + if (mask_underflow) { + float_exception_flags |= float_flag_underflow; + *dxc = 0x10; + r->exp += fmt->exp_re_bias; + } + } else if (r_flags & float_flag_inexact) { + float_exception_flags |= float_flag_inexact; + if (mask_inexact) { + bool saved_r_sign, saved_r_precise_sign; + + /* + * Check whether remainder was truncated (rounded + * toward zero) or incremented. + */ + saved_r_sign = r->sign; + saved_r_precise_sign = r_precise.sign; + r->sign = false; + r_precise.sign = false; + if (parts64_compare(r, &r_precise, status, true) < + float_relation_equal) { + *dxc = 0x8; + } else { + *dxc = 0xc; + } + r->sign = saved_r_sign; + r_precise.sign = saved_r_precise_sign; + } + } + } + *cc = 0; + } else if (n->exp > (1 << (fmt->exp_size - 1)) - 1) { + n->exp -= fmt->exp_re_bias; + *cc = r->cls == float_class_zero ? 1 : 3; + } else { + *cc = r->cls == float_class_zero ? 0 : 2; + } + + /* Adjust signs of zero results */ + if (r->cls == float_class_zero) { + r->sign = a->sign; + } + if (n->cls == float_class_zero) { + n->sign = a->sign ^ b->sign; + } + + status->float_exception_flags = float_exception_flags; + } +} + +#define DEFINE_S390_DIVIDE_TO_INTEGER(floatN) \ +static void floatN ## _s390_divide_to_integer(floatN a, floatN b, \ + int final_quotient_rounding_mode, bool mask_underflow, bool mask_inexact, \ + floatN *r, floatN *n, uint32_t *cc, int *dxc, float_status *status) \ +{ \ + FloatParts64 pa = floatN ## _unpack_canonical(a, status); \ + FloatParts64 pb = floatN ## _unpack_canonical(b, status); \ + FloatParts64 pr, pn; \ + parts_s390_divide_to_integer(&pa, &pb, final_quotient_rounding_mode, \ + mask_underflow, mask_inexact, \ + &floatN ## _params, \ + &pr, &pn, cc, dxc, status); \ + *r = floatN ## _round_pack_canonical(&pr, status); \ + *n = floatN ## _round_pack_canonical(&pn, status); \ +} + +DEFINE_S390_DIVIDE_TO_INTEGER(float32) +DEFINE_S390_DIVIDE_TO_INTEGER(float64) + void HELPER(dib)(CPUS390XState *env, uint32_t r1, uint32_t r2, uint32_t r3, uint32_t m4, uint32_t bits) { -- 2.43.0 Use softfloat-parts.h so that we can more naturally perform the required operations witha single rounding step. This happens to also simplify the NaN detection step. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- target/arm/tcg/vec_helper.c | 77 +++++++++++++++++++------------------ 1 file changed, 40 insertions(+), 37 deletions(-) diff --git a/target/arm/tcg/vec_helper.c b/target/arm/tcg/vec_helper.c index XXXXXXX..XXXXXXX 100644 --- a/target/arm/tcg/vec_helper.c +++ b/target/arm/tcg/vec_helper.c @@ -XXX,XX +XXX,XX @@ #include "helper.h" #include "tcg/tcg-gvec-desc.h" #include "fpu/softfloat.h" +#include "fpu/softfloat-parts.h" #include "qemu/int128.h" #include "crypto/clmul.h" #include "vec_internal.h" @@ -XXX,XX +XXX,XX @@ float32 bfdotadd(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst) float32 bfdotadd_ebf(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst, float_status *fpst_odd) { + /* Unpack two BFloat16 into two Float32, trivially. */ float32 s1r = e1 << 16; float32 s1c = e1 & 0xffff0000u; float32 s2r = e2 << 16; float32 s2c = e2 & 0xffff0000u; float32 t32; + /* + * Compare f16_dotadd() in sme_helper.c, but here we have + * bfloat16 inputs. In particular that means that we do not + * want the FPCR.FZ16 flush semantics, so we use the normal + * float_status for the input handling here. + */ + FloatParts64 p1r = float32_unpack_canonical(s1r, fpst); + FloatParts64 p1c = float32_unpack_canonical(s1c, fpst); + FloatParts64 p2r = float32_unpack_canonical(s2r, fpst); + FloatParts64 p2c = float32_unpack_canonical(s2c, fpst); + + int all_mask = (float_cmask(p1r.cls) | float_cmask(p1c.cls) | + float_cmask(p1r.cls) | float_cmask(p1c.cls)); + /* C.f. FPProcessNaNs4 */ - if (float32_is_any_nan(s1r) || float32_is_any_nan(s1c) || - float32_is_any_nan(s2r) || float32_is_any_nan(s2c)) { - if (float32_is_signaling_nan(s1r, fpst)) { - t32 = s1r; - } else if (float32_is_signaling_nan(s1c, fpst)) { - t32 = s1c; - } else if (float32_is_signaling_nan(s2r, fpst)) { - t32 = s2r; - } else if (float32_is_signaling_nan(s2c, fpst)) { - t32 = s2c; - } else if (float32_is_any_nan(s1r)) { - t32 = s1r; - } else if (float32_is_any_nan(s1c)) { - t32 = s1c; - } else if (float32_is_any_nan(s2r)) { - t32 = s2r; + if (unlikely(all_mask & float_cmask_anynan)) { + if (unlikely(all_mask & float_cmask_snan)) { + if (p1r.cls == float_class_snan) { + t32 = s1r; + } else if (p1c.cls == float_class_snan) { + t32 = s1c; + } else if (p2r.cls == float_class_snan) { + t32 = s2r; + } else { + t32 = s2c; + } } else { - t32 = s2c; + if (p1r.cls == float_class_qnan) { + t32 = s1r; + } else if (p1c.cls == float_class_qnan) { + t32 = s1c; + } else if (p2r.cls == float_class_qnan) { + t32 = s2r; + } else { + t32 = s2c; + } } /* * FPConvertNaN(FPProcessNaN(t32)) will be done as part * of the final addition below. */ } else { - /* - * Compare f16_dotadd() in sme_helper.c, but here we have - * bfloat16 inputs. In particular that means that we do not - * want the FPCR.FZ16 flush semantics, so we use the normal - * float_status for the input handling here. - */ - float64 e1r = float32_to_float64(s1r, fpst); - float64 e1c = float32_to_float64(s1c, fpst); - float64 e2r = float32_to_float64(s2r, fpst); - float64 e2c = float32_to_float64(s2c, fpst); - float64 t64; - /* * The ARM pseudocode function FPDot performs both multiplies - * and the add with a single rounding operation. Emulate this - * by performing the first multiply in round-to-odd, then doing - * the second multiply as fused multiply-add, and rounding to - * float32 all in one step. + * and the add with a single rounding operation. */ - t64 = float64_mul(e1r, e2r, fpst_odd); - t64 = float64r32_muladd(e1c, e2c, t64, 0, fpst); + FloatParts64 tmp = parts64_mul(&p1r, &p2r, fpst); + tmp = parts64_muladd(&p1c, &p2c, &tmp, 0, fpst); - /* This conversion is exact, because we've already rounded. */ - t32 = float64_to_float32(t64, fpst); + t32 = float32_round_pack_canonical(&tmp, fpst); } /* The final accumulation step is not fused. */ -- 2.43.0 This argument is no longer used. Reviewed-by: Philippe Mathieu-Daudé <philmd@linaro.org> Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- target/arm/tcg/vec_internal.h | 12 +++------ target/arm/tcg/sme_helper.c | 6 ++--- target/arm/tcg/vec_helper.c | 50 +++++++++++++++-------------------- 3 files changed, 29 insertions(+), 39 deletions(-) diff --git a/target/arm/tcg/vec_internal.h b/target/arm/tcg/vec_internal.h index XXXXXXX..XXXXXXX 100644 --- a/target/arm/tcg/vec_internal.h +++ b/target/arm/tcg/vec_internal.h @@ -XXX,XX +XXX,XX @@ float32 bfdotadd(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst); * @sum: addend * @e1, @e2: multiplicand vectors * @fpst: floating-point status to use - * @fpst_odd: floating-point status to use for round-to-odd operations * * BFloat16 2-way dot product of @e1 & @e2, accumulating with @sum. * The @e1 and @e2 operands correspond to the 32-bit source vector @@ -XXX,XX +XXX,XX @@ float32 bfdotadd(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst); * Corresponds to the ARM pseudocode function BFDotAdd, specialized * for the FPCR.EBF == 1 case. */ -float32 bfdotadd_ebf(float32 sum, uint32_t e1, uint32_t e2, - float_status *fpst, float_status *fpst_odd); +float32 bfdotadd_ebf(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst); /** * is_ebf: * @env: CPU state * @statusp: pointer to floating point status to fill in - * @oddstatusp: pointer to floating point status to fill in for round-to-odd * * Determine whether a BFDotAdd operation should use FPCR.EBF = 0 - * or FPCR.EBF = 1 semantics. On return, has initialized *statusp - * and *oddstatusp to suitable float_status arguments to use with either - * bfdotadd() or bfdotadd_ebf(). + * or FPCR.EBF = 1 semantics. On return, has initialized *statusp as suitable + * for float_status arguments to either bfdotadd() or bfdotadd_ebf(). * Returns true for EBF = 1, false for EBF = 0. (The caller should use this * to decide whether to call bfdotadd() or bfdotadd_ebf().) */ -bool is_ebf(CPUARMState *env, float_status *statusp, float_status *oddstatusp); +bool is_ebf(CPUARMState *env, float_status *statusp); /* * Negate as for FPCR.AH=1 -- do not negate NaNs. diff --git a/target/arm/tcg/sme_helper.c b/target/arm/tcg/sme_helper.c index XXXXXXX..XXXXXXX 100644 --- a/target/arm/tcg/sme_helper.c +++ b/target/arm/tcg/sme_helper.c @@ -XXX,XX +XXX,XX @@ static void do_bfmopa_w(void *vza, void *vzn, void *vzm, uint32_t desc, uint32_t negx, bool ah_neg) { intptr_t row, col, oprsz = simd_maxsz(desc); - float_status fpst, fpst_odd; + float_status fpst; - if (is_ebf(env, &fpst, &fpst_odd)) { + if (is_ebf(env, &fpst)) { for (row = 0; row < oprsz; ) { uint16_t prow = pn[H2(row >> 4)]; do { @@ -XXX,XX +XXX,XX @@ static void do_bfmopa_w(void *vza, void *vzn, void *vzm, uint32_t m = *(uint32_t *)(vzm + H1_4(col)); m = f16mop_adj_pair(m, pcol, 0); - *a = bfdotadd_ebf(*a, n, m, &fpst, &fpst_odd); + *a = bfdotadd_ebf(*a, n, m, &fpst); } col += 4; pcol >>= 4; diff --git a/target/arm/tcg/vec_helper.c b/target/arm/tcg/vec_helper.c index XXXXXXX..XXXXXXX 100644 --- a/target/arm/tcg/vec_helper.c +++ b/target/arm/tcg/vec_helper.c @@ -XXX,XX +XXX,XX @@ DO_MMLA_B(gvec_usmmla_b, do_usmmla_b) * BFloat16 Dot Product */ -bool is_ebf(CPUARMState *env, float_status *statusp, float_status *oddstatusp) +bool is_ebf(CPUARMState *env, float_status *statusp) { /* * For BFDOT, BFMMLA, etc, the behaviour depends on FPCR.EBF. @@ -XXX,XX +XXX,XX @@ bool is_ebf(CPUARMState *env, float_status *statusp, float_status *oddstatusp) *statusp = env->vfp.fp_status[is_a64(env) ? FPST_A64 : FPST_A32]; set_default_nan_mode(true, statusp); - if (ebf) { - /* EBF=1 needs to do a step with round-to-odd semantics */ - *oddstatusp = *statusp; - set_float_rounding_mode(float_round_to_odd, oddstatusp); - } else { + if (!ebf) { set_flush_to_zero(true, statusp); set_flush_inputs_to_zero(true, statusp); set_float_rounding_mode(float_round_to_odd_inf, statusp); @@ -XXX,XX +XXX,XX @@ float32 bfdotadd(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst) return t1; } -float32 bfdotadd_ebf(float32 sum, uint32_t e1, uint32_t e2, - float_status *fpst, float_status *fpst_odd) +float32 bfdotadd_ebf(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst) { /* Unpack two BFloat16 into two Float32, trivially. */ float32 s1r = e1 << 16; @@ -XXX,XX +XXX,XX @@ void HELPER(gvec_bfdot)(void *vd, void *vn, void *vm, void *va, intptr_t i, opr_sz = simd_oprsz(desc); float32 *d = vd, *a = va; uint32_t *n = vn, *m = vm; - float_status fpst, fpst_odd; + float_status fpst; - if (is_ebf(env, &fpst, &fpst_odd)) { + if (is_ebf(env, &fpst)) { for (i = 0; i < opr_sz / 4; ++i) { - d[i] = bfdotadd_ebf(a[i], n[i], m[i], &fpst, &fpst_odd); + d[i] = bfdotadd_ebf(a[i], n[i], m[i], &fpst); } } else { for (i = 0; i < opr_sz / 4; ++i) { @@ -XXX,XX +XXX,XX @@ void HELPER(gvec_bfdot_idx)(void *vd, void *vn, void *vm, intptr_t eltspersegment = MIN(16 / 4, elements); float32 *d = vd, *a = va; uint32_t *n = vn, *m = vm; - float_status fpst, fpst_odd; + float_status fpst; - if (is_ebf(env, &fpst, &fpst_odd)) { + if (is_ebf(env, &fpst)) { for (i = 0; i < elements; i += eltspersegment) { uint32_t m_idx = m[i + H4(index)]; for (j = i; j < i + eltspersegment; j++) { - d[j] = bfdotadd_ebf(a[j], n[j], m_idx, &fpst, &fpst_odd); + d[j] = bfdotadd_ebf(a[j], n[j], m_idx, &fpst); } } } else { @@ -XXX,XX +XXX,XX @@ void HELPER(sme2_bfvdot_idx)(void *vd, void *vn, void *vm, uint16_t *n0 = vn; uint16_t *n1 = vn + sizeof(ARMVectorReg); uint32_t *m = vm; - float_status fpst, fpst_odd; + float_status fpst; - if (is_ebf(env, &fpst, &fpst_odd)) { + if (is_ebf(env, &fpst)) { for (i = 0; i < elements; i += eltspersegment) { uint32_t m_idx = m[i + H4(idx)]; for (j = 0; j < eltspersegment; j++) { uint32_t nn = (n0[H2(2 * (i + j) + sel)]) | (n1[H2(2 * (i + j) + sel)] << 16); - d[i + H4(j)] = bfdotadd_ebf(a[i + H4(j)], nn, m_idx, - &fpst, &fpst_odd); + d[i + H4(j)] = bfdotadd_ebf(a[i + H4(j)], nn, m_idx, &fpst); } } } else { @@ -XXX,XX +XXX,XX @@ void HELPER(gvec_bfmmla)(void *vd, void *vn, void *vm, void *va, intptr_t s, opr_sz = simd_oprsz(desc); float32 *d = vd, *a = va; uint32_t *n = vn, *m = vm; - float_status fpst, fpst_odd; + float_status fpst; - if (is_ebf(env, &fpst, &fpst_odd)) { + if (is_ebf(env, &fpst)) { for (s = 0; s < opr_sz / 4; s += 4) { float32 sum00, sum01, sum10, sum11; @@ -XXX,XX +XXX,XX @@ void HELPER(gvec_bfmmla)(void *vd, void *vn, void *vm, void *va, * i j i k j k */ sum00 = a[s + H4(0 + 0)]; - sum00 = bfdotadd_ebf(sum00, n[s + H4(0 + 0)], m[s + H4(0 + 0)], &fpst, &fpst_odd); - sum00 = bfdotadd_ebf(sum00, n[s + H4(0 + 1)], m[s + H4(0 + 1)], &fpst, &fpst_odd); + sum00 = bfdotadd_ebf(sum00, n[s + H4(0 + 0)], m[s + H4(0 + 0)], &fpst); + sum00 = bfdotadd_ebf(sum00, n[s + H4(0 + 1)], m[s + H4(0 + 1)], &fpst); sum01 = a[s + H4(0 + 1)]; - sum01 = bfdotadd_ebf(sum01, n[s + H4(0 + 0)], m[s + H4(2 + 0)], &fpst, &fpst_odd); - sum01 = bfdotadd_ebf(sum01, n[s + H4(0 + 1)], m[s + H4(2 + 1)], &fpst, &fpst_odd); + sum01 = bfdotadd_ebf(sum01, n[s + H4(0 + 0)], m[s + H4(2 + 0)], &fpst); + sum01 = bfdotadd_ebf(sum01, n[s + H4(0 + 1)], m[s + H4(2 + 1)], &fpst); sum10 = a[s + H4(2 + 0)]; - sum10 = bfdotadd_ebf(sum10, n[s + H4(2 + 0)], m[s + H4(0 + 0)], &fpst, &fpst_odd); - sum10 = bfdotadd_ebf(sum10, n[s + H4(2 + 1)], m[s + H4(0 + 1)], &fpst, &fpst_odd); + sum10 = bfdotadd_ebf(sum10, n[s + H4(2 + 0)], m[s + H4(0 + 0)], &fpst); + sum10 = bfdotadd_ebf(sum10, n[s + H4(2 + 1)], m[s + H4(0 + 1)], &fpst); sum11 = a[s + H4(2 + 1)]; - sum11 = bfdotadd_ebf(sum11, n[s + H4(2 + 0)], m[s + H4(2 + 0)], &fpst, &fpst_odd); - sum11 = bfdotadd_ebf(sum11, n[s + H4(2 + 1)], m[s + H4(2 + 1)], &fpst, &fpst_odd); + sum11 = bfdotadd_ebf(sum11, n[s + H4(2 + 0)], m[s + H4(2 + 0)], &fpst); + sum11 = bfdotadd_ebf(sum11, n[s + H4(2 + 1)], m[s + H4(2 + 1)], &fpst); d[s + H4(0 + 0)] = sum00; d[s + H4(0 + 1)] = sum01; -- 2.43.0 Use softfloat-parts.h so that we can more naturally perform the required operations witha single rounding step. This happens to also simplify the NaN detection step. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- target/arm/tcg/sme_helper.c | 96 ++++++++++++++++--------------------- 1 file changed, 40 insertions(+), 56 deletions(-) diff --git a/target/arm/tcg/sme_helper.c b/target/arm/tcg/sme_helper.c index XXXXXXX..XXXXXXX 100644 --- a/target/arm/tcg/sme_helper.c +++ b/target/arm/tcg/sme_helper.c @@ -XXX,XX +XXX,XX @@ #include "accel/tcg/helper-retaddr.h" #include "qemu/int128.h" #include "fpu/softfloat.h" +#include "fpu/softfloat-parts.h" #include "vec_internal.h" #include "sve_ldst_internal.h" @@ -XXX,XX +XXX,XX @@ static inline uint32_t bf16mop_ah_neg_adj_pair(uint32_t pair, uint32_t pg) } static float32 f16_dotadd(float32 sum, uint32_t e1, uint32_t e2, - float_status *s_f16, float_status *s_std, - float_status *s_odd) + float_status *s_f16, float_status *s_std) { /* - * We need three different float_status for different parts of this + * We need two different float_status for different parts of this * operation: * - the input conversion of the float16 values must use the * f16-specific float_status, so that the FPCR.FZ16 control is applied * - operations on float32 including the final accumulation must use * the normal float_status, so that FPCR.FZ is applied - * - we have pre-set-up copy of s_std which is set to round-to-odd, - * for the multiply (see below) */ float16 h1r = e1 & 0xffff; float16 h1c = e1 >> 16; @@ -XXX,XX +XXX,XX @@ static float32 f16_dotadd(float32 sum, uint32_t e1, uint32_t e2, float16 h2c = e2 >> 16; float32 t32; + FloatParts64 p1r = float16_unpack_canonical(h1r, s_f16); + FloatParts64 p1c = float16_unpack_canonical(h1c, s_f16); + FloatParts64 p2r = float16_unpack_canonical(h2r, s_f16); + FloatParts64 p2c = float16_unpack_canonical(h2c, s_f16); + + int all_mask = (float_cmask(p1r.cls) | float_cmask(p1c.cls) | + float_cmask(p1r.cls) | float_cmask(p1c.cls)); + /* C.f. FPProcessNaNs4 */ - if (float16_is_any_nan(h1r) || float16_is_any_nan(h1c) || - float16_is_any_nan(h2r) || float16_is_any_nan(h2c)) { + if (unlikely(all_mask & float_cmask_anynan)) { float16 t16; - if (float16_is_signaling_nan(h1r, s_f16)) { - t16 = h1r; - } else if (float16_is_signaling_nan(h1c, s_f16)) { - t16 = h1c; - } else if (float16_is_signaling_nan(h2r, s_f16)) { - t16 = h2r; - } else if (float16_is_signaling_nan(h2c, s_f16)) { - t16 = h2c; - } else if (float16_is_any_nan(h1r)) { - t16 = h1r; - } else if (float16_is_any_nan(h1c)) { - t16 = h1c; - } else if (float16_is_any_nan(h2r)) { - t16 = h2r; + if (unlikely(all_mask & float_cmask_snan)) { + if (p1r.cls == float_class_snan) { + t16 = h1r; + } else if (p1c.cls == float_class_snan) { + t16 = h1c; + } else if (p2r.cls == float_class_snan) { + t16 = h2r; + } else { + t16 = h2c; + } } else { - t16 = h2c; + if (p1r.cls == float_class_qnan) { + t16 = h1r; + } else if (p1c.cls == float_class_qnan) { + t16 = h1c; + } else if (p2r.cls == float_class_qnan) { + t16 = h2r; + } else { + t16 = h2c; + } } t32 = float16_to_float32(t16, true, s_f16); } else { - float64 e1r = float16_to_float64(h1r, true, s_f16); - float64 e1c = float16_to_float64(h1c, true, s_f16); - float64 e2r = float16_to_float64(h2r, true, s_f16); - float64 e2c = float16_to_float64(h2c, true, s_f16); - float64 t64; - /* * The ARM pseudocode function FPDot performs both multiplies - * and the add with a single rounding operation. Emulate this - * by performing the first multiply in round-to-odd, then doing - * the second multiply as fused multiply-add, and rounding to - * float32 all in one step. + * and the add with a single rounding operation. */ - t64 = float64_mul(e1r, e2r, s_odd); - t64 = float64r32_muladd(e1c, e2c, t64, 0, s_std); + FloatParts64 tmp = parts64_mul(&p1r, &p2r, s_std); + tmp = parts64_muladd(&p1c, &p2c, &tmp, 0, s_std); - /* This conversion is exact, because we've already rounded. */ - t32 = float64_to_float32(t64, s_std); + t32 = float32_round_pack_canonical(&tmp, s_std); } /* The final accumulation step is not fused. */ @@ -XXX,XX +XXX,XX @@ static void do_fmopa_w_h(void *vza, void *vzn, void *vzm, uint16_t *pn, uint32_t negx, bool ah_neg) { intptr_t row, col, oprsz = simd_maxsz(desc); - float_status fpst_odd = env->vfp.fp_status[FPST_ZA]; - - set_float_rounding_mode(float_round_to_odd, &fpst_odd); for (row = 0; row < oprsz; ) { uint16_t prow = pn[H2(row >> 4)]; @@ -XXX,XX +XXX,XX @@ static void do_fmopa_w_h(void *vza, void *vzn, void *vzm, uint16_t *pn, m = f16mop_adj_pair(m, pcol, 0); *a = f16_dotadd(*a, n, m, &env->vfp.fp_status[FPST_ZA_F16], - &env->vfp.fp_status[FPST_ZA], - &fpst_odd); + &env->vfp.fp_status[FPST_ZA]); } col += 4; pcol >>= 4; @@ -XXX,XX +XXX,XX @@ void HELPER(sme2_fdot_h)(void *vd, void *vn, void *vm, void *va, bool za = extract32(desc, SIMD_DATA_SHIFT, 1); float_status *fpst_std = &env->vfp.fp_status[za ? FPST_ZA : FPST_A64]; float_status *fpst_f16 = &env->vfp.fp_status[za ? FPST_ZA_F16 : FPST_A64_F16]; - float_status fpst_odd = *fpst_std; float32 *d = vd, *a = va; uint32_t *n = vn, *m = vm; - set_float_rounding_mode(float_round_to_odd, &fpst_odd); - for (i = 0; i < oprsz / sizeof(float32); ++i) { d[H4(i)] = f16_dotadd(a[H4(i)], n[H4(i)], m[H4(i)], - fpst_f16, fpst_std, &fpst_odd); + fpst_f16, fpst_std); } } @@ -XXX,XX +XXX,XX @@ void HELPER(sme2_fdot_idx_h)(void *vd, void *vn, void *vm, void *va, bool za = extract32(desc, SIMD_DATA_SHIFT + 2, 1); float_status *fpst_std = &env->vfp.fp_status[za ? FPST_ZA : FPST_A64]; float_status *fpst_f16 = &env->vfp.fp_status[za ? FPST_ZA_F16 : FPST_A64_F16]; - float_status fpst_odd = *fpst_std; float32 *d = vd, *a = va; uint32_t *n = vn, *m = (uint32_t *)vm + H4(idx); - set_float_rounding_mode(float_round_to_odd, &fpst_odd); - for (i = 0; i < elements; i += eltspersegment) { uint32_t mm = m[i]; for (j = 0; j < eltspersegment; ++j) { d[H4(i + j)] = f16_dotadd(a[H4(i + j)], n[H4(i + j)], mm, - fpst_f16, fpst_std, &fpst_odd); + fpst_f16, fpst_std); } } } @@ -XXX,XX +XXX,XX @@ void HELPER(sme2_fvdot_idx_h)(void *vd, void *vn, void *vm, void *va, intptr_t eltspersegment = MIN(4, elements); int idx = extract32(desc, SIMD_DATA_SHIFT, 2); int sel = extract32(desc, SIMD_DATA_SHIFT + 2, 1); - float_status fpst_odd, *fpst_std, *fpst_f16; float32 *d = vd, *a = va; uint16_t *n0 = vn; uint16_t *n1 = vn + sizeof(ARMVectorReg); uint32_t *m = (uint32_t *)vm + H4(idx); - fpst_std = &env->vfp.fp_status[FPST_ZA]; - fpst_f16 = &env->vfp.fp_status[FPST_ZA_F16]; - fpst_odd = *fpst_std; - set_float_rounding_mode(float_round_to_odd, &fpst_odd); - for (i = 0; i < elements; i += eltspersegment) { uint32_t mm = m[i]; for (j = 0; j < eltspersegment; ++j) { uint32_t nn = (n0[H2(2 * (i + j) + sel)]) | (n1[H2(2 * (i + j) + sel)] << 16); d[i + H4(j)] = f16_dotadd(a[i + H4(j)], nn, mm, - fpst_f16, fpst_std, &fpst_odd); + &env->vfp.fp_status[FPST_ZA_F16], + &env->vfp.fp_status[FPST_ZA]); } } } -- 2.43.0 Notice the conflict between saturation and rebias_overflow, which should never happen. Otherwise, saturate to INT32_MAX and allow the usual overflow to max, infinity, or dnan. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat-parts.c.inc | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ static void partsN(uncanon_normal)(FloatPartsN *p, float_status *s, g_assert_not_reached(); } - exp = p->exp + fmt->exp_bias; + /* Because exp_bias is positive, we can only overflow past INT_MAX. */ + if (sadd32_overflow(p->exp, fmt->exp_bias, &exp)) { + /* + * rebias_overflow wants to compute a modulo exponent, which + * conflicts with saturation. That said, saturation can only + * happen with scalbn, which is not a PowerPC operation. + */ + assert(!s->rebias_overflow); + exp = INT32_MAX; + } + if (likely(exp > 0)) { if (p->frac_lo & round_mask) { flags |= float_flag_inexact; -- 2.43.0 Avoid exponent overflow as well as checking that we don't lose information with opposing scaling. Use it in partsN(scalbn) and partsN(round_to_int_normal). Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat.c | 17 +++++++++++++++++ fpu/softfloat-parts.c.inc | 5 ++--- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/fpu/softfloat.c b/fpu/softfloat.c index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -XXX,XX +XXX,XX @@ static float128 QEMU_FLATTEN float128_pack_raw(const FloatParts128 *p) *----------------------------------------------------------------------------*/ #include "softfloat-specialize.c.inc" +static int32_t exp_scalbn(int32_t exp, int32_t scale) +{ + /* + * Catch chains of scaling which lose information. + * In particular, if the exponent has been saturated, + * do not allow it to become unsaturated. + */ + if (unlikely(exp == INT32_MAX)) { + assert(scale >= 0); + } else if (unlikely(exp == INT32_MIN)) { + assert(scale <= 0); + } else { + exp = sadd32_saturate(exp, scale); + } + return exp; +} + /* * Helper functions for softfloat-parts.c.inc, per-size operations. */ diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ static bool partsN(round_to_int_normal)(FloatPartsN *a, FloatRoundMode rmode, uint64_t frac_lsb, frac_lsbm1, rnd_even_mask, rnd_mask, inc; int shift_adj; - scale = MIN(MAX(scale, -0x10000), 0x10000); - a->exp += scale; + a->exp = exp_scalbn(a->exp, scale); if (a->exp < 0) { bool one; @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(scalbn)(const FloatPartsN *a, int n, float_status *s) case float_class_normal: { FloatPartsN r = *a; - r.exp += MIN(MAX(n, -0x10000), 0x10000); + r.exp = exp_scalbn(r.exp, n); return r; } default: -- 2.43.0 Fixes for scalbn and muladd. r~ Richard Henderson (12): fpu: Return struct from parts{64,128}_scalbn fpu: Reorganize partsN(muladd) fpu: Return struct from parts{64,128}_muladd fpu: Hoist nan check in partsN_addsub fpu: Simplify 0 +/- N case in parts_addsub fpu: Use parts64_round_to_int in parts_s390_divide_to_integer target/s390x: Move float{32,64}_s390_divide_to_integer target/arm: Use FloatParts64 in bfdotadd_ebf target/arm: Drop oddstatus from is_ebf and bfdotadd_ebf target/arm: Use FloatParts64 in f16_dotadd fpu: Saturate the exponent in uncanon_normal fpu: Introduce exp_scalbn include/fpu/softfloat-parts.h | 12 ++ include/fpu/softfloat.h | 11 -- target/arm/tcg/vec_internal.h | 12 +- fpu/softfloat.c | 236 ++++++++-------------------------- target/arm/tcg/sme_helper.c | 102 +++++++-------- target/arm/tcg/vec_helper.c | 127 +++++++++--------- target/s390x/tcg/fpu_helper.c | 135 +++++++++++++++++++ fpu/softfloat-parts.c.inc | 236 ++++++++++++++++------------------ 8 files changed, 424 insertions(+), 447 deletions(-) -- 2.43.0 At the same time, export. Reviewed-by: Philippe Mathieu-Daudé <philmd@linaro.org> Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- include/fpu/softfloat-parts.h | 3 +++ fpu/softfloat.c | 18 +++++++++--------- fpu/softfloat-parts.c.inc | 14 ++++++++------ 3 files changed, 20 insertions(+), 15 deletions(-) diff --git a/include/fpu/softfloat-parts.h b/include/fpu/softfloat-parts.h index XXXXXXX..XXXXXXX 100644 --- a/include/fpu/softfloat-parts.h +++ b/include/fpu/softfloat-parts.h @@ -XXX,XX +XXX,XX @@ FloatParts128 parts128_round_to_int(const FloatParts128 *a, FloatParts64 parts64_round_to_fmt(const FloatParts64 *p, float_status *s, const FloatFmt *fmt); +FloatParts64 parts64_scalbn(const FloatParts64 *a, int n, float_status *s); +FloatParts128 parts128_scalbn(const FloatParts128 *a, int n, float_status *s); + #endif diff --git a/fpu/softfloat.c b/fpu/softfloat.c index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -XXX,XX +XXX,XX @@ float16 float16_muladd_scalbn(float16 a, float16 b, float16 c, /* Before rounding, scale. */ if (scale) { - parts64_scalbn(pr, scale, status); + *pr = parts64_scalbn(pr, scale, status); } parts64_uncanon(pr, status, &float16_params, false); /* After rounding, apply negate result, especially for -0.0. */ @@ -XXX,XX +XXX,XX @@ float32_muladd_scalbn(float32 a, float32 b, float32 c, /* Before rounding, scale. */ if (scale) { - parts64_scalbn(pr, scale, status); + *pr = parts64_scalbn(pr, scale, status); } parts64_uncanon(pr, status, &float32_params, false); /* After rounding, apply negate result, especially for -0.0. */ @@ -XXX,XX +XXX,XX @@ float64_muladd_scalbn(float64 a, float64 b, float64 c, /* Before rounding, scale. */ if (scale) { - parts64_scalbn(pr, scale, status); + *pr = parts64_scalbn(pr, scale, status); } parts64_uncanon(pr, status, &float64_params, false); /* After rounding, apply negate result, especially for -0.0. */ @@ -XXX,XX +XXX,XX @@ float16 float16_scalbn(float16 a, int n, float_status *status) { FloatParts64 p = float16_unpack_canonical(a, status); - parts64_scalbn(&p, n, status); + p = parts64_scalbn(&p, n, status); return float16_round_pack_canonical(&p, status); } @@ -XXX,XX +XXX,XX @@ float32 float32_scalbn(float32 a, int n, float_status *status) { FloatParts64 p = float32_unpack_canonical(a, status); - parts64_scalbn(&p, n, status); + p = parts64_scalbn(&p, n, status); return float32_round_pack_canonical(&p, status); } @@ -XXX,XX +XXX,XX @@ float64 float64_scalbn(float64 a, int n, float_status *status) { FloatParts64 p = float64_unpack_canonical(a, status); - parts64_scalbn(&p, n, status); + p = parts64_scalbn(&p, n, status); return float64_round_pack_canonical(&p, status); } @@ -XXX,XX +XXX,XX @@ bfloat16 bfloat16_scalbn(bfloat16 a, int n, float_status *status) { FloatParts64 p = bfloat16_unpack_canonical(a, status); - parts64_scalbn(&p, n, status); + p = parts64_scalbn(&p, n, status); return bfloat16_round_pack_canonical(&p, status); } @@ -XXX,XX +XXX,XX @@ float128 float128_scalbn(float128 a, int n, float_status *status) { FloatParts128 p = float128_unpack_canonical(a, status); - parts128_scalbn(&p, n, status); + p = parts128_scalbn(&p, n, status); return float128_round_pack_canonical(&p, status); } @@ -XXX,XX +XXX,XX @@ floatx80 floatx80_scalbn(floatx80 a, int n, float_status *status) if (!floatx80_unpack_canonical(&p, a, status)) { return floatx80_default_nan(status); } - parts128_scalbn(&p, n, status); + p = parts128_scalbn(&p, n, status); return floatx80_round_pack_canonical(&p, status); } diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ FloatRelation partsN(compare)(const FloatPartsN *a, const FloatPartsN *b, /* * Multiply A by 2 raised to the power N. */ -static void partsN(scalbn)(FloatPartsN *a, int n, float_status *s) +FloatPartsN partsN(scalbn)(const FloatPartsN *a, int n, float_status *s) { switch (a->cls) { case float_class_snan: case float_class_qnan: - *a = partsN(return_nan)(a, s); - break; + return partsN(return_nan)(a, s); case float_class_zero: case float_class_inf: - break; + return *a; case float_class_denormal: float_raise(float_flag_input_denormal_used, s); /* fall through */ case float_class_normal: - a->exp += MIN(MAX(n, -0x10000), 0x10000); - break; + { + FloatPartsN r = *a; + r.exp += MIN(MAX(n, -0x10000), 0x10000); + return r; + } default: g_assert_not_reached(); } -- 2.43.0 Check the likely case of normal product and normal or zero addend first; shift NaN and infinity detection down; end with zero product + addend. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat-parts.c.inc | 156 +++++++++++++++++--------------------- 1 file changed, 71 insertions(+), 85 deletions(-) diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, FloatPartsN *c, int flags, float_status *s) { - int ab_mask, abc_mask; - FloatPartsW p_widen, c_widen; + int ab_mask = float_cmask(a->cls) | float_cmask(b->cls); + int c_mask = float_cmask(c->cls); + int abc_mask = ab_mask | c_mask; + bool c_sign = c->sign ^ !!(flags & float_muladd_negate_c); + bool p_sign = a->sign ^ b->sign ^ !!(flags & float_muladd_negate_product); - ab_mask = float_cmask(a->cls) | float_cmask(b->cls); - abc_mask = float_cmask(c->cls) | ab_mask; + /* + * The "likely" case is A and B normal, so that the product is normal, + * and C normal or zero so that the result is normal. + */ + int likely_mask = ab_mask | (c_mask & ~float_cmask_zero); + if (likely(cmask_is_only_normals(likely_mask))) { + record_denormals_used(abc_mask, s); + + /* Perform the multiplication step. */ + FloatPartsW p_widen = { .sign = p_sign, .exp = a->exp + b->exp + 1 }; + fracN(mulw)(&p_widen, a, b); + if (!(p_widen.frac_hi & DECOMPOSED_IMPLICIT_BIT)) { + fracW(add)(&p_widen, &p_widen, &p_widen); + p_widen.exp -= 1; + } + + /* Perform the addition step. */ + if (!(c_mask & float_cmask_zero)) { + /* Zero-extend C to less significant bits. */ + FloatPartsW c_widen = { .sign = c_sign, .exp = c->exp }; + fracN(widen)(&c_widen, c); + + if (p_sign == c_sign) { + partsW(add_normal)(&p_widen, &c_widen); + } else if (!partsW(sub_normal)(&p_widen, &c_widen)) { + goto return_sub_zero; + } + } + + /* Narrow with sticky bit, for proper rounding later. */ + fracN(truncjam)(a, &p_widen); + a->sign = p_widen.sign; + a->exp = p_widen.exp; + return a; + } /* * It is implementation-defined whether the cases of (0,inf,qnan) @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, return a; } - if (flags & float_muladd_negate_c) { - c->sign ^= 1; + if (unlikely(ab_mask == float_cmask_infzero)) { + /* Inf * Zero == NaN */ + float_raise(float_flag_invalid | float_flag_invalid_imz, s); + goto d_nan; } - /* Compute the sign of the product into A. */ - a->sign ^= b->sign; - if (flags & float_muladd_negate_product) { - a->sign ^= 1; - } - - if (unlikely(!cmask_is_only_normals(ab_mask))) { - if (unlikely(ab_mask == float_cmask_infzero)) { - float_raise(float_flag_invalid | float_flag_invalid_imz, s); + if (unlikely(ab_mask & float_cmask_inf)) { + if ((c_mask & float_cmask_inf) && p_sign != c_sign) { + /* Inf - Inf == NaN */ + float_raise(float_flag_invalid | float_flag_invalid_isi, s); goto d_nan; } - - if (ab_mask & float_cmask_inf) { - if (c->cls == float_class_inf && a->sign != c->sign) { - float_raise(float_flag_invalid | float_flag_invalid_isi, s); - goto d_nan; - } - goto return_inf; - } - - g_assert(ab_mask & float_cmask_zero); - if (is_anynorm(c->cls)) { - *a = *c; - goto finish_sign; - } - if (c->cls == float_class_zero) { - if (flags & float_muladd_suppress_add_product_zero) { - a->sign = c->sign; - } else if (a->sign != c->sign) { - goto return_sub_zero; - } - goto return_zero; - } - g_assert(c->cls == float_class_inf); + /* Inf + C == Inf */ + record_denormals_used(abc_mask, s); + a->sign = p_sign; + a->cls = float_class_inf; + return a; } - - if (unlikely(c->cls == float_class_inf)) { - a->sign = c->sign; - goto return_inf; - } - - /* Perform the multiplication step. */ - p_widen.sign = a->sign; - p_widen.exp = a->exp + b->exp + 1; - fracN(mulw)(&p_widen, a, b); - if (!(p_widen.frac_hi & DECOMPOSED_IMPLICIT_BIT)) { - fracW(add)(&p_widen, &p_widen, &p_widen); - p_widen.exp -= 1; - } - - /* Perform the addition step. */ - if (c->cls != float_class_zero) { - /* Zero-extend C to less significant bits. */ - fracN(widen)(&c_widen, c); - c_widen.exp = c->exp; - - if (a->sign == c->sign) { - partsW(add_normal)(&p_widen, &c_widen); - } else if (!partsW(sub_normal)(&p_widen, &c_widen)) { - goto return_sub_zero; - } - } - - /* Narrow with sticky bit, for proper rounding later. */ - fracN(truncjam)(a, &p_widen); - a->sign = p_widen.sign; - a->exp = p_widen.exp; - - finish_sign: - /* - * All result types except for "return the default NaN - * because this is an Invalid Operation" go through here; - * this matches the set of cases where we consumed a - * denormal input. - */ record_denormals_used(abc_mask, s); - return a; + + /* Only remaining cases are zero product or inf addend. */ + assert((ab_mask & float_cmask_zero) | (c_mask & float_cmask_inf)); + + /* + * P + Inf == Inf, or + * 0 + C == C, + * except for 0 - 0, which needs special rounding, + * except for when we want to suppress this addition step. + */ + if (!(c_mask & float_cmask_zero) + || p_sign == c_sign + || (flags & float_muladd_suppress_add_product_zero)) { + c->sign = c_sign; + return c; + } return_sub_zero: + /* 0 - 0 == -0 for round_down, +0 otherwise. */ a->sign = s->float_rounding_mode == float_round_down; - return_zero: a->cls = float_class_zero; - goto finish_sign; - - return_inf: - a->cls = float_class_inf; - goto finish_sign; + return a; d_nan: *a = partsN(default_nan)(s); -- 2.43.0 At the same time, export. Reviewed-by: Philippe Mathieu-Daudé <philmd@linaro.org> Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- include/fpu/softfloat-parts.h | 9 ++++ fpu/softfloat.c | 83 +++++++++++++++++------------------ fpu/softfloat-parts.c.inc | 38 ++++++++-------- 3 files changed, 69 insertions(+), 61 deletions(-) diff --git a/include/fpu/softfloat-parts.h b/include/fpu/softfloat-parts.h index XXXXXXX..XXXXXXX 100644 --- a/include/fpu/softfloat-parts.h +++ b/include/fpu/softfloat-parts.h @@ -XXX,XX +XXX,XX @@ FloatParts64 parts64_mul(const FloatParts64 *a, const FloatParts64 *b, FloatParts128 parts128_mul(const FloatParts128 *a, const FloatParts128 *b, float_status *s); +FloatParts64 parts64_muladd(const FloatParts64 *a, + const FloatParts64 *b, + const FloatParts64 *c, + int flags, float_status *s); +FloatParts128 parts128_muladd(const FloatParts128 *a, + const FloatParts128 *b, + const FloatParts128 *c, + int flags, float_status *s); + FloatParts64 parts64_round_to_int(const FloatParts64 *a, FloatRoundMode rmode, int scale, float_status *s, diff --git a/fpu/softfloat.c b/fpu/softfloat.c index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -XXX,XX +XXX,XX @@ float16 float16_muladd_scalbn(float16 a, float16 b, float16 c, FloatParts64 pa = float16_unpack_canonical(a, status); FloatParts64 pb = float16_unpack_canonical(b, status); FloatParts64 pc = float16_unpack_canonical(c, status); - FloatParts64 *pr = parts64_muladd(&pa, &pb, &pc, flags, status); + FloatParts64 pr = parts64_muladd(&pa, &pb, &pc, flags, status); /* Before rounding, scale. */ if (scale) { - *pr = parts64_scalbn(pr, scale, status); + pr = parts64_scalbn(&pr, scale, status); } - parts64_uncanon(pr, status, &float16_params, false); + parts64_uncanon(&pr, status, &float16_params, false); /* After rounding, apply negate result, especially for -0.0. */ - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return pack_raw64(pr, &float16_params); + return pack_raw64(&pr, &float16_params); } float16 float16_muladd(float16 a, float16 b, float16 c, @@ -XXX,XX +XXX,XX @@ float32_muladd_scalbn(float32 a, float32 b, float32 c, FloatParts64 pa = float32_unpack_canonical(a, status); FloatParts64 pb = float32_unpack_canonical(b, status); FloatParts64 pc = float32_unpack_canonical(c, status); - FloatParts64 *pr = parts64_muladd(&pa, &pb, &pc, flags, status); + FloatParts64 pr = parts64_muladd(&pa, &pb, &pc, flags, status); /* Before rounding, scale. */ if (scale) { - *pr = parts64_scalbn(pr, scale, status); + pr = parts64_scalbn(&pr, scale, status); } - parts64_uncanon(pr, status, &float32_params, false); + parts64_uncanon(&pr, status, &float32_params, false); /* After rounding, apply negate result, especially for -0.0. */ - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return pack_raw64(pr, &float32_params); + return pack_raw64(&pr, &float32_params); } float64 QEMU_SOFTFLOAT_ATTR @@ -XXX,XX +XXX,XX @@ float64_muladd_scalbn(float64 a, float64 b, float64 c, FloatParts64 pa = float64_unpack_canonical(a, status); FloatParts64 pb = float64_unpack_canonical(b, status); FloatParts64 pc = float64_unpack_canonical(c, status); - FloatParts64 *pr = parts64_muladd(&pa, &pb, &pc, flags, status); + FloatParts64 pr = parts64_muladd(&pa, &pb, &pc, flags, status); /* Before rounding, scale. */ if (scale) { - *pr = parts64_scalbn(pr, scale, status); + pr = parts64_scalbn(&pr, scale, status); } - parts64_uncanon(pr, status, &float64_params, false); + parts64_uncanon(&pr, status, &float64_params, false); /* After rounding, apply negate result, especially for -0.0. */ - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return pack_raw64(pr, &float64_params); + return pack_raw64(&pr, &float64_params); } static bool force_soft_fma; @@ -XXX,XX +XXX,XX @@ float64 float64r32_muladd(float64 a, float64 b, float64 c, FloatParts64 pa = float64_unpack_canonical(a, status); FloatParts64 pb = float64_unpack_canonical(b, status); FloatParts64 pc = float64_unpack_canonical(c, status); - FloatParts64 *pr = parts64_muladd(&pa, &pb, &pc, flags, status); + FloatParts64 pr = parts64_muladd(&pa, &pb, &pc, flags, status); /* Round before applying negate result. */ - parts64_uncanon(pr, status, &float32_params, false); - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + parts64_uncanon(&pr, status, &float32_params, false); + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return float64r32_pack_raw(pr); + return float64r32_pack_raw(&pr); } bfloat16 bfloat16_muladd(bfloat16 a, bfloat16 b, bfloat16 c, @@ -XXX,XX +XXX,XX @@ bfloat16 bfloat16_muladd(bfloat16 a, bfloat16 b, bfloat16 c, FloatParts64 pa = bfloat16_unpack_canonical(a, status); FloatParts64 pb = bfloat16_unpack_canonical(b, status); FloatParts64 pc = bfloat16_unpack_canonical(c, status); - FloatParts64 *pr = parts64_muladd(&pa, &pb, &pc, flags, status); + FloatParts64 pr = parts64_muladd(&pa, &pb, &pc, flags, status); /* Round before applying negate result. */ - parts64_uncanon(pr, status, &bfloat16_params, false); - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + parts64_uncanon(&pr, status, &bfloat16_params, false); + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return pack_raw64(pr, &bfloat16_params); + return pack_raw64(&pr, &bfloat16_params); } float128 float128_muladd(float128 a, float128 b, float128 c, @@ -XXX,XX +XXX,XX @@ float128 float128_muladd(float128 a, float128 b, float128 c, FloatParts128 pa = float128_unpack_canonical(a, status); FloatParts128 pb = float128_unpack_canonical(b, status); FloatParts128 pc = float128_unpack_canonical(c, status); - FloatParts128 *pr = parts128_muladd(&pa, &pb, &pc, flags, status); + FloatParts128 pr = parts128_muladd(&pa, &pb, &pc, flags, status); /* Round before applying negate result. */ - parts128_uncanon(pr, status, &float128_params, false); - if ((flags & float_muladd_negate_result) && !is_nan(pr->cls)) { - pr->sign ^= 1; + parts128_uncanon(&pr, status, &float128_params, false); + if ((flags & float_muladd_negate_result) && !is_nan(pr.cls)) { + pr.sign ^= 1; } - return float128_pack_raw(pr); + return float128_pack_raw(&pr); } /* @@ -XXX,XX +XXX,XX @@ float32 float32_exp2(float32 a, float_status *status) rp = float64_unpack_canonical(float64_one, status); for (int i = 0; i < 15; i++) { tp = float64_unpack_canonical(float32_exp2_coefficients[i], status); - rp = *parts64_muladd(&tp, &xnp, &rp, 0, status); + rp = parts64_muladd(&tp, &xnp, &rp, 0, status); xnp = parts64_mul(&xnp, &xp, status); } @@ -XXX,XX +XXX,XX @@ static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, n->sign = a->sign ^ b->sign; *cc = 0; } else { - FloatParts64 *q, q_buf, *r_precise, r_precise_buf; + FloatParts64 *q, q_buf, r_precise; int float_exception_flags = 0; bool is_q_smallish; uint32_t r_flags; @@ -XXX,XX +XXX,XX @@ static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, 0, fmt->frac_size); /* Compute precise remainder */ - r_precise_buf = *b; - r_precise = parts64_muladd(&r_precise_buf, n, a, + r_precise = parts64_muladd(b, n, a, float_muladd_negate_product, status); /* Round remainder to the target format */ - *r = *r_precise; + *r = r_precise; status->float_exception_flags = 0; *r = parts64_round_to_fmt(r, status, fmt); r_flags = status->float_exception_flags; @@ -XXX,XX +XXX,XX @@ static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, * toward zero) or incremented. */ saved_r_sign = r->sign; - saved_r_precise_sign = r_precise->sign; + saved_r_precise_sign = r_precise.sign; r->sign = false; - r_precise->sign = false; - if (parts64_compare(r, r_precise, status, true) < + r_precise.sign = false; + if (parts64_compare(r, &r_precise, status, true) < float_relation_equal) { *dxc = 0x8; } else { *dxc = 0xc; } r->sign = saved_r_sign; - r_precise->sign = saved_r_precise_sign; + r_precise.sign = saved_r_precise_sign; } } } diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(mul)(const FloatPartsN *a, const FloatPartsN *b, * Requires A and C extracted into a double-sized structure to provide the * extra space for the widening multiply. */ -static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, - FloatPartsN *c, - int flags, float_status *s) +FloatPartsN partsN(muladd)(const FloatPartsN *a, const FloatPartsN *b, + const FloatPartsN *c, int flags, float_status *s) { int ab_mask = float_cmask(a->cls) | float_cmask(b->cls); int c_mask = float_cmask(c->cls); @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, } /* Narrow with sticky bit, for proper rounding later. */ - fracN(truncjam)(a, &p_widen); - a->sign = p_widen.sign; - a->exp = p_widen.exp; - return a; + FloatPartsN r = { + .sign = p_widen.sign, + .exp = p_widen.exp, + .cls = float_class_normal, + }; + fracN(truncjam)(&r, &p_widen); + return r; } /* @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, * off to the target-specific pick-a-NaN routine. */ if (unlikely(abc_mask & float_cmask_anynan)) { - *a = partsN(pick_nan_muladd)(a, b, c, s, ab_mask, abc_mask); - return a; + return partsN(pick_nan_muladd)(a, b, c, s, ab_mask, abc_mask); } if (unlikely(ab_mask == float_cmask_infzero)) { @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, } /* Inf + C == Inf */ record_denormals_used(abc_mask, s); - a->sign = p_sign; - a->cls = float_class_inf; - return a; + return (FloatPartsN){ .sign = p_sign, .cls = float_class_inf }; } record_denormals_used(abc_mask, s); @@ -XXX,XX +XXX,XX @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b, if (!(c_mask & float_cmask_zero) || p_sign == c_sign || (flags & float_muladd_suppress_add_product_zero)) { - c->sign = c_sign; - return c; + FloatPartsN r = *c; + r.sign = c_sign; + return r; } return_sub_zero: /* 0 - 0 == -0 for round_down, +0 otherwise. */ - a->sign = s->float_rounding_mode == float_round_down; - a->cls = float_class_zero; - return a; + return (FloatPartsN){ + .sign = s->float_rounding_mode == float_round_down, + .cls = float_class_zero + }; d_nan: - *a = partsN(default_nan)(s); - return a; + return partsN(default_nan)(s); } /* -- 2.43.0 The nan test had been down below because it was unlikely. But if we have to have one anyway because of denormals, we might as well take care of them right away. Reviewed-by: Philippe Mathieu-Daudé <philmd@linaro.org> Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat-parts.c.inc | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, { int ab_mask = float_cmask(a_orig->cls) | float_cmask(b_orig->cls); + if (unlikely(ab_mask & float_cmask_anynan)) { + return partsN(pick_nan)(a_orig, b_orig, s); + } + /* * For addition and subtraction, we will consume an * input denormal unless the other input is a NaN. */ - if (!(ab_mask & float_cmask_anynan)) { - record_denormals_used(ab_mask, s); - } + record_denormals_used(ab_mask, s); FloatPartsN a = *a_orig; FloatPartsN b = *b_orig; @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, return a; } - if (unlikely(ab_mask & float_cmask_anynan)) { - goto p_nan; - } - if (ab_mask & float_cmask_inf) { if (a.cls != float_class_inf) { /* N - Inf */ @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, return a; } - if (unlikely(ab_mask & float_cmask_anynan)) { - goto p_nan; - } - if (ab_mask & float_cmask_inf) { a.cls = float_class_inf; return a; @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, g_assert(a.cls == float_class_zero); g_assert(is_anynorm(b.cls)); return b; - - p_nan: - return partsN(pick_nan)(a_orig, b_orig, s); } /* -- 2.43.0 Consolidate the tests for zero and anynorm. Add comments for a few cases. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat-parts.c.inc | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, } if (ab_mask == float_cmask_zero) { + /* 0 - 0 */ a.sign = s->float_rounding_mode == float_round_down; return a; } @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(addsub)(const FloatPartsN *a_orig, } if (ab_mask == float_cmask_zero) { + /* 0 + 0 */ return a; } if (ab_mask & float_cmask_inf) { + /* N + Inf or Inf + N */ a.cls = float_class_inf; return a; } } - if (b.cls == float_class_zero) { - g_assert(is_anynorm(a.cls)); - return a; - } - - g_assert(a.cls == float_class_zero); - g_assert(is_anynorm(b.cls)); - return b; + /* 0 +/- N or N +/- 0 */ + assert((ab_mask & float_cmask_zero) && (ab_mask & float_cmask_anynorm)); + return b.cls == float_class_zero ? a : b; } /* -- 2.43.0 We will not expose parts_round_to_int_normal. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat.c | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/fpu/softfloat.c b/fpu/softfloat.c index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -XXX,XX +XXX,XX @@ static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, * Rounding of partial quotient may be inexact. This is the whole point * of distinguishing partial quotients, so ignore the exception. */ - *n = *q; - parts64_round_to_int_normal(n, - is_q_smallish - ? final_quotient_rounding_mode - : float_round_to_zero, - 0, fmt->frac_size); + *n = parts64_round_to_int(q, + is_q_smallish + ? final_quotient_rounding_mode + : float_round_to_zero, + 0, status, fmt); /* Compute precise remainder */ r_precise = parts64_muladd(b, n, a, -- 2.43.0 Now that we've exposed enough infrastructure, this can be implemented in the backend that needs it. Reviewed-by: Philippe Mathieu-Daudé <philmd@linaro.org> Reviewed-by: Ilya Leoshkevich <iii@linux.ibm.com> Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- include/fpu/softfloat.h | 11 --- fpu/softfloat.c | 137 ---------------------------------- target/s390x/tcg/fpu_helper.c | 135 +++++++++++++++++++++++++++++++++ 3 files changed, 135 insertions(+), 148 deletions(-) diff --git a/include/fpu/softfloat.h b/include/fpu/softfloat.h index XXXXXXX..XXXXXXX 100644 --- a/include/fpu/softfloat.h +++ b/include/fpu/softfloat.h @@ -XXX,XX +XXX,XX @@ static inline bool float128_unordered_quiet(float128 a, float128 b, *----------------------------------------------------------------------------*/ float128 float128_default_nan(float_status *status); -#define DECLARE_S390_DIVIDE_TO_INTEGER(floatN) \ -void floatN ## _s390_divide_to_integer(floatN a, floatN b, \ - int final_quotient_rounding_mode, \ - bool mask_underflow, bool mask_inexact, \ - floatN *r, floatN *n, \ - uint32_t *cc, int *dxc, \ - float_status *status) -DECLARE_S390_DIVIDE_TO_INTEGER(float32); -DECLARE_S390_DIVIDE_TO_INTEGER(float64); - - #endif /* SOFTFLOAT_H */ diff --git a/fpu/softfloat.c b/fpu/softfloat.c index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -XXX,XX +XXX,XX @@ floatx80 floatx80_round(floatx80 a, float_status *status) return floatx80_round_pack_canonical(&p, status); } -static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, - int final_quotient_rounding_mode, - bool mask_underflow, bool mask_inexact, - const FloatFmt *fmt, - FloatParts64 *r, FloatParts64 *n, - uint32_t *cc, int *dxc, - float_status *status) -{ - /* POp table "Results: DIVIDE TO INTEGER (Part 1 of 2)" */ - if ((float_cmask(a->cls) | float_cmask(b->cls)) & float_cmask_anynan) { - *r = parts64_pick_nan(a, b, status); - *n = *r; - *cc = 1; - } else if (a->cls == float_class_inf || b->cls == float_class_zero) { - *r = parts64_default_nan(status); - *n = *r; - *cc = 1; - status->float_exception_flags |= float_flag_invalid; - } else if (b->cls == float_class_inf) { - *r = *a; - n->cls = float_class_zero; - n->sign = a->sign ^ b->sign; - *cc = 0; - } else { - FloatParts64 *q, q_buf, r_precise; - int float_exception_flags = 0; - bool is_q_smallish; - uint32_t r_flags; - - /* Compute precise quotient */ - q_buf = parts64_div(a, b, status); - q = &q_buf; - - /* - * Check whether two closest integers can be precisely represented, - * i.e., all their bits fit into the fractional part. - */ - is_q_smallish = q->exp < (fmt->frac_size + 1); - - /* - * Final quotient is rounded using final-quotient-rounding method, and - * partial quotient is rounded toward zero. - * - * Rounding of partial quotient may be inexact. This is the whole point - * of distinguishing partial quotients, so ignore the exception. - */ - *n = parts64_round_to_int(q, - is_q_smallish - ? final_quotient_rounding_mode - : float_round_to_zero, - 0, status, fmt); - - /* Compute precise remainder */ - r_precise = parts64_muladd(b, n, a, - float_muladd_negate_product, status); - - /* Round remainder to the target format */ - *r = r_precise; - status->float_exception_flags = 0; - *r = parts64_round_to_fmt(r, status, fmt); - r_flags = status->float_exception_flags; - - /* POp table "Results: DIVIDE TO INTEGER (Part 2 of 2)" */ - if (is_q_smallish) { - if (r->cls != float_class_zero) { - if (r->exp < 2 - (1 << (fmt->exp_size - 1))) { - if (mask_underflow) { - float_exception_flags |= float_flag_underflow; - *dxc = 0x10; - r->exp += fmt->exp_re_bias; - } - } else if (r_flags & float_flag_inexact) { - float_exception_flags |= float_flag_inexact; - if (mask_inexact) { - bool saved_r_sign, saved_r_precise_sign; - - /* - * Check whether remainder was truncated (rounded - * toward zero) or incremented. - */ - saved_r_sign = r->sign; - saved_r_precise_sign = r_precise.sign; - r->sign = false; - r_precise.sign = false; - if (parts64_compare(r, &r_precise, status, true) < - float_relation_equal) { - *dxc = 0x8; - } else { - *dxc = 0xc; - } - r->sign = saved_r_sign; - r_precise.sign = saved_r_precise_sign; - } - } - } - *cc = 0; - } else if (n->exp > (1 << (fmt->exp_size - 1)) - 1) { - n->exp -= fmt->exp_re_bias; - *cc = r->cls == float_class_zero ? 1 : 3; - } else { - *cc = r->cls == float_class_zero ? 0 : 2; - } - - /* Adjust signs of zero results */ - if (r->cls == float_class_zero) { - r->sign = a->sign; - } - if (n->cls == float_class_zero) { - n->sign = a->sign ^ b->sign; - } - - status->float_exception_flags = float_exception_flags; - } -} - -#define DEFINE_S390_DIVIDE_TO_INTEGER(floatN) \ -void floatN ## _s390_divide_to_integer(floatN a, floatN b, \ - int final_quotient_rounding_mode, \ - bool mask_underflow, bool mask_inexact, \ - floatN *r, floatN *n, \ - uint32_t *cc, int *dxc, \ - float_status *status) \ -{ \ - FloatParts64 pa = floatN ## _unpack_canonical(a, status); \ - FloatParts64 pb = floatN ## _unpack_canonical(b, status); \ - FloatParts64 pr, pn; \ - parts_s390_divide_to_integer(&pa, &pb, final_quotient_rounding_mode, \ - mask_underflow, mask_inexact, \ - &floatN ## _params, \ - &pr, &pn, cc, dxc, status); \ - *r = floatN ## _round_pack_canonical(&pr, status); \ - *n = floatN ## _round_pack_canonical(&pn, status); \ -} - -DEFINE_S390_DIVIDE_TO_INTEGER(float32) -DEFINE_S390_DIVIDE_TO_INTEGER(float64) - static void __attribute__((constructor)) softfloat_init(void) { union_float64 ua, ub, uc, ur; diff --git a/target/s390x/tcg/fpu_helper.c b/target/s390x/tcg/fpu_helper.c index XXXXXXX..XXXXXXX 100644 --- a/target/s390x/tcg/fpu_helper.c +++ b/target/s390x/tcg/fpu_helper.c @@ -XXX,XX +XXX,XX @@ #include "tcg_s390x.h" #include "exec/helper-proto.h" #include "fpu/softfloat.h" +#include "fpu/softfloat-parts.h" /* #define DEBUG_HELPER */ #ifdef DEBUG_HELPER @@ -XXX,XX +XXX,XX @@ Int128 HELPER(dxb)(CPUS390XState *env, Int128 a, Int128 b) return RET128(ret); } +static void parts_s390_divide_to_integer(FloatParts64 *a, FloatParts64 *b, + int final_quotient_rounding_mode, + bool mask_underflow, bool mask_inexact, + const FloatFmt *fmt, + FloatParts64 *r, FloatParts64 *n, + uint32_t *cc, int *dxc, + float_status *status) +{ + /* POp table "Results: DIVIDE TO INTEGER (Part 1 of 2)" */ + if ((float_cmask(a->cls) | float_cmask(b->cls)) & float_cmask_anynan) { + *r = parts64_pick_nan(a, b, status); + *n = *r; + *cc = 1; + } else if (a->cls == float_class_inf || b->cls == float_class_zero) { + *r = parts64_default_nan(status); + *n = *r; + *cc = 1; + status->float_exception_flags |= float_flag_invalid; + } else if (b->cls == float_class_inf) { + *r = *a; + n->cls = float_class_zero; + n->sign = a->sign ^ b->sign; + *cc = 0; + } else { + FloatParts64 *q, q_buf, r_precise; + int float_exception_flags = 0; + bool is_q_smallish; + uint32_t r_flags; + + /* Compute precise quotient */ + q_buf = parts64_div(a, b, status); + q = &q_buf; + + /* + * Check whether two closest integers can be precisely represented, + * i.e., all their bits fit into the fractional part. + */ + is_q_smallish = q->exp < (fmt->frac_size + 1); + + /* + * Final quotient is rounded using final-quotient-rounding method, and + * partial quotient is rounded toward zero. + * + * Rounding of partial quotient may be inexact. This is the whole point + * of distinguishing partial quotients, so ignore the exception. + */ + *n = parts64_round_to_int(q, + is_q_smallish + ? final_quotient_rounding_mode + : float_round_to_zero, + 0, status, fmt); + + /* Compute precise remainder */ + r_precise = parts64_muladd(b, n, a, + float_muladd_negate_product, status); + + /* Round remainder to the target format */ + *r = r_precise; + status->float_exception_flags = 0; + *r = parts64_round_to_fmt(r, status, fmt); + r_flags = status->float_exception_flags; + + /* POp table "Results: DIVIDE TO INTEGER (Part 2 of 2)" */ + if (is_q_smallish) { + if (r->cls != float_class_zero) { + if (r->exp < 2 - (1 << (fmt->exp_size - 1))) { + if (mask_underflow) { + float_exception_flags |= float_flag_underflow; + *dxc = 0x10; + r->exp += fmt->exp_re_bias; + } + } else if (r_flags & float_flag_inexact) { + float_exception_flags |= float_flag_inexact; + if (mask_inexact) { + bool saved_r_sign, saved_r_precise_sign; + + /* + * Check whether remainder was truncated (rounded + * toward zero) or incremented. + */ + saved_r_sign = r->sign; + saved_r_precise_sign = r_precise.sign; + r->sign = false; + r_precise.sign = false; + if (parts64_compare(r, &r_precise, status, true) < + float_relation_equal) { + *dxc = 0x8; + } else { + *dxc = 0xc; + } + r->sign = saved_r_sign; + r_precise.sign = saved_r_precise_sign; + } + } + } + *cc = 0; + } else if (n->exp > (1 << (fmt->exp_size - 1)) - 1) { + n->exp -= fmt->exp_re_bias; + *cc = r->cls == float_class_zero ? 1 : 3; + } else { + *cc = r->cls == float_class_zero ? 0 : 2; + } + + /* Adjust signs of zero results */ + if (r->cls == float_class_zero) { + r->sign = a->sign; + } + if (n->cls == float_class_zero) { + n->sign = a->sign ^ b->sign; + } + + status->float_exception_flags = float_exception_flags; + } +} + +#define DEFINE_S390_DIVIDE_TO_INTEGER(floatN) \ +static void floatN ## _s390_divide_to_integer(floatN a, floatN b, \ + int final_quotient_rounding_mode, bool mask_underflow, bool mask_inexact, \ + floatN *r, floatN *n, uint32_t *cc, int *dxc, float_status *status) \ +{ \ + FloatParts64 pa = floatN ## _unpack_canonical(a, status); \ + FloatParts64 pb = floatN ## _unpack_canonical(b, status); \ + FloatParts64 pr, pn; \ + parts_s390_divide_to_integer(&pa, &pb, final_quotient_rounding_mode, \ + mask_underflow, mask_inexact, \ + &floatN ## _params, \ + &pr, &pn, cc, dxc, status); \ + *r = floatN ## _round_pack_canonical(&pr, status); \ + *n = floatN ## _round_pack_canonical(&pn, status); \ +} + +DEFINE_S390_DIVIDE_TO_INTEGER(float32) +DEFINE_S390_DIVIDE_TO_INTEGER(float64) + void HELPER(dib)(CPUS390XState *env, uint32_t r1, uint32_t r2, uint32_t r3, uint32_t m4, uint32_t bits) { -- 2.43.0 Use softfloat-parts.h so that we can more naturally perform the required operations witha single rounding step. This happens to also simplify the NaN detection step. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- target/arm/tcg/vec_helper.c | 77 +++++++++++++++++++------------------ 1 file changed, 40 insertions(+), 37 deletions(-) diff --git a/target/arm/tcg/vec_helper.c b/target/arm/tcg/vec_helper.c index XXXXXXX..XXXXXXX 100644 --- a/target/arm/tcg/vec_helper.c +++ b/target/arm/tcg/vec_helper.c @@ -XXX,XX +XXX,XX @@ #include "helper.h" #include "tcg/tcg-gvec-desc.h" #include "fpu/softfloat.h" +#include "fpu/softfloat-parts.h" #include "qemu/int128.h" #include "crypto/clmul.h" #include "vec_internal.h" @@ -XXX,XX +XXX,XX @@ float32 bfdotadd(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst) float32 bfdotadd_ebf(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst, float_status *fpst_odd) { + /* Unpack two BFloat16 into two Float32, trivially. */ float32 s1r = e1 << 16; float32 s1c = e1 & 0xffff0000u; float32 s2r = e2 << 16; float32 s2c = e2 & 0xffff0000u; float32 t32; + /* + * Compare f16_dotadd() in sme_helper.c, but here we have + * bfloat16 inputs. In particular that means that we do not + * want the FPCR.FZ16 flush semantics, so we use the normal + * float_status for the input handling here. + */ + FloatParts64 p1r = float32_unpack_canonical(s1r, fpst); + FloatParts64 p1c = float32_unpack_canonical(s1c, fpst); + FloatParts64 p2r = float32_unpack_canonical(s2r, fpst); + FloatParts64 p2c = float32_unpack_canonical(s2c, fpst); + + int all_mask = (float_cmask(p1r.cls) | float_cmask(p1c.cls) | + float_cmask(p1r.cls) | float_cmask(p1c.cls)); + /* C.f. FPProcessNaNs4 */ - if (float32_is_any_nan(s1r) || float32_is_any_nan(s1c) || - float32_is_any_nan(s2r) || float32_is_any_nan(s2c)) { - if (float32_is_signaling_nan(s1r, fpst)) { - t32 = s1r; - } else if (float32_is_signaling_nan(s1c, fpst)) { - t32 = s1c; - } else if (float32_is_signaling_nan(s2r, fpst)) { - t32 = s2r; - } else if (float32_is_signaling_nan(s2c, fpst)) { - t32 = s2c; - } else if (float32_is_any_nan(s1r)) { - t32 = s1r; - } else if (float32_is_any_nan(s1c)) { - t32 = s1c; - } else if (float32_is_any_nan(s2r)) { - t32 = s2r; + if (unlikely(all_mask & float_cmask_anynan)) { + if (unlikely(all_mask & float_cmask_snan)) { + if (p1r.cls == float_class_snan) { + t32 = s1r; + } else if (p1c.cls == float_class_snan) { + t32 = s1c; + } else if (p2r.cls == float_class_snan) { + t32 = s2r; + } else { + t32 = s2c; + } } else { - t32 = s2c; + if (p1r.cls == float_class_qnan) { + t32 = s1r; + } else if (p1c.cls == float_class_qnan) { + t32 = s1c; + } else if (p2r.cls == float_class_qnan) { + t32 = s2r; + } else { + t32 = s2c; + } } /* * FPConvertNaN(FPProcessNaN(t32)) will be done as part * of the final addition below. */ } else { - /* - * Compare f16_dotadd() in sme_helper.c, but here we have - * bfloat16 inputs. In particular that means that we do not - * want the FPCR.FZ16 flush semantics, so we use the normal - * float_status for the input handling here. - */ - float64 e1r = float32_to_float64(s1r, fpst); - float64 e1c = float32_to_float64(s1c, fpst); - float64 e2r = float32_to_float64(s2r, fpst); - float64 e2c = float32_to_float64(s2c, fpst); - float64 t64; - /* * The ARM pseudocode function FPDot performs both multiplies - * and the add with a single rounding operation. Emulate this - * by performing the first multiply in round-to-odd, then doing - * the second multiply as fused multiply-add, and rounding to - * float32 all in one step. + * and the add with a single rounding operation. */ - t64 = float64_mul(e1r, e2r, fpst_odd); - t64 = float64r32_muladd(e1c, e2c, t64, 0, fpst); + FloatParts64 tmp = parts64_mul(&p1r, &p2r, fpst); + tmp = parts64_muladd(&p1c, &p2c, &tmp, 0, fpst); - /* This conversion is exact, because we've already rounded. */ - t32 = float64_to_float32(t64, fpst); + t32 = float32_round_pack_canonical(&tmp, fpst); } /* The final accumulation step is not fused. */ -- 2.43.0 This argument is no longer used. Reviewed-by: Philippe Mathieu-Daudé <philmd@linaro.org> Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- target/arm/tcg/vec_internal.h | 12 +++------ target/arm/tcg/sme_helper.c | 6 ++--- target/arm/tcg/vec_helper.c | 50 +++++++++++++++-------------------- 3 files changed, 29 insertions(+), 39 deletions(-) diff --git a/target/arm/tcg/vec_internal.h b/target/arm/tcg/vec_internal.h index XXXXXXX..XXXXXXX 100644 --- a/target/arm/tcg/vec_internal.h +++ b/target/arm/tcg/vec_internal.h @@ -XXX,XX +XXX,XX @@ float32 bfdotadd(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst); * @sum: addend * @e1, @e2: multiplicand vectors * @fpst: floating-point status to use - * @fpst_odd: floating-point status to use for round-to-odd operations * * BFloat16 2-way dot product of @e1 & @e2, accumulating with @sum. * The @e1 and @e2 operands correspond to the 32-bit source vector @@ -XXX,XX +XXX,XX @@ float32 bfdotadd(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst); * Corresponds to the ARM pseudocode function BFDotAdd, specialized * for the FPCR.EBF == 1 case. */ -float32 bfdotadd_ebf(float32 sum, uint32_t e1, uint32_t e2, - float_status *fpst, float_status *fpst_odd); +float32 bfdotadd_ebf(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst); /** * is_ebf: * @env: CPU state * @statusp: pointer to floating point status to fill in - * @oddstatusp: pointer to floating point status to fill in for round-to-odd * * Determine whether a BFDotAdd operation should use FPCR.EBF = 0 - * or FPCR.EBF = 1 semantics. On return, has initialized *statusp - * and *oddstatusp to suitable float_status arguments to use with either - * bfdotadd() or bfdotadd_ebf(). + * or FPCR.EBF = 1 semantics. On return, has initialized *statusp as suitable + * for float_status arguments to either bfdotadd() or bfdotadd_ebf(). * Returns true for EBF = 1, false for EBF = 0. (The caller should use this * to decide whether to call bfdotadd() or bfdotadd_ebf().) */ -bool is_ebf(CPUARMState *env, float_status *statusp, float_status *oddstatusp); +bool is_ebf(CPUARMState *env, float_status *statusp); /* * Negate as for FPCR.AH=1 -- do not negate NaNs. diff --git a/target/arm/tcg/sme_helper.c b/target/arm/tcg/sme_helper.c index XXXXXXX..XXXXXXX 100644 --- a/target/arm/tcg/sme_helper.c +++ b/target/arm/tcg/sme_helper.c @@ -XXX,XX +XXX,XX @@ static void do_bfmopa_w(void *vza, void *vzn, void *vzm, uint32_t desc, uint32_t negx, bool ah_neg) { intptr_t row, col, oprsz = simd_maxsz(desc); - float_status fpst, fpst_odd; + float_status fpst; - if (is_ebf(env, &fpst, &fpst_odd)) { + if (is_ebf(env, &fpst)) { for (row = 0; row < oprsz; ) { uint16_t prow = pn[H2(row >> 4)]; do { @@ -XXX,XX +XXX,XX @@ static void do_bfmopa_w(void *vza, void *vzn, void *vzm, uint32_t m = *(uint32_t *)(vzm + H1_4(col)); m = f16mop_adj_pair(m, pcol, 0); - *a = bfdotadd_ebf(*a, n, m, &fpst, &fpst_odd); + *a = bfdotadd_ebf(*a, n, m, &fpst); } col += 4; pcol >>= 4; diff --git a/target/arm/tcg/vec_helper.c b/target/arm/tcg/vec_helper.c index XXXXXXX..XXXXXXX 100644 --- a/target/arm/tcg/vec_helper.c +++ b/target/arm/tcg/vec_helper.c @@ -XXX,XX +XXX,XX @@ DO_MMLA_B(gvec_usmmla_b, do_usmmla_b) * BFloat16 Dot Product */ -bool is_ebf(CPUARMState *env, float_status *statusp, float_status *oddstatusp) +bool is_ebf(CPUARMState *env, float_status *statusp) { /* * For BFDOT, BFMMLA, etc, the behaviour depends on FPCR.EBF. @@ -XXX,XX +XXX,XX @@ bool is_ebf(CPUARMState *env, float_status *statusp, float_status *oddstatusp) *statusp = env->vfp.fp_status[is_a64(env) ? FPST_A64 : FPST_A32]; set_default_nan_mode(true, statusp); - if (ebf) { - /* EBF=1 needs to do a step with round-to-odd semantics */ - *oddstatusp = *statusp; - set_float_rounding_mode(float_round_to_odd, oddstatusp); - } else { + if (!ebf) { set_flush_to_zero(true, statusp); set_flush_inputs_to_zero(true, statusp); set_float_rounding_mode(float_round_to_odd_inf, statusp); @@ -XXX,XX +XXX,XX @@ float32 bfdotadd(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst) return t1; } -float32 bfdotadd_ebf(float32 sum, uint32_t e1, uint32_t e2, - float_status *fpst, float_status *fpst_odd) +float32 bfdotadd_ebf(float32 sum, uint32_t e1, uint32_t e2, float_status *fpst) { /* Unpack two BFloat16 into two Float32, trivially. */ float32 s1r = e1 << 16; @@ -XXX,XX +XXX,XX @@ void HELPER(gvec_bfdot)(void *vd, void *vn, void *vm, void *va, intptr_t i, opr_sz = simd_oprsz(desc); float32 *d = vd, *a = va; uint32_t *n = vn, *m = vm; - float_status fpst, fpst_odd; + float_status fpst; - if (is_ebf(env, &fpst, &fpst_odd)) { + if (is_ebf(env, &fpst)) { for (i = 0; i < opr_sz / 4; ++i) { - d[i] = bfdotadd_ebf(a[i], n[i], m[i], &fpst, &fpst_odd); + d[i] = bfdotadd_ebf(a[i], n[i], m[i], &fpst); } } else { for (i = 0; i < opr_sz / 4; ++i) { @@ -XXX,XX +XXX,XX @@ void HELPER(gvec_bfdot_idx)(void *vd, void *vn, void *vm, intptr_t eltspersegment = MIN(16 / 4, elements); float32 *d = vd, *a = va; uint32_t *n = vn, *m = vm; - float_status fpst, fpst_odd; + float_status fpst; - if (is_ebf(env, &fpst, &fpst_odd)) { + if (is_ebf(env, &fpst)) { for (i = 0; i < elements; i += eltspersegment) { uint32_t m_idx = m[i + H4(index)]; for (j = i; j < i + eltspersegment; j++) { - d[j] = bfdotadd_ebf(a[j], n[j], m_idx, &fpst, &fpst_odd); + d[j] = bfdotadd_ebf(a[j], n[j], m_idx, &fpst); } } } else { @@ -XXX,XX +XXX,XX @@ void HELPER(sme2_bfvdot_idx)(void *vd, void *vn, void *vm, uint16_t *n0 = vn; uint16_t *n1 = vn + sizeof(ARMVectorReg); uint32_t *m = vm; - float_status fpst, fpst_odd; + float_status fpst; - if (is_ebf(env, &fpst, &fpst_odd)) { + if (is_ebf(env, &fpst)) { for (i = 0; i < elements; i += eltspersegment) { uint32_t m_idx = m[i + H4(idx)]; for (j = 0; j < eltspersegment; j++) { uint32_t nn = (n0[H2(2 * (i + j) + sel)]) | (n1[H2(2 * (i + j) + sel)] << 16); - d[i + H4(j)] = bfdotadd_ebf(a[i + H4(j)], nn, m_idx, - &fpst, &fpst_odd); + d[i + H4(j)] = bfdotadd_ebf(a[i + H4(j)], nn, m_idx, &fpst); } } } else { @@ -XXX,XX +XXX,XX @@ void HELPER(gvec_bfmmla)(void *vd, void *vn, void *vm, void *va, intptr_t s, opr_sz = simd_oprsz(desc); float32 *d = vd, *a = va; uint32_t *n = vn, *m = vm; - float_status fpst, fpst_odd; + float_status fpst; - if (is_ebf(env, &fpst, &fpst_odd)) { + if (is_ebf(env, &fpst)) { for (s = 0; s < opr_sz / 4; s += 4) { float32 sum00, sum01, sum10, sum11; @@ -XXX,XX +XXX,XX @@ void HELPER(gvec_bfmmla)(void *vd, void *vn, void *vm, void *va, * i j i k j k */ sum00 = a[s + H4(0 + 0)]; - sum00 = bfdotadd_ebf(sum00, n[s + H4(0 + 0)], m[s + H4(0 + 0)], &fpst, &fpst_odd); - sum00 = bfdotadd_ebf(sum00, n[s + H4(0 + 1)], m[s + H4(0 + 1)], &fpst, &fpst_odd); + sum00 = bfdotadd_ebf(sum00, n[s + H4(0 + 0)], m[s + H4(0 + 0)], &fpst); + sum00 = bfdotadd_ebf(sum00, n[s + H4(0 + 1)], m[s + H4(0 + 1)], &fpst); sum01 = a[s + H4(0 + 1)]; - sum01 = bfdotadd_ebf(sum01, n[s + H4(0 + 0)], m[s + H4(2 + 0)], &fpst, &fpst_odd); - sum01 = bfdotadd_ebf(sum01, n[s + H4(0 + 1)], m[s + H4(2 + 1)], &fpst, &fpst_odd); + sum01 = bfdotadd_ebf(sum01, n[s + H4(0 + 0)], m[s + H4(2 + 0)], &fpst); + sum01 = bfdotadd_ebf(sum01, n[s + H4(0 + 1)], m[s + H4(2 + 1)], &fpst); sum10 = a[s + H4(2 + 0)]; - sum10 = bfdotadd_ebf(sum10, n[s + H4(2 + 0)], m[s + H4(0 + 0)], &fpst, &fpst_odd); - sum10 = bfdotadd_ebf(sum10, n[s + H4(2 + 1)], m[s + H4(0 + 1)], &fpst, &fpst_odd); + sum10 = bfdotadd_ebf(sum10, n[s + H4(2 + 0)], m[s + H4(0 + 0)], &fpst); + sum10 = bfdotadd_ebf(sum10, n[s + H4(2 + 1)], m[s + H4(0 + 1)], &fpst); sum11 = a[s + H4(2 + 1)]; - sum11 = bfdotadd_ebf(sum11, n[s + H4(2 + 0)], m[s + H4(2 + 0)], &fpst, &fpst_odd); - sum11 = bfdotadd_ebf(sum11, n[s + H4(2 + 1)], m[s + H4(2 + 1)], &fpst, &fpst_odd); + sum11 = bfdotadd_ebf(sum11, n[s + H4(2 + 0)], m[s + H4(2 + 0)], &fpst); + sum11 = bfdotadd_ebf(sum11, n[s + H4(2 + 1)], m[s + H4(2 + 1)], &fpst); d[s + H4(0 + 0)] = sum00; d[s + H4(0 + 1)] = sum01; -- 2.43.0 Use softfloat-parts.h so that we can more naturally perform the required operations witha single rounding step. This happens to also simplify the NaN detection step. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- target/arm/tcg/sme_helper.c | 96 ++++++++++++++++--------------------- 1 file changed, 40 insertions(+), 56 deletions(-) diff --git a/target/arm/tcg/sme_helper.c b/target/arm/tcg/sme_helper.c index XXXXXXX..XXXXXXX 100644 --- a/target/arm/tcg/sme_helper.c +++ b/target/arm/tcg/sme_helper.c @@ -XXX,XX +XXX,XX @@ #include "accel/tcg/helper-retaddr.h" #include "qemu/int128.h" #include "fpu/softfloat.h" +#include "fpu/softfloat-parts.h" #include "vec_internal.h" #include "sve_ldst_internal.h" @@ -XXX,XX +XXX,XX @@ static inline uint32_t bf16mop_ah_neg_adj_pair(uint32_t pair, uint32_t pg) } static float32 f16_dotadd(float32 sum, uint32_t e1, uint32_t e2, - float_status *s_f16, float_status *s_std, - float_status *s_odd) + float_status *s_f16, float_status *s_std) { /* - * We need three different float_status for different parts of this + * We need two different float_status for different parts of this * operation: * - the input conversion of the float16 values must use the * f16-specific float_status, so that the FPCR.FZ16 control is applied * - operations on float32 including the final accumulation must use * the normal float_status, so that FPCR.FZ is applied - * - we have pre-set-up copy of s_std which is set to round-to-odd, - * for the multiply (see below) */ float16 h1r = e1 & 0xffff; float16 h1c = e1 >> 16; @@ -XXX,XX +XXX,XX @@ static float32 f16_dotadd(float32 sum, uint32_t e1, uint32_t e2, float16 h2c = e2 >> 16; float32 t32; + FloatParts64 p1r = float16_unpack_canonical(h1r, s_f16); + FloatParts64 p1c = float16_unpack_canonical(h1c, s_f16); + FloatParts64 p2r = float16_unpack_canonical(h2r, s_f16); + FloatParts64 p2c = float16_unpack_canonical(h2c, s_f16); + + int all_mask = (float_cmask(p1r.cls) | float_cmask(p1c.cls) | + float_cmask(p1r.cls) | float_cmask(p1c.cls)); + /* C.f. FPProcessNaNs4 */ - if (float16_is_any_nan(h1r) || float16_is_any_nan(h1c) || - float16_is_any_nan(h2r) || float16_is_any_nan(h2c)) { + if (unlikely(all_mask & float_cmask_anynan)) { float16 t16; - if (float16_is_signaling_nan(h1r, s_f16)) { - t16 = h1r; - } else if (float16_is_signaling_nan(h1c, s_f16)) { - t16 = h1c; - } else if (float16_is_signaling_nan(h2r, s_f16)) { - t16 = h2r; - } else if (float16_is_signaling_nan(h2c, s_f16)) { - t16 = h2c; - } else if (float16_is_any_nan(h1r)) { - t16 = h1r; - } else if (float16_is_any_nan(h1c)) { - t16 = h1c; - } else if (float16_is_any_nan(h2r)) { - t16 = h2r; + if (unlikely(all_mask & float_cmask_snan)) { + if (p1r.cls == float_class_snan) { + t16 = h1r; + } else if (p1c.cls == float_class_snan) { + t16 = h1c; + } else if (p2r.cls == float_class_snan) { + t16 = h2r; + } else { + t16 = h2c; + } } else { - t16 = h2c; + if (p1r.cls == float_class_qnan) { + t16 = h1r; + } else if (p1c.cls == float_class_qnan) { + t16 = h1c; + } else if (p2r.cls == float_class_qnan) { + t16 = h2r; + } else { + t16 = h2c; + } } t32 = float16_to_float32(t16, true, s_f16); } else { - float64 e1r = float16_to_float64(h1r, true, s_f16); - float64 e1c = float16_to_float64(h1c, true, s_f16); - float64 e2r = float16_to_float64(h2r, true, s_f16); - float64 e2c = float16_to_float64(h2c, true, s_f16); - float64 t64; - /* * The ARM pseudocode function FPDot performs both multiplies - * and the add with a single rounding operation. Emulate this - * by performing the first multiply in round-to-odd, then doing - * the second multiply as fused multiply-add, and rounding to - * float32 all in one step. + * and the add with a single rounding operation. */ - t64 = float64_mul(e1r, e2r, s_odd); - t64 = float64r32_muladd(e1c, e2c, t64, 0, s_std); + FloatParts64 tmp = parts64_mul(&p1r, &p2r, s_std); + tmp = parts64_muladd(&p1c, &p2c, &tmp, 0, s_std); - /* This conversion is exact, because we've already rounded. */ - t32 = float64_to_float32(t64, s_std); + t32 = float32_round_pack_canonical(&tmp, s_std); } /* The final accumulation step is not fused. */ @@ -XXX,XX +XXX,XX @@ static void do_fmopa_w_h(void *vza, void *vzn, void *vzm, uint16_t *pn, uint32_t negx, bool ah_neg) { intptr_t row, col, oprsz = simd_maxsz(desc); - float_status fpst_odd = env->vfp.fp_status[FPST_ZA]; - - set_float_rounding_mode(float_round_to_odd, &fpst_odd); for (row = 0; row < oprsz; ) { uint16_t prow = pn[H2(row >> 4)]; @@ -XXX,XX +XXX,XX @@ static void do_fmopa_w_h(void *vza, void *vzn, void *vzm, uint16_t *pn, m = f16mop_adj_pair(m, pcol, 0); *a = f16_dotadd(*a, n, m, &env->vfp.fp_status[FPST_ZA_F16], - &env->vfp.fp_status[FPST_ZA], - &fpst_odd); + &env->vfp.fp_status[FPST_ZA]); } col += 4; pcol >>= 4; @@ -XXX,XX +XXX,XX @@ void HELPER(sme2_fdot_h)(void *vd, void *vn, void *vm, void *va, bool za = extract32(desc, SIMD_DATA_SHIFT, 1); float_status *fpst_std = &env->vfp.fp_status[za ? FPST_ZA : FPST_A64]; float_status *fpst_f16 = &env->vfp.fp_status[za ? FPST_ZA_F16 : FPST_A64_F16]; - float_status fpst_odd = *fpst_std; float32 *d = vd, *a = va; uint32_t *n = vn, *m = vm; - set_float_rounding_mode(float_round_to_odd, &fpst_odd); - for (i = 0; i < oprsz / sizeof(float32); ++i) { d[H4(i)] = f16_dotadd(a[H4(i)], n[H4(i)], m[H4(i)], - fpst_f16, fpst_std, &fpst_odd); + fpst_f16, fpst_std); } } @@ -XXX,XX +XXX,XX @@ void HELPER(sme2_fdot_idx_h)(void *vd, void *vn, void *vm, void *va, bool za = extract32(desc, SIMD_DATA_SHIFT + 2, 1); float_status *fpst_std = &env->vfp.fp_status[za ? FPST_ZA : FPST_A64]; float_status *fpst_f16 = &env->vfp.fp_status[za ? FPST_ZA_F16 : FPST_A64_F16]; - float_status fpst_odd = *fpst_std; float32 *d = vd, *a = va; uint32_t *n = vn, *m = (uint32_t *)vm + H4(idx); - set_float_rounding_mode(float_round_to_odd, &fpst_odd); - for (i = 0; i < elements; i += eltspersegment) { uint32_t mm = m[i]; for (j = 0; j < eltspersegment; ++j) { d[H4(i + j)] = f16_dotadd(a[H4(i + j)], n[H4(i + j)], mm, - fpst_f16, fpst_std, &fpst_odd); + fpst_f16, fpst_std); } } } @@ -XXX,XX +XXX,XX @@ void HELPER(sme2_fvdot_idx_h)(void *vd, void *vn, void *vm, void *va, intptr_t eltspersegment = MIN(4, elements); int idx = extract32(desc, SIMD_DATA_SHIFT, 2); int sel = extract32(desc, SIMD_DATA_SHIFT + 2, 1); - float_status fpst_odd, *fpst_std, *fpst_f16; float32 *d = vd, *a = va; uint16_t *n0 = vn; uint16_t *n1 = vn + sizeof(ARMVectorReg); uint32_t *m = (uint32_t *)vm + H4(idx); - fpst_std = &env->vfp.fp_status[FPST_ZA]; - fpst_f16 = &env->vfp.fp_status[FPST_ZA_F16]; - fpst_odd = *fpst_std; - set_float_rounding_mode(float_round_to_odd, &fpst_odd); - for (i = 0; i < elements; i += eltspersegment) { uint32_t mm = m[i]; for (j = 0; j < eltspersegment; ++j) { uint32_t nn = (n0[H2(2 * (i + j) + sel)]) | (n1[H2(2 * (i + j) + sel)] << 16); d[i + H4(j)] = f16_dotadd(a[i + H4(j)], nn, mm, - fpst_f16, fpst_std, &fpst_odd); + &env->vfp.fp_status[FPST_ZA_F16], + &env->vfp.fp_status[FPST_ZA]); } } } -- 2.43.0 Notice the conflict between saturation and rebias_overflow, which should never happen. Otherwise, saturate to INT32_MAX and allow the usual overflow to max, infinity, or dnan. Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat-parts.c.inc | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ static void partsN(uncanon_normal)(FloatPartsN *p, float_status *s, g_assert_not_reached(); } - exp = p->exp + fmt->exp_bias; + /* Because exp_bias is positive, we can only overflow past INT_MAX. */ + if (sadd32_overflow(p->exp, fmt->exp_bias, &exp)) { + /* + * rebias_overflow wants to compute a modulo exponent, which + * conflicts with saturation. That said, saturation can only + * happen with scalbn, which is not a PowerPC operation. + */ + assert(!s->rebias_overflow); + exp = INT32_MAX; + } + if (likely(exp > 0)) { if (p->frac_lo & round_mask) { flags |= float_flag_inexact; -- 2.43.0 Avoid exponent overflow as well as checking that we don't lose information with opposing scaling. Use it in partsN(scalbn) and partsN(round_to_int_normal). Signed-off-by: Richard Henderson <richard.henderson@linaro.org> --- fpu/softfloat.c | 17 +++++++++++++++++ fpu/softfloat-parts.c.inc | 5 ++--- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/fpu/softfloat.c b/fpu/softfloat.c index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -XXX,XX +XXX,XX @@ static float128 QEMU_FLATTEN float128_pack_raw(const FloatParts128 *p) *----------------------------------------------------------------------------*/ #include "softfloat-specialize.c.inc" +static int32_t exp_scalbn(int32_t exp, int32_t scale) +{ + /* + * Catch chains of scaling which lose information. + * In particular, if the exponent has been saturated, + * do not allow it to become unsaturated. + */ + if (unlikely(exp == INT32_MAX)) { + assert(scale >= 0); + } else if (unlikely(exp == INT32_MIN)) { + assert(scale <= 0); + } else { + exp = sadd32_saturate(exp, scale); + } + return exp; +} + /* * Helper functions for softfloat-parts.c.inc, per-size operations. */ diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index XXXXXXX..XXXXXXX 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -XXX,XX +XXX,XX @@ static bool partsN(round_to_int_normal)(FloatPartsN *a, FloatRoundMode rmode, uint64_t frac_lsb, frac_lsbm1, rnd_even_mask, rnd_mask, inc; int shift_adj; - scale = MIN(MAX(scale, -0x10000), 0x10000); - a->exp += scale; + a->exp = exp_scalbn(a->exp, scale); if (a->exp < 0) { bool one; @@ -XXX,XX +XXX,XX @@ FloatPartsN partsN(scalbn)(const FloatPartsN *a, int n, float_status *s) case float_class_normal: { FloatPartsN r = *a; - r.exp += MIN(MAX(n, -0x10000), 0x10000); + r.exp = exp_scalbn(r.exp, n); return r; } default: -- 2.43.0