1/* SPDX-License-Identifier: MIT */
2/* origin: musl src/math/fma.c. Ported to generic Rust algorithm in 2025, TG. */
3
4use super::support::{DInt, FpResult, HInt, IntTy, Round, Status};
5use super::{CastFrom, CastInto, Float, Int, MinInt};
6
7/// Fused multiply add (f64)
8///
9/// Computes `(x*y)+z`, rounded as one ternary operation (i.e. calculated with infinite precision).
10#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
11pub fn fma(x: f64, y: f64, z: f64) -> f64 {
12 fma_round(x, y, z, Round::Nearest).val
13}
14
15/// Fused multiply add (f128)
16///
17/// Computes `(x*y)+z`, rounded as one ternary operation (i.e. calculated with infinite precision).
18#[cfg(f128_enabled)]
19#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
20pub fn fmaf128(x: f128, y: f128, z: f128) -> f128 {
21 fma_round(x, y, z, Round::Nearest).val
22}
23
24/// Fused multiply-add that works when there is not a larger float size available. Computes
25/// `(x * y) + z`.
26pub fn fma_round<F>(x: F, y: F, z: F, _round: Round) -> FpResult<F>
27where
28 F: Float,
29 F: CastFrom<F::SignedInt>,
30 F: CastFrom<i8>,
31 F::Int: HInt,
32 u32: CastInto<F::Int>,
33{
34 let one = IntTy::<F>::ONE;
35 let zero = IntTy::<F>::ZERO;
36
37 // Normalize such that the top of the mantissa is zero and we have a guard bit.
38 let nx = Norm::from_float(x);
39 let ny = Norm::from_float(y);
40 let nz = Norm::from_float(z);
41
42 if nx.is_zero_nan_inf() || ny.is_zero_nan_inf() {
43 // Value will overflow, defer to non-fused operations.
44 return FpResult::ok(x * y + z);
45 }
46
47 if nz.is_zero_nan_inf() {
48 if nz.is_zero() {
49 // Empty add component means we only need to multiply.
50 return FpResult::ok(x * y);
51 }
52 // `z` is NaN or infinity, which sets the result.
53 return FpResult::ok(z);
54 }
55
56 // multiply: r = x * y
57 let zhi: F::Int;
58 let zlo: F::Int;
59 let (mut rlo, mut rhi) = nx.m.widen_mul(ny.m).lo_hi();
60
61 // Exponent result of multiplication
62 let mut e: i32 = nx.e + ny.e;
63 // Needed shift to align `z` to the multiplication result
64 let mut d: i32 = nz.e - e;
65 let sbits = F::BITS as i32;
66
67 // Scale `z`. Shift `z <<= kz`, `r >>= kr`, so `kz+kr == d`, set `e = e+kr` (== ez-kz)
68 if d > 0 {
69 // The magnitude of `z` is larger than `x * y`
70 if d < sbits {
71 // Maximum shift of one `F::BITS` means shifted `z` will fit into `2 * F::BITS`. Shift
72 // it into `(zhi, zlo)`. No exponent adjustment necessary.
73 zlo = nz.m << d;
74 zhi = nz.m >> (sbits - d);
75 } else {
76 // Shift larger than `sbits`, `z` only needs the top half `zhi`. Place it there (acts
77 // as a shift by `sbits`).
78 zlo = zero;
79 zhi = nz.m;
80 d -= sbits;
81
82 // `z`'s exponent is large enough that it now needs to be taken into account.
83 e = nz.e - sbits;
84
85 if d == 0 {
86 // Exactly `sbits`, nothing to do
87 } else if d < sbits {
88 // Remaining shift fits within `sbits`. Leave `z` in place, shift `x * y`
89 rlo = (rhi << (sbits - d)) | (rlo >> d);
90 // Set the sticky bit
91 rlo |= IntTy::<F>::from((rlo << (sbits - d)) != zero);
92 rhi = rhi >> d;
93 } else {
94 // `z`'s magnitude is enough that `x * y` is irrelevant. It was nonzero, so set
95 // the sticky bit.
96 rlo = one;
97 rhi = zero;
98 }
99 }
100 } else {
101 // `z`'s magnitude once shifted fits entirely within `zlo`
102 zhi = zero;
103 d = -d;
104 if d == 0 {
105 // No shift needed
106 zlo = nz.m;
107 } else if d < sbits {
108 // Shift s.t. `nz.m` fits into `zlo`
109 let sticky = IntTy::<F>::from((nz.m << (sbits - d)) != zero);
110 zlo = (nz.m >> d) | sticky;
111 } else {
112 // Would be entirely shifted out, only set the sticky bit
113 zlo = one;
114 }
115 }
116
117 /* addition */
118
119 let mut neg = nx.neg ^ ny.neg;
120 let samesign: bool = !neg ^ nz.neg;
121 let mut rhi_nonzero = true;
122
123 if samesign {
124 // r += z
125 rlo = rlo.wrapping_add(zlo);
126 rhi += zhi + IntTy::<F>::from(rlo < zlo);
127 } else {
128 // r -= z
129 let (res, borrow) = rlo.overflowing_sub(zlo);
130 rlo = res;
131 rhi = rhi.wrapping_sub(zhi.wrapping_add(IntTy::<F>::from(borrow)));
132 if (rhi >> (F::BITS - 1)) != zero {
133 rlo = rlo.signed().wrapping_neg().unsigned();
134 rhi = rhi.signed().wrapping_neg().unsigned() - IntTy::<F>::from(rlo != zero);
135 neg = !neg;
136 }
137 rhi_nonzero = rhi != zero;
138 }
139
140 /* Construct result */
141
142 // Shift result into `rhi`, left-aligned. Last bit is sticky
143 if rhi_nonzero {
144 // `d` > 0, need to shift both `rhi` and `rlo` into result
145 e += sbits;
146 d = rhi.leading_zeros() as i32 - 1;
147 rhi = (rhi << d) | (rlo >> (sbits - d));
148 // Update sticky
149 rhi |= IntTy::<F>::from((rlo << d) != zero);
150 } else if rlo != zero {
151 // `rhi` is zero, `rlo` is the entire result and needs to be shifted
152 d = rlo.leading_zeros() as i32 - 1;
153 if d < 0 {
154 // Shift and set sticky
155 rhi = (rlo >> 1) | (rlo & one);
156 } else {
157 rhi = rlo << d;
158 }
159 } else {
160 // exact +/- 0.0
161 return FpResult::ok(x * y + z);
162 }
163
164 e -= d;
165
166 // Use int->float conversion to populate the significand.
167 // i is in [1 << (BITS - 2), (1 << (BITS - 1)) - 1]
168 let mut i: F::SignedInt = rhi.signed();
169
170 if neg {
171 i = -i;
172 }
173
174 // `|r|` is in `[0x1p62,0x1p63]` for `f64`
175 let mut r: F = F::cast_from_lossy(i);
176
177 /* Account for subnormal and rounding */
178
179 // Unbiased exponent for the maximum value of `r`
180 let max_pow = F::BITS - 1 + F::EXP_BIAS;
181
182 let mut status = Status::OK;
183
184 if e < -(max_pow as i32 - 2) {
185 // Result is subnormal before rounding
186 if e == -(max_pow as i32 - 1) {
187 let mut c = F::from_parts(false, max_pow, zero);
188 if neg {
189 c = -c;
190 }
191
192 if r == c {
193 // Min normal after rounding,
194 status.set_underflow(true);
195 r = F::MIN_POSITIVE_NORMAL.copysign(r);
196 return FpResult::new(r, status);
197 }
198
199 if (rhi << (F::SIG_BITS + 1)) != zero {
200 // Account for truncated bits. One bit will be lost in the `scalbn` call, add
201 // another top bit to avoid double rounding if inexact.
202 let iu: F::Int = (rhi >> 1) | (rhi & one) | (one << (F::BITS - 2));
203 i = iu.signed();
204
205 if neg {
206 i = -i;
207 }
208
209 r = F::cast_from_lossy(i);
210
211 // Remove the top bit
212 r = F::cast_from(2i8) * r - c;
213 status.set_underflow(true);
214 }
215 } else {
216 // Only round once when scaled
217 d = F::EXP_BITS as i32 - 1;
218 let sticky = IntTy::<F>::from(rhi << (F::BITS as i32 - d) != zero);
219 i = (((rhi >> d) | sticky) << d).signed();
220
221 if neg {
222 i = -i;
223 }
224
225 r = F::cast_from_lossy(i);
226 }
227 }
228
229 // Use our exponent to scale the final value.
230 FpResult::new(super::generic::scalbn(r, e), status)
231}
232
233/// Representation of `F` that has handled subnormals.
234#[derive(Clone, Copy, Debug)]
235struct Norm<F: Float> {
236 /// Normalized significand with one guard bit, unsigned.
237 m: F::Int,
238 /// Exponent of the mantissa such that `m * 2^e = x`. Accounts for the shift in the mantissa
239 /// and the guard bit; that is, 1.0 will normalize as `m = 1 << 53` and `e = -53`.
240 e: i32,
241 neg: bool,
242}
243
244impl<F: Float> Norm<F> {
245 /// Unbias the exponent and account for the mantissa's precision, including the guard bit.
246 const EXP_UNBIAS: u32 = F::EXP_BIAS + F::SIG_BITS + 1;
247
248 /// Values greater than this had a saturated exponent (infinity or NaN), OR were zero and we
249 /// adjusted the exponent such that it exceeds this threashold.
250 const ZERO_INF_NAN: u32 = F::EXP_SAT - Self::EXP_UNBIAS;
251
252 fn from_float(x: F) -> Self {
253 let mut ix = x.to_bits();
254 let mut e = x.ex() as i32;
255 let neg = x.is_sign_negative();
256 if e == 0 {
257 // Normalize subnormals by multiplication
258 let scale_i = F::BITS - 1;
259 let scale_f = F::from_parts(false, scale_i + F::EXP_BIAS, F::Int::ZERO);
260 let scaled = x * scale_f;
261 ix = scaled.to_bits();
262 e = scaled.ex() as i32;
263 e = if e == 0 {
264 // If the exponent is still zero, the input was zero. Artifically set this value
265 // such that the final `e` will exceed `ZERO_INF_NAN`.
266 1 << F::EXP_BITS
267 } else {
268 // Otherwise, account for the scaling we just did.
269 e - scale_i as i32
270 };
271 }
272
273 e -= Self::EXP_UNBIAS as i32;
274
275 // Absolute value, set the implicit bit, and shift to create a guard bit
276 ix &= F::SIG_MASK;
277 ix |= F::IMPLICIT_BIT;
278 ix <<= 1;
279
280 Self { m: ix, e, neg }
281 }
282
283 /// True if the value was zero, infinity, or NaN.
284 fn is_zero_nan_inf(self) -> bool {
285 self.e >= Self::ZERO_INF_NAN as i32
286 }
287
288 /// The only value we have
289 fn is_zero(self) -> bool {
290 // The only exponent that strictly exceeds this value is our sentinel value for zero.
291 self.e > Self::ZERO_INF_NAN as i32
292 }
293}
294
295#[cfg(test)]
296mod tests {
297 use super::*;
298
299 /// Test the generic `fma_round` algorithm for a given float.
300 fn spec_test<F>()
301 where
302 F: Float,
303 F: CastFrom<F::SignedInt>,
304 F: CastFrom<i8>,
305 F::Int: HInt,
306 u32: CastInto<F::Int>,
307 {
308 let x = F::from_bits(F::Int::ONE);
309 let y = F::from_bits(F::Int::ONE);
310 let z = F::ZERO;
311
312 let fma = |x, y, z| fma_round(x, y, z, Round::Nearest).val;
313
314 // 754-2020 says "When the exact result of (a × b) + c is non-zero yet the result of
315 // fusedMultiplyAdd is zero because of rounding, the zero result takes the sign of the
316 // exact result"
317 assert_biteq!(fma(x, y, z), F::ZERO);
318 assert_biteq!(fma(x, -y, z), F::NEG_ZERO);
319 assert_biteq!(fma(-x, y, z), F::NEG_ZERO);
320 assert_biteq!(fma(-x, -y, z), F::ZERO);
321 }
322
323 #[test]
324 fn spec_test_f32() {
325 spec_test::<f32>();
326 }
327
328 #[test]
329 fn spec_test_f64() {
330 spec_test::<f64>();
331
332 let expect_underflow = [
333 (
334 hf64!("0x1.0p-1070"),
335 hf64!("0x1.0p-1070"),
336 hf64!("0x1.ffffffffffffp-1023"),
337 hf64!("0x0.ffffffffffff8p-1022"),
338 ),
339 (
340 // FIXME: we raise underflow but this should only be inexact (based on C and
341 // `rustc_apfloat`).
342 hf64!("0x1.0p-1070"),
343 hf64!("0x1.0p-1070"),
344 hf64!("-0x1.0p-1022"),
345 hf64!("-0x1.0p-1022"),
346 ),
347 ];
348
349 for (x, y, z, res) in expect_underflow {
350 let FpResult { val, status } = fma_round(x, y, z, Round::Nearest);
351 assert_biteq!(val, res);
352 assert_eq!(status, Status::UNDERFLOW);
353 }
354 }
355
356 #[test]
357 #[cfg(f128_enabled)]
358 fn spec_test_f128() {
359 spec_test::<f128>();
360 }
361
362 #[test]
363 fn fma_segfault() {
364 // These two inputs cause fma to segfault on release due to overflow:
365 assert_eq!(
366 fma(
367 -0.0000000000000002220446049250313,
368 -0.0000000000000002220446049250313,
369 -0.0000000000000002220446049250313
370 ),
371 -0.00000000000000022204460492503126,
372 );
373
374 let result = fma(-0.992, -0.992, -0.992);
375 //force rounding to storage format on x87 to prevent superious errors.
376 #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
377 let result = force_eval!(result);
378 assert_eq!(result, -0.007936000000000007,);
379 }
380
381 #[test]
382 fn fma_sbb() {
383 assert_eq!(fma(-(1.0 - f64::EPSILON), f64::MIN, f64::MIN), -3991680619069439e277);
384 }
385
386 #[test]
387 fn fma_underflow() {
388 assert_eq!(fma(1.1102230246251565e-16, -9.812526705433188e-305, 1.0894e-320), 0.0,);
389 }
390}
391