Fix FMA sign handling in float_bvt for infinity and zero results

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.
+//
+// The infinity-result tests use __CPROVER_fmaf directly to bypass the
+// library wrapper's exception-raising assertions for inf inputs. The
+// zero-result tests use ordinary fmaf with finite inputs, which the
+// wrapper passes through unchanged.
+
+#include <assert.h>
+#include <math.h>
+
+#ifndef _MSC_VER
+#  include <fenv.h>
+#endif
+
+extern int __CPROVER_rounding_mode;
+float __CPROVER_fmaf(float, float, float);
+
+// 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 = __CPROVER_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 = __CPROVER_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 = __CPROVER_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 = __CPROVER_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;
+}

regression/cbmc/Float-fma-sign/test.desc

@@ -0,0 +1,10 @@
+CORE broken-cprover-smt-backend no-new-smt
+main.c
+--floatbv
+^EXIT=0$
+^SIGNAL=0$
+^VERIFICATION SUCCESSFUL$
+--
+^warning: ignoring
+--
+FMA sign handling for infinity and zero results.

regression/cbmc/Float-fma-sign/test_smt.desc

@@ -0,0 +1,10 @@
+CORE broken-cprover-smt-backend 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).

src/ansi-c/library/math.c

@@ -3518,12 +3518,12 @@ double fma(double x, double y, double z)
 float fmaf(float x, float y, float z)
 {
   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
+  else if(isinf(z)) // this is an over-approximation
   {
     feraiseexcept(FE_INVALID);
   }
@@ -3546,12 +3546,12 @@ float fmaf(float x, float y, float z)
 long double fmal(long double x, long double y, long double z)
 {
   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
+  else if(isinf(z)) // this is an over-approximation
   {
     feraiseexcept(FE_INVALID);
   }

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,
+    infinity_sign,
+    if_exprt(result.zero, zero_sign, add_sub_sign));
+
   return rounder(result, rm, spec);
 }