Fix subnormal division precision loss in float_utilst and float_bvt

regression/cbmc/Float-div-subnormal/main.c

@@ -0,0 +1,47 @@
+// Test division with a subnormal dividend.
+//
+// Bug: float_utilst and float_bvt used an insufficient division width,
+// causing 1-ULP errors when the unpacked dividend had many leading zeros
+// (which happens for subnormal operands).  After the rounder's
+// normalization shift left-shifts the quotient to place the leading 1 at
+// the top of the buffer, the bits below the round position were all zero,
+// so the `have_remainder` sticky bit was effectively lost: a true
+// round-up (round=1, sticky=1) was mistaken for a tie (round=1, sticky=0)
+// and broken to even.  The error only surfaces when the quotient lands
+// back in the normal range and needs rounding, i.e. for subnormal/subnormal
+// division; a subnormal divided by a normal yields an (exact) subnormal or
+// zero and does not trigger it.
+
+#include <assert.h>
+#include <float.h>
+
+int main()
+{
+  // Smallest positive subnormal divided by three times the smallest
+  // subnormal: mathematically 1/3, whose IEEE-754 round-to-nearest-even
+  // value is 0x1.555556p-2.  The buggy build returned 0x1.555554p-2
+  // (1 ULP low) because the sticky bit was lost.
+  float a, b;
+  __CPROVER_assume(a == 0x1p-149f);
+  __CPROVER_assume(b == 0x3p-149f);
+  float r = a / b;
+  assert(r == 0x1.555556p-2f);
+
+  // A second subnormal/subnormal case, 7/4 = 1.75 exactly.
+  float c, d;
+  __CPROVER_assume(c == 0x7p-149f);
+  __CPROVER_assume(d == 0x4p-149f);
+  float r2 = c / d;
+  assert(r2 == 0x1.cp+0f);
+
+  // Exact subnormal/normal sanity check: FLT_MIN is the smallest positive
+  // normal float; FLT_MIN * 0.5f is subnormal, and dividing it by 2 stays
+  // exactly representable.
+  float e, f;
+  __CPROVER_assume(e == FLT_MIN * 0.5f);
+  __CPROVER_assume(f == 2.0f);
+  float r3 = e / f;
+  assert(r3 == FLT_MIN * 0.25f);
+
+  return 0;
+}

regression/cbmc/Float-div-subnormal/test.desc

@@ -0,0 +1,10 @@
+CORE no-new-smt
+main.c
+
+^EXIT=0$
+^SIGNAL=0$
+^VERIFICATION SUCCESSFUL$
+--
+^warning: ignoring
+--
+Division with subnormal dividend.

regression/cbmc/Float-div-subnormal/test_smt.desc

@@ -0,0 +1,10 @@
+CORE no-new-smt
+main.c
+--smt2 --z3
+^EXIT=0$
+^SIGNAL=0$
+^VERIFICATION SUCCESSFUL$
+--
+^warning: ignoring
+--
+Division with subnormal dividend (SMT path).

src/solvers/floatbv/float_bv.cpp

@@ -757,11 +757,24 @@ exprt float_bvt::div(
 
   std::size_t fraction_width=
     to_unsignedbv_type(unpacked1.fraction.type()).get_width();
-  std::size_t div_width=fraction_width*2+1;
+
+  // Pre-normalize the dividend.  A subnormal dividend has leading zeros in
+  // its fraction (the hidden bit is 0); dividing it directly would give a
+  // quotient with too few significant bits, so the `have_remainder` sticky
+  // bit is lost during the rounder's normalization shift and ties round the
+  // wrong way (a 1-ULP error).  Left-aligning the fraction and decreasing the
+  // exponent by the same amount preserves the value but guarantees a
+  // normalized dividend, so the minimal division width suffices and
+  // normal/normal division pays no extra cost.
+  exprt dividend_fraction = unpacked1.fraction;
+  exprt dividend_exponent = unpacked1.exponent;
+  normalization_shift(dividend_fraction, dividend_exponent);
+
+  std::size_t div_width = fraction_width * 2 + 1;
 
   // pad fraction1 with zeros
   const concatenation_exprt fraction1(
-    unpacked1.fraction,
+    dividend_fraction,
     from_integer(0, unsignedbv_typet(div_width - fraction_width)),
     unsignedbv_typet(div_width));
 
@@ -785,12 +798,12 @@ exprt float_bvt::div(
     concatenation_exprt(
       result.fraction, have_remainder, unsignedbv_typet(div_width+1));
 
-  // We will subtract the exponents;
-  // to account for overflow, we add a bit.
-  const typecast_exprt exponent1(
-    unpacked1.exponent, signedbv_typet(spec.e + 1));
+  // We will subtract the exponents; allow two extra bits so the more
+  // negative exponent of a freshly normalized (formerly subnormal) dividend
+  // cannot overflow the subtraction.
+  const typecast_exprt exponent1(dividend_exponent, signedbv_typet(spec.e + 2));
   const typecast_exprt exponent2(
-    unpacked2.exponent, signedbv_typet(spec.e + 1));
+    unpacked2.exponent, signedbv_typet(spec.e + 2));
 
   // subtract exponents
   const minus_exprt added_exponent(exponent1, exponent2);

src/solvers/floatbv/float_utils.cpp

@@ -615,10 +615,22 @@ bvt float_utilst::div(const bvt &src1, const bvt &src2)
   const unbiased_floatt unpacked1=unpack(src1);
   const unbiased_floatt unpacked2=unpack(src2);
 
-  std::size_t div_width=unpacked1.fraction.size()*2+1;
+  // Pre-normalize the dividend.  A subnormal dividend has leading zeros in
+  // its fraction (the hidden bit is 0); dividing it directly would give a
+  // quotient with too few significant bits, so the `have_remainder` sticky
+  // bit is lost during the rounder's normalization shift and ties round the
+  // wrong way (a 1-ULP error).  Left-aligning the fraction and decreasing the
+  // exponent by the same amount preserves the value but guarantees a
+  // normalized dividend, so the minimal division width suffices and
+  // normal/normal division pays no extra cost.
+  bvt dividend_fraction = unpacked1.fraction;
+  bvt dividend_exponent = unpacked1.exponent;
+  normalization_shift(dividend_fraction, dividend_exponent);
+
+  std::size_t div_width = dividend_fraction.size() * 2 + 1;
 
   // pad fraction1 with zeros
-  bvt fraction1=unpacked1.fraction;
+  bvt fraction1 = dividend_fraction;
   fraction1.reserve(div_width);
   while(fraction1.size()<div_width)
     fraction1.insert(fraction1.begin(), const_literal(false));
@@ -643,10 +655,12 @@ bvt float_utilst::div(const bvt &src1, const bvt &src2)
   // We will subtract the exponents;
   // to account for overflow, we add a bit.
   // we add a second bit for the adjust by extra fraction bits
-  const bvt exponent1=
-    bv_utils.sign_extension(unpacked1.exponent, unpacked1.exponent.size()+2);
-  const bvt exponent2=
-    bv_utils.sign_extension(unpacked2.exponent, unpacked2.exponent.size()+2);
+  // The dividend exponent may have grown during normalization; extend both
+  // exponents to a common width.
+  const std::size_t exp_width =
+    std::max(dividend_exponent.size(), unpacked2.exponent.size()) + 2;
+  const bvt exponent1 = bv_utils.sign_extension(dividend_exponent, exp_width);
+  const bvt exponent2 = bv_utils.sign_extension(unpacked2.exponent, exp_width);
 
   // subtract exponents
   bvt added_exponent=bv_utils.sub(exponent1, exponent2);