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.
+//
+// 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 <assert.h>
+#include <math.h>
+
+#ifndef _MSC_VER
+#  include <fenv.h>
+#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;
+}

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.

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).

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 <math.h> 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);
   }

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);
 }