Soft-path FMA emulation: incorrect results under directed rounding modes (and two RNE-mode edge cases)

[SolAstrius]

Following up on a question raised in #246: does the odd-rounding soft FMA path in fpu_lib.h stay correct under non-RNE rounding modes?

I fuzzed the soft paths of fpu_fma32 and fpu_fma64_raw (builtin FMA compiled out) against two independent references — AArch64 hardware FMA and glibc's soft fma/fmaf on x86-64, which agree with each other on every counterexample. 5M exponent-biased random triples per mode:

Mode	fma32 finite mismatches	fma64 finite mismatches	fma64 inf/NaN mismatches
RNE	0	6	312,444
RUP	188,408	11,099	156,646
RDN	380,579	11,139	302,509
RTZ	380,255	132,652	0

Four distinct failure classes, one minimal counterexample each (bit patterns; soft = current soft path, ref = correctly rounded):

1. fma32, directed modes — TwoSum error term invalid outside RNE, odd-round nudges the wrong way:
   RUP a=2f2cfc98 b=b45a0d25 c=c33f5ab8 → soft=c33f5ab7 ref=c33f5ab8
2. fma64, directed modes — same, plus the Dekker product splitting also assumes RNE:
   RUP a=0e111f2c376c3b7c b=b09ab33b286b11f2 c=d1ae852da52f2f73 → soft=…2f72 ref=…2f73
3. fma64, RNE, subnormal-range result — the "u_l is representable" hypothesis of the odd-rounding paper is violated when the product's error term underflows:
   RNE a=1e492e4a5e074268 b=20636e8053091041 c=80043de497288d74 → soft=80043de41cd50b65 ref=…0b64
4. fma64, product overflow, all modes — mul = inf makes fpu_mul_error64 compute inf − inf = NaN, which propagates to the result instead of ±inf:
   RNE a=69cdf850bf334cfd b=7ca8649ebebeef4e c=988d1e217a30a612 → soft=NaN ref=+inf

Classes 3 and 4 affect RNE too, i.e. any current build that takes the soft path (e.g. baseline x86-64 without -mfma). fma32 is immune to overflow (the f64 product of two f32s can't overflow) and appears fully correct in RNE.

Reproducers below. fma_rm_test.c is standalone (replicates the soft paths verbatim); real_check.c compiles against the actual src/util/fpu_lib.h (-I src/util) and confirms each counterexample above on the real code. Verified on macOS/AArch64 (clang, -O1 -ffp-model=strict) and Linux/x86-64 (gcc 15, -O1 -frounding-math).

fma_rm_test.c — standalone fuzzer

// Brute-force RVVM soft-FMA emulation vs hardware FMA under non-RNE rounding.
// Replicates fpu_fma32 / fpu_fma64_raw soft paths from src/util/fpu_lib.h.
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <fenv.h>
#include <stdlib.h>

#pragma STDC FENV_ACCESS ON

static inline uint32_t f2u(float f)  { uint32_t u; memcpy(&u, &f, 4); return u; }
static inline float    u2f(uint32_t u) { float f; memcpy(&f, &u, 4); return f; }
static inline uint64_t d2u(double d) { uint64_t u; memcpy(&u, &d, 8); return u; }
static inline double   u2d(uint64_t u) { double d; memcpy(&d, &u, 8); return d; }

static double add_error64(double fs, double fa, double fb)
{
    volatile double fi = fs - fa;
    volatile double t1 = fb - fi;
    volatile double t2 = fs - fi;
    volatile double t3 = fa - t2;
    return t1 + t3;
}

static double mul_error64(double fm, double a, double b)
{
    volatile double ha = u2d(d2u(a) & ~0x7FFFFFFULL);
    volatile double la = a - ha;
    volatile double hb = u2d(d2u(b) & ~0x7FFFFFFULL);
    volatile double lb = b - hb;
    volatile double p1 = ha * hb - fm;
    volatile double p2 = p1 + ha * lb;
    volatile double p3 = p2 + la * hb;
    return p3 + la * lb;
}

static double odd_round64(double val, double err)
{
    uint64_t uv = d2u(val), ue = d2u(err);
    if (ue && !(uv & 1)) {
        return u2d((ue >> 63) ? uv - 1 : uv + 1);
    }
    return val;
}

// fpu_fma32 soft path
static float soft_fma32(float a, float b, float c)
{
    volatile double mul = (double)a * (double)b; // exact
    volatile double add = (double)c;
    volatile double sum = mul + add;
    double err = add_error64(sum, mul, add);
    double res = odd_round64(sum, err);
    return (float)res;
}

// fpu_fma64_raw soft path
static double soft_fma64(double a, double b, double c)
{
    volatile double mul = a * b;
    double e_m = mul_error64(mul, a, b);
    volatile double sum = mul + c;
    double e_s = add_error64(sum, mul, c);
    volatile double e_f = e_m + e_s;
    double err = odd_round64(e_f, add_error64(e_f, e_s, e_m));
    volatile double ret = sum + err;
    return ret;
}

static uint64_t rng_state = 0x243F6A8885A308D3ULL;
static uint64_t rng(void)
{
    rng_state ^= rng_state << 13; rng_state ^= rng_state >> 7; rng_state ^= rng_state << 17;
    return rng_state;
}

// random float with biased exponent to provoke cancellation/overlap
static float rand_f32(void)
{
    uint32_t u = (uint32_t)rng();
    if (rng() & 1) { // squeeze exponent into a narrow band around a random center
        int32_t center = 64 + (int32_t)(rng() % 128);
        int32_t e = center + (int32_t)(rng() % 48) - 24;
        if (e < 0) e = 0; if (e > 254) e = 254;
        u = (u & 0x807FFFFFU) | ((uint32_t)e << 23);
    }
    return u2f(u);
}

static double rand_f64(void)
{
    uint64_t u = rng();
    if (rng() & 1) {
        int64_t center = 512 + (int64_t)(rng() % 1024);
        int64_t e = center + (int64_t)(rng() % 120) - 60;
        if (e < 0) e = 0; if (e > 2046) e = 2046;
        u = (u & 0x800FFFFFFFFFFFFFULL) | ((uint64_t)e << 52);
    }
    return u2d(u);
}

static const int modes[] = { FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO };
static const char* mode_names[] = { "RNE", "RUP", "RDN", "RTZ" };

int main(int argc, char** argv)
{
    long iters = argc > 1 ? atol(argv[1]) : 10000000;
    long fail32 = 0, fail64 = 0;
    int shown32 = 0, shown64 = 0;

    for (int m = 0; m < 4; m++) {
        fesetround(modes[m]);
        long f32c = 0, f64c = 0, f32inf = 0, f64inf = 0;
        for (long i = 0; i < iters; i++) {
            float a = rand_f32(), b = rand_f32(), c = rand_f32();
            if (!isfinite(a) || !isfinite(b) || !isfinite(c)) continue;
            float got = soft_fma32(a, b, c);
            float ref = fmaf(a, b, c);
            if (f2u(got) != f2u(ref) && !(isnan(got) && isnan(ref))) {
                if (!isfinite(got) || !isfinite(ref)) { f32inf++; goto skip32; }
                f32c++;
                if (shown32 < 3) {
                    printf("[%s] fma32 MISMATCH a=%08x b=%08x c=%08x got=%08x ref=%08x  (a=%g b=%g c=%g)\n",
                           mode_names[m], f2u(a), f2u(b), f2u(c), f2u(got), f2u(ref), a, b, c);
                    shown32++;
                }
            }
            skip32:;
            double da = rand_f64(), db = rand_f64(), dc = rand_f64();
            if (!isfinite(da) || !isfinite(db) || !isfinite(dc)) continue;
            double dgot = soft_fma64(da, db, dc);
            double dref = fma(da, db, dc);
            if (d2u(dgot) != d2u(dref) && !(isnan(dgot) && isnan(dref))) {
                if (!isfinite(dgot) || !isfinite(dref)) { f64inf++; continue; }
                f64c++;
                if (shown64 < 3) {
                    printf("[%s] fma64 MISMATCH a=%016llx b=%016llx c=%016llx got=%016llx ref=%016llx\n",
                           mode_names[m], (unsigned long long)d2u(da), (unsigned long long)d2u(db),
                           (unsigned long long)d2u(dc), (unsigned long long)d2u(dgot), (unsigned long long)d2u(dref));
                    shown64++;
                }
            }
        }
        printf("[%s] %ld iters: fma32 finite-fails=%ld (inf/nan-fails=%ld) fma64 finite-fails=%ld (inf/nan-fails=%ld)\n", mode_names[m], iters, f32c, f32inf, f64c, f64inf);
        fail32 += f32c; fail64 += f64c;
        shown32 = shown64 = 0;
    }
    fesetround(FE_TONEAREST);
    printf("TOTAL: fma32=%ld fma64=%ld\n", fail32, fail64);
    return (fail32 || fail64) ? 1 : 0;
}

real_check.c — counterexample checker against the real fpu_lib.h

// Verify mismatches against the real RVVM fpu_lib.h soft path (builtin FMA forced off)
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <fenv.h>
#define FPU_LIB_DISABLE_BUILTIN_FMA 1
#include "fpu_lib.h"
#undef FPU_LIB_OPTIMAL_BUILTIN_FMA
// re-include won't work; instead call the soft sub-path directly:
static fpu_f32_t soft32(fpu_f32_t a, fpu_f32_t b, fpu_f32_t c)
{
    fpu_f64_t mul = fpu_mul64(fpu_fcvt_f32_to_f64(a), fpu_fcvt_f32_to_f64(b));
    fpu_f64_t add = fpu_fcvt_f32_to_f64(c);
    fpu_f64_t sum = fpu_add64(mul, add);
    fpu_f64_t err = fpu_add_error64(sum, mul, add);
    fpu_f64_t res = fpu_odd_round64(sum, err);
    return fpu_fcvt_f64_to_f32(res);
}
static fpu_f64_t soft64(fpu_f64_t a, fpu_f64_t b, fpu_f64_t c)
{
    fpu_f64_t mul = fpu_mul64(a, b);
    fpu_f64_t e_m = fpu_mul_error64(mul, a, b);
    fpu_f64_t sum = fpu_add64(mul, c);
    fpu_f64_t e_s = fpu_add_error64(sum, mul, c);
    fpu_f64_t e_f = fpu_add64(e_m, e_s);
    fpu_f64_t err = fpu_odd_round64(e_f, fpu_add_error64(e_f, e_s, e_m));
    return fpu_add64(sum, err);
}
static void case32(int mode, const char* mn, uint32_t ua, uint32_t ub, uint32_t uc)
{
    fesetround(mode);
    fpu_f32_t r = soft32(fpu_bit_u32_to_f32(ua), fpu_bit_u32_to_f32(ub), fpu_bit_u32_to_f32(uc));
    float a, b, c; memcpy(&a,&ua,4); memcpy(&b,&ub,4); memcpy(&c,&uc,4);
    float ref = fmaf(a, b, c);
    uint32_t ur; memcpy(&ur, &ref, 4);
    printf("[%s] fma32 a=%08x b=%08x c=%08x rvvm=%08x hw=%08x %s\n", mn, ua, ub, uc,
           fpu_bit_f32_to_u32(r), ur, fpu_bit_f32_to_u32(r)==ur ? "OK" : "MISMATCH");
}
static void case64(int mode, const char* mn, uint64_t ua, uint64_t ub, uint64_t uc)
{
    fesetround(mode);
    fpu_f64_t r = soft64(fpu_bit_u64_to_f64(ua), fpu_bit_u64_to_f64(ub), fpu_bit_u64_to_f64(uc));
    double a, b, c; memcpy(&a,&ua,8); memcpy(&b,&ub,8); memcpy(&c,&uc,8);
    double ref = fma(a, b, c);
    uint64_t ur; memcpy(&ur, &ref, 8);
    printf("[%s] fma64 a=%016llx b=%016llx c=%016llx rvvm=%016llx hw=%016llx %s\n", mn, ua, ub, uc,
           (unsigned long long)fpu_bit_f64_to_u64(r), (unsigned long long)ur,
           fpu_bit_f64_to_u64(r)==ur ? "OK" : "MISMATCH");
}
int main(void)
{
    case32(FE_UPWARD,   "RUP", 0x2f2cfc98, 0xb45a0d25, 0xc33f5ab8);
    case32(FE_DOWNWARD, "RDN", 0x8a552b83, 0xa6ab6216, 0xa3f9d530);
    case32(FE_TOWARDZERO,"RTZ", 0x249a113c, 0x8409344a, 0x877b1ecb);
    case64(FE_UPWARD,   "RUP", 0x0e111f2c376c3b7cULL, 0xb09ab33b286b11f2ULL, 0xd1ae852da52f2f73ULL);
    case64(FE_TOWARDZERO,"RTZ", 0x21855c5cb02ac126ULL, 0x89691a36d54af81eULL, 0x28c853cd1a11f2baULL);
    // RNE finite fail (subnormal region)
    case64(FE_TONEAREST,"RNE", 0x1e492e4a5e074268ULL, 0x20636e8053091041ULL, 0x80043de497288d74ULL);
    // RNE overflow: NaN instead of inf
    case64(FE_TONEAREST,"RNE", 0x69cdf850bf334cfdULL, 0x7ca8649ebebeef4eULL, 0x988d1e217a30a612ULL);
    return 0;
}

[purplesyringa]

Related: fpu_odd_round64 performs all logic exclusively on the bitwise representation of the value, meaning that it may return nonsensical results on NaN, +-Infinity, and around 0. Specifically I'm worried that uv - 1 can turn 0 into NaN, Infinity into a real value, and change the mantissa of the NaN to be non-canonical, and similarly uv + 1 can do the opposite.