From 2b07fc2cd7aefe642a44b307c75a19878aac8568 Mon Sep 17 00:00:00 2001 From: Michael Tautschnig Date: Mon, 27 Apr 2026 19:12:07 +0000 Subject: [PATCH] Fix subnormal division precision loss in float_utilst and float_bvt When dividing a subnormal by a normal, the subnormal's fraction has leading zeros (no hidden bit), causing the quotient to have fewer significant bits than needed for correct rounding (1-ULP errors). Fix: increase div_width by spec.f extra bits to compensate for the worst-case leading zeros in a subnormal fraction. Co-authored-by: Kiro --- regression/cbmc/Float-div-subnormal/main.c | 32 +++++++++++++++++++ regression/cbmc/Float-div-subnormal/test.desc | 10 ++++++ .../cbmc/Float-div-subnormal/test_smt.desc | 10 ++++++ src/solvers/floatbv/float_bv.cpp | 5 ++- src/solvers/floatbv/float_utils.cpp | 5 ++- 5 files changed, 60 insertions(+), 2 deletions(-) create mode 100644 regression/cbmc/Float-div-subnormal/main.c create mode 100644 regression/cbmc/Float-div-subnormal/test.desc create mode 100644 regression/cbmc/Float-div-subnormal/test_smt.desc diff --git a/regression/cbmc/Float-div-subnormal/main.c b/regression/cbmc/Float-div-subnormal/main.c new file mode 100644 index 00000000000..1c7b02e41a5 --- /dev/null +++ b/regression/cbmc/Float-div-subnormal/main.c @@ -0,0 +1,32 @@ +// Test division with subnormal dividend. +// +// Bug: float_utilst and float_bvt had insufficient division width, +// causing 1-ULP errors when the dividend is subnormal (leading zeros +// in the fraction reduce the effective precision of the quotient). + +#include +#include +#include + +int main() +{ + // A subnormal divided by a normal. + // FLT_MIN / 2 is the smallest normal; values below are subnormal. + float a, b; + __CPROVER_assume(a == FLT_MIN * 0.5f); // subnormal + __CPROVER_assume(b == 2.0f); + + float r = a / b; + // The result should be exactly FLT_MIN * 0.25, also subnormal. + assert(r == FLT_MIN * 0.25f); + + // Another case: subnormal / large normal + float c, d; + __CPROVER_assume(c == FLT_MIN * 0.75f); + __CPROVER_assume(d == 3.0f); + + float r2 = c / d; + assert(r2 == 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..3bc31034cfb --- /dev/null +++ b/regression/cbmc/Float-div-subnormal/test.desc @@ -0,0 +1,10 @@ +CORE no-new-smt +main.c +--floatbv +^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..77db19ffb41 100644 --- a/src/solvers/floatbv/float_bv.cpp +++ b/src/solvers/floatbv/float_bv.cpp @@ -757,7 +757,10 @@ 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; + // Division width: we need enough bits for the quotient to have + // full precision even when the dividend is subnormal. A subnormal + // has up to f leading zeros in the fraction, so we add f extra bits. + std::size_t div_width = fraction_width * 2 + 1 + spec.f; // pad fraction1 with zeros const concatenation_exprt fraction1( diff --git a/src/solvers/floatbv/float_utils.cpp b/src/solvers/floatbv/float_utils.cpp index 13d2ccded79..cfb2d8e5eaf 100644 --- a/src/solvers/floatbv/float_utils.cpp +++ b/src/solvers/floatbv/float_utils.cpp @@ -615,7 +615,10 @@ 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; + // Division width: we need enough bits for the quotient to have + // full precision even when the dividend is subnormal. A subnormal + // has up to f leading zeros in the fraction, so we add f extra bits. + std::size_t div_width = unpacked1.fraction.size() * 2 + 1 + spec.f; // pad fraction1 with zeros bvt fraction1=unpacked1.fraction;