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(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(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(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(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}
172