Add round_to_integral to float_bvt, rewrite in float_utilst

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 <assert.h>
+#include <float.h>
+#include <math.h>
+
+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;
+}

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.

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.

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,160 @@ exprt float_bvt::conversion(
   }
 }
 
+/// Round a floating-point value to integral.
+///
+/// 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 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 logic.
+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 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);
+
+  // ±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_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),
+    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
+
+  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);
+
+  // 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 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};
+  return if_exprt(or_exprt(is_special, exp_ge_f), src, result);
+}
+
 exprt float_bvt::isnormal(
   const exprt &src,
   const ieee_float_spect &spec)

src/solvers/floatbv/float_bv.h

@@ -93,6 +93,27 @@ 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 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.
+  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