libm/math/
fma.rs

Help
1/* SPDX-License-Identifier: MIT */
2/* origin: musl src/math/fma.c, fmaf.c Ported to generic Rust algorithm in 2025, TG. */
3
4use super::generic;
5use crate::support::Round;
6
7// Placeholder so we can have `fmaf16` in the `Float` trait.
8#[allow(unused)]
9#[cfg(f16_enabled)]
10#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
11pub(crate) fn fmaf16(_x: f16, _y: f16, _z: f16) -> f16 {
12    unimplemented!()
13}
14
15/// Floating multiply add (f32)
16///
17/// Computes `(x*y)+z`, rounded as one ternary operation (i.e. calculated with infinite precision).
18#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
19pub fn fmaf(x: f32, y: f32, z: f32) -> f32 {
20    select_implementation! {
21        name: fmaf,
22        use_arch: any(
23            all(target_arch = "aarch64", target_feature = "neon"),
24            target_feature = "sse2",
25        ),
26        args: x, y, z,
27    }
28
29    generic::fma_wide_round(x, y, z, Round::Nearest).val
30}
31
32/// Fused multiply add (f64)
33///
34/// Computes `(x*y)+z`, rounded as one ternary operation (i.e. calculated with infinite precision).
35#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
36pub fn fma(x: f64, y: f64, z: f64) -> f64 {
37    select_implementation! {
38        name: fma,
39        use_arch: any(
40            all(target_arch = "aarch64", target_feature = "neon"),
41            target_feature = "sse2",
42        ),
43        args: x, y, z,
44    }
45
46    generic::fma_round(x, y, z, Round::Nearest).val
47}
48
49/// Fused multiply add (f128)
50///
51/// Computes `(x*y)+z`, rounded as one ternary operation (i.e. calculated with infinite precision).
52#[cfg(f128_enabled)]
53#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
54pub fn fmaf128(x: f128, y: f128, z: f128) -> f128 {
55    generic::fma_round(x, y, z, Round::Nearest).val
56}
57
58#[cfg(test)]
59mod tests {
60    use super::*;
61    use crate::support::{CastFrom, CastInto, Float, FpResult, HInt, MinInt, Round, Status};
62
63    /// Test the generic `fma_round` algorithm for a given float.
64    fn spec_test<F>(f: impl Fn(F, F, F) -> F)
65    where
66        F: Float,
67        F: CastFrom<F::SignedInt>,
68        F: CastFrom<i8>,
69        F::Int: HInt,
70        u32: CastInto<F::Int>,
71    {
72        let x = F::from_bits(F::Int::ONE);
73        let y = F::from_bits(F::Int::ONE);
74        let z = F::ZERO;
75
76        // 754-2020 says "When the exact result of (a × b) + c is non-zero yet the result of
77        // fusedMultiplyAdd is zero because of rounding, the zero result takes the sign of the
78        // exact result"
79        assert_biteq!(f(x, y, z), F::ZERO);
80        assert_biteq!(f(x, -y, z), F::NEG_ZERO);
81        assert_biteq!(f(-x, y, z), F::NEG_ZERO);
82        assert_biteq!(f(-x, -y, z), F::ZERO);
83    }
84
85    #[test]
86    fn spec_test_f32() {
87        spec_test::<f32>(fmaf);
88
89        // Also do a small check that the non-widening version works for f32 (this should ideally
90        // get tested some more).
91        spec_test::<f32>(|x, y, z| generic::fma_round(x, y, z, Round::Nearest).val);
92    }
93
94    #[test]
95    fn spec_test_f64() {
96        spec_test::<f64>(fma);
97
98        let expect_underflow = [
99            (
100                hf64!("0x1.0p-1070"),
101                hf64!("0x1.0p-1070"),
102                hf64!("0x1.ffffffffffffp-1023"),
103                hf64!("0x0.ffffffffffff8p-1022"),
104            ),
105            (
106                // FIXME: we raise underflow but this should only be inexact (based on C and
107                // `rustc_apfloat`).
108                hf64!("0x1.0p-1070"),
109                hf64!("0x1.0p-1070"),
110                hf64!("-0x1.0p-1022"),
111                hf64!("-0x1.0p-1022"),
112            ),
113        ];
114
115        for (x, y, z, res) in expect_underflow {
116            let FpResult { val, status } = generic::fma_round(x, y, z, Round::Nearest);
117            assert_biteq!(val, res);
118            assert_eq!(status, Status::UNDERFLOW);
119        }
120    }
121
122    #[test]
123    #[cfg(f128_enabled)]
124    fn spec_test_f128() {
125        spec_test::<f128>(fmaf128);
126    }
127
128    #[test]
129    fn issue_263() {
130        let a = f32::from_bits(1266679807);
131        let b = f32::from_bits(1300234242);
132        let c = f32::from_bits(1115553792);
133        let expected = f32::from_bits(1501560833);
134        assert_eq!(fmaf(a, b, c), expected);
135    }
136
137    #[test]
138    fn fma_segfault() {
139        // These two inputs cause fma to segfault on release due to overflow:
140        assert_eq!(
141            fma(
142                -0.0000000000000002220446049250313,
143                -0.0000000000000002220446049250313,
144                -0.0000000000000002220446049250313
145            ),
146            -0.00000000000000022204460492503126,
147        );
148
149        let result = fma(-0.992, -0.992, -0.992);
150        //force rounding to storage format on x87 to prevent superious errors.
151        #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
152        let result = force_eval!(result);
153        assert_eq!(result, -0.007936000000000007,);
154    }
155
156    #[test]
157    fn fma_sbb() {
158        assert_eq!(
159            fma(-(1.0 - f64::EPSILON), f64::MIN, f64::MIN),
160            -3991680619069439e277
161        );
162    }
163
164    #[test]
165    fn fma_underflow() {
166        assert_eq!(
167            fma(1.1102230246251565e-16, -9.812526705433188e-305, 1.0894e-320),
168            0.0,
169        );
170    }
171}
libm/math/fma.rs

libm/math/
fma.rs
