malachite_base/num/arithmetic/
mul_add_mul.rs

1// Copyright © 2026 Mikhail Hogrefe
2//
3// This file is part of Malachite.
4//
5// Malachite is free software: you can redistribute it and/or modify it under the terms of the GNU
6// Lesser General Public License (LGPL) as published by the Free Software Foundation; either version
7// 3 of the License, or (at your option) any later version. See <https://www.gnu.org/licenses/>.
8
9use crate::num::arithmetic::traits::{MulAddMul, MulAddMulAssign, UnsignedAbs};
10use crate::num::basic::signeds::PrimitiveSigned;
11use crate::num::basic::unsigneds::PrimitiveUnsigned;
12use crate::num::conversion::traits::WrappingFrom;
13
14// Where the exact value of $xy \pm zw$ sits relative to the type's range. The products are formed
15// at double width, so no intermediate overflow can hide a result that would have fit.
16pub(crate) enum Wide<T> {
17    Fits(T),
18    Above,
19    Below,
20}
21
22// The exact value of $xy + zw$, or of $xy - zw$ when `sub` is set, for an unsigned type.
23pub(crate) fn mul_add_mul_wide_unsigned<T: PrimitiveUnsigned>(
24    x: T,
25    y: T,
26    z: T,
27    w: T,
28    sub: bool,
29) -> Wide<T> {
30    let (p_1, p_0) = T::x_mul_y_to_zz(x, y);
31    let (q_1, q_0) = T::x_mul_y_to_zz(z, w);
32    let (r_1, r_0) = if sub {
33        if (p_1, p_0) < (q_1, q_0) {
34            return Wide::Below;
35        }
36        T::xx_sub_yy_to_zz(p_1, p_0, q_1, q_0)
37    } else {
38        let (r_1, r_0) = T::xx_add_yy_to_zz(p_1, p_0, q_1, q_0);
39        // Two products can carry out of double width, since $(2^W-1)^2 + (2^W-1)^2 \geq 2^{2W}$,
40        // and `xx_add_yy_to_zz` wraps without saying so. The sum wrapped iff it is now smaller than
41        // one of the addends.
42        if (r_1, r_0) < (p_1, p_0) {
43            return Wide::Above;
44        }
45        (r_1, r_0)
46    };
47    if r_1 == T::ZERO {
48        Wide::Fits(r_0)
49    } else {
50        Wide::Above
51    }
52}
53
54// The exact value of $xy + zw$, or of $xy - zw$ when `sub` is set, for a signed type. The two
55// products are combined as a sign and a double-width magnitude, so the sign of an out-of-range
56// result is known and the saturating variants can pick the right bound.
57pub(crate) fn mul_add_mul_wide_signed<
58    U: PrimitiveUnsigned,
59    T: PrimitiveSigned + UnsignedAbs<Output = U> + WrappingFrom<U>,
60>(
61    x: T,
62    y: T,
63    z: T,
64    w: T,
65    sub: bool,
66) -> Wide<T> {
67    let p_neg = (x < T::ZERO) != (y < T::ZERO);
68    let q_neg = ((z < T::ZERO) != (w < T::ZERO)) != sub;
69    let (p_1, p_0) = U::x_mul_y_to_zz(x.unsigned_abs(), y.unsigned_abs());
70    let (q_1, q_0) = U::x_mul_y_to_zz(z.unsigned_abs(), w.unsigned_abs());
71    let (neg, r_1, r_0) = if p_neg == q_neg {
72        let (r_1, r_0) = U::xx_add_yy_to_zz(p_1, p_0, q_1, q_0);
73        (p_neg, r_1, r_0)
74    } else if (p_1, p_0) >= (q_1, q_0) {
75        let (r_1, r_0) = U::xx_sub_yy_to_zz(p_1, p_0, q_1, q_0);
76        (p_neg, r_1, r_0)
77    } else {
78        let (r_1, r_0) = U::xx_sub_yy_to_zz(q_1, q_0, p_1, p_0);
79        (q_neg, r_1, r_0)
80    };
81    if r_1 != U::ZERO {
82        return if neg { Wide::Below } else { Wide::Above };
83    }
84    if neg {
85        // The negative bound has one more magnitude than the positive one.
86        if r_0 <= T::MIN.unsigned_abs() {
87            Wide::Fits(T::wrapping_from(r_0).wrapping_neg())
88        } else {
89            Wide::Below
90        }
91    } else if r_0 <= T::MAX.unsigned_abs() {
92        Wide::Fits(T::wrapping_from(r_0))
93    } else {
94        Wide::Above
95    }
96}
97
98macro_rules! impl_mul_add_mul_primitive_int {
99    ($t:ident) => {
100        impl MulAddMul for $t {
101            type Output = $t;
102
103            /// Adds the products of two pairs of numbers.
104            ///
105            /// $f(x, y, z, w) = xy + zw$.
106            ///
107            /// Both products and their sum wrap on overflow, as they do for
108            /// [`add_mul`](super::traits::AddMul::add_mul).
109            ///
110            /// # Worst-case complexity
111            /// Constant time and additional memory.
112            ///
113            /// # Examples
114            /// See [here](super::mul_add_mul#mul_add_mul).
115            #[inline]
116            fn mul_add_mul(self, y: $t, z: $t, w: $t) -> $t {
117                self.wrapping_mul(y).wrapping_add(z.wrapping_mul(w))
118            }
119        }
120
121        impl MulAddMulAssign for $t {
122            /// Adds the products of two pairs of numbers, in place.
123            ///
124            /// $x \gets xy + zw$.
125            ///
126            /// Both products and their sum wrap on overflow, as they do for
127            /// [`add_mul`](super::traits::AddMul::add_mul).
128            ///
129            /// # Worst-case complexity
130            /// Constant time and additional memory.
131            ///
132            /// # Examples
133            /// See [here](super::mul_add_mul#mul_add_mul_assign).
134            #[inline]
135            fn mul_add_mul_assign(&mut self, y: $t, z: $t, w: $t) {
136                *self = self.wrapping_mul(y).wrapping_add(z.wrapping_mul(w));
137            }
138        }
139    };
140}
141apply_to_primitive_ints!(impl_mul_add_mul_primitive_int);
malachite_base/num/arithmetic/mul_add_mul.rs

malachite_base/num/arithmetic/
mul_add_mul.rs
