diff --git a/regression/cbmc/Float-div-subnormal/main.c b/regression/cbmc/Float-div-subnormal/main.c new file mode 100644 index 00000000000..b33709fb944 --- /dev/null +++ b/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 +#include + +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; +} diff --git a/regression/cbmc/Float-div-subnormal/test.desc b/regression/cbmc/Float-div-subnormal/test.desc new file mode 100644 index 00000000000..07192cbf394 --- /dev/null +++ b/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. diff --git a/regression/cbmc/Float-div-subnormal/test_smt.desc b/regression/cbmc/Float-div-subnormal/test_smt.desc new file mode 100644 index 00000000000..7e43bb5601d --- /dev/null +++ b/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). diff --git a/src/solvers/floatbv/float_bv.cpp b/src/solvers/floatbv/float_bv.cpp index ee657069fe6..9e057e7846e 100644 --- a/src/solvers/floatbv/float_bv.cpp +++ b/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); diff --git a/src/solvers/floatbv/float_utils.cpp b/src/solvers/floatbv/float_utils.cpp index ba72bb09b21..db09b384a27 100644 --- a/src/solvers/floatbv/float_utils.cpp +++ b/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()