From 9c7f6e425797829499111d3a8b456e2378446193 Mon Sep 17 00:00:00 2001 From: Sol Astrius Phoenix Date: Sat, 18 Jul 2026 00:08:30 +0300 Subject: [PATCH 01/10] fpu_lib: round to integral values on the encoding, not via biasing adds fpu_round_fXX_internal biased the value by +/-0.5 (or +/-1.0) with a real FP add and let the caller's int cast truncate. The add is itself a rounded op, which broke it two ways: it raised a spurious INEXACT for the discarded fraction (leaking NX alongside e.g. the NV of an out-of-range fcvt), and it mis-rounded whenever the biasing add was inexact or exact on the wrong side -- rne(prev(0.5)) returned 1 because prev(0.5) + 0.5 rounds up to 1.0, and rne(2.5) returned 3 because 2.5 + 0.5 is exactly 3.0, i.e. the trick implements ties-away, not ties-to-even. Decide on the fraction bits of the encoding instead: no host FP arithmetic, so no flag leaks, no double rounding, and no dependence on the host rounding mode (which USE_SOFT_FENV hosts don't have). RMM is now a first-class case -- only DYN falls back to the tracked mode, so an explicit rmm request is honored instead of being swapped for the dynamic mode. Verified bit-exact against host rint()/round() in all four modes over the full fractional exponent range of both widths (2.0e9 cases). Signed-off-by: Sol Astrius Phoenix --- src/util/fpu_lib.c | 99 +++++++++++++++++++++++++++++++++++----------- 1 file changed, 77 insertions(+), 22 deletions(-) diff --git a/src/util/fpu_lib.c b/src/util/fpu_lib.c index c82f0451e..15c50d8ed 100644 --- a/src/util/fpu_lib.c +++ b/src/util/fpu_lib.c @@ -405,54 +405,109 @@ slow_path uint32_t fpu_fclass64(fpu_f64_t d) return ret; } +/* + * Round to an integral value directly on the encoding: no host FP arithmetic, so + * no spurious exception flags and no double rounding. Only DYN falls back to the + * tracked mode; an explicit RMM request is honored as-is. + */ slow_path fpu_f32_t fpu_round_f32_internal(fpu_f32_t f, uint32_t mode) { - uint32_t u = fpu_bit_f32_to_u32(f); - uint32_t s = u & FPU_LIB_FP32_SIGNEDFP_MASK; - if (unlikely(mode > FPU_LIB_ROUND_UP)) { + const uint32_t u = fpu_bit_f32_to_u32(f); + const uint32_t s = u & FPU_LIB_FP32_SIGNEDFP_MASK; + const int32_t e = fpu_exponent32(f); + bool away = false; + if (unlikely(mode > FPU_LIB_ROUND_MM)) { mode = fpu_get_rounding_mode(); } + if (e >= 23 || !(u << 1)) { + return f; // Already integral: |f| >= 2^23, +/-0, inf, NaN + } + if (e < 0) { + // |f| in (0, 1) rounds to +/-0 or +/-1; the 0.5 tie goes to even == 0 + switch (mode) { + case FPU_LIB_ROUND_NE: + away = (e == -1) && (u & FPU_LIB_FP32_MANTISSA_MASK); + break; + case FPU_LIB_ROUND_MM: + away = (e == -1); + break; + case FPU_LIB_ROUND_DN: + away = !!s; + break; + case FPU_LIB_ROUND_UP: + away = !s; + break; + } + return fpu_bit_u32_to_f32(s | (away ? 0x3F800000U : 0)); + } + // |f| in [1, 2^23): split the mantissa into integer part and fraction bits + const uint32_t frac = u & (FPU_LIB_FP32_MANTISSA_MASK >> e); + const uint32_t half = 0x00400000U >> e; + const uint32_t step = 0x00800000U >> e; // 1.0 at this exponent; carries into a new binade switch (mode) { case FPU_LIB_ROUND_NE: + away = frac > half || (frac == half && (u & step)); + break; case FPU_LIB_ROUND_MM: - return fpu_add32(f, fpu_bit_u32_to_f32(0x3F000000U | s)); + away = frac >= half; + break; case FPU_LIB_ROUND_DN: - if (s && fpu_is_fractional32(f)) { - return fpu_sub32(f, fpu_bit_u32_to_f32(0x3F800000U)); - } + away = s && frac; break; case FPU_LIB_ROUND_UP: - if (!s && fpu_is_fractional32(f)) { - return fpu_add32(f, fpu_bit_u32_to_f32(0x3F800000U)); - } + away = !s && frac; break; } - return f; + return fpu_bit_u32_to_f32((u - frac) + (away ? step : 0)); } slow_path fpu_f64_t fpu_round_f64_internal(fpu_f64_t d, uint32_t mode) { - uint64_t u = fpu_bit_f64_to_u64(d); - uint64_t s = u & FPU_LIB_FP64_SIGNEDFP_MASK; - if (unlikely(mode > FPU_LIB_ROUND_UP)) { + const uint64_t u = fpu_bit_f64_to_u64(d); + const uint64_t s = u & FPU_LIB_FP64_SIGNEDFP_MASK; + const int32_t e = fpu_exponent64(d); + bool away = false; + if (unlikely(mode > FPU_LIB_ROUND_MM)) { mode = fpu_get_rounding_mode(); } + if (e >= 52 || !(u << 1)) { + return d; // Already integral: |d| >= 2^52, +/-0, inf, NaN + } + if (e < 0) { + switch (mode) { + case FPU_LIB_ROUND_NE: + away = (e == -1) && (u & FPU_LIB_FP64_MANTISSA_MASK); + break; + case FPU_LIB_ROUND_MM: + away = (e == -1); + break; + case FPU_LIB_ROUND_DN: + away = !!s; + break; + case FPU_LIB_ROUND_UP: + away = !s; + break; + } + return fpu_bit_u64_to_f64(s | (away ? 0x3FF0000000000000ULL : 0)); + } + const uint64_t frac = u & (FPU_LIB_FP64_MANTISSA_MASK >> e); + const uint64_t half = 0x0008000000000000ULL >> e; + const uint64_t step = 0x0010000000000000ULL >> e; switch (mode) { case FPU_LIB_ROUND_NE: + away = frac > half || (frac == half && (u & step)); + break; case FPU_LIB_ROUND_MM: - return fpu_add64(d, fpu_bit_u64_to_f64(0x3FE0000000000000ULL | s)); + away = frac >= half; + break; case FPU_LIB_ROUND_DN: - if (s && fpu_is_fractional64(d)) { - return fpu_sub64(d, fpu_bit_u64_to_f64(0x3FF0000000000000ULL)); - } + away = s && frac; break; case FPU_LIB_ROUND_UP: - if (!s && fpu_is_fractional64(d)) { - return fpu_add64(d, fpu_bit_u64_to_f64(0x3FF0000000000000ULL)); - } + away = !s && frac; break; } - return d; + return fpu_bit_u64_to_f64((u - frac) + (away ? step : 0)); } #if defined(USE_SOFT_FPU_SQRT) From 5f27427c2c1ce1d98ff92138d14f78ddf7343050 Mon Sep 17 00:00:00 2001 From: Sol Astrius Phoenix Date: Sat, 18 Jul 2026 00:12:28 +0300 Subject: [PATCH 02/10] riscv_fpu: honor the static rounding-mode field Arithmetic and conversion ops only ever ran in the mode the frm CSR left on the host FPU; a static rm field (fadd.s ...,rtz etc.) was silently ignored. Apply a differing static mode around the op, gated on the op actually being rounding-capable: funct3 is an rm field only there, while on fsgnj/fmin/fmax/fcmp/fclass/fmv it encodes the operation itself and must not drive the host mode. Ops that take rm as an explicit argument (fcvt to integer, fround) consume the field directly and need no wrap. A static rmm field runs in RNE for now, which differs from roundTiesToAway only on exact halfway ties; those are covered once the exact RMM fixups land (#242). The dynamic-RMM preparation is gated to rm == DYN so it no longer fires under a static field. Signed-off-by: Sol Astrius Phoenix --- src/cpu/riscv_fpu.c | 48 ++++++++++++++++++++++++++++++++++++++++++--- src/cpu/riscv_fpu.h | 7 +++++++ 2 files changed, 52 insertions(+), 3 deletions(-) diff --git a/src/cpu/riscv_fpu.c b/src/cpu/riscv_fpu.c index 689a0af34..b065c6fab 100644 --- a/src/cpu/riscv_fpu.c +++ b/src/cpu/riscv_fpu.c @@ -93,7 +93,31 @@ static void riscv_prepare_rmm(rvvm_hart_t* vm, const uint32_t insn, const size_t } } -slow_path void riscv_emulate_f_opc_op(rvvm_hart_t* vm, const uint32_t insn) +// funct3 is an rm field only on the rounding-capable OP-FP ops; on +// fsgnj/fmin/fmax/fcmp/fclass/fmv it encodes the operation itself +static forceinline bool riscv_f_op_is_rounding(const uint32_t insn) +{ + switch (insn & 0xFE000000UL) { + case 0x00000000UL: // fadd.s + case 0x02000000UL: // fadd.d + case 0x08000000UL: // fsub.s + case 0x0A000000UL: // fsub.d + case 0x10000000UL: // fmul.s + case 0x12000000UL: // fmul.d + case 0x18000000UL: // fdiv.s + case 0x1A000000UL: // fdiv.d + case 0x58000000UL: // fsqrt.s + case 0x5A000000UL: // fsqrt.d + case 0x40000000UL: // fcvt.s.d, fround.s (Zfa) + case 0x42000000UL: // fcvt.d.s, fround.d (Zfa) + case 0xD0000000UL: // fcvt.s.w[u]/l[u] + case 0xD2000000UL: // fcvt.d.w[u]/l[u] + return true; + } + return false; +} + +static slow_path void riscv_emulate_f_opc_op_impl(rvvm_hart_t* vm, const uint32_t insn) { const size_t rds = bit_ext_u32(insn, 7, 5); const uint32_t rm = bit_ext_u32(insn, 12, 3); @@ -102,8 +126,8 @@ slow_path void riscv_emulate_f_opc_op(rvvm_hart_t* vm, const uint32_t insn) if (likely(riscv_fpu_is_enabled(vm))) { - if (unlikely(vm->csr.fcsr >> 5 == 0x04)) { - // Handle RMM rounding + if (unlikely(vm->csr.fcsr >> 5 == 0x04) && rm == 0x07) { + // Handle dynamic RMM rounding; a static rm field overrides frm riscv_prepare_rmm(vm, insn, rs1, rs2); } @@ -406,4 +430,22 @@ slow_path void riscv_emulate_f_opc_op(rvvm_hart_t* vm, const uint32_t insn) riscv_illegal_insn(vm, insn); } +slow_path void riscv_emulate_f_opc_op(rvvm_hart_t* vm, const uint32_t insn) +{ + const uint32_t rm = bit_ext_u32(insn, 12, 3); + // A static rm field on a rounding-capable op overrides the frm-tracked host + // mode around the op; ops taking rm as an argument consume the field directly + if (unlikely(rm != 0x07) && riscv_fpu_rm_is_valid(rm) && riscv_f_op_is_rounding(insn)) { + const uint32_t prev = fpu_get_rounding_mode(); + const uint32_t need = riscv_fpu_host_rm(rm); + if (need != prev) { + fpu_set_rounding_mode(need); + riscv_emulate_f_opc_op_impl(vm, insn); + fpu_set_rounding_mode(prev); + return; + } + } + riscv_emulate_f_opc_op_impl(vm, insn); +} + #endif diff --git a/src/cpu/riscv_fpu.h b/src/cpu/riscv_fpu.h index 0ad5354fb..734367832 100644 --- a/src/cpu/riscv_fpu.h +++ b/src/cpu/riscv_fpu.h @@ -26,6 +26,13 @@ static forceinline bool riscv_fpu_rm_is_valid(uint32_t rm) return rm > 1; } +// Host mode implementing a given rm: RMM has no host equivalent and runs in RNE, +// which differs only on exact halfway ties +static forceinline uint32_t riscv_fpu_host_rm(uint32_t rm) +{ + return (rm == FPU_LIB_ROUND_MM) ? FPU_LIB_ROUND_NE : rm; +} + // Bit-precise register read (raw low 32 bits, no NaN-box check) -- for fmv.x.w static forceinline fpu_f32_t riscv_view_s(rvvm_hart_t* vm, size_t reg) { From 8cf6ee167786ceeb2331e4280454b370c83342e7 Mon Sep 17 00:00:00 2001 From: Sol Astrius Phoenix Date: Sat, 18 Jul 2026 00:14:10 +0300 Subject: [PATCH 03/10] riscv_fpu: honor the rounding mode in the FMA ops fpu_fma rounds in the host mode, so a static rm field on fmadd/fmsub/fnmsub/fnmadd was silently computed in whatever mode frm left set. Route the FMA ops through riscv_fma32/64, which apply a differing static mode around the op using the same riscv_fpu_host_rm rule as the OP-FP dispatch. As there, a static rmm field runs in RNE until the exact ties-away fixups (#242) extend to FMA. Signed-off-by: Sol Astrius Phoenix --- src/cpu/riscv_fpu.h | 88 ++++++++++++++++++++++++++++++++------------- 1 file changed, 64 insertions(+), 24 deletions(-) diff --git a/src/cpu/riscv_fpu.h b/src/cpu/riscv_fpu.h index 734367832..a9d5f1638 100644 --- a/src/cpu/riscv_fpu.h +++ b/src/cpu/riscv_fpu.h @@ -133,6 +133,38 @@ static forceinline void riscv_emulate_f_opc_store(rvvm_hart_t* vm, const uint32_ riscv_illegal_insn(vm, insn); } +// FMA in the op's rounding mode: a static rm field that differs from the +// frm-tracked host mode is applied around the op, as in the OP-FP dispatch +static forceinline fpu_f32_t riscv_fma32(uint32_t rm, fpu_f32_t a, fpu_f32_t b, fpu_f32_t c) +{ + if (unlikely(rm != 0x07)) { + const uint32_t prev = fpu_get_rounding_mode(); + const uint32_t need = riscv_fpu_host_rm(rm); + if (need != prev) { + fpu_set_rounding_mode(need); + const fpu_f32_t r = fpu_fma32(a, b, c); + fpu_set_rounding_mode(prev); + return r; + } + } + return fpu_fma32(a, b, c); +} + +static forceinline fpu_f64_t riscv_fma64(uint32_t rm, fpu_f64_t a, fpu_f64_t b, fpu_f64_t c) +{ + if (unlikely(rm != 0x07)) { + const uint32_t prev = fpu_get_rounding_mode(); + const uint32_t need = riscv_fpu_host_rm(rm); + if (need != prev) { + fpu_set_rounding_mode(need); + const fpu_f64_t r = fpu_fma64(a, b, c); + fpu_set_rounding_mode(prev); + return r; + } + } + return fpu_fma64(a, b, c); +} + static forceinline void riscv_emulate_f_fmadd(rvvm_hart_t* vm, const uint32_t insn) { const size_t rds = bit_ext_u32(insn, 7, 5); @@ -145,15 +177,17 @@ static forceinline void riscv_emulate_f_fmadd(rvvm_hart_t* vm, const uint32_t in switch (bit_ext_u32(insn, 25, 2)) { case 0x0: // fmadd.s riscv_emit_s(vm, rds, - fpu_fma32(riscv_read_s(vm, rs1), // - riscv_read_s(vm, rs2), // - riscv_read_s(vm, rs3))); + riscv_fma32(rm, // + riscv_read_s(vm, rs1), // + riscv_read_s(vm, rs2), // + riscv_read_s(vm, rs3))); return; case 0x1: // fmadd.d riscv_emit_d(vm, rds, - fpu_fma64(riscv_view_d(vm, rs1), // - riscv_view_d(vm, rs2), // - riscv_view_d(vm, rs3))); + riscv_fma64(rm, // + riscv_view_d(vm, rs1), // + riscv_view_d(vm, rs2), // + riscv_view_d(vm, rs3))); return; } } @@ -173,15 +207,17 @@ static forceinline void riscv_emulate_f_fmsub(rvvm_hart_t* vm, const uint32_t in switch (bit_ext_u32(insn, 25, 2)) { case 0x0: // fmsub.s riscv_emit_s(vm, rds, - fpu_fma32(riscv_read_s(vm, rs1), // - riscv_read_s(vm, rs2), // - fpu_neg32(riscv_read_s(vm, rs3)))); + riscv_fma32(rm, // + riscv_read_s(vm, rs1), // + riscv_read_s(vm, rs2), // + fpu_neg32(riscv_read_s(vm, rs3)))); return; case 0x1: // fmsub.d riscv_emit_d(vm, rds, - fpu_fma64(riscv_view_d(vm, rs1), // - riscv_view_d(vm, rs2), // - fpu_neg64(riscv_view_d(vm, rs3)))); + riscv_fma64(rm, // + riscv_view_d(vm, rs1), // + riscv_view_d(vm, rs2), // + fpu_neg64(riscv_view_d(vm, rs3)))); return; } } @@ -201,15 +237,17 @@ static forceinline void riscv_emulate_f_fnmsub(rvvm_hart_t* vm, const uint32_t i switch (bit_ext_u32(insn, 25, 2)) { case 0x0: // fnmsub.s riscv_emit_s(vm, rds, - fpu_fma32(fpu_neg32(riscv_read_s(vm, rs1)), // - riscv_read_s(vm, rs2), // - riscv_read_s(vm, rs3))); + riscv_fma32(rm, // + fpu_neg32(riscv_read_s(vm, rs1)), // + riscv_read_s(vm, rs2), // + riscv_read_s(vm, rs3))); return; case 0x1: // fnmsub.d riscv_emit_d(vm, rds, - fpu_fma64(fpu_neg64(riscv_view_d(vm, rs1)), // - riscv_view_d(vm, rs2), // - riscv_view_d(vm, rs3))); + riscv_fma64(rm, // + fpu_neg64(riscv_view_d(vm, rs1)), // + riscv_view_d(vm, rs2), // + riscv_view_d(vm, rs3))); return; } } @@ -230,15 +268,17 @@ static forceinline void riscv_emulate_f_fnmadd(rvvm_hart_t* vm, const uint32_t i case 0x0: // fnmadd.s = -(rs1*rs2) - rs3; negate operands so the single // rounding sees the correctly-signed result (directed modes) riscv_emit_s(vm, rds, - fpu_fma32(fpu_neg32(riscv_read_s(vm, rs1)), // - riscv_read_s(vm, rs2), // - fpu_neg32(riscv_read_s(vm, rs3)))); + riscv_fma32(rm, // + fpu_neg32(riscv_read_s(vm, rs1)), // + riscv_read_s(vm, rs2), // + fpu_neg32(riscv_read_s(vm, rs3)))); return; case 0x1: // fnmadd.d riscv_emit_d(vm, rds, - fpu_fma64(fpu_neg64(riscv_view_d(vm, rs1)), // - riscv_view_d(vm, rs2), // - fpu_neg64(riscv_view_d(vm, rs3)))); + riscv_fma64(rm, // + fpu_neg64(riscv_view_d(vm, rs1)), // + riscv_view_d(vm, rs2), // + fpu_neg64(riscv_view_d(vm, rs3)))); return; } } From b17bfc5fcbe0f2fe58863641c35e870f73006285 Mon Sep 17 00:00:00 2001 From: Sol Astrius Phoenix Date: Sat, 18 Jul 2026 00:21:38 +0300 Subject: [PATCH 04/10] fpu_lib: set the FMA underflow flag by IEEE after-rounding tininess IEEE (and RISC-V) underflow means tiny after rounding: the result, rounded with an unbounded exponent range, lies below the minimum normal. Hosts that detect tininess before rounding (e.g. aarch64) disagree exactly when an FMA result lands on +/- the minimum normal, flagging a spurious UF for a result that rounded up out of the tiny range. On that boundary, recompute the after-rounding verdict and force the flag: - f32: the exact product widens into f64, the sum rounds once at 53 bits, and the scaled conversion once at 24; 53 >= 2*24 + 2 makes the double rounding innocuous, so the scaled magnitude compared against 1.0 is the unbounded-exponent verdict. - f64: no wider type exists, so rescale the op by 2^52 into the normal range, where the single fused rounding is already unbounded- equivalent (safety bounds are derived in a comment; an overflowing a*2^52 for b == 0 yields NaN, which correctly reads as not-tiny). The verdict merges over a flag snapshot from before the op, so a sticky UF from earlier ops is preserved and nothing the recompute raises can leak. fpu_fma64 is split into a bare fpu_fma64_raw plus the flag handling so the f64 fixup can reuse the fused op without re-entering it. On aarch64, fma(2^-62*(1-2^-23), 2^-64*(1+2^-23), 0) previously returned the correct 2^-126 but flagged NX|UF instead of NX; the f64 analog and the negative side misbehaved the same way. Signed-off-by: Sol Astrius Phoenix --- src/util/fpu_lib.h | 84 ++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 78 insertions(+), 6 deletions(-) diff --git a/src/util/fpu_lib.h b/src/util/fpu_lib.h index e636436af..f235e4939 100644 --- a/src/util/fpu_lib.h +++ b/src/util/fpu_lib.h @@ -311,6 +311,9 @@ fpu_f64_t fpu_sqrt64_soft_internal(fpu_f64_t d); #define FPU_LIB_FP32_NEGATIVE_ZERO 0x80000000U #define FPU_LIB_FP64_NEGATIVE_ZERO 0x8000000000000000ULL +#define FPU_LIB_FP32_MINIMUM_NORM 0x00800000U +#define FPU_LIB_FP64_MINIMUM_NORM 0x0010000000000000ULL + static forceinline uint16_t fpu_bit_f16_to_u16(fpu_f16_t f) { #if defined(USE_SOFT_FPU_ENCAP) @@ -1005,6 +1008,32 @@ static forceinline fpu_f64_t fpu_div64(fpu_f64_t a, fpu_f64_t b) return div; } +/* + * IEEE (and RISC-V) underflow means tiny *after* rounding: the result, rounded to + * full precision with an unbounded exponent range, lies below the minimum normal. + * Some hosts (e.g. aarch64) detect tininess before rounding instead; the two + * disagree exactly when the result is +/- the minimum normal, so on that boundary + * the after-rounding verdict is recomputed and the UF flag forced to match. + * old_exceptions preserves a UF that was already sticky before the op. + */ +static forceinline void fpu_fma32_fixup_uf(fpu_f32_t a, fpu_f32_t b, fpu_f32_t c, uint32_t old_exceptions) +{ + uint32_t exceptions = fpu_get_exceptions(); + // The product is exact in f64, the sum rounds once at 53 bits, the scaled + // conversion once at 24: 53 >= 2*24 + 2 makes the double rounding innocuous, + // so |rn| is the result of unbounded-exponent rounding, tiny iff below 1.0 + fpu_f64_t sum = fpu_add64(fpu_mul64(fpu_fcvt_f32_to_f64(a), fpu_fcvt_f32_to_f64(b)), + fpu_fcvt_f32_to_f64(c)); + fpu_f32_t rn = fpu_fcvt_f64_to_f32(fpu_mul64(sum, fpu_bit_u64_to_f64(0x47D0000000000000ULL))); // 2^126 + bool tiny = (fpu_bit_f32_to_u32(rn) & FPU_LIB_FP32_NOSIGNED_MASK) < 0x3F800000U; + if (tiny) { + exceptions |= FPU_LIB_FLAG_UF; + } else { + exceptions = (exceptions & ~FPU_LIB_FLAG_UF) | (old_exceptions & FPU_LIB_FLAG_UF); + } + fpu_set_exceptions(exceptions); +} + static forceinline func_opt_size fpu_f32_t fpu_fma32(fpu_f32_t a, fpu_f32_t b, fpu_f32_t c) { uint32_t old_exceptions = fpu_get_exceptions(); @@ -1024,6 +1053,10 @@ static forceinline func_opt_size fpu_f32_t fpu_fma32(fpu_f32_t a, fpu_f32_t b, f #endif #endif + if (unlikely((fpu_bit_f32_to_u32(ret) & FPU_LIB_FP32_NOSIGNED_MASK) == FPU_LIB_FP32_MINIMUM_NORM)) { + fpu_fma32_fixup_uf(a, b, c, old_exceptions); + } + uint32_t exceptions = fpu_get_exceptions(); if (invalid) { fpu_raise_invalid(); @@ -1034,13 +1067,11 @@ static forceinline func_opt_size fpu_f32_t fpu_fma32(fpu_f32_t a, fpu_f32_t b, f return ret; } -static forceinline fpu_f64_t fpu_fma64(fpu_f64_t a, fpu_f64_t b, fpu_f64_t c) +// The bare fused op: no soft NV check, no flag fixups +static forceinline fpu_f64_t fpu_fma64_raw(fpu_f64_t a, fpu_f64_t b, fpu_f64_t c) { - uint32_t old_exceptions = fpu_get_exceptions(); - bool invalid = fpu_fma64_invalid_soft(a, b, c); - #if defined(FPU_LIB_OPTIMAL_BUILTIN_FMA) - fpu_f64_t ret = fpu_wrap_f64(__builtin_fma(fpu_raw_f64(a), fpu_raw_f64(b), fpu_raw_f64(c))); + return fpu_wrap_f64(__builtin_fma(fpu_raw_f64(a), fpu_raw_f64(b), fpu_raw_f64(c))); #else fpu_f64_t mul = fpu_mul64(a, b); fpu_f64_t e_m = fpu_mul_error64(mul, a, b); @@ -1048,8 +1079,49 @@ static forceinline fpu_f64_t fpu_fma64(fpu_f64_t a, fpu_f64_t b, fpu_f64_t 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)); - fpu_f64_t ret = fpu_add64(sum, err); + return fpu_add64(sum, err); #endif +} + +/* + * f64 counterpart of fpu_fma32_fixup_uf: no wider type exists, so rescale the op + * by 2^52 into the normal range, where the single fused rounding is already + * unbounded-equivalent, and test the scaled result against 2^-970. + * + * The scaling is safe: the result rounds to +/-2^-1022, so at least one of a*b, c + * has a set bit at 2^-1022 or below (else the sum is either 0 or larger). If it is + * in c, |c| <= 2^-969 and by triangle inequality |a*b| <= 2^-968; if it is in the + * exact product (up to 106 bits wide), |a*b| <= 2^-916 and |c| <= 2^-915. Either + * way |c|*2^52 is tiny and, for b != 0, |b| >= 2^-1074 gives |a| <= 2^159, so + * a*2^52 cannot overflow. If b == 0 then a is unbounded and a*2^52 may overflow to + * inf, making rs NaN -- which correctly reads as "not tiny", since a result of + * exactly +/-2^-1022 from b == 0 is the exact value of c, and exact never + * underflows. + */ +static forceinline void fpu_fma64_fixup_uf(fpu_f64_t a, fpu_f64_t b, fpu_f64_t c, uint32_t old_exceptions) +{ + uint32_t exceptions = fpu_get_exceptions(); + fpu_f64_t scale = fpu_bit_u64_to_f64(0x4330000000000000ULL); // 2^52 + fpu_f64_t rs = fpu_fma64_raw(fpu_mul64(a, scale), b, fpu_mul64(c, scale)); + bool tiny = (fpu_bit_f64_to_u64(rs) & FPU_LIB_FP64_NOSIGNED_MASK) < 0x0350000000000000ULL; // 2^-970 + if (tiny) { + exceptions |= FPU_LIB_FLAG_UF; + } else { + exceptions = (exceptions & ~FPU_LIB_FLAG_UF) | (old_exceptions & FPU_LIB_FLAG_UF); + } + fpu_set_exceptions(exceptions); +} + +static forceinline fpu_f64_t fpu_fma64(fpu_f64_t a, fpu_f64_t b, fpu_f64_t c) +{ + uint32_t old_exceptions = fpu_get_exceptions(); + bool invalid = fpu_fma64_invalid_soft(a, b, c); + + fpu_f64_t ret = fpu_fma64_raw(a, b, c); + + if (unlikely((fpu_bit_f64_to_u64(ret) & FPU_LIB_FP64_NOSIGNED_MASK) == FPU_LIB_FP64_MINIMUM_NORM)) { + fpu_fma64_fixup_uf(a, b, c, old_exceptions); + } uint32_t exceptions = fpu_get_exceptions(); if (invalid) { From 1587da8db0e6e5eace36a0d486f34c9ba1f20713 Mon Sep 17 00:00:00 2001 From: Sol Astrius Phoenix Date: Mon, 20 Jul 2026 17:43:00 +0300 Subject: [PATCH 05/10] fpu_lib: clarify the integral-rounding edge cases Review feedback on #246: spell out the 0.5 ties-to-even case, the e == 0 odd-integer test through the biased-exponent LSB, and the mantissa-to- exponent carry on rounding away; use shifted constants for half/step and annotate the 1.0 literals. No functional change (re-verified bit-exact against host rint()/round(), 2.0e9 cases). Signed-off-by: Sol Astrius Phoenix --- src/util/fpu_lib.c | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/util/fpu_lib.c b/src/util/fpu_lib.c index 15c50d8ed..270cbb7c5 100644 --- a/src/util/fpu_lib.c +++ b/src/util/fpu_lib.c @@ -423,9 +423,11 @@ slow_path fpu_f32_t fpu_round_f32_internal(fpu_f32_t f, uint32_t mode) return f; // Already integral: |f| >= 2^23, +/-0, inf, NaN } if (e < 0) { - // |f| in (0, 1) rounds to +/-0 or +/-1; the 0.5 tie goes to even == 0 + // |f| in (0, 1) rounds to +/-0 or +/-1 switch (mode) { case FPU_LIB_ROUND_NE: + // |f| in (0.5, 1) is nearest to 1; exactly 0.5 is a tie and goes + // to the even neighbour, which is 0 away = (e == -1) && (u & FPU_LIB_FP32_MANTISSA_MASK); break; case FPU_LIB_ROUND_MM: @@ -438,14 +440,17 @@ slow_path fpu_f32_t fpu_round_f32_internal(fpu_f32_t f, uint32_t mode) away = !s; break; } - return fpu_bit_u32_to_f32(s | (away ? 0x3F800000U : 0)); + return fpu_bit_u32_to_f32(s | (away ? 0x3F800000U : 0)); // +/-1.0 or +/-0.0 } // |f| in [1, 2^23): split the mantissa into integer part and fraction bits const uint32_t frac = u & (FPU_LIB_FP32_MANTISSA_MASK >> e); - const uint32_t half = 0x00400000U >> e; - const uint32_t step = 0x00800000U >> e; // 1.0 at this exponent; carries into a new binade + const uint32_t half = (1U << 22) >> e; + const uint32_t step = (1U << 23) >> e; // 1.0 at this exponent switch (mode) { case FPU_LIB_ROUND_NE: + // On a tie, round up iff the integer part is odd. For e == 0 the step + // bit is the exponent LSB rather than a mantissa bit, but the biased + // exponent of [1, 2) is odd (0x7F/0x3FF), matching its odd integer 1. away = frac > half || (frac == half && (u & step)); break; case FPU_LIB_ROUND_MM: @@ -458,6 +463,10 @@ slow_path fpu_f32_t fpu_round_f32_internal(fpu_f32_t f, uint32_t mode) away = !s && frac; break; } + // Adding step may carry from the mantissa into the exponent field: that only + // happens when the truncated mantissa wraps to zero, i.e. when rounding away + // lands exactly on the next power of two, where the carry is the intended + // encoding (the same trick as the classic nextafter bit-increment). return fpu_bit_u32_to_f32((u - frac) + (away ? step : 0)); } @@ -488,11 +497,11 @@ slow_path fpu_f64_t fpu_round_f64_internal(fpu_f64_t d, uint32_t mode) away = !s; break; } - return fpu_bit_u64_to_f64(s | (away ? 0x3FF0000000000000ULL : 0)); + return fpu_bit_u64_to_f64(s | (away ? 0x3FF0000000000000ULL : 0)); // +/-1.0 or +/-0.0 } const uint64_t frac = u & (FPU_LIB_FP64_MANTISSA_MASK >> e); - const uint64_t half = 0x0008000000000000ULL >> e; - const uint64_t step = 0x0010000000000000ULL >> e; + const uint64_t half = (1ULL << 51) >> e; + const uint64_t step = (1ULL << 52) >> e; // 1.0 at this exponent switch (mode) { case FPU_LIB_ROUND_NE: away = frac > half || (frac == half && (u & step)); From b3a3f3ae6ff6c1a4a954a83c2281859806e9cac1 Mon Sep 17 00:00:00 2001 From: Sol Astrius Phoenix Date: Mon, 20 Jul 2026 17:43:04 +0300 Subject: [PATCH 06/10] riscv_fpu: share static-rm handling between the OP-FP and FMA dispatches Review feedback on #246: hoist the static-rm application out of the per-op FMA helpers into riscv_fpu_static_rm_enter/leave, used identically by the OP-FP wrapper and the four FMA emulate functions -- one code path, less forceinline bloat. Both paths reject rm 5/6 before entering (the FMA functions via their existing riscv_fpu_rm_is_valid gate). Gate the RMM preparation on the effective mode instead of only the dynamic one: a static rmm field now behaves exactly like frm == RMM rather than approximating RNE, so the two select the same (interim) RMM path until the exact fixups land. Rename the predicate to riscv_f_op_is_implicitly_rounding to distinguish the ops that consume rm as an explicit argument. Signed-off-by: Sol Astrius Phoenix --- src/cpu/riscv_fpu.c | 31 +++++------ src/cpu/riscv_fpu.h | 132 +++++++++++++++++++++++--------------------- 2 files changed, 83 insertions(+), 80 deletions(-) diff --git a/src/cpu/riscv_fpu.c b/src/cpu/riscv_fpu.c index b065c6fab..6026e40a3 100644 --- a/src/cpu/riscv_fpu.c +++ b/src/cpu/riscv_fpu.c @@ -94,8 +94,10 @@ static void riscv_prepare_rmm(rvvm_hart_t* vm, const uint32_t insn, const size_t } // funct3 is an rm field only on the rounding-capable OP-FP ops; on -// fsgnj/fmin/fmax/fcmp/fclass/fmv it encodes the operation itself -static forceinline bool riscv_f_op_is_rounding(const uint32_t insn) +// fsgnj/fmin/fmax/fcmp/fclass/fmv it encodes the operation itself. "Implicitly" +// rounding: ops that take rm as an argument (fcvt to integer, fround) consume +// the field themselves and are not listed here. +static forceinline bool riscv_f_op_is_implicitly_rounding(const uint32_t insn) { switch (insn & 0xFE000000UL) { case 0x00000000UL: // fadd.s @@ -126,8 +128,9 @@ static slow_path void riscv_emulate_f_opc_op_impl(rvvm_hart_t* vm, const uint32_ if (likely(riscv_fpu_is_enabled(vm))) { - if (unlikely(vm->csr.fcsr >> 5 == 0x04) && rm == 0x07) { - // Handle dynamic RMM rounding; a static rm field overrides frm + if (unlikely(((rm == 0x07) ? vm->csr.fcsr >> 5 : rm) == 0x04)) { + // Handle RMM rounding in the effective mode: a static rmm field + // behaves exactly like frm == RMM riscv_prepare_rmm(vm, insn, rs1, rs2); } @@ -433,19 +436,15 @@ static slow_path void riscv_emulate_f_opc_op_impl(rvvm_hart_t* vm, const uint32_ slow_path void riscv_emulate_f_opc_op(rvvm_hart_t* vm, const uint32_t insn) { const uint32_t rm = bit_ext_u32(insn, 12, 3); - // A static rm field on a rounding-capable op overrides the frm-tracked host - // mode around the op; ops taking rm as an argument consume the field directly - if (unlikely(rm != 0x07) && riscv_fpu_rm_is_valid(rm) && riscv_f_op_is_rounding(insn)) { - const uint32_t prev = fpu_get_rounding_mode(); - const uint32_t need = riscv_fpu_host_rm(rm); - if (need != prev) { - fpu_set_rounding_mode(need); - riscv_emulate_f_opc_op_impl(vm, insn); - fpu_set_rounding_mode(prev); - return; - } + // A static rm field on an implicitly rounding op overrides the frm-tracked + // host mode around the op + if (unlikely(rm != 0x07) && riscv_fpu_rm_is_valid(rm) && riscv_f_op_is_implicitly_rounding(insn)) { + const uint32_t prev_rm = riscv_fpu_static_rm_enter(rm); + riscv_emulate_f_opc_op_impl(vm, insn); + riscv_fpu_static_rm_leave(prev_rm); + } else { + riscv_emulate_f_opc_op_impl(vm, insn); } - riscv_emulate_f_opc_op_impl(vm, insn); } #endif diff --git a/src/cpu/riscv_fpu.h b/src/cpu/riscv_fpu.h index a9d5f1638..2756c6adf 100644 --- a/src/cpu/riscv_fpu.h +++ b/src/cpu/riscv_fpu.h @@ -33,6 +33,30 @@ static forceinline uint32_t riscv_fpu_host_rm(uint32_t rm) return (rm == FPU_LIB_ROUND_MM) ? FPU_LIB_ROUND_NE : rm; } +/* + * A rounding-capable op runs under its effective mode: the frm-tracked host mode + * when rm == DYN, or the static rm field applied around the op. Enter returns the + * mode to restore afterwards via leave, or DYN when nothing was changed. + */ +static forceinline uint32_t riscv_fpu_static_rm_enter(uint32_t rm) +{ + if (unlikely(rm != 0x07)) { + const uint32_t prev = fpu_get_rounding_mode(); + if (riscv_fpu_host_rm(rm) != prev) { + fpu_set_rounding_mode(riscv_fpu_host_rm(rm)); + return prev; + } + } + return 0x07; +} + +static forceinline void riscv_fpu_static_rm_leave(uint32_t prev) +{ + if (unlikely(prev != 0x07)) { + fpu_set_rounding_mode(prev); + } +} + // Bit-precise register read (raw low 32 bits, no NaN-box check) -- for fmv.x.w static forceinline fpu_f32_t riscv_view_s(rvvm_hart_t* vm, size_t reg) { @@ -133,38 +157,6 @@ static forceinline void riscv_emulate_f_opc_store(rvvm_hart_t* vm, const uint32_ riscv_illegal_insn(vm, insn); } -// FMA in the op's rounding mode: a static rm field that differs from the -// frm-tracked host mode is applied around the op, as in the OP-FP dispatch -static forceinline fpu_f32_t riscv_fma32(uint32_t rm, fpu_f32_t a, fpu_f32_t b, fpu_f32_t c) -{ - if (unlikely(rm != 0x07)) { - const uint32_t prev = fpu_get_rounding_mode(); - const uint32_t need = riscv_fpu_host_rm(rm); - if (need != prev) { - fpu_set_rounding_mode(need); - const fpu_f32_t r = fpu_fma32(a, b, c); - fpu_set_rounding_mode(prev); - return r; - } - } - return fpu_fma32(a, b, c); -} - -static forceinline fpu_f64_t riscv_fma64(uint32_t rm, fpu_f64_t a, fpu_f64_t b, fpu_f64_t c) -{ - if (unlikely(rm != 0x07)) { - const uint32_t prev = fpu_get_rounding_mode(); - const uint32_t need = riscv_fpu_host_rm(rm); - if (need != prev) { - fpu_set_rounding_mode(need); - const fpu_f64_t r = fpu_fma64(a, b, c); - fpu_set_rounding_mode(prev); - return r; - } - } - return fpu_fma64(a, b, c); -} - static forceinline void riscv_emulate_f_fmadd(rvvm_hart_t* vm, const uint32_t insn) { const size_t rds = bit_ext_u32(insn, 7, 5); @@ -174,22 +166,25 @@ static forceinline void riscv_emulate_f_fmadd(rvvm_hart_t* vm, const uint32_t in const size_t rs3 = insn >> 27; if (likely(riscv_fpu_is_enabled(vm) && riscv_fpu_rm_is_valid(rm))) { + // A static rm field overrides the frm-tracked host mode, as in the OP-FP dispatch + const uint32_t prev_rm = riscv_fpu_static_rm_enter(rm); switch (bit_ext_u32(insn, 25, 2)) { case 0x0: // fmadd.s riscv_emit_s(vm, rds, - riscv_fma32(rm, // - riscv_read_s(vm, rs1), // - riscv_read_s(vm, rs2), // - riscv_read_s(vm, rs3))); + fpu_fma32(riscv_read_s(vm, rs1), // + riscv_read_s(vm, rs2), // + riscv_read_s(vm, rs3))); + riscv_fpu_static_rm_leave(prev_rm); return; case 0x1: // fmadd.d riscv_emit_d(vm, rds, - riscv_fma64(rm, // - riscv_view_d(vm, rs1), // - riscv_view_d(vm, rs2), // - riscv_view_d(vm, rs3))); + fpu_fma64(riscv_view_d(vm, rs1), // + riscv_view_d(vm, rs2), // + riscv_view_d(vm, rs3))); + riscv_fpu_static_rm_leave(prev_rm); return; } + riscv_fpu_static_rm_leave(prev_rm); } riscv_illegal_insn(vm, insn); @@ -204,22 +199,25 @@ static forceinline void riscv_emulate_f_fmsub(rvvm_hart_t* vm, const uint32_t in const size_t rs3 = insn >> 27; if (likely(riscv_fpu_is_enabled(vm) && riscv_fpu_rm_is_valid(rm))) { + // A static rm field overrides the frm-tracked host mode, as in the OP-FP dispatch + const uint32_t prev_rm = riscv_fpu_static_rm_enter(rm); switch (bit_ext_u32(insn, 25, 2)) { case 0x0: // fmsub.s riscv_emit_s(vm, rds, - riscv_fma32(rm, // - riscv_read_s(vm, rs1), // - riscv_read_s(vm, rs2), // - fpu_neg32(riscv_read_s(vm, rs3)))); + fpu_fma32(riscv_read_s(vm, rs1), // + riscv_read_s(vm, rs2), // + fpu_neg32(riscv_read_s(vm, rs3)))); + riscv_fpu_static_rm_leave(prev_rm); return; case 0x1: // fmsub.d riscv_emit_d(vm, rds, - riscv_fma64(rm, // - riscv_view_d(vm, rs1), // - riscv_view_d(vm, rs2), // - fpu_neg64(riscv_view_d(vm, rs3)))); + fpu_fma64(riscv_view_d(vm, rs1), // + riscv_view_d(vm, rs2), // + fpu_neg64(riscv_view_d(vm, rs3)))); + riscv_fpu_static_rm_leave(prev_rm); return; } + riscv_fpu_static_rm_leave(prev_rm); } riscv_illegal_insn(vm, insn); @@ -234,22 +232,25 @@ static forceinline void riscv_emulate_f_fnmsub(rvvm_hart_t* vm, const uint32_t i const size_t rs3 = insn >> 27; if (likely(riscv_fpu_is_enabled(vm) && riscv_fpu_rm_is_valid(rm))) { + // A static rm field overrides the frm-tracked host mode, as in the OP-FP dispatch + const uint32_t prev_rm = riscv_fpu_static_rm_enter(rm); switch (bit_ext_u32(insn, 25, 2)) { case 0x0: // fnmsub.s riscv_emit_s(vm, rds, - riscv_fma32(rm, // - fpu_neg32(riscv_read_s(vm, rs1)), // - riscv_read_s(vm, rs2), // - riscv_read_s(vm, rs3))); + fpu_fma32(fpu_neg32(riscv_read_s(vm, rs1)), // + riscv_read_s(vm, rs2), // + riscv_read_s(vm, rs3))); + riscv_fpu_static_rm_leave(prev_rm); return; case 0x1: // fnmsub.d riscv_emit_d(vm, rds, - riscv_fma64(rm, // - fpu_neg64(riscv_view_d(vm, rs1)), // - riscv_view_d(vm, rs2), // - riscv_view_d(vm, rs3))); + fpu_fma64(fpu_neg64(riscv_view_d(vm, rs1)), // + riscv_view_d(vm, rs2), // + riscv_view_d(vm, rs3))); + riscv_fpu_static_rm_leave(prev_rm); return; } + riscv_fpu_static_rm_leave(prev_rm); } riscv_illegal_insn(vm, insn); @@ -264,23 +265,26 @@ static forceinline void riscv_emulate_f_fnmadd(rvvm_hart_t* vm, const uint32_t i const size_t rs3 = insn >> 27; if (likely(riscv_fpu_is_enabled(vm) && riscv_fpu_rm_is_valid(rm))) { + // A static rm field overrides the frm-tracked host mode, as in the OP-FP dispatch + const uint32_t prev_rm = riscv_fpu_static_rm_enter(rm); switch (bit_ext_u32(insn, 25, 2)) { case 0x0: // fnmadd.s = -(rs1*rs2) - rs3; negate operands so the single // rounding sees the correctly-signed result (directed modes) riscv_emit_s(vm, rds, - riscv_fma32(rm, // - fpu_neg32(riscv_read_s(vm, rs1)), // - riscv_read_s(vm, rs2), // - fpu_neg32(riscv_read_s(vm, rs3)))); + fpu_fma32(fpu_neg32(riscv_read_s(vm, rs1)), // + riscv_read_s(vm, rs2), // + fpu_neg32(riscv_read_s(vm, rs3)))); + riscv_fpu_static_rm_leave(prev_rm); return; case 0x1: // fnmadd.d riscv_emit_d(vm, rds, - riscv_fma64(rm, // - fpu_neg64(riscv_view_d(vm, rs1)), // - riscv_view_d(vm, rs2), // - fpu_neg64(riscv_view_d(vm, rs3)))); + fpu_fma64(fpu_neg64(riscv_view_d(vm, rs1)), // + riscv_view_d(vm, rs2), // + fpu_neg64(riscv_view_d(vm, rs3)))); + riscv_fpu_static_rm_leave(prev_rm); return; } + riscv_fpu_static_rm_leave(prev_rm); } riscv_illegal_insn(vm, insn); From c2379db76a7d33aa6a52216c618d1f7b0e195f1b Mon Sep 17 00:00:00 2001 From: Sol Astrius Phoenix Date: Mon, 20 Jul 2026 17:43:08 +0300 Subject: [PATCH 07/10] fpu_lib: odd-round the FMA32 underflow recompute The boundary recompute rounded the widened sum once at 53 bits and then converted at 24; when the 53-bit rounding lands exactly on a 24-bit tie midpoint the conversion resolves the tie blindly and the tininess verdict flips. Concretely, an fma32 with significand product 2^47 + t (t < 2^18, quantum 2^-198, negative) and c = 2^-126 has an exact value just below the midpoint 2^-126 - 2^-151: the delivered result is 2^-126 but the unbounded rounding is below the minimum normal, so UF must be raised -- the plain recompute concluded not-tiny and cleared it. Odd-round the sum against its exact TwoSum error before scaling, as fpu_fma32 itself does: odd rounding at 53 bits composes correctly with any narrower final rounding (53 >= 24 + 2), so the scaled conversion yields the true unbounded-exponent result. Also tighten the f64 bounds comment (|a| <= 2^158, per review). Signed-off-by: Sol Astrius Phoenix --- src/util/fpu_lib.h | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/util/fpu_lib.h b/src/util/fpu_lib.h index f235e4939..c3203db7a 100644 --- a/src/util/fpu_lib.h +++ b/src/util/fpu_lib.h @@ -1019,12 +1019,15 @@ static forceinline fpu_f64_t fpu_div64(fpu_f64_t a, fpu_f64_t b) static forceinline void fpu_fma32_fixup_uf(fpu_f32_t a, fpu_f32_t b, fpu_f32_t c, uint32_t old_exceptions) { uint32_t exceptions = fpu_get_exceptions(); - // The product is exact in f64, the sum rounds once at 53 bits, the scaled - // conversion once at 24: 53 >= 2*24 + 2 makes the double rounding innocuous, - // so |rn| is the result of unbounded-exponent rounding, tiny iff below 1.0 - fpu_f64_t sum = fpu_add64(fpu_mul64(fpu_fcvt_f32_to_f64(a), fpu_fcvt_f32_to_f64(b)), - fpu_fcvt_f32_to_f64(c)); - fpu_f32_t rn = fpu_fcvt_f64_to_f32(fpu_mul64(sum, fpu_bit_u64_to_f64(0x47D0000000000000ULL))); // 2^126 + // The product is exact in f64 and the sum is odd-rounded at 53 bits, so the + // scaled 24-bit conversion cannot double-round across a tie (53 >= 24 + 2): + // |rn| is the unbounded-exponent rounding of the exact result, tiny iff the + // 2^126-scaled magnitude stays below 1.0 + 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 res = fpu_odd_round64(sum, fpu_add_error64(sum, mul, add)); + fpu_f32_t rn = fpu_fcvt_f64_to_f32(fpu_mul64(res, fpu_bit_u64_to_f64(0x47D0000000000000ULL))); // 2^126 bool tiny = (fpu_bit_f32_to_u32(rn) & FPU_LIB_FP32_NOSIGNED_MASK) < 0x3F800000U; if (tiny) { exceptions |= FPU_LIB_FLAG_UF; @@ -1092,7 +1095,7 @@ static forceinline fpu_f64_t fpu_fma64_raw(fpu_f64_t a, fpu_f64_t b, fpu_f64_t c * has a set bit at 2^-1022 or below (else the sum is either 0 or larger). If it is * in c, |c| <= 2^-969 and by triangle inequality |a*b| <= 2^-968; if it is in the * exact product (up to 106 bits wide), |a*b| <= 2^-916 and |c| <= 2^-915. Either - * way |c|*2^52 is tiny and, for b != 0, |b| >= 2^-1074 gives |a| <= 2^159, so + * way |c|*2^52 is tiny and, for b != 0, |b| >= 2^-1074 gives |a| <= 2^158, so * a*2^52 cannot overflow. If b == 0 then a is unbounded and a*2^52 may overflow to * inf, making rs NaN -- which correctly reads as "not tiny", since a result of * exactly +/-2^-1022 from b == 0 is the exact value of c, and exact never From de1b6f655878eafc3d04d465222e53d845a6b3f6 Mon Sep 17 00:00:00 2001 From: Sol Astrius Phoenix Date: Mon, 20 Jul 2026 17:51:30 +0300 Subject: [PATCH 08/10] fpu_lib: clang-format the touched declarations Signed-off-by: Sol Astrius Phoenix --- src/util/fpu_lib.c | 16 ++++++++-------- src/util/fpu_lib.h | 22 +++++++++++----------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/util/fpu_lib.c b/src/util/fpu_lib.c index 270cbb7c5..9d39191bf 100644 --- a/src/util/fpu_lib.c +++ b/src/util/fpu_lib.c @@ -412,10 +412,10 @@ slow_path uint32_t fpu_fclass64(fpu_f64_t d) */ slow_path fpu_f32_t fpu_round_f32_internal(fpu_f32_t f, uint32_t mode) { - const uint32_t u = fpu_bit_f32_to_u32(f); - const uint32_t s = u & FPU_LIB_FP32_SIGNEDFP_MASK; - const int32_t e = fpu_exponent32(f); - bool away = false; + const uint32_t u = fpu_bit_f32_to_u32(f); + const uint32_t s = u & FPU_LIB_FP32_SIGNEDFP_MASK; + const int32_t e = fpu_exponent32(f); + bool away = false; if (unlikely(mode > FPU_LIB_ROUND_MM)) { mode = fpu_get_rounding_mode(); } @@ -472,10 +472,10 @@ slow_path fpu_f32_t fpu_round_f32_internal(fpu_f32_t f, uint32_t mode) slow_path fpu_f64_t fpu_round_f64_internal(fpu_f64_t d, uint32_t mode) { - const uint64_t u = fpu_bit_f64_to_u64(d); - const uint64_t s = u & FPU_LIB_FP64_SIGNEDFP_MASK; - const int32_t e = fpu_exponent64(d); - bool away = false; + const uint64_t u = fpu_bit_f64_to_u64(d); + const uint64_t s = u & FPU_LIB_FP64_SIGNEDFP_MASK; + const int32_t e = fpu_exponent64(d); + bool away = false; if (unlikely(mode > FPU_LIB_ROUND_MM)) { mode = fpu_get_rounding_mode(); } diff --git a/src/util/fpu_lib.h b/src/util/fpu_lib.h index c3203db7a..f75898261 100644 --- a/src/util/fpu_lib.h +++ b/src/util/fpu_lib.h @@ -1023,12 +1023,12 @@ static forceinline void fpu_fma32_fixup_uf(fpu_f32_t a, fpu_f32_t b, fpu_f32_t c // scaled 24-bit conversion cannot double-round across a tie (53 >= 24 + 2): // |rn| is the unbounded-exponent rounding of the exact result, tiny iff the // 2^126-scaled magnitude stays below 1.0 - 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 res = fpu_odd_round64(sum, fpu_add_error64(sum, mul, add)); - fpu_f32_t rn = fpu_fcvt_f64_to_f32(fpu_mul64(res, fpu_bit_u64_to_f64(0x47D0000000000000ULL))); // 2^126 - bool tiny = (fpu_bit_f32_to_u32(rn) & FPU_LIB_FP32_NOSIGNED_MASK) < 0x3F800000U; + 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 res = fpu_odd_round64(sum, fpu_add_error64(sum, mul, add)); + fpu_f32_t rn = fpu_fcvt_f64_to_f32(fpu_mul64(res, fpu_bit_u64_to_f64(0x47D0000000000000ULL))); // 2^126 + bool tiny = (fpu_bit_f32_to_u32(rn) & FPU_LIB_FP32_NOSIGNED_MASK) < 0x3F800000U; if (tiny) { exceptions |= FPU_LIB_FLAG_UF; } else { @@ -1103,10 +1103,10 @@ static forceinline fpu_f64_t fpu_fma64_raw(fpu_f64_t a, fpu_f64_t b, fpu_f64_t c */ static forceinline void fpu_fma64_fixup_uf(fpu_f64_t a, fpu_f64_t b, fpu_f64_t c, uint32_t old_exceptions) { - uint32_t exceptions = fpu_get_exceptions(); - fpu_f64_t scale = fpu_bit_u64_to_f64(0x4330000000000000ULL); // 2^52 - fpu_f64_t rs = fpu_fma64_raw(fpu_mul64(a, scale), b, fpu_mul64(c, scale)); - bool tiny = (fpu_bit_f64_to_u64(rs) & FPU_LIB_FP64_NOSIGNED_MASK) < 0x0350000000000000ULL; // 2^-970 + uint32_t exceptions = fpu_get_exceptions(); + fpu_f64_t scale = fpu_bit_u64_to_f64(0x4330000000000000ULL); // 2^52 + fpu_f64_t rs = fpu_fma64_raw(fpu_mul64(a, scale), b, fpu_mul64(c, scale)); + bool tiny = (fpu_bit_f64_to_u64(rs) & FPU_LIB_FP64_NOSIGNED_MASK) < 0x0350000000000000ULL; // 2^-970 if (tiny) { exceptions |= FPU_LIB_FLAG_UF; } else { @@ -1118,7 +1118,7 @@ static forceinline void fpu_fma64_fixup_uf(fpu_f64_t a, fpu_f64_t b, fpu_f64_t c static forceinline fpu_f64_t fpu_fma64(fpu_f64_t a, fpu_f64_t b, fpu_f64_t c) { uint32_t old_exceptions = fpu_get_exceptions(); - bool invalid = fpu_fma64_invalid_soft(a, b, c); + bool invalid = fpu_fma64_invalid_soft(a, b, c); fpu_f64_t ret = fpu_fma64_raw(a, b, c); From bd053f08e4efcb723884671a53ceabb9be403b49 Mon Sep 17 00:00:00 2001 From: Sol Astrius Phoenix Date: Mon, 20 Jul 2026 19:58:04 +0300 Subject: [PATCH 09/10] fpu_lib: decide FMA32 tininess by an f64 threshold compare Per review: instead of reconstructing the unbounded 24-bit rounding through odd rounding, scaling and a conversion, compare the 53-bit sum directly against the magnitude below which that rounding falls under 2^-126 -- the to-nearest midpoint 2^-126 - 2^-151 (exclusive; an exact tie there resolves up), 2^-126 - 2^-150 inclusive when rounding away from zero, and 2^-126 exclusive when rounding toward zero. The sum alone is ambiguous only when it lands exactly on the to-nearest threshold, where the exact TwoSum error breaks the tie -- and round-to- nearest is precisely where 2Sum is exact, so the directed-rounding 2Sum caveat is never in play: in the directed modes the sum errs in the same direction as the rounding, making the plain compare already conclusive. Covers the same below-midpoint counterexample as before, plus new vectors for the exact midpoint (err == 0), the away-inclusive threshold under RNE and RUP, and the from-above RDN case. Signed-off-by: Sol Astrius Phoenix --- src/util/fpu_lib.h | 35 ++++++++++++++++++++++++++++------- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/src/util/fpu_lib.h b/src/util/fpu_lib.h index f75898261..d877be1cf 100644 --- a/src/util/fpu_lib.h +++ b/src/util/fpu_lib.h @@ -1019,16 +1019,37 @@ static forceinline fpu_f64_t fpu_div64(fpu_f64_t a, fpu_f64_t b) static forceinline void fpu_fma32_fixup_uf(fpu_f32_t a, fpu_f32_t b, fpu_f32_t c, uint32_t old_exceptions) { uint32_t exceptions = fpu_get_exceptions(); - // The product is exact in f64 and the sum is odd-rounded at 53 bits, so the - // scaled 24-bit conversion cannot double-round across a tie (53 >= 24 + 2): - // |rn| is the unbounded-exponent rounding of the exact result, tiny iff the - // 2^126-scaled magnitude stays below 1.0 + // The product is exact in f64, so sum is the exact result rounded once at 53 + // bits; tininess is decided by comparing it in f64 against the mode-dependent + // magnitude below which the unbounded 24-bit rounding falls under 2^-126: + // - to-nearest: the midpoint 2^-126 - 2^-151, exclusive (both NE and MM + // resolve an exact-midpoint tie up to 2^-126). Only here can sum sit exactly + // on the threshold with the true result on either side; the exact TwoSum + // error breaks the tie (2Sum is exact in round-to-nearest). + // - rounding away from zero: the largest 24-bit value below, 2^-126 - 2^-150, + // inclusive. sum is rounded away too, so sum <= T already implies the exact + // result is <= T. + // - rounding toward zero: 2^-126 itself, exclusive; likewise sum < T implies + // the exact result is < T. 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 res = fpu_odd_round64(sum, fpu_add_error64(sum, mul, add)); - fpu_f32_t rn = fpu_fcvt_f64_to_f32(fpu_mul64(res, fpu_bit_u64_to_f64(0x47D0000000000000ULL))); // 2^126 - bool tiny = (fpu_bit_f32_to_u32(rn) & FPU_LIB_FP32_NOSIGNED_MASK) < 0x3F800000U; + uint64_t bits = fpu_bit_f64_to_u64(sum); + uint64_t mag = bits & FPU_LIB_FP64_NOSIGNED_MASK; + uint32_t mode = fpu_get_rounding_mode(); + bool away = (mode == FPU_LIB_ROUND_UP) == !(bits >> 63); + bool tiny; + if (mode == FPU_LIB_ROUND_NE || mode == FPU_LIB_ROUND_MM) { + tiny = mag < 0x380FFFFFF0000000ULL; // 2^-126 - 2^-151 + if (mag == 0x380FFFFFF0000000ULL) { + uint64_t err = fpu_bit_f64_to_u64(fpu_add_error64(sum, mul, add)); + tiny = (err << 1) != 0 && ((err ^ bits) >> 63); // exact result below the midpoint + } + } else if ((mode == FPU_LIB_ROUND_UP || mode == FPU_LIB_ROUND_DN) && away) { + tiny = mag <= 0x380FFFFFE0000000ULL; // 2^-126 - 2^-150 + } else { + tiny = mag < 0x3810000000000000ULL; // 2^-126 + } if (tiny) { exceptions |= FPU_LIB_FLAG_UF; } else { From 4c1c294afadf91f093dc5ca6fe6f7e637a8c4923 Mon Sep 17 00:00:00 2001 From: Sol Astrius Phoenix Date: Mon, 20 Jul 2026 20:25:05 +0300 Subject: [PATCH 10/10] riscv_fpu: keep the RMM preparation away from non-arithmetic ops Gating the RMM preparation on the effective mode let it fire for ops outside its add/sub/mul/div cases, where the default branch derived the rounding direction from rs1 reinterpreted as f64 and left that mode on the host with nothing restoring it: a static-rmm fcvt would poison the tracked mode for every later dynamic-rm op in the same frm block (F-fcvt.l.s/s.w/wu.s regressions). Make the default case a no-op -- sqrt has no exact ties to synthesize, and ops taking rm as an argument handle RMM natively in the integral rounder -- and have the static-rm enter always report a mode to restore, so a mode change made by the op itself (the RMM preparation) is undone by leave. Net effect beyond fixing the regression: fcvt and fsqrt under RMM (both static and dynamic) now round exactly instead of inheriting the directed hack; F goes 16->21 and D 45->51 on the ACT4 suite with no new failures. Signed-off-by: Sol Astrius Phoenix --- src/cpu/riscv_fpu.c | 5 +++-- src/cpu/riscv_fpu.h | 8 ++++---- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/cpu/riscv_fpu.c b/src/cpu/riscv_fpu.c index 6026e40a3..aa2daa341 100644 --- a/src/cpu/riscv_fpu.c +++ b/src/cpu/riscv_fpu.c @@ -81,8 +81,9 @@ static void riscv_prepare_rmm(rvvm_hart_t* vm, const uint32_t insn, const size_t neg = fpu_signbit64(riscv_view_d(vm, rs1)) != fpu_signbit64(riscv_view_d(vm, rs2)); break; default: - neg = fpu_signbit64(riscv_view_d(vm, rs1)); - break; + // Only add/sub/mul/div need the directed synthesis: sqrt has no exact + // ties, and ops taking rm as an argument handle RMM natively + return; } // Round to positive/negative infinity based on the result sign diff --git a/src/cpu/riscv_fpu.h b/src/cpu/riscv_fpu.h index 2756c6adf..a2379dd1e 100644 --- a/src/cpu/riscv_fpu.h +++ b/src/cpu/riscv_fpu.h @@ -41,11 +41,11 @@ static forceinline uint32_t riscv_fpu_host_rm(uint32_t rm) static forceinline uint32_t riscv_fpu_static_rm_enter(uint32_t rm) { if (unlikely(rm != 0x07)) { + // Always report a mode to restore: the op itself may change it further + // (the RMM preparation), and leave must undo that too const uint32_t prev = fpu_get_rounding_mode(); - if (riscv_fpu_host_rm(rm) != prev) { - fpu_set_rounding_mode(riscv_fpu_host_rm(rm)); - return prev; - } + fpu_set_rounding_mode(riscv_fpu_host_rm(rm)); + return prev; } return 0x07; }