riscv_fpu: static/FMA rounding modes, exact integral rounding, FMA underflow tininess

[SolAstrius]

Summary

Reworked resubmission of the remaining #239 commits, addressing @purplesyringa's review. Four commits, each standalone on staging; no dependency on #242 (RMM), and the two deliberately overlap nowhere except the shared riscv_fpu_host_rm helper.

1. fpu_lib: round to integral values on the encoding, not via biasing adds

Replaces the NACK'd INEXACT patch. The review showed the old f + 0.5 trick was broken beyond the flag leak (this comment): the biasing add itself rounds, so rne(prev(0.5)) returned 1 and rne(2.5) returned 3 (ties-away, not ties-to-even). fpu_round_fXX_internal now decides on the fraction bits of the encoding: no host FP arithmetic, so no spurious flags, no double rounding, and no dependence on a host rounding mode (relevant for the future USE_SOFT_FENV story). RMM is a first-class case — only DYN falls back to the tracked mode, so fcvt/fround with a static rmm field round correctly.

Verified bit-exact against host rint()/round() in all four modes over the full fractional exponent range of both widths (2.0e9 cases, zero mismatches).

2. riscv_fpu: honor the static rounding-mode field

Reworked per review: the mode change is now gated per-instruction (riscv_f_op_is_rounding) — funct3 is an rm field only on rounding-capable ops, so fsgnj/fcmp/fclass/fmv never trigger it. Ops that consume rm as an argument (fcvt-to-int, fround) need no wrap. A static rmm field runs in RNE for now (differs only on exact ties; wired to the exact fixups after #242 lands). The dynamic-RMM preparation is gated to rm == DYN so a static field is no longer clobbered by it.

3. riscv_fpu: honor the rounding mode in the FMA ops

The fnmadd operand-negation part already landed via #241; this is the rm plumbing, now sharing riscv_fpu_host_rm with the OP-FP dispatch instead of duplicating the effective-mode logic (review ask).

4. fpu_lib: set the FMA underflow flag by IEEE after-rounding tininess

Reworked per review, and moved into fpu_fma32/64 proper; fpu_fma32's own rounding is unchanged (see the review thread for why the direct downcast mis-rounds):

• the NX fast-path gate is gone (it conflated "set by this op" with "ever set");
• the verdict merges over a flag snapshot from before the op, so a sticky UF from earlier ops survives;
• the fixup runs while the effective mode is still set (it lives inside the op now);
• f64 scales by 2^52 instead of 2^200, with the operand-bound proof from the review included as a comment; the b == 0 overflow-to-NaN case reads correctly as not-tiny;
• f32 decides tininess by comparing the widened 53-bit sum against the mode's threshold in f64, with the exact TwoSum error breaking the tie when the sum lands exactly on the to-nearest threshold.

On aarch64 (before-rounding tininess), fma(2^-62*(1-2^-23), 2^-64*(1+2^-23), 0) returned the correct 2^-126 but flagged NX|UF instead of NX; targeted tests cover both boundary directions, an exact boundary, a genuinely-tiny-at-boundary case, and sticky-UF preservation — all pass with the fix, three fail without it.

Conformance (riscv-arch-test ACT4, Spike reference)

build	F	D	I	M	total
staging	13/82	24/114	51/51	13/13	101/260
this PR	21/82	51/114	51/51	13/13	136/260

No regressions. The remaining F/D failures are, by sampled diagnosis, exactly the two intentionally-excluded categories:

• RMM vectors — the old riscv_prepare_rmm directed hack (e.g. it reads a NaN-boxed f32 through the f64 default: path and rounds fsqrt the wrong way). Deleted and replaced by riscv_fpu: RMM (roundTiesToAway) for the scalar OP-FP ops #242.
• NaN result canonicalization — arithmetic results emit host NaN payloads where the spec wants the canonical NaN (fadd.s → 0x7fd825fe vs 0x7fc00000). That's the separate flag-gated patch we agreed on; needs a decision on compile-time vs runtime flag first.

A local combined build (this PR + #242 + the canonicalization patch) reaches 244/260, and every one of its 16 residual failures is an FMA test failing on an frm=rmm tie sub-case — that's the one #239 commit not yet resubmitted anywhere ("round the FMA family to nearest, ties-away under RMM"); it needs #242's fixup machinery, so it goes on top of #242 once that merges, closing the gap to the 260/260 the old chain reported.

cc @purplesyringa — draft until reviewed, to keep GitHub from merging things on its own again.

[purplesyringa]

I noticed that the first commit (the one about float-to-integer rounding) kept

    if (unlikely(fpu_is_zero64_soft(d))) {
        return 0;
    }

in fpu_round_f64_to_u32 (and similar methods). I'm not sure what role that if used to play, but judging from unlikely, maybe it's not needed anymore? @LekKit do you remember why that was needed?

[purplesyringa]

Just throwing ideas out there: I think fp-to-int roundings can be optimized further. I don't expect you to do this here if you don't want to, I'm mostly just glad it's sound, but you can try this if you have time.

The old FP-based code tried to force correct roundings with formulas like f + 0.5, which have edge cases, because you're trying to fit complex rules in one formula. However, you can force the FPU to do the correct rounding with a different method. Let's say 0 <= f < 2^52 - 1 is a floating-point value. Then f + 2^52 will have an exponent of exactly 52, so the mantissa field is exactly round(f). Since addition uses the rounding mode from fenv, you can effectively control which rounding mode this method implements. See also.

This method introduces a few subtleties. For instance, there are some edge cases around f ~= 2^52, e.g. for f = 2^52 - 0.5, f + 2^52 can round up to 2^53, which has an empty mantissa field, but you can fix this by computing bitcast(f + 2^52) - bitcast(2^52) instead of masking bits. Point is, you'll still need some special casing for e >= 52, but you don't need a switch on the runtime type, you can just use the host mode for everything and avoid manually implementing FP logic.

Of course, RMM is the elephant in the room. That will still require some sort of fixup.

[purplesyringa]

My original review mentioned that computing sum is already exact in the 32-bit case, so handling err/res/soft_fenv_... is unnecessary and you can directly downcast sum to f32 and obtain the correct result. The PR mentions that you fixed fpu_fma32, but I don't see that? Did you run into some problems while trying that?

[SolAstrius]

I tried exactly that and it mis-rounds — sum isn't exact (the product is, but product + c can span more bits than f64 holds), and downcasting the 53-bit rounding erases 24-bit tie-breaks:

// cc -O2 tie_probe.c -lm && ./a.out
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
volatile float a, b, c;
int main(void)
{
    uint32_t ua = 0x6BE58000, ub = 0x3A123B00, uc = 0x117D0E3E, r1, r2;
    memcpy((void*)&a, &ua, 4); memcpy((void*)&b, &ub, 4); memcpy((void*)&c, &uc, 4);
    float f1 = fmaf(a, b, c), f2 = (float)((double)a * b + c);
    memcpy(&r1, &f1, 4); memcpy(&r2, &f2, 4);
    printf("fmaf     %08x\ndowncast %08x\n", r1, r2);
}

# arm64 and x86_64, identical:
fmaf     668317e5
downcast 668317e4

a*b sits exactly on a 24-bit tie and c breaks it upward; the intermediate rounding loses that. (2p+2 doesn't apply — the addends aren't 24-bit values.) So fpu_fma32's own rounding keeps the odd-round; the PR text claiming I "fixed" it was wrong, rewording.

[purplesyringa]

Oh, now I see, (float)(a * b + c) has double rounding, one after + and one due to (float), thanks. I'll re-review the code taking this into account in a bit.

[purplesyringa]

Hopefully you can explain this to me.

If I understand it correctly, RISC-V says UF corresponds to the value after rounding for all operations. In particular, I'm assuming this includes fcvt.s.d, which is implemented via fpu_fcvt_f32_to_f64.

Since f32 denormals easily fit in f64, the only function in fpu_fma32 that can set UF is fpu_fcvt_f32_to_f64. So if fcvt.s.d sets UF correctly, then fpu_fcvt_f32_to_f64 should be sufficient here would FMA-specific logic, and this function could be much, much simpler (I think you wouldn't need to handle exceptions at all, actually, but I'm not certain about that). This would certainly remove fpu_fma32_fixup_uf at least.

And if fpu_fcvt_f32_to_f64 is broken, then that's where the fix should live, not here, I think. This also raises the topic of whether any other arithmetic operations need similar adjustments -- e.g. what specifically makes FMA so different from multiplication that only FMA needs special handling?

Could you explain where my logic fails?

[SolAstrius]

The premise fails: fcvt.s.d doesn't set UF correctly on aarch64 either. (You mean the narrowing fpu_fcvt_f64_to_f32 — the widening one is exact.)

// cc -O2 uf_probe.c -lm && ./a.out
// Every case rounds up from below onto exactly +min-normal:
// after-rounding tininess says NX only.
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <fenv.h>
volatile float fa, fb, fr;
volatile double da, db, dr;
static void flags(const char* op, uint64_t bits)
{
    int e = fetestexcept(FE_ALL_EXCEPT);
    printf("%-8s %016llx  %s%s\n", op, (unsigned long long)bits,
           (e & FE_INEXACT) ? "NX " : "", (e & FE_UNDERFLOW) ? "UF" : "");
    feclearexcept(FE_ALL_EXCEPT);
}
int main(void)
{
    uint32_t a32 = 0x207FFFFE, b32 = 0x1F800001, r32; // 2^-62*(1-2^-23), 2^-64*(1+2^-23)
    uint64_t a64 = 0x200FFFFFFFFFFFFE, b64 = 0x1FF0000000000001, r64;
    memcpy((void*)&fa, &a32, 4); memcpy((void*)&fb, &b32, 4);
    memcpy((void*)&da, &a64, 8); memcpy((void*)&db, &b64, 8);
    feclearexcept(FE_ALL_EXCEPT);
    fr = fa * fb;         memcpy(&r32, (void*)&fr, 4); flags("mul.s", r32);
    fr = fmaf(fa, fb, 0); memcpy(&r32, (void*)&fr, 4); flags("fma.s", r32);
    dr = da * db;         memcpy(&r64, (void*)&dr, 8); flags("mul.d", r64);
    dr = fma(da, db, 0);  memcpy(&r64, (void*)&dr, 8); flags("fma.d", r64);
    dr = ldexp(1.0 - ldexp(1.0, -46), -126); feclearexcept(FE_ALL_EXCEPT);
    fr = (float)dr;       memcpy(&r32, (void*)&fr, 4); flags("cvt.s.d", r32);
}

# x86_64:                        # arm64 (M1):
mul.s    0000000000800000  NX    mul.s    0000000000800000  NX UF
fma.s    0000000000800000  NX    fma.s    0000000000800000  NX UF
mul.d    0010000000000000  NX    mul.d    0010000000000000  NX UF
fma.d    0010000000000000  NX    fma.d    0010000000000000  NX UF
cvt.s.d  0000000000800000  NX    cvt.s.d  0000000000800000  NX UF

ARM detects tininess before rounding on every instruction, x86 after. So nothing makes FMA special — mul/div/cvt all have the same bug; only add/sub are safe, since a sum near the boundary is always exact (shared quantum when both operands are small, Sterbenz's lemma on cancellation).

The right end state is fixing this once for all of them, agreed — but doing it correctly needs a flags snapshot from before the op, or a not-tiny verdict can't clear the bogus UF without also clearing a sticky one. fma already pays for that snapshot at entry; the hot ops don't. So: FMA here, cvt/mul/div/sqrt together as a follow-up.

[purplesyringa]

I think this about concludes my review for now. The logic I didn't comment on looks good.

[LekKit]

I noticed that the first commit (the one about float-to-integer rounding) kept

    if (unlikely(fpu_is_zero64_soft(d))) {
        return 0;
    }

in fpu_round_f64_to_u32 (and similar methods). I'm not sure what role that if used to play, but judging from unlikely, maybe it's not needed anymore? @LekKit do you remember why that was needed?

AFAIK the original rounding routine was written in a way that did not account for literal zero (which is already rounded anyways). Now this seems redundant, but i would instead verify if it actually produces expected result.

[purplesyringa]

LGTM on everything except the fpu_fma32 UF handling, I still don't understand why we have some complicated logic here.

Say, consider the example in the commit message. We have an exact result of 2^-126 - 2^-151 = 2^-127 * (1 + (1 - 2^-24)). Round-to-nearest with an unbounded exponent produces 2^-127 * (1 + (1 - 2^23)), which is below the smallest denormal, and bounded rounding produces 2^-126, which is not tiny, so yes, setting tiny from the value already rounded to f32 is incorrect, as the commit description correctly mentions.

But we never used sum rounded to f32 for calculating tiny, we first multiplied it by 2^126, so there was never a situation where the exponent in this rounding would get bounded. And even if it did, we can just compare sum to 2^-126 while staying in f64 -- in fact, that should be faster and cleaner.

So I don't see why you're piling on more logic instead of removing it. See also: the two unresolved review comments on this topic, in particular I'm wondering if other operations, like addition or f64-to-f32 downcasting, should also get special handling for UF, and if so, if 32-bit FMA UF handling would look simpler if it was based on that logic.

[SolAstrius]

In b3a3f3a I have introduced a regression; fixing it.

== b17bfc5 (fpu_lib: set the FMA underflow flag by IEEE a)
  F-fcvt.l.s-00 PASS
  F-fcvt.s.w-00 PASS
  F-fcvt.wu.s-00 PASS
== 1587da8 (fpu_lib: clarify the integral-rounding edge c)
  F-fcvt.l.s-00 PASS
  F-fcvt.s.w-00 PASS
  F-fcvt.wu.s-00 PASS
== b3a3f3a (riscv_fpu: share static-rm handling between t)
  F-fcvt.l.s-00 FAIL
  F-fcvt.s.w-00 FAIL
  F-fcvt.wu.s-00 FAIL
== bd053f0 (fpu_lib: decide FMA32 tininess by an f64 thre)
  F-fcvt.l.s-00 FAIL
  F-fcvt.s.w-00 FAIL
  F-fcvt.wu.s-00 FAIL

[SolAstrius]

Switched to the f64 compare in bd053f0; the per-mode thresholds are in the comment, the rest is in the two threads. One correction to the walk-through: 2^-126 - 2^-151 in the commit message is the midpoint, not the exact result (ambiguous wording, mine). The failing input sits just below it — e.g. 2^-126 - 2^-151 - 2^-197, whose sum collapses onto the midpoint exactly — which is why comparing sum alone needs the TwoSum tie-break there, and only there. Done committing on this one.