diff --git a/src/cpu/riscv_fpu.c b/src/cpu/riscv_fpu.c index 689a0af34..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 @@ -93,7 +94,33 @@ 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. "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 + 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 +129,9 @@ 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(((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); } @@ -406,4 +434,18 @@ 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 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); + } +} + #endif diff --git a/src/cpu/riscv_fpu.h b/src/cpu/riscv_fpu.h index 0ad5354fb..a2379dd1e 100644 --- a/src/cpu/riscv_fpu.h +++ b/src/cpu/riscv_fpu.h @@ -26,6 +26,37 @@ 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; +} + +/* + * 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)) { + // 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(); + 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) { @@ -135,20 +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, 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, 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); @@ -163,20 +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, 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, 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); @@ -191,20 +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, 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, 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); @@ -219,6 +265,8 @@ 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) @@ -226,14 +274,17 @@ static forceinline void riscv_emulate_f_fnmadd(rvvm_hart_t* vm, const uint32_t i 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, 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); diff --git a/src/util/fpu_lib.c b/src/util/fpu_lib.c index c82f0451e..9d39191bf 100644 --- a/src/util/fpu_lib.c +++ b/src/util/fpu_lib.c @@ -405,54 +405,118 @@ 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 + 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: + 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)); // +/-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 = (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: - 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; + // 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)); } 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)); // +/-1.0 or +/-0.0 + } + const uint64_t frac = u & (FPU_LIB_FP64_MANTISSA_MASK >> 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)); + 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) diff --git a/src/util/fpu_lib.h b/src/util/fpu_lib.h index e636436af..d877be1cf 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,56 @@ 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, 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); + 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 { + 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 +1077,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 +1091,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 +1103,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^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 + * 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) {