Skip to content

Commit 6df204e

Browse files
authored
Merge pull request #8884 from tautschnig/fma
Add fused multiply-add (FMA) support with single rounding
2 parents 54feda7 + 6756b84 commit 6df204e

39 files changed

Lines changed: 768 additions & 90 deletions
Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
#include <assert.h>
2+
#include <math.h>
3+
4+
int main()
5+
{
6+
// Basic FMA
7+
assert(fmaf(2.0f, 3.0f, 4.0f) == 10.0f);
8+
assert(fmaf(-3.0f, 2.0f, 10.0f) == 4.0f);
9+
10+
// FMA vs mul+add: FMA preserves precision lost by separate mul+add.
11+
// 0x1.fffffep+23 * 0x1.000002p+0 + 1.0:
12+
// mul rounds product to 0x1p+24, then +1 = 0x1p+24 (unchanged)
13+
// FMA keeps exact product, +1 = 0x1.000002p+24
14+
assert(0x1.fffffep+23f * 0x1.000002p+0f + 1.0f == 0x1p+24f);
15+
assert(fmaf(0x1.fffffep+23f, 0x1.000002p+0f, 1.0f) == 0x1.000002p+24f);
16+
17+
// FMA for remainder: fma(-n, y, x) computes x - n*y with single rounding
18+
float x = 0x1.d55556p+0f;
19+
float y = 0x1.555556p-2f;
20+
assert(x - 5.0f * y == 0x1.55555p-3f); // mul+sub: two roundings
21+
assert(fmaf(-5.0f, y, x) == 0x1.555554p-3f); // FMA: one rounding
22+
}
Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
// IEEE 754 FMA special cases
2+
#include <assert.h>
3+
#include <math.h>
4+
5+
int main()
6+
{
7+
// fma(NaN, y, z) = NaN
8+
assert(isnan(fmaf(NAN, 1.0f, 0.0f)));
9+
10+
// fma(x, NaN, z) = NaN
11+
assert(isnan(fmaf(1.0f, NAN, 0.0f)));
12+
13+
// fma(x, y, NaN) = NaN
14+
assert(isnan(fmaf(1.0f, 1.0f, NAN)));
15+
16+
// fma(0, inf, z) = NaN (0 * inf is undefined)
17+
assert(isnan(fmaf(0.0f, INFINITY, 1.0f)));
18+
19+
// fma(inf, 0, z) = NaN
20+
assert(isnan(fmaf(INFINITY, 0.0f, 1.0f)));
21+
22+
// fma(inf, x, -inf) = NaN when inf*x = +inf (inf + (-inf) is undefined)
23+
assert(isnan(fmaf(INFINITY, 1.0f, -INFINITY)));
24+
25+
// fma(inf, x, z) = +inf for finite z
26+
assert(isinf(fmaf(INFINITY, 2.0f, 3.0f)));
27+
28+
// fma(x, y, inf) = +inf for finite x*y
29+
assert(isinf(fmaf(2.0f, 3.0f, INFINITY)));
30+
31+
return 0;
32+
}
Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
CORE smt-backend
2+
precision.c
3+
--smt2 --outfile -
4+
float_bv\.floatbv_fma
5+
^EXIT=0$
6+
^SIGNAL=0$
7+
--
8+
^warning: ignoring
9+
--
10+
Verify that float_bvt generates valid SMT2 for FMA expressions
11+
(non-FPA bitvector encoding).
Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
CORE smt-backend
2+
precision.c
3+
--smt2 --fpa --outfile -
4+
fp\.fma
5+
^EXIT=0$
6+
^SIGNAL=0$
7+
--
8+
^warning: ignoring
9+
--
10+
Verify that the SMT2 FPA backend emits fp.fma (not a bitvector encoding).
Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
CORE smt-backend
2+
special_cases.c
3+
--smt2 --fpa --no-built-in-assertions
4+
^EXIT=0$
5+
^SIGNAL=0$
6+
^VERIFICATION SUCCESSFUL$
7+
--
8+
^warning: ignoring
9+
--
10+
Run FMA special cases through an external SMT solver with FPA theory.
11+
Verifies NaN, infinity, and 0*inf handling end-to-end.
Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
CORE smt-backend
2+
special_cases.c
3+
--smt2 --fpa --outfile -
4+
fp\.fma
5+
^EXIT=0$
6+
^SIGNAL=0$
7+
--
8+
^warning: ignoring
9+
--
10+
Verify that FMA special cases (NaN, inf) emit fp.fma via FPA backend.
Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
CORE
2+
precision.c
3+
4+
^EXIT=0$
5+
^SIGNAL=0$
6+
^VERIFICATION SUCCESSFUL$
7+
--
8+
^warning: ignoring
Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
CORE
2+
special_cases.c
3+
--floatbv --no-built-in-assertions
4+
^EXIT=0$
5+
^SIGNAL=0$
6+
^VERIFICATION SUCCESSFUL$
7+
--
8+
^warning: ignoring
9+
--
10+
IEEE 754 FMA special cases: NaN propagation, 0*inf, inf+(-inf), inf results.
Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
// fmaf should compute x*y+z with a single rounding.
2+
// The C library model used to do x*y then +z (two roundings).
3+
// Example: fmaf(1+eps, 1+eps, -(1+2*eps)) where eps = 2^-23
4+
// Exact result: eps^2 = 2^-46 > 0
5+
// Double rounding: (1+eps)*(1+eps) rounds to 1+2*eps, then +c = 0
6+
// True FMA: should give eps^2 which is > 0
7+
#include <assert.h>
8+
#include <math.h>
9+
10+
int main()
11+
{
12+
float a = 1.0f + 1.1920929e-7f; // 1 + 2^-23
13+
float c = -(1.0f + 2.3841858e-7f); // -(1 + 2^-22)
14+
float r = fmaf(a, a, c);
15+
assert(r > 0.0f);
16+
return 0;
17+
}
Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
CORE no-new-smt
2+
main.c
3+
--floatbv --no-built-in-assertions
4+
^EXIT=0$
5+
^SIGNAL=0$
6+
^VERIFICATION SUCCESSFUL$
7+
--
8+
fmaf with single rounding via __CPROVER_fmaf built-in.

0 commit comments

Comments
 (0)