diff --git a/regression/cbmc/Float-fma-sign/main.c b/regression/cbmc/Float-fma-sign/main.c new file mode 100644 index 00000000000..9a0df12c44d --- /dev/null +++ b/regression/cbmc/Float-fma-sign/main.c @@ -0,0 +1,120 @@ +// Test that FMA correctly computes the sign of infinity and zero results. +// +// Bug: float_bvt::fma always used the add/sub sign as the result sign, +// ignoring the special-case rules for infinity and zero. As a result: +// - For infinity results, the sign could be wrong when the product is +// infinite but the addend is finite (or vice versa), because the +// product sign and addend sign differ. +// - For exact-zero results (cancellation), the signed-zero conventions +// for the rounding mode were not respected. In particular, under +// round-to-minus-inf (FE_DOWNWARD), +0 + (-0) must yield -0, but the +// pre-fix code could produce +0. +// +// All cases below call plain `fmaf` (not `__CPROVER_fmaf`) so that the +// library wrapper is exercised too: with the tightened wrapper that +// only raises FE_INVALID for actual inf - inf, none of these inputs +// trigger a spurious exception. This also gives the test free +// coverage of the wrapper's flag handling per IEEE-754. + +#include +#include + +#ifndef _MSC_VER +# include +#endif + +extern int __CPROVER_rounding_mode; + +// Force an exact-zero FMA result by exact cancellation between the +// product and the addend, then check the resulting signbit. +void testFmaZero(int mode, float x, float y, float z, int expected_sign) +{ +#ifndef _MSC_VER + int error = fesetround(mode); + assert(error == 0); + + float r = fmaf(x, y, z); + assert(r == 0.0f); + assert(signbit(r) == expected_sign); +#endif + + return; +} + +int main() +{ + float a, b, c; + + // Test 1: (+inf) * (+2) + 42 = +inf (sign from product) + __CPROVER_assume(isinf(a) && a > 0); + __CPROVER_assume(b == 2.0f); + __CPROVER_assume(c == 42.0f); + + __CPROVER_rounding_mode = 0; + float r1 = fmaf(a, b, c); + assert(isinf(r1) && r1 > 0); + + // Test 2: (+inf) * (-3) + 100 = -inf (sign from product) + float d, e, f; + __CPROVER_assume(isinf(d) && d > 0); + __CPROVER_assume(e == -3.0f); + __CPROVER_assume(f == 100.0f); + + __CPROVER_rounding_mode = 0; + float r2 = fmaf(d, e, f); + assert(isinf(r2) && r2 < 0); + + // Test 3: 2 * 3 + (+inf) = +inf (sign from addend) + float g, h, i; + __CPROVER_assume(g == 2.0f); + __CPROVER_assume(h == 3.0f); + __CPROVER_assume(isinf(i) && i > 0); + + __CPROVER_rounding_mode = 0; + float r3 = fmaf(g, h, i); + assert(isinf(r3) && r3 > 0); + + // Test 4: 2 * 3 + (-inf) = -inf (sign from addend) + float j, k, l; + __CPROVER_assume(j == 2.0f); + __CPROVER_assume(k == 3.0f); + __CPROVER_assume(isinf(l) && l < 0); + + __CPROVER_rounding_mode = 0; + float r4 = fmaf(j, k, l); + assert(isinf(r4) && r4 < 0); + +#ifndef _MSC_VER + // Zero-result sign tests, exercising the new zero_sign computation in + // float_bvt::fma (the OR-of-signs term under FE_DOWNWARD vs the + // AND-of-signs term in other modes). + + // Exact cancellation: 1 * 1 + (-1) = 0 + // FE_TONEAREST yields +0 (AND of signs: +,- -> +). + testFmaZero(FE_TONEAREST, 1.0f, 1.0f, -1.0f, 0); + // FE_DOWNWARD yields -0 (OR of signs: +,- -> -). + testFmaZero(FE_DOWNWARD, 1.0f, 1.0f, -1.0f, 1); + + // Exact cancellation with negated product: (-1) * 1 + 1 = 0 + // FE_TONEAREST yields +0 (AND of signs: -,+ -> +). + testFmaZero(FE_TONEAREST, -1.0f, 1.0f, 1.0f, 0); + // FE_DOWNWARD yields -0 (OR of signs: -,+ -> -). + testFmaZero(FE_DOWNWARD, -1.0f, 1.0f, 1.0f, 1); + + // Both contributing signs negative: (-0) * 1 + (-0) = -0 + // The product's sign is the XOR of the factor signs, so (-0) * 1 has + // sign -, and addend has sign -; AND yields -, OR yields -. + testFmaZero(FE_TONEAREST, -0.0f, 1.0f, -0.0f, 1); + testFmaZero(FE_DOWNWARD, -0.0f, 1.0f, -0.0f, 1); + + // Both contributing signs positive: 0 * 1 + 0 = +0 in every mode. + testFmaZero(FE_TONEAREST, 0.0f, 1.0f, 0.0f, 0); + testFmaZero(FE_DOWNWARD, 0.0f, 1.0f, 0.0f, 0); + + // Mixed: (+0) * 1 + (-0) = +0 under FE_TONEAREST, -0 under FE_DOWNWARD. + testFmaZero(FE_TONEAREST, 0.0f, 1.0f, -0.0f, 0); + testFmaZero(FE_DOWNWARD, 0.0f, 1.0f, -0.0f, 1); +#endif + + return 0; +} diff --git a/regression/cbmc/Float-fma-sign/test.desc b/regression/cbmc/Float-fma-sign/test.desc new file mode 100644 index 00000000000..e1fd0d03593 --- /dev/null +++ b/regression/cbmc/Float-fma-sign/test.desc @@ -0,0 +1,10 @@ +CORE broken-cprover-smt-backend thorough-paths no-new-smt +main.c +--floatbv +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +FMA sign handling for infinity and zero results. diff --git a/regression/cbmc/Float-fma-sign/test_smt.desc b/regression/cbmc/Float-fma-sign/test_smt.desc new file mode 100644 index 00000000000..2624e55fc58 --- /dev/null +++ b/regression/cbmc/Float-fma-sign/test_smt.desc @@ -0,0 +1,10 @@ +CORE broken-cprover-smt-backend thorough-paths no-new-smt +main.c +--smt2 --z3 +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +FMA sign handling for infinity and zero results (SMT path). diff --git a/src/ansi-c/library/math.c b/src/ansi-c/library/math.c index 32b6afa9fa1..988edc4aa77 100644 --- a/src/ansi-c/library/math.c +++ b/src/ansi-c/library/math.c @@ -3488,14 +3488,22 @@ long double powl(long double x, long double y) double fma(double x, double y, double z) { - // IEEE 754: raise FE_INVALID for 0*inf or inf+(-inf) + // IEEE 754: raise FE_INVALID for 0*inf or inf*0 (regardless of z). if( (isinf(x) || isinf(y)) && (fpclassify(x) == FP_ZERO || fpclassify(y) == FP_ZERO)) { feraiseexcept(FE_INVALID); } - else if(isinf(z)) // this is an over-approximation + // IEEE 754: raise FE_INVALID for the inf - inf subcase, i.e. when + // the product x*y is infinite (one factor is infinite and neither + // is zero -- the 0*inf case is handled above) and the addend z is + // infinite with the opposite sign. Plain finite_nonzero*finite_nonzero + // + (+/-inf), or inf*finite_nonzero + same-signed-inf, is exact and + // raises no flag. + else if( + (isinf(x) || isinf(y)) && isinf(z) && + ((signbit(x) != 0) != (signbit(y) != 0)) != (signbit(z) != 0)) { feraiseexcept(FE_INVALID); } @@ -3517,13 +3525,21 @@ double fma(double x, double y, double z) float fmaf(float x, float y, float z) { + // IEEE 754: raise FE_INVALID for 0*inf or inf*0 (regardless of z). + // Use the type-generic isinf() macro: some headers (e.g. + // macOS-14 arm64) only expose the C99 generic form, not the typed + // isinff()/isinfl() variants. if( - (isinff(x) || isinff(y)) && + (isinf(x) || isinf(y)) && (fpclassify(x) == FP_ZERO || fpclassify(y) == FP_ZERO)) { feraiseexcept(FE_INVALID); } - else if(isinff(z)) // this is an over-approximation + // IEEE 754: raise FE_INVALID for the inf - inf subcase; see fma() + // above for the rationale. + else if( + (isinf(x) || isinf(y)) && isinf(z) && + ((signbit(x) != 0) != (signbit(y) != 0)) != (signbit(z) != 0)) { feraiseexcept(FE_INVALID); } @@ -3545,13 +3561,20 @@ float fmaf(float x, float y, float z) long double fmal(long double x, long double y, long double z) { + // IEEE 754: raise FE_INVALID for 0*inf or inf*0 (regardless of z). + // Use the type-generic isinf() macro for portability; see fmaf() + // above. if( - (isinfl(x) || isinfl(y)) && + (isinf(x) || isinf(y)) && (fpclassify(x) == FP_ZERO || fpclassify(y) == FP_ZERO)) { feraiseexcept(FE_INVALID); } - else if(isinfl(z)) // this is an over-approximation + // IEEE 754: raise FE_INVALID for the inf - inf subcase; see fma() + // above for the rationale. + else if( + (isinf(x) || isinf(y)) && isinf(z) && + ((signbit(x) != 0) != (signbit(y) != 0)) != (signbit(z) != 0)) { feraiseexcept(FE_INVALID); } diff --git a/src/solvers/floatbv/float_bv.cpp b/src/solvers/floatbv/float_bv.cpp index ee657069fe6..d7f053801a6 100644 --- a/src/solvers/floatbv/float_bv.cpp +++ b/src/solvers/floatbv/float_bv.cpp @@ -932,7 +932,6 @@ exprt float_bvt::fma( // Sign exprt add_sub_sign = notequal_exprt( if_exprt(c_bigger, unpacked_add.sign, prod_sign), fraction_sign); - result.sign = add_sub_sign; // NaN exprt prod_inf = or_exprt(unpacked_lhs.infinity, unpacked_rhs.infinity); @@ -950,6 +949,38 @@ exprt float_bvt::fma( result.infinity = and_exprt(not_exprt(result.NaN), or_exprt(prod_inf, unpacked_add.infinity)); + // Zero? + // Note that: + // 1. The zero flag isn't used apart from in divide and is only set on + // unpack. + // 2. The wide-precision exact product means a non-zero true result has at + // least one bit set in result.fraction; subnormals can't round to zero, + // so this test is sound here. + // 3. The rules for sign are different for zero. + result.zero = and_exprt{ + not_exprt{or_exprt{result.infinity, result.NaN}}, + equal_exprt{result.fraction, from_integer(0, result.fraction.type())}}; + + // Sign. Keep in sync with the analogous block in float_bvt::add_sub. + // For an infinity result, use the product sign if the product itself is + // infinite, otherwise the addend sign. For a zero result, follow the + // signed-zero conventions: round-to-minus-inf yields a negative zero unless + // both contributing signs are positive; all other modes yield a positive + // zero unless both contributing signs are negative. Otherwise use the + // standard add/sub sign computed above. + exprt infinity_sign = if_exprt{prod_inf, prod_sign, unpacked_add.sign}; + + const rounding_mode_bitst rounding_mode_bits{rm}; + exprt zero_sign = if_exprt{ + rounding_mode_bits.round_to_minus_inf, + or_exprt{prod_sign, unpacked_add.sign}, + and_exprt{prod_sign, unpacked_add.sign}}; + + result.sign = if_exprt{ + result.infinity, + std::move(infinity_sign), + if_exprt{result.zero, std::move(zero_sign), add_sub_sign}}; + return rounder(result, rm, spec); }