From e99c43735a413e2566b4ad36eeda6dd061a9b939 Mon Sep 17 00:00:00 2001 From: Richard Henderson Date: Sat, 7 Nov 2020 11:19:32 -0800 Subject: [PATCH] softfloat: Move the binary point to the msb MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rather than point the binary point at msb-1, put it at the msb. Use uadd64_overflow to detect when addition overflows instead of DECOMPOSED_OVERFLOW_BIT. This reduces the number of special cases within the code, such as shifting an int64_t either left or right during conversion. Reviewed-by: Alex Bennée Signed-off-by: Richard Henderson --- fpu/softfloat.c | 169 +++++++++++++++++++----------------------------- 1 file changed, 66 insertions(+), 103 deletions(-) diff --git a/fpu/softfloat.c b/fpu/softfloat.c index 67cfa0fd82..cd777743f1 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -503,9 +503,8 @@ typedef struct { bool sign; } FloatParts; -#define DECOMPOSED_BINARY_POINT (64 - 2) +#define DECOMPOSED_BINARY_POINT 63 #define DECOMPOSED_IMPLICIT_BIT (1ull << DECOMPOSED_BINARY_POINT) -#define DECOMPOSED_OVERFLOW_BIT (DECOMPOSED_IMPLICIT_BIT << 1) /* Structure holding all of the relevant parameters for a format. * exp_size: the size of the exponent field @@ -657,7 +656,7 @@ static FloatParts sf_canonicalize(FloatParts part, const FloatFmt *parm, part.cls = float_class_zero; part.frac = 0; } else { - int shift = clz64(part.frac) - 1; + int shift = clz64(part.frac); part.cls = float_class_normal; part.exp = parm->frac_shift - parm->exp_bias - shift + 1; part.frac <<= shift; @@ -727,9 +726,8 @@ static FloatParts round_canonical(FloatParts p, float_status *s, if (likely(exp > 0)) { if (frac & round_mask) { flags |= float_flag_inexact; - frac += inc; - if (frac & DECOMPOSED_OVERFLOW_BIT) { - frac >>= 1; + if (uadd64_overflow(frac, inc, &frac)) { + frac = (frac >> 1) | DECOMPOSED_IMPLICIT_BIT; exp++; } } @@ -758,9 +756,12 @@ static FloatParts round_canonical(FloatParts p, float_status *s, p.cls = float_class_zero; goto do_zero; } else { - bool is_tiny = s->tininess_before_rounding - || (exp < 0) - || !((frac + inc) & DECOMPOSED_OVERFLOW_BIT); + bool is_tiny = s->tininess_before_rounding || (exp < 0); + + if (!is_tiny) { + uint64_t discard; + is_tiny = !uadd64_overflow(frac, inc, &discard); + } shift64RightJamming(frac, 1 - exp, &frac); if (frac & round_mask) { @@ -985,7 +986,7 @@ static FloatParts addsub_floats(FloatParts a, FloatParts b, bool subtract, a.cls = float_class_zero; a.sign = s->float_rounding_mode == float_round_down; } else { - int shift = clz64(a.frac) - 1; + int shift = clz64(a.frac); a.frac = a.frac << shift; a.exp = a.exp - shift; a.sign = a_sign; @@ -1022,9 +1023,10 @@ static FloatParts addsub_floats(FloatParts a, FloatParts b, bool subtract, shift64RightJamming(a.frac, b.exp - a.exp, &a.frac); a.exp = b.exp; } - a.frac += b.frac; - if (a.frac & DECOMPOSED_OVERFLOW_BIT) { + + if (uadd64_overflow(a.frac, b.frac, &a.frac)) { shift64RightJamming(a.frac, 1, &a.frac); + a.frac |= DECOMPOSED_IMPLICIT_BIT; a.exp += 1; } return a; @@ -1219,16 +1221,17 @@ static FloatParts mul_floats(FloatParts a, FloatParts b, float_status *s) int exp = a.exp + b.exp; mul64To128(a.frac, b.frac, &hi, &lo); - shift128RightJamming(hi, lo, DECOMPOSED_BINARY_POINT, &hi, &lo); - if (lo & DECOMPOSED_OVERFLOW_BIT) { - shift64RightJamming(lo, 1, &lo); + if (hi & DECOMPOSED_IMPLICIT_BIT) { exp += 1; + } else { + hi <<= 1; } + hi |= (lo != 0); /* Re-use a */ a.exp = exp; a.sign = sign; - a.frac = lo; + a.frac = hi; return a; } /* handle all the NaN cases */ @@ -1411,56 +1414,41 @@ static FloatParts muladd_floats(FloatParts a, FloatParts b, FloatParts c, p_exp = a.exp + b.exp; - /* Multiply of 2 62-bit numbers produces a (2*62) == 124-bit - * result. - */ mul64To128(a.frac, b.frac, &hi, &lo); - /* binary point now at bit 124 */ - /* check for overflow */ - if (hi & (1ULL << (DECOMPOSED_BINARY_POINT * 2 + 1 - 64))) { - shift128RightJamming(hi, lo, 1, &hi, &lo); + /* Renormalize to the msb. */ + if (hi & DECOMPOSED_IMPLICIT_BIT) { p_exp += 1; + } else { + shortShift128Left(hi, lo, 1, &hi, &lo); } /* + add/sub */ - if (c.cls == float_class_zero) { - /* move binary point back to 62 */ - shift128RightJamming(hi, lo, DECOMPOSED_BINARY_POINT, &hi, &lo); - } else { + if (c.cls != float_class_zero) { int exp_diff = p_exp - c.exp; if (p_sign == c.sign) { /* Addition */ if (exp_diff <= 0) { - shift128RightJamming(hi, lo, - DECOMPOSED_BINARY_POINT - exp_diff, - &hi, &lo); - lo += c.frac; + shift64RightJamming(hi, -exp_diff, &hi); p_exp = c.exp; + if (uadd64_overflow(hi, c.frac, &hi)) { + shift64RightJamming(hi, 1, &hi); + hi |= DECOMPOSED_IMPLICIT_BIT; + p_exp += 1; + } } else { - uint64_t c_hi, c_lo; - /* shift c to the same binary point as the product (124) */ - c_hi = c.frac >> 2; - c_lo = 0; - shift128RightJamming(c_hi, c_lo, - exp_diff, - &c_hi, &c_lo); - add128(hi, lo, c_hi, c_lo, &hi, &lo); - /* move binary point back to 62 */ - shift128RightJamming(hi, lo, DECOMPOSED_BINARY_POINT, &hi, &lo); + uint64_t c_hi, c_lo, over; + shift128RightJamming(c.frac, 0, exp_diff, &c_hi, &c_lo); + add192(0, hi, lo, 0, c_hi, c_lo, &over, &hi, &lo); + if (over) { + shift64RightJamming(hi, 1, &hi); + hi |= DECOMPOSED_IMPLICIT_BIT; + p_exp += 1; + } } - - if (lo & DECOMPOSED_OVERFLOW_BIT) { - shift64RightJamming(lo, 1, &lo); - p_exp += 1; - } - } else { /* Subtraction */ - uint64_t c_hi, c_lo; - /* make C binary point match product at bit 124 */ - c_hi = c.frac >> 2; - c_lo = 0; + uint64_t c_hi = c.frac, c_lo = 0; if (exp_diff <= 0) { shift128RightJamming(hi, lo, -exp_diff, &hi, &lo); @@ -1495,20 +1483,15 @@ static FloatParts muladd_floats(FloatParts a, FloatParts b, FloatParts c, /* Normalizing to a binary point of 124 is the correct adjust for the exponent. However since we're shifting, we might as well put the binary point back - at 62 where we really want it. Therefore shift as + at 63 where we really want it. Therefore shift as if we're leaving 1 bit at the top of the word, but adjust the exponent as if we're leaving 3 bits. */ - shift -= 1; - if (shift >= 64) { - lo = lo << (shift - 64); - } else { - hi = (hi << shift) | (lo >> (64 - shift)); - lo = hi | ((lo << shift) != 0); - } - p_exp -= shift - 2; + shift128Left(hi, lo, shift, &hi, &lo); + p_exp -= shift; } } } + hi |= (lo != 0); if (flags & float_muladd_halve_result) { p_exp -= 1; @@ -1518,7 +1501,7 @@ static FloatParts muladd_floats(FloatParts a, FloatParts b, FloatParts c, a.cls = float_class_normal; a.sign = p_sign ^ sign_flip; a.exp = p_exp; - a.frac = lo; + a.frac = hi; return a; } @@ -1742,25 +1725,17 @@ static FloatParts div_floats(FloatParts a, FloatParts b, float_status *s) * exponent to match. * * The udiv_qrnnd algorithm that we're using requires normalization, - * i.e. the msb of the denominator must be set. Since we know that - * DECOMPOSED_BINARY_POINT is msb-1, the inputs must be shifted left - * by one (more), and the remainder must be shifted right by one. + * i.e. the msb of the denominator must be set, which is already true. */ if (a.frac < b.frac) { exp -= 1; - shift128Left(0, a.frac, DECOMPOSED_BINARY_POINT + 2, &n1, &n0); - } else { shift128Left(0, a.frac, DECOMPOSED_BINARY_POINT + 1, &n1, &n0); + } else { + shift128Left(0, a.frac, DECOMPOSED_BINARY_POINT, &n1, &n0); } - q = udiv_qrnnd(&r, n1, n0, b.frac << 1); + q = udiv_qrnnd(&r, n1, n0, b.frac); - /* - * Set lsb if there is a remainder, to set inexact. - * As mentioned above, to find the actual value of the remainder we - * would need to shift right, but (1) we are only concerned about - * non-zero-ness, and (2) the remainder will always be even because - * both inputs to the division primitive are even. - */ + /* Set lsb if there is a remainder, to set inexact. */ a.frac = q | (r != 0); a.sign = sign; a.exp = exp; @@ -2135,12 +2110,12 @@ static FloatParts round_to_int(FloatParts a, FloatRoundMode rmode, if (a.frac & rnd_mask) { s->float_exception_flags |= float_flag_inexact; - a.frac += inc; - a.frac &= ~rnd_mask; - if (a.frac & DECOMPOSED_OVERFLOW_BIT) { + if (uadd64_overflow(a.frac, inc, &a.frac)) { a.frac >>= 1; + a.frac |= DECOMPOSED_IMPLICIT_BIT; a.exp++; } + a.frac &= ~rnd_mask; } } break; @@ -2213,10 +2188,8 @@ static int64_t round_to_int_and_pack(FloatParts in, FloatRoundMode rmode, case float_class_zero: return 0; case float_class_normal: - if (p.exp < DECOMPOSED_BINARY_POINT) { + if (p.exp <= DECOMPOSED_BINARY_POINT) { r = p.frac >> (DECOMPOSED_BINARY_POINT - p.exp); - } else if (p.exp - DECOMPOSED_BINARY_POINT < 2) { - r = p.frac << (p.exp - DECOMPOSED_BINARY_POINT); } else { r = UINT64_MAX; } @@ -2498,10 +2471,8 @@ static uint64_t round_to_uint_and_pack(FloatParts in, FloatRoundMode rmode, return 0; } - if (p.exp < DECOMPOSED_BINARY_POINT) { + if (p.exp <= DECOMPOSED_BINARY_POINT) { r = p.frac >> (DECOMPOSED_BINARY_POINT - p.exp); - } else if (p.exp - DECOMPOSED_BINARY_POINT < 2) { - r = p.frac << (p.exp - DECOMPOSED_BINARY_POINT); } else { s->float_exception_flags = orig_flags | float_flag_invalid; return max; @@ -2765,11 +2736,11 @@ static FloatParts int_to_float(int64_t a, int scale, float_status *status) f = -f; r.sign = true; } - shift = clz64(f) - 1; + shift = clz64(f); scale = MIN(MAX(scale, -0x10000), 0x10000); r.exp = DECOMPOSED_BINARY_POINT - shift + scale; - r.frac = (shift < 0 ? DECOMPOSED_IMPLICIT_BIT : f << shift); + r.frac = f << shift; } return r; @@ -2920,21 +2891,16 @@ bfloat16 int16_to_bfloat16(int16_t a, float_status *status) static FloatParts uint_to_float(uint64_t a, int scale, float_status *status) { FloatParts r = { .sign = false }; + int shift; if (a == 0) { r.cls = float_class_zero; } else { scale = MIN(MAX(scale, -0x10000), 0x10000); + shift = clz64(a); r.cls = float_class_normal; - if ((int64_t)a < 0) { - r.exp = DECOMPOSED_BINARY_POINT + 1 + scale; - shift64RightJamming(a, 1, &a); - r.frac = a; - } else { - int shift = clz64(a) - 1; - r.exp = DECOMPOSED_BINARY_POINT - shift + scale; - r.frac = a << shift; - } + r.exp = DECOMPOSED_BINARY_POINT - shift + scale; + r.frac = a << shift; } return r; @@ -3475,12 +3441,9 @@ static FloatParts sqrt_float(FloatParts a, float_status *s, const FloatFmt *p) /* We need two overflow bits at the top. Adding room for that is a * right shift. If the exponent is odd, we can discard the low bit * by multiplying the fraction by 2; that's a left shift. Combine - * those and we shift right if the exponent is even. + * those and we shift right by 1 if the exponent is odd, otherwise 2. */ - a_frac = a.frac; - if (!(a.exp & 1)) { - a_frac >>= 1; - } + a_frac = a.frac >> (2 - (a.exp & 1)); a.exp >>= 1; /* Bit-by-bit computation of sqrt. */ @@ -3488,10 +3451,10 @@ static FloatParts sqrt_float(FloatParts a, float_status *s, const FloatFmt *p) s_frac = 0; /* Iterate from implicit bit down to the 3 extra bits to compute a - * properly rounded result. Remember we've inserted one more bit - * at the top, so these positions are one less. + * properly rounded result. Remember we've inserted two more bits + * at the top, so these positions are two less. */ - bit = DECOMPOSED_BINARY_POINT - 1; + bit = DECOMPOSED_BINARY_POINT - 2; last_bit = MAX(p->frac_shift - 4, 0); do { uint64_t q = 1ULL << bit; @@ -3507,7 +3470,7 @@ static FloatParts sqrt_float(FloatParts a, float_status *s, const FloatFmt *p) /* Undo the right shift done above. If there is any remaining * fraction, the result is inexact. Set the sticky bit. */ - a.frac = (r_frac << 1) + (a_frac != 0); + a.frac = (r_frac << 2) + (a_frac != 0); return a; }