From 208a5e8c02d7a1295004505b0ac0b1e9174c4ab1 Mon Sep 17 00:00:00 2001 From: Michael Tautschnig Date: Tue, 28 Apr 2026 11:46:45 +0000 Subject: [PATCH 1/2] Add round_to_integral to float_bvt, rewrite in float_utilst Add float_bvt::round_to_integral (a new expression-rewrite for floatbv_round_to_integral_exprt) and rewrite float_utilst::round_to_integral in the same shape. Both encodings replace the previous add-magic-subtract-magic algorithm with a direct bitvector approach: for each possible biased exponent E in [bias, bias+f-1], mask off the trailing fractional bits of the packed representation and apply the rounding-mode-specific increment. Values with |x| < 1 are handled separately; special values and values already integer short-circuit to src. User-visible fix: round_to_integral now works on the SMT non-FPA path. Previously, 'cbmc --smt2 --z3' on any rintf/rint program tripped UNEXPECTEDCASE("TODO floatbv_round_to_integral without FPA") in smt2_conv.cpp. The new dispatch routes through the existing bvfp_set / convert_floatbv define-fun cache, so multiple call sites share a single SMT model definition rather than duplicating the float_bvt rewrite at each call site. The new encoding uses an O(f^2) expression-tree (an f-deep ITE cascade where each branch contains an f-bit adder), versus O(f) for the magic-number formulation. For double (f=52) that is roughly a 25x growth in node count. The two-place duplication between float_bv.cpp and float_utils.cpp is intentional and matches the rest of these files' structure; comments mark them as needing to be kept in lockstep. A loop-free barrel-shifter formulation is plausible as a follow-up. The regression test exercises rintf and rint across all five rounding modes, special values (NaN, +/-Inf, +/-0), denormals, negative ties, and an already-integer input above the loop bound (exp_ge_f early-out). The SMT variant fails on develop with EXIT=134 (UNEXPECTEDCASE) and passes here. Co-authored-by: Kiro --- .../cbmc/Float-round-to-integral/main.c | 165 ++++++++++++++++++ .../cbmc/Float-round-to-integral/test.desc | 12 ++ .../Float-round-to-integral/test_smt.desc | 12 ++ src/solvers/floatbv/float_bv.cpp | 134 ++++++++++++++ src/solvers/floatbv/float_bv.h | 18 ++ src/solvers/floatbv/float_utils.cpp | 123 ++++++++++--- src/solvers/smt2/smt2_conv.cpp | 9 +- 7 files changed, 452 insertions(+), 21 deletions(-) create mode 100644 regression/cbmc/Float-round-to-integral/main.c create mode 100644 regression/cbmc/Float-round-to-integral/test.desc create mode 100644 regression/cbmc/Float-round-to-integral/test_smt.desc diff --git a/regression/cbmc/Float-round-to-integral/main.c b/regression/cbmc/Float-round-to-integral/main.c new file mode 100644 index 00000000000..b98b038a452 --- /dev/null +++ b/regression/cbmc/Float-round-to-integral/main.c @@ -0,0 +1,165 @@ +// Test that round-to-integral (rintf, rint) works correctly under various +// rounding modes, special-value inputs, denormals, and at the boundaries +// where the f-deep ITE cascade dispatches. +// +// Without the rewrite, `cbmc --smt2 --z3` (the non-FPA SMT path, exercised +// by test_smt.desc) hits an UNEXPECTEDCASE invariant in +// smt2_conv.cpp::convert_floatbv_round_to_integral, so this regression +// pins down the user-visible fix on that path. + +#include +#include +#include + +extern int __CPROVER_rounding_mode; + +int main() +{ + // Rounding-mode encoding (matches ieee_floatt::rounding_modet): + // 0 = ROUND_TO_EVEN + // 1 = ROUND_TO_MINUS_INF + // 2 = ROUND_TO_PLUS_INF + // 3 = ROUND_TO_ZERO + // 4 = ROUND_TO_AWAY + + // ROUND_TO_EVEN: ties round to even. + __CPROVER_rounding_mode = 0; + float a = 2.5f; + float r1 = rintf(a); + assert(r1 == 2.0f); + + __CPROVER_rounding_mode = 0; + float b = 3.5f; + float r2 = rintf(b); + assert(r2 == 4.0f); + + // Negative ties under ROUND_TO_EVEN. + __CPROVER_rounding_mode = 0; + float bn1 = -2.5f; + float rn1 = rintf(bn1); + assert(rn1 == -2.0f); + + __CPROVER_rounding_mode = 0; + float bn2 = -3.5f; + float rn2 = rintf(bn2); + assert(rn2 == -4.0f); + + // ROUND_TO_PLUS_INF (ceil). + __CPROVER_rounding_mode = 2; + float c = 2.3f; + float r3 = rintf(c); + assert(r3 == 3.0f); + + // ROUND_TO_MINUS_INF (floor). + __CPROVER_rounding_mode = 1; + float d = 2.7f; + float r4 = rintf(d); + assert(r4 == 2.0f); + + // ROUND_TO_ZERO (truncate). + __CPROVER_rounding_mode = 3; + float e = -2.7f; + float r5 = rintf(e); + assert(r5 == -2.0f); + + // ROUND_TO_AWAY: 0.5 rounds away from zero. + __CPROVER_rounding_mode = 4; + float aw1 = 0.5f; + float r_aw1 = rintf(aw1); + assert(r_aw1 == 1.0f); + + __CPROVER_rounding_mode = 4; + float aw2 = -0.5f; + float r_aw2 = rintf(aw2); + assert(r_aw2 == -1.0f); + + __CPROVER_rounding_mode = 4; + float aw3 = 2.5f; + float r_aw3 = rintf(aw3); + assert(r_aw3 == 3.0f); + + // Already-integral input (small). + __CPROVER_rounding_mode = 0; + float f = 42.0f; + float r6 = rintf(f); + assert(r6 == 42.0f); + + // Already-integral input with unbiased exponent >= f for float (f=23): + // 2^24 = 16777216 has biased exp 1023+24 = 1047, well past the loop + // upper bound. Exercises the exp_ge_f early-out. + __CPROVER_rounding_mode = 0; + float big_int = 16777216.0f; + float r_big = rintf(big_int); + assert(r_big == 16777216.0f); + + // Special values: should pass through unchanged via is_special. + // NaN. + __CPROVER_rounding_mode = 0; + float nan_in; + __CPROVER_assume(isnan(nan_in)); + float r_nan = rintf(nan_in); + assert(isnan(r_nan)); + + // +Infinity. + __CPROVER_rounding_mode = 0; + float inf_pos; + __CPROVER_assume(inf_pos == INFINITY); + float r_ip = rintf(inf_pos); + assert(isinf(r_ip) && r_ip > 0); + + // -Infinity. + __CPROVER_rounding_mode = 0; + float inf_neg; + __CPROVER_assume(inf_neg == -INFINITY); + float r_in = rintf(inf_neg); + assert(isinf(r_in) && r_in < 0); + + // +0.0. + __CPROVER_rounding_mode = 0; + float pz = 0.0f; + float r_pz = rintf(pz); + assert(r_pz == 0.0f); + + // -0.0 should preserve sign under any mode. + __CPROVER_rounding_mode = 0; + float nz = -0.0f; + float r_nz = rintf(nz); + assert(r_nz == 0.0f); // both ±0 compare equal to 0 + assert(signbit(r_nz)); + + // Denormals: rintf of a positive denormal under ROUND_TO_EVEN -> +0. + __CPROVER_rounding_mode = 0; + float den = FLT_MIN / 2.0f; + float r_den = rintf(den); + assert(r_den == 0.0f); + + // Denormal under ROUND_TO_PLUS_INF (ceil) -> +1. + __CPROVER_rounding_mode = 2; + float den2 = FLT_MIN / 2.0f; + float r_den2 = rintf(den2); + assert(r_den2 == 1.0f); + + // double / rint coverage. Exercises a different fraction width (f=52). + __CPROVER_rounding_mode = 0; + double da = 2.5; + double dr1 = rint(da); + assert(dr1 == 2.0); + + __CPROVER_rounding_mode = 0; + double db = 3.5; + double dr2 = rint(db); + assert(dr2 == 4.0); + + __CPROVER_rounding_mode = 1; + double dc = -2.7; + double dr3 = rint(dc); + assert(dr3 == -3.0); + + // double already-integral past the loop bound (2^53 = 9007199254740992). + __CPROVER_rounding_mode = 0; + double dbig = 9007199254740992.0; + double dr_big = rint(dbig); + assert(dr_big == 9007199254740992.0); + + return 0; +} diff --git a/regression/cbmc/Float-round-to-integral/test.desc b/regression/cbmc/Float-round-to-integral/test.desc new file mode 100644 index 00000000000..cfde2775906 --- /dev/null +++ b/regression/cbmc/Float-round-to-integral/test.desc @@ -0,0 +1,12 @@ +CORE no-new-smt +main.c +--floatbv +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +Round-to-integral (rintf, rint) under all rounding modes, with coverage +of special values (NaN, ±Inf, ±0), denormals, negative ties, and +already-integral inputs that hit the exp_ge_f early-out. diff --git a/regression/cbmc/Float-round-to-integral/test_smt.desc b/regression/cbmc/Float-round-to-integral/test_smt.desc new file mode 100644 index 00000000000..58784182a0d --- /dev/null +++ b/regression/cbmc/Float-round-to-integral/test_smt.desc @@ -0,0 +1,12 @@ +CORE smt-backend no-new-smt +main.c +--smt2 --z3 +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +Round-to-integral on the non-FPA SMT path. Without the rewrite, +smt2_conv hits UNEXPECTEDCASE("TODO floatbv_round_to_integral without +FPA"), so this test pins down the user-visible fix. diff --git a/src/solvers/floatbv/float_bv.cpp b/src/solvers/floatbv/float_bv.cpp index ee657069fe6..a1c447b6246 100644 --- a/src/solvers/floatbv/float_bv.cpp +++ b/src/solvers/floatbv/float_bv.cpp @@ -163,6 +163,12 @@ exprt float_bvt::convert(const exprt &expr) const const auto &op = to_unary_expr(expr).op(); return isnormal(op, get_spec(op)); } + else if(expr.id() == ID_floatbv_round_to_integral) + { + const auto &rti_expr = to_floatbv_round_to_integral_expr(expr); + return round_to_integral( + rti_expr.op(), rti_expr.rounding_mode(), get_spec(expr)); + } else if(expr.id()==ID_lt) { const auto &rel_expr = to_binary_relation_expr(expr); @@ -496,6 +502,134 @@ exprt float_bvt::conversion( } } +/// Round a floating-point value to integral. +/// +/// Strategy: enumerate the possible biased exponents that have a non-trivial +/// fractional part (`bias..bias+f-1`) and, for each one, mask off the bottom +/// `drop = f - unbiased_exp` bits of the packed representation and +/// conditionally increment by `2^drop`. Values with `|x| < 1` are handled +/// separately (result is ±0 or ±1 depending on rounding mode and sign). +/// Special values (NaN, ±Inf, ±0) and values whose unbiased exponent already +/// meets or exceeds the fraction width short-circuit to `src`. +/// +/// The encoding builds an `f`-deep ITE cascade where each branch contains an +/// `f`-bit adder, giving an O(f^2) expression-tree size. For `double` (f=52) +/// that is roughly a 25× growth in node count over the prior +/// add-magic-subtract-magic formulation; the trade is correctness on the SMT +/// non-FPA path (which previously hit `UNEXPECTEDCASE`) and a more directly +/// auditable algorithm. A loop-free formulation using a single +/// barrel-shifted mask plus a variable-position increment is plausible as a +/// follow-up. +/// +/// `float_utilst::round_to_integral` is the bvt-level mirror of this code; +/// keep the two implementations in lockstep when changing the rounding +/// tables. +exprt float_bvt::round_to_integral( + const exprt &src, + const exprt &rm, + const ieee_float_spect &spec) const +{ + const unbiased_floatt unpacked = unpack(src, spec); + const exprt is_special = + or_exprt(unpacked.zero, or_exprt(unpacked.NaN, unpacked.infinity)); + const exprt exp_ge_f = binary_relation_exprt( + unpacked.exponent, ID_ge, from_integer(spec.f, unpacked.exponent.type())); + + const floatbv_typet float_type = spec.to_type(); + const std::size_t width = spec.width(); + const unsignedbv_typet uint_type{width}; + const exprt src_uint = extractbits_exprt{src, 0, uint_type}; + + const rounding_mode_bitst rounding_mode_bits(rm); + const exprt biased_exp = get_exponent(src, spec); + + // ±0 and ±1 as unsigned bitvectors + ieee_floatt pz{spec, ieee_floatt::ROUND_TO_ZERO, 0}; + ieee_floatt nz{spec, ieee_floatt::ROUND_TO_ZERO, 0}; + nz.set_sign(true); + ieee_floatt p1{spec, ieee_floatt::ROUND_TO_ZERO, 1}; + ieee_floatt n1{spec, ieee_floatt::ROUND_TO_ZERO, -1}; + + const exprt signed_zero_uint = if_exprt( + sign_bit(src), + from_integer(nz.pack(), uint_type), + from_integer(pz.pack(), uint_type)); + const exprt signed_one_uint = if_exprt( + sign_bit(src), + from_integer(n1.pack(), uint_type), + from_integer(p1.pack(), uint_type)); + + // For |x| < 1 (biased exponent < bias): result is ±0 or ±1. + // |x| >= 0.5 iff biased exponent >= bias-1 (for normals). + // |x| == 0.5 iff biased exponent == bias-1 AND fraction == 0. + const exprt bias_m1 = from_integer(spec.bias() - 1, unsignedbv_typet(spec.e)); + const exprt exp_ge_bias_m1 = + binary_relation_exprt(biased_exp, ID_ge, bias_m1); + const exprt exp_eq_bias_m1 = equal_exprt(biased_exp, bias_m1); + const exprt frac_zero = fraction_all_zeros(src, spec); + const exprt abs_eq_half = and_exprt(exp_eq_bias_m1, frac_zero); + const exprt abs_gt_half = and_exprt(exp_ge_bias_m1, not_exprt(abs_eq_half)); + + // clang-format off + const exprt round_up = if_exprt( + rounding_mode_bits.round_to_even, abs_gt_half, + if_exprt(rounding_mode_bits.round_to_away, + or_exprt(abs_gt_half, abs_eq_half), + if_exprt(rounding_mode_bits.round_to_plus_inf, + not_exprt(unpacked.sign), // any positive non-zero rounds up + if_exprt(rounding_mode_bits.round_to_minus_inf, + unpacked.sign, // any negative non-zero rounds down (to -1) + false_exprt())))); + // clang-format on + + exprt result_uint = if_exprt(round_up, signed_one_uint, signed_zero_uint); + + // For each unbiased exponent 0..f-1 (biased: bias..bias+f-1): + for(std::size_t eu = 0; eu < static_cast(spec.f); eu++) + { + const mp_integer be = eu + spec.bias(); + const std::size_t drop = spec.f - eu; + + mp_integer mask_val = power(2, width) - power(2, drop); + const exprt masked = + bitand_exprt(src_uint, from_integer(mask_val, uint_type)); + + const exprt rbit = extractbit_exprt(src_uint, drop - 1); + exprt sticky = false_exprt(); + for(std::size_t i = 0; i + 1 < drop; i++) + sticky = or_exprt(sticky, extractbit_exprt(src_uint, i)); + const exprt lsb = extractbit_exprt(src_uint, drop); + + // clang-format off + const exprt inc = if_exprt( + rounding_mode_bits.round_to_even, + and_exprt(rbit, or_exprt(lsb, sticky)), + if_exprt(rounding_mode_bits.round_to_away, rbit, + if_exprt(rounding_mode_bits.round_to_plus_inf, + and_exprt(not_exprt(unpacked.sign), or_exprt(rbit, sticky)), + if_exprt(rounding_mode_bits.round_to_minus_inf, + and_exprt(unpacked.sign, or_exprt(rbit, sticky)), + false_exprt())))); + // clang-format on + + // Increment: add 1 at position `drop`. The add runs on the full packed + // representation, so for the largest in-loop biased exponent + // (`bias + f - 1`) a fraction-bit carry can propagate into the exponent + // bits and bump the biased exponent to `bias + f`. That is intentional + // and bounded — `bias + f` stays well below the NaN/Inf range + // (`2*bias + 1`), so no saturation is needed. + const exprt inc_val = from_integer(power(2, drop), uint_type); + const exprt branch = if_exprt(inc, plus_exprt(masked, inc_val), masked); + const exprt match = + equal_exprt(biased_exp, from_integer(be, unsignedbv_typet(spec.e))); + result_uint = if_exprt(match, branch, result_uint); + } + + // Bitwise reinterpret unsigned as float via extractbits + const exprt result = extractbits_exprt{result_uint, 0, float_type}; + return if_exprt(or_exprt(is_special, exp_ge_f), src, result); +} + exprt float_bvt::isnormal( const exprt &src, const ieee_float_spect &spec) diff --git a/src/solvers/floatbv/float_bv.h b/src/solvers/floatbv/float_bv.h index be1cb853809..1ef7e92d404 100644 --- a/src/solvers/floatbv/float_bv.h +++ b/src/solvers/floatbv/float_bv.h @@ -93,6 +93,24 @@ class float_bvt const ieee_float_spect &src_spec, const ieee_float_spect &dest_spec) const; + /// Round a floating-point value to integral, IEEE 754 semantics. + /// \param src: packed bit-vector encoding of an IEEE float of the given + /// `spec`. + /// \param rm: rounding mode encoded as `ieee_floatt::rounding_modet`. + /// \param spec: target IEEE format. + /// \return Packed bit-vector encoding of `roundToIntegral(rm, src)`. + /// Special values (NaN, ±Inf, ±0) and values whose unbiased exponent + /// already meets or exceeds the fraction width short-circuit to `src`. + /// \note The encoding uses an f-deep ITE cascade where each branch + /// contains an f-bit adder, giving an O(f^2) expression-tree size. + /// `float_utilst::round_to_integral` is the bvt-level mirror; keep + /// the two implementations in lockstep when changing the rounding + /// tables. + exprt round_to_integral( + const exprt &src, + const exprt &rm, + const ieee_float_spect &spec) const; + // relations enum class relt { LT, LE, EQ, GT, GE }; static exprt diff --git a/src/solvers/floatbv/float_utils.cpp b/src/solvers/floatbv/float_utils.cpp index ba72bb09b21..8d9677e1100 100644 --- a/src/solvers/floatbv/float_utils.cpp +++ b/src/solvers/floatbv/float_utils.cpp @@ -152,34 +152,117 @@ bvt float_utilst::build_constant(const ieee_float_valuet &src) return pack(bias(result)); } +/// Round a floating-point value to integral, IEEE 754 semantics. +/// +/// Bvt-level mirror of `float_bvt::round_to_integral` (see the Doxygen on +/// that method for the full algorithmic description). The two +/// implementations are line-for-line translations between the `exprt` API +/// and the `literalt`/`bvt` API; keep them in lockstep when changing the +/// rounding tables. +/// +/// The encoding builds an `f`-deep `bv_utils.select` cascade where each +/// branch contains an `f`-bit adder, giving an O(f^2) circuit size. For +/// `double` (f=52) that is roughly a 25× growth in node count over the +/// prior add-magic-subtract-magic formulation. A loop-free formulation +/// using a single barrel-shifted mask plus a variable-position increment +/// is plausible as a follow-up. bvt float_utilst::round_to_integral(const bvt &src) { PRECONDITION(src.size() == spec.width()); - // Zero? NaN? Infinity? - auto unpacked = unpack(src); - auto is_special = prop.lor({unpacked.zero, unpacked.NaN, unpacked.infinity}); - - // add 2^f, where f is the number of fraction bits, - // by adding f to the exponent - auto magic_number = ieee_floatt{ - spec, ieee_floatt::rounding_modet::ROUND_TO_ZERO, power(2, spec.f)}; - - auto magic_number_bv = build_constant(magic_number); - - // abs(x) >= magic_number? If so, then there is no fractional part. - literalt ge_magic_number = relation(abs(src), relt::GE, magic_number_bv); - - magic_number_bv.back() = src.back(); // copy sign bit + const unbiased_floatt unpacked = unpack(src); + const literalt is_special = + prop.lor({unpacked.zero, unpacked.NaN, unpacked.infinity}); + + // If unbiased exponent >= f, the number is already integral. + const bvt f_const = bv_utils.build_constant(spec.f, unpacked.exponent.size()); + const literalt exp_ge_f = + !bv_utils.signed_less_than(unpacked.exponent, f_const); + + const bvt biased_exp = get_exponent(src); + + // ±0 and ±1 as packed bitvectors + ieee_floatt pz{spec, ieee_floatt::ROUND_TO_ZERO, 0}; + ieee_floatt nz{spec, ieee_floatt::ROUND_TO_ZERO, 0}; + nz.set_sign(true); + ieee_floatt p1{spec, ieee_floatt::ROUND_TO_ZERO, 1}; + ieee_floatt n1{spec, ieee_floatt::ROUND_TO_ZERO, -1}; + + const bvt signed_zero = + bv_utils.select(sign_bit(src), build_constant(nz), build_constant(pz)); + const bvt signed_one = + bv_utils.select(sign_bit(src), build_constant(n1), build_constant(p1)); + + // For |x| < 1 (biased exponent < bias): result is ±0 or ±1. + // |x| >= 0.5 iff biased exponent >= bias-1. + // |x| == 0.5 iff biased exponent == bias-1 AND fraction == 0. + const bvt bias_m1 = bv_utils.build_constant(spec.bias() - 1, spec.e); + const literalt exp_ge_bm1 = !bv_utils.unsigned_less_than(biased_exp, bias_m1); + const literalt exp_eq_bm1 = bv_utils.equal(biased_exp, bias_m1); + const literalt frac_zero = fraction_all_zeros(src); + const literalt abs_eq_half = prop.land(exp_eq_bm1, frac_zero); + const literalt abs_gt_half = prop.land(exp_ge_bm1, !abs_eq_half); - auto tmp1 = add_sub(src, magic_number_bv, false); + // clang-format off + const literalt round_up = prop.lselect( + rounding_mode_bits.round_to_even, abs_gt_half, + prop.lselect(rounding_mode_bits.round_to_away, + prop.lor(abs_gt_half, abs_eq_half), + prop.lselect(rounding_mode_bits.round_to_plus_inf, + !unpacked.sign, + prop.lselect(rounding_mode_bits.round_to_minus_inf, + unpacked.sign, + const_literal(false))))); + // clang-format on - auto tmp2 = add_sub(tmp1, magic_number_bv, true); + bvt result = bv_utils.select(round_up, signed_one, signed_zero); - // restore the original sign bit - tmp2.back() = src.back(); + // For each unbiased exponent 0..f-1 (biased: bias..bias+f-1): + for(std::size_t eu = 0; eu < static_cast(spec.f); eu++) + { + const mp_integer be = eu + spec.bias(); + const std::size_t drop = spec.f - eu; + + // Mask: clear bottom 'drop' bits + bvt masked = src; + for(std::size_t i = 0; i < drop; i++) + masked[i] = const_literal(false); + + // Round bit, sticky bit, least kept bit + const literalt rbit = src[drop - 1]; + literalt sticky = const_literal(false); + for(std::size_t i = 0; i + 1 < drop; i++) + sticky = prop.lor(sticky, src[i]); + const literalt lsb = src[drop]; + + // clang-format off + const literalt inc = prop.lselect( + rounding_mode_bits.round_to_even, + prop.land(rbit, prop.lor(lsb, sticky)), + prop.lselect(rounding_mode_bits.round_to_away, rbit, + prop.lselect(rounding_mode_bits.round_to_plus_inf, + prop.land(!unpacked.sign, prop.lor(rbit, sticky)), + prop.lselect(rounding_mode_bits.round_to_minus_inf, + prop.land(unpacked.sign, prop.lor(rbit, sticky)), + const_literal(false))))); + // clang-format on + + // Increment: add 1 at position `drop`. The add runs on the full packed + // representation, so for the largest in-loop biased exponent + // (`bias + f - 1`) a fraction-bit carry can propagate into the exponent + // bits and bump the biased exponent to `bias + f`. That is intentional + // and bounded — `bias + f` stays well below the NaN/Inf range + // (`2*bias + 1`), so no saturation is needed. + const bvt inc_val = bv_utils.build_constant(power(2, drop), spec.width()); + const bvt incremented = bv_utils.add(masked, inc_val); + const bvt branch = bv_utils.select(inc, incremented, masked); + + const bvt be_const = bv_utils.build_constant(be, spec.e); + const literalt match = bv_utils.equal(biased_exp, be_const); + result = bv_utils.select(match, branch, result); + } - return bv_utils.select(prop.lor(is_special, ge_magic_number), src, tmp2); + return bv_utils.select(prop.lor(is_special, exp_ge_f), src, result); } bvt float_utilst::conversion( diff --git a/src/solvers/smt2/smt2_conv.cpp b/src/solvers/smt2/smt2_conv.cpp index e7b9893a882..3d0095c9edd 100644 --- a/src/solvers/smt2/smt2_conv.cpp +++ b/src/solvers/smt2/smt2_conv.cpp @@ -3526,7 +3526,13 @@ void smt2_convt::convert_floatbv_round_to_integral( out << ")"; } else - UNEXPECTEDCASE("TODO floatbv_round_to_integral without FPA"); + { + // Routed via the bvfp_set / define-fun cache (see the + // ID_floatbv_round_to_integral entry in convert_expr's id-list near the + // bvfp_set block) so that multiple call sites share a single SMT model + // definition rather than duplicating the float_bvt rewrite. + convert_floatbv(expr); + } } void smt2_convt::convert_struct(const struct_exprt &expr) @@ -5707,6 +5713,7 @@ void smt2_convt::find_symbols(const exprt &expr) expr.id() == ID_floatbv_div || expr.id() == ID_floatbv_fma || expr.id() == ID_floatbv_typecast || + expr.id() == ID_floatbv_round_to_integral || expr.id() == ID_ieee_float_equal || expr.id() == ID_ieee_float_notequal || ((expr.id() == ID_lt || From b44ec4f7a2955d61167e59a1b25085f8ee8abe7b Mon Sep 17 00:00:00 2001 From: Michael Tautschnig Date: Fri, 29 May 2026 15:23:16 +0000 Subject: [PATCH 2/2] round_to_integral: barrel-shifter formulation in float_bvt and float_utilst MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the f-deep ITE cascade introduced in the previous commit with a single barrel-shifter pass: drop = f - max(0, unbiased_exp) (0..f) masked = (src >> drop) << drop rbit = bit (drop - 1) of src lsb = bit (drop) of src sticky = (src << (width + 1 - drop)) != 0 inc_val = 1 << drop result = masked + (inc ? inc_val : 0) The per-rounding-mode increment decision (round_to_even / _away / _plus_inf / _minus_inf / _zero) is unchanged from the loop version. The |x| < 1 and exp_ge_f / is_special branches are likewise unchanged. The fraction-bit carry into the exponent at the largest relevant biased exponent is also unchanged and remains intentional and bounded. This brings the round_to_integral encoding from O(f^2) expression nodes (an f-deep ITE cascade where each branch contains an f-bit adder) down to O(1) at the float_bvt / float_utilst level — a constant number of variable-amount shifts plus a single adder, all at width spec.width(). The bit-blasting / SMT layer expands the variable-amount shifts into O(log f)-depth circuits, so the asymptotic gate count is O(width * log f) rather than O(f^2). float_bvt::round_to_integral and float_utilst::round_to_integral are mirror translations of the same algorithm between the exprt and literalt/bvt APIs; comments mark them as needing to be kept in lockstep, as before. Co-authored-by: Kiro --- src/solvers/floatbv/float_bv.cpp | 142 ++++++++++++++++------------ src/solvers/floatbv/float_bv.h | 7 +- src/solvers/floatbv/float_utils.cpp | 139 ++++++++++++++++----------- 3 files changed, 174 insertions(+), 114 deletions(-) diff --git a/src/solvers/floatbv/float_bv.cpp b/src/solvers/floatbv/float_bv.cpp index a1c447b6246..47262d947b8 100644 --- a/src/solvers/floatbv/float_bv.cpp +++ b/src/solvers/floatbv/float_bv.cpp @@ -504,26 +504,26 @@ exprt float_bvt::conversion( /// Round a floating-point value to integral. /// -/// Strategy: enumerate the possible biased exponents that have a non-trivial -/// fractional part (`bias..bias+f-1`) and, for each one, mask off the bottom -/// `drop = f - unbiased_exp` bits of the packed representation and -/// conditionally increment by `2^drop`. Values with `|x| < 1` are handled -/// separately (result is ±0 or ±1 depending on rounding mode and sign). -/// Special values (NaN, ±Inf, ±0) and values whose unbiased exponent already +/// Strategy: compute `drop = max(0, f - unbiased_exp)` — the number of +/// trailing fractional bits to discard — and form the rounded result as +/// `(src >> drop) << drop` plus a conditional `1 << drop` increment +/// determined by the rounding mode and the round/sticky/lsb bits at +/// variable position `drop`. Values with `|x| < 1` are handled separately +/// (result is ±0 or ±1 depending on rounding mode and sign). Special +/// values (NaN, ±Inf, ±0) and values whose unbiased exponent already /// meets or exceeds the fraction width short-circuit to `src`. /// -/// The encoding builds an `f`-deep ITE cascade where each branch contains an -/// `f`-bit adder, giving an O(f^2) expression-tree size. For `double` (f=52) -/// that is roughly a 25× growth in node count over the prior -/// add-magic-subtract-magic formulation; the trade is correctness on the SMT -/// non-FPA path (which previously hit `UNEXPECTEDCASE`) and a more directly -/// auditable algorithm. A loop-free formulation using a single -/// barrel-shifted mask plus a variable-position increment is plausible as a -/// follow-up. +/// The encoding builds an O(1)-deep expression tree at the float_bvt +/// level — only constant-depth shifts, masks, and a single adder, all +/// at width `spec.width()`. The bit-blasting / SMT-translation layer +/// expands the variable-amount shifts into O(log f)-depth circuits. +/// This replaces an earlier formulation that built an f-deep ITE +/// cascade (one branch per possible biased exponent), which was +/// O(f^2) in expression size — roughly 25× the node count for double. /// -/// `float_utilst::round_to_integral` is the bvt-level mirror of this code; -/// keep the two implementations in lockstep when changing the rounding -/// tables. +/// `float_utilst::round_to_integral` is the bvt-level mirror of this +/// code; keep the two implementations in lockstep when changing the +/// rounding logic. exprt float_bvt::round_to_integral( const exprt &src, const exprt &rm, @@ -539,6 +539,8 @@ exprt float_bvt::round_to_integral( const std::size_t width = spec.width(); const unsignedbv_typet uint_type{width}; const exprt src_uint = extractbits_exprt{src, 0, uint_type}; + const exprt zero_uint = from_integer(0, uint_type); + const exprt one_uint = from_integer(1, uint_type); const rounding_mode_bitst rounding_mode_bits(rm); const exprt biased_exp = get_exponent(src, spec); @@ -571,7 +573,7 @@ exprt float_bvt::round_to_integral( const exprt abs_gt_half = and_exprt(exp_ge_bias_m1, not_exprt(abs_eq_half)); // clang-format off - const exprt round_up = if_exprt( + const exprt round_up_lt1 = if_exprt( rounding_mode_bits.round_to_even, abs_gt_half, if_exprt(rounding_mode_bits.round_to_away, or_exprt(abs_gt_half, abs_eq_half), @@ -582,48 +584,72 @@ exprt float_bvt::round_to_integral( false_exprt())))); // clang-format on - exprt result_uint = if_exprt(round_up, signed_one_uint, signed_zero_uint); + const exprt result_lt1 = + if_exprt(round_up_lt1, signed_one_uint, signed_zero_uint); + + // |x| >= 1 branch: barrel-shifter formulation. + // + // drop = f - max(0, unbiased_exp). For unbiased_exp >= f, the result is + // overridden by the `exp_ge_f -> src` early-out below; for unbiased_exp + // < 0, overridden by the |x| < 1 branch above. So in the path where this + // matters, drop ∈ [1, f]. + const signedbv_typet exp_type = + signedbv_typet{to_signedbv_type(unpacked.exponent.type()).get_width()}; + const exprt unbiased_exp_clamped_lo = if_exprt( + binary_relation_exprt(unpacked.exponent, ID_lt, from_integer(0, exp_type)), + from_integer(0, exp_type), + unpacked.exponent); + const exprt drop_signed = + minus_exprt(from_integer(spec.f, exp_type), unbiased_exp_clamped_lo); + const exprt drop = typecast_exprt(drop_signed, uint_type); + + // masked = (src >> drop) << drop + const exprt shifted_right = lshr_exprt{src_uint, drop}; + const exprt masked = shl_exprt{shifted_right, drop}; + + // rbit = bit (drop - 1) of src; obtained by shifting right by drop-1 + // and reading the LSB. When drop == 0 (only on the exp_ge_f path) the + // value is irrelevant, since the final select picks `src`. + const exprt drop_m1 = minus_exprt(drop, one_uint); + const exprt shifted_to_rbit = lshr_exprt{src_uint, drop_m1}; + const exprt rbit = + equal_exprt(bitand_exprt(shifted_to_rbit, one_uint), one_uint); + + // lsb = bit (drop) of src + const exprt lsb = + equal_exprt(bitand_exprt(shifted_right, one_uint), one_uint); + + // sticky = OR of bits 0..(drop-2) ≡ (src AND ((1 << (drop-1)) - 1)) != 0 + const exprt sticky_mask = minus_exprt(shl_exprt{one_uint, drop_m1}, one_uint); + const exprt sticky = + notequal_exprt(bitand_exprt(src_uint, sticky_mask), zero_uint); - // For each unbiased exponent 0..f-1 (biased: bias..bias+f-1): - for(std::size_t eu = 0; eu < static_cast(spec.f); eu++) - { - const mp_integer be = eu + spec.bias(); - const std::size_t drop = spec.f - eu; - - mp_integer mask_val = power(2, width) - power(2, drop); - const exprt masked = - bitand_exprt(src_uint, from_integer(mask_val, uint_type)); - - const exprt rbit = extractbit_exprt(src_uint, drop - 1); - exprt sticky = false_exprt(); - for(std::size_t i = 0; i + 1 < drop; i++) - sticky = or_exprt(sticky, extractbit_exprt(src_uint, i)); - const exprt lsb = extractbit_exprt(src_uint, drop); + // clang-format off + const exprt inc = if_exprt( + rounding_mode_bits.round_to_even, + and_exprt(rbit, or_exprt(lsb, sticky)), + if_exprt(rounding_mode_bits.round_to_away, rbit, + if_exprt(rounding_mode_bits.round_to_plus_inf, + and_exprt(not_exprt(unpacked.sign), or_exprt(rbit, sticky)), + if_exprt(rounding_mode_bits.round_to_minus_inf, + and_exprt(unpacked.sign, or_exprt(rbit, sticky)), + false_exprt())))); + // clang-format on - // clang-format off - const exprt inc = if_exprt( - rounding_mode_bits.round_to_even, - and_exprt(rbit, or_exprt(lsb, sticky)), - if_exprt(rounding_mode_bits.round_to_away, rbit, - if_exprt(rounding_mode_bits.round_to_plus_inf, - and_exprt(not_exprt(unpacked.sign), or_exprt(rbit, sticky)), - if_exprt(rounding_mode_bits.round_to_minus_inf, - and_exprt(unpacked.sign, or_exprt(rbit, sticky)), - false_exprt())))); - // clang-format on - - // Increment: add 1 at position `drop`. The add runs on the full packed - // representation, so for the largest in-loop biased exponent - // (`bias + f - 1`) a fraction-bit carry can propagate into the exponent - // bits and bump the biased exponent to `bias + f`. That is intentional - // and bounded — `bias + f` stays well below the NaN/Inf range - // (`2*bias + 1`), so no saturation is needed. - const exprt inc_val = from_integer(power(2, drop), uint_type); - const exprt branch = if_exprt(inc, plus_exprt(masked, inc_val), masked); - const exprt match = - equal_exprt(biased_exp, from_integer(be, unsignedbv_typet(spec.e))); - result_uint = if_exprt(match, branch, result_uint); - } + // Increment by `1 << drop`. The add runs on the full packed + // representation, so a fraction-bit carry can propagate into the + // exponent bits, bumping the biased exponent by one. That is + // intentional and bounded: at the largest relevant biased exponent + // (`bias + f - 1`), the result reaches `bias + f`, which still stays + // well below the NaN/Inf range (`2*bias + 1`), so no saturation is + // needed. + const exprt inc_val = shl_exprt{one_uint, drop}; + const exprt result_ge1 = if_exprt(inc, plus_exprt(masked, inc_val), masked); + + // Select between the two branches by exponent magnitude. + const exprt bias_e = from_integer(spec.bias(), unsignedbv_typet(spec.e)); + const exprt abs_lt_1 = binary_relation_exprt(biased_exp, ID_lt, bias_e); + const exprt result_uint = if_exprt(abs_lt_1, result_lt1, result_ge1); // Bitwise reinterpret unsigned as float via extractbits const exprt result = extractbits_exprt{result_uint, 0, float_type}; diff --git a/src/solvers/floatbv/float_bv.h b/src/solvers/floatbv/float_bv.h index 1ef7e92d404..4de8e22ae5e 100644 --- a/src/solvers/floatbv/float_bv.h +++ b/src/solvers/floatbv/float_bv.h @@ -101,8 +101,11 @@ class float_bvt /// \return Packed bit-vector encoding of `roundToIntegral(rm, src)`. /// Special values (NaN, ±Inf, ±0) and values whose unbiased exponent /// already meets or exceeds the fraction width short-circuit to `src`. - /// \note The encoding uses an f-deep ITE cascade where each branch - /// contains an f-bit adder, giving an O(f^2) expression-tree size. + /// \note The encoding is a barrel-shifter pass: a constant-depth + /// expression tree (a handful of variable-amount shifts plus a + /// single adder, all at width `spec.width()`). The bit-blasting / + /// SMT layer expands the variable-amount shifts into + /// O(log f)-depth circuits. /// `float_utilst::round_to_integral` is the bvt-level mirror; keep /// the two implementations in lockstep when changing the rounding /// tables. diff --git a/src/solvers/floatbv/float_utils.cpp b/src/solvers/floatbv/float_utils.cpp index 8d9677e1100..52f4d903ff2 100644 --- a/src/solvers/floatbv/float_utils.cpp +++ b/src/solvers/floatbv/float_utils.cpp @@ -158,14 +158,15 @@ bvt float_utilst::build_constant(const ieee_float_valuet &src) /// that method for the full algorithmic description). The two /// implementations are line-for-line translations between the `exprt` API /// and the `literalt`/`bvt` API; keep them in lockstep when changing the -/// rounding tables. +/// rounding logic. /// -/// The encoding builds an `f`-deep `bv_utils.select` cascade where each -/// branch contains an `f`-bit adder, giving an O(f^2) circuit size. For -/// `double` (f=52) that is roughly a 25× growth in node count over the -/// prior add-magic-subtract-magic formulation. A loop-free formulation -/// using a single barrel-shifted mask plus a variable-position increment -/// is plausible as a follow-up. +/// The encoding builds an O(1)-deep computation at the bvt level — a +/// constant number of variable-amount shifts, masks, and a single adder, +/// all at width `spec.width()`. The bit-blasting layer expands the +/// variable-amount shifts into O(log f)-depth circuits. This replaces an +/// earlier formulation that built an f-deep `bv_utils.select` cascade +/// (one branch per possible biased exponent), which was O(f^2) in +/// circuit size. bvt float_utilst::round_to_integral(const bvt &src) { PRECONDITION(src.size() == spec.width()); @@ -175,7 +176,8 @@ bvt float_utilst::round_to_integral(const bvt &src) prop.lor({unpacked.zero, unpacked.NaN, unpacked.infinity}); // If unbiased exponent >= f, the number is already integral. - const bvt f_const = bv_utils.build_constant(spec.f, unpacked.exponent.size()); + const std::size_t exp_width = unpacked.exponent.size(); + const bvt f_const = bv_utils.build_constant(spec.f, exp_width); const literalt exp_ge_f = !bv_utils.signed_less_than(unpacked.exponent, f_const); @@ -204,7 +206,7 @@ bvt float_utilst::round_to_integral(const bvt &src) const literalt abs_gt_half = prop.land(exp_ge_bm1, !abs_eq_half); // clang-format off - const literalt round_up = prop.lselect( + const literalt round_up_lt1 = prop.lselect( rounding_mode_bits.round_to_even, abs_gt_half, prop.lselect(rounding_mode_bits.round_to_away, prop.lor(abs_gt_half, abs_eq_half), @@ -215,52 +217,81 @@ bvt float_utilst::round_to_integral(const bvt &src) const_literal(false))))); // clang-format on - bvt result = bv_utils.select(round_up, signed_one, signed_zero); + const bvt result_lt1 = bv_utils.select(round_up_lt1, signed_one, signed_zero); + + // |x| >= 1 branch: barrel-shifter formulation. + // + // drop = f - max(0, unbiased_exp). For unbiased_exp >= f, the result is + // overridden by `exp_ge_f -> src`; for unbiased_exp < 0, overridden by + // the |x| < 1 branch above. So in the path where this matters, + // drop ∈ [1, f]. + const std::size_t width = spec.width(); + const literalt unbiased_neg = bv_utils.signed_less_than( + unpacked.exponent, bv_utils.build_constant(0, exp_width)); + const bvt unbiased_clamped_lo = bv_utils.select( + unbiased_neg, bv_utils.build_constant(0, exp_width), unpacked.exponent); + const bvt drop_exp_width = bv_utils.sub( + bv_utils.build_constant(spec.f, exp_width), unbiased_clamped_lo); + const bvt drop = bv_utils.zero_extension(drop_exp_width, width); + + const bvt one_bv = bv_utils.build_constant(1, width); + + // masked = (src >> drop) << drop + const bvt shifted_right = + bv_utils.shift(src, bv_utilst::shiftt::SHIFT_LRIGHT, drop); + const bvt masked = + bv_utils.shift(shifted_right, bv_utilst::shiftt::SHIFT_LEFT, drop); + + // rbit = bit (drop - 1) of src; obtained by shifting right by drop-1 + // and reading the LSB. When drop == 0 (only on the exp_ge_f path) the + // value is irrelevant since the final select picks `src`. + const bvt drop_m1 = bv_utils.sub(drop, one_bv); + const bvt shifted_to_rbit = + bv_utils.shift(src, bv_utilst::shiftt::SHIFT_LRIGHT, drop_m1); + const literalt rbit = shifted_to_rbit[0]; + + // lsb = bit (drop) of src + const literalt lsb = shifted_right[0]; + + // sticky = OR of bits 0..(drop-2) of src. + // Push those bits into the top of a width-bit value via a left shift + // by `width + 1 - drop`, then test for non-zero. The round bit + // (position drop-1) is shifted out; bits at and above drop are also + // out. + const bvt width_p1 = bv_utils.build_constant(width + 1, width); + const bvt sticky_shift_amount = bv_utils.sub(width_p1, drop); + const bvt sticky_pushed = + bv_utils.shift(src, bv_utilst::shiftt::SHIFT_LEFT, sticky_shift_amount); + const literalt sticky = !bv_utils.is_zero(sticky_pushed); - // For each unbiased exponent 0..f-1 (biased: bias..bias+f-1): - for(std::size_t eu = 0; eu < static_cast(spec.f); eu++) - { - const mp_integer be = eu + spec.bias(); - const std::size_t drop = spec.f - eu; - - // Mask: clear bottom 'drop' bits - bvt masked = src; - for(std::size_t i = 0; i < drop; i++) - masked[i] = const_literal(false); - - // Round bit, sticky bit, least kept bit - const literalt rbit = src[drop - 1]; - literalt sticky = const_literal(false); - for(std::size_t i = 0; i + 1 < drop; i++) - sticky = prop.lor(sticky, src[i]); - const literalt lsb = src[drop]; - - // clang-format off - const literalt inc = prop.lselect( - rounding_mode_bits.round_to_even, - prop.land(rbit, prop.lor(lsb, sticky)), - prop.lselect(rounding_mode_bits.round_to_away, rbit, - prop.lselect(rounding_mode_bits.round_to_plus_inf, - prop.land(!unpacked.sign, prop.lor(rbit, sticky)), - prop.lselect(rounding_mode_bits.round_to_minus_inf, - prop.land(unpacked.sign, prop.lor(rbit, sticky)), - const_literal(false))))); - // clang-format on - - // Increment: add 1 at position `drop`. The add runs on the full packed - // representation, so for the largest in-loop biased exponent - // (`bias + f - 1`) a fraction-bit carry can propagate into the exponent - // bits and bump the biased exponent to `bias + f`. That is intentional - // and bounded — `bias + f` stays well below the NaN/Inf range - // (`2*bias + 1`), so no saturation is needed. - const bvt inc_val = bv_utils.build_constant(power(2, drop), spec.width()); - const bvt incremented = bv_utils.add(masked, inc_val); - const bvt branch = bv_utils.select(inc, incremented, masked); - - const bvt be_const = bv_utils.build_constant(be, spec.e); - const literalt match = bv_utils.equal(biased_exp, be_const); - result = bv_utils.select(match, branch, result); - } + // clang-format off + const literalt inc = prop.lselect( + rounding_mode_bits.round_to_even, + prop.land(rbit, prop.lor(lsb, sticky)), + prop.lselect(rounding_mode_bits.round_to_away, rbit, + prop.lselect(rounding_mode_bits.round_to_plus_inf, + prop.land(!unpacked.sign, prop.lor(rbit, sticky)), + prop.lselect(rounding_mode_bits.round_to_minus_inf, + prop.land(unpacked.sign, prop.lor(rbit, sticky)), + const_literal(false))))); + // clang-format on + + // Increment by `1 << drop`. The add runs on the full packed + // representation, so a fraction-bit carry can propagate into the + // exponent bits, bumping the biased exponent by one. That is + // intentional and bounded: at the largest relevant biased exponent + // (`bias + f - 1`), the result reaches `bias + f`, which still stays + // well below the NaN/Inf range (`2*bias + 1`), so no saturation is + // needed. + const bvt inc_val = + bv_utils.shift(one_bv, bv_utilst::shiftt::SHIFT_LEFT, drop); + const bvt incremented = bv_utils.add(masked, inc_val); + const bvt result_ge1 = bv_utils.select(inc, incremented, masked); + + // Select between the two branches by exponent magnitude. + const bvt bias_e = bv_utils.build_constant(spec.bias(), spec.e); + const literalt abs_lt_1 = bv_utils.unsigned_less_than(biased_exp, bias_e); + const bvt result = bv_utils.select(abs_lt_1, result_lt1, result_ge1); return bv_utils.select(prop.lor(is_special, exp_ge_f), src, result); }