crypto_bigint/modular/
mul.rs

1use super::reduction::montgomery_reduction;
2use crate::{CtLt, Limb, Odd, Uint, UintRef, WideWord, Word, word};
3
4/// Computes an "almost" Montgomery multiplication of `x` and `y` into `out`, that is
5/// `out = x * y * 2^(-n*W) mod m + am` assuming `k = -1/m mod 2^W`,
6/// where `W` is the bit size of the limb, and `n * W` is the full bit size of the integer.
7///
8/// NOTE: `out` is assumed to be pre-zeroized.
9///
10/// Unlike the standard Montgomery multiplication, we are reducing the final result only if
11/// it overflows `2^(n*W)`, not when it overflows `m`.
12///
13/// This means that this function does not assume `x` and `y` are reduced `mod m`,
14/// and the result will be correct `mod m`, but potentially greater than `m`,
15/// and smaller than `2^(n * W) + m`.
16/// See "Efficient Software Implementations of Modular Exponentiation" by S. Gueron for details
17/// (<https://eprint.iacr.org/2011/239.pdf>).
18///
19/// This function exhibits certain properties which were discovered via randomized tests,
20/// but (to my knowledge at this moment) have not been proven formally.
21/// Hereinafter we denote `f(x) = floor(x / m)`, that is `f` is the number of subtractions
22/// of the modulus required to fully reduce `x`.
23///
24/// 1. In general, if `f(x) = k` and `f(y) = n`, then `f(AMM(x, y)) <= min(k, n) + 1`.
25///    That is the "reduction error" grows with every operation,
26///    but is determined by the argument with the lower error.
27/// 2. To retrieve the number from Montgomery form we MM it by 1. In this case `f(AMM(x, 1)) = 0`,
28///    that is the result is always fully reduced regardless of `f(x)`.
29/// 3. `f(AMM(x, x)) <= 1` regardless of `f(x)`. That is, squaring resets the error to at most 1.
30#[inline]
31pub(crate) fn almost_montgomery_mul(
32    x: &UintRef,
33    y: &UintRef,
34    out: &mut UintRef,
35    modulus: &UintRef,
36    mod_neg_inv: Limb,
37) {
38    let overflow = montgomery_multiply_inner(
39        x.as_limbs(),
40        y.as_limbs(),
41        out.as_mut_limbs(),
42        modulus.as_limbs(),
43        mod_neg_inv,
44    );
45    let overflow = word::choice_from_lsb(overflow.0);
46    out.conditional_borrowing_sub_assign(modulus, overflow);
47}
48
49/// Ensure the output of an "almost" Montgomery multiplication is properly reduced.
50///
51/// Using the properties of `almost_montgomery_mul()` (see its documentation):
52/// - We have an incoming `x` which is fully reduced (`floor(x / modulus) = 0`).
53/// - We build an array of `powers` which are produced by multiplying the previous power by
54///   `x`, so for each power `floor(power / modulus) <= 1`.
55/// - Then we take turns squaring the accumulator `z` (bringing `floor(z / modulus)` to 1
56///   regardless of the previous reduction level) and multiplying by a power of `x`
57///   (bringing `floor(z / modulus)` to at most 2).
58/// - Then we either exit the loop, or square again, which brings `floor(z / modulus)` back
59///   to 1.
60///
61/// We now need to reduce `z` at most twice to bring it within `[0, modulus)`.
62pub(crate) fn almost_montgomery_reduce(z: &mut UintRef, modulus: &UintRef) {
63    z.conditional_borrowing_sub_assign(modulus, !z.ct_lt(modulus));
64    z.conditional_borrowing_sub_assign(modulus, !z.ct_lt(modulus));
65}
66
67/// Based on Algorithm 14.36 in Handbook of Applied Cryptography
68/// <https://cacr.uwaterloo.ca/hac/about/chap14.pdf>
69///
70/// Multiply `x` and `y` in Montgomery form, producing `x•y•R^-1 mod modulus + a•modulus`.
71///
72/// This algorithm roughly corresponds to the Finely Integrated Operand Scanning (FIOS)
73/// method of "Analyzing and Comparing Montgomery Multiplication Algorithms" by Koc et al
74/// <https://www.microsoft.com/en-us/research/wp-content/uploads/1996/01/j37acmon.pdf>
75/// but using wide words to track the intermediate products and carry.
76///
77/// The final conditional subtraction of the modulus to produce a result in the range
78/// `[0, modulus)` is not performed here, and must be performed by the caller. In some
79/// cases this may be deferred, as demonstrated by the `almost_montgomery_mul` method.
80#[inline(always)]
81#[allow(clippy::cast_possible_truncation)]
82pub const fn montgomery_multiply_inner(
83    x: &[Limb],
84    y: &[Limb],
85    out: &mut [Limb],
86    modulus: &[Limb],
87    mod_neg_inv: Limb,
88) -> Limb {
89    let nlimbs = modulus.len();
90    assert!(nlimbs == x.len() && nlimbs == y.len() && nlimbs == out.len());
91
92    let mut meta_carry = 0;
93
94    let mut i = 0;
95    while i < nlimbs {
96        let xi = x[i];
97        // Write B = 2^Word::BITS. Then axy <= (B - 1)^2 + (B - 1) = B^2 - B.
98        let axy = (xi.0 as WideWord) * (y[0].0 as WideWord) + out[0].0 as WideWord;
99        let u = (axy as Word).wrapping_mul(mod_neg_inv.0);
100
101        // Add the low halves separately: lo <= 2B - 2 fits in WideWord. This keeps
102        // the carry in arithmetic form without a double-width overflow flag.
103        let um = (u as WideWord) * (modulus[0].0 as WideWord);
104        let lo = (axy as Word as WideWord) + (um as Word as WideWord);
105        // floor((axy + um) / B) <= 2B - 2, including any carry from the low halves.
106        let mut carry = (axy >> Word::BITS) + (um >> Word::BITS) + (lo >> Word::BITS);
107
108        let mut j = 1;
109        while j < nlimbs {
110            // axy <= B^2 - B.
111            let axy = (xi.0 as WideWord) * (y[j].0 as WideWord) + out[j].0 as WideWord;
112            // umc <= (B - 1)^2 + (2B - 2) = B^2 - 1.
113            let umc = (u as WideWord) * (modulus[j].0 as WideWord) + carry;
114            let lo = (axy as Word as WideWord) + (umc as Word as WideWord);
115            out[j - 1] = Limb(lo as Word);
116            // floor((axy + umc) / B) <= floor((2B^2 - B - 1) / B) = 2B - 2.
117            carry = (axy >> Word::BITS) + (umc >> Word::BITS) + (lo >> Word::BITS);
118            j += 1;
119        }
120
121        // meta_carry <= 1, so this sum is at most 2B - 1 and still fits in WideWord.
122        carry += meta_carry;
123        (out[nlimbs - 1], meta_carry) = (Limb(carry as Word), carry >> Word::BITS);
124
125        i += 1;
126    }
127
128    Limb(meta_carry as Word)
129}
130
131/// Computes the Montgomery product of `a` and `b` modulo `modulus`, where
132/// `a` and `b` are in Montgomery form.
133pub(crate) const fn mul_montgomery_form<const LIMBS: usize>(
134    a: &Uint<LIMBS>,
135    b: &Uint<LIMBS>,
136    modulus: &Odd<Uint<LIMBS>>,
137    mod_neg_inv: Limb,
138) -> Uint<LIMBS> {
139    let mut out = Uint::<LIMBS>::ZERO;
140    let carry = montgomery_multiply_inner(
141        &a.limbs,
142        &b.limbs,
143        &mut out.limbs,
144        &modulus.as_ref().limbs,
145        mod_neg_inv,
146    );
147    out.try_sub_with_carry(carry, modulus.as_ref()).0
148}
149
150/// Computes the Montgomery squaring of `a` modulo `modulus` where
151/// `a` is in Montgomery form.
152pub(crate) const fn square_montgomery_form<const LIMBS: usize>(
153    a: &Uint<LIMBS>,
154    modulus: &Odd<Uint<LIMBS>>,
155    mod_neg_inv: Limb,
156) -> Uint<LIMBS> {
157    // One case in which it appears to be more efficient to perform a wide squaring and reduction.
158    if LIMBS == 4 {
159        let lower_upper = a.widening_square();
160        return montgomery_reduction(&lower_upper, modulus, mod_neg_inv);
161    }
162
163    let mut out = Uint::<LIMBS>::ZERO;
164    let carry = montgomery_multiply_inner(
165        &a.limbs,
166        &a.limbs,
167        &mut out.limbs,
168        &modulus.as_ref().limbs,
169        mod_neg_inv,
170    );
171    out.try_sub_with_carry(carry, modulus.as_ref()).0
172}
173
174/// Computes a repeated Montgomery squaring of `a` modulo `modulus` where
175/// `a` is in Montgomery form.
176///
177/// This method is variable time in `n`.
178///
179/// Assumes `a` is fully reduced, i.e. `a < modulus`.
180#[inline(always)]
181pub(crate) const fn square_repeat_montgomery_form<const LIMBS: usize>(
182    a: &Uint<LIMBS>,
183    n: u32,
184    modulus: &Odd<Uint<LIMBS>>,
185    mod_neg_inv: Limb,
186) -> Uint<LIMBS> {
187    if n == 0 {
188        return *a;
189    }
190    if n == 1 {
191        return square_montgomery_form(a, modulus, mod_neg_inv);
192    }
193
194    let mut i = 0;
195    let mut out = *a;
196    let mut base;
197
198    while i != n {
199        (base, out) = (out, Uint::ZERO);
200        let carry = montgomery_multiply_inner(
201            &base.limbs,
202            &base.limbs,
203            &mut out.limbs,
204            &modulus.as_ref().limbs,
205            mod_neg_inv,
206        );
207        // Because `base < modulus` here, the raw product is below `2 * modulus` and a single
208        // conditional subtraction suffices.
209        out = out.try_sub_with_carry(carry, modulus.as_ref()).0;
210        i += 1;
211    }
212
213    out
214}
215
216#[cfg(test)]
217mod tests {
218    use super::montgomery_multiply_inner;
219    use crate::{Limb, Uint, Word, modular::FixedMontyParams};
220    use num_bigint::BigUint;
221
222    fn to_biguint(limbs: &[Limb]) -> BigUint {
223        limbs.iter().rev().fold(BigUint::from(0u8), |value, limb| {
224            (value << Word::BITS) + BigUint::from(limb.0)
225        })
226    }
227
228    fn check_product<const N: usize>(x: [Limb; N], y: [Limb; N], modulus: [Limb; N]) {
229        let modulus = Uint::new(modulus).to_odd().unwrap();
230        let params = FixedMontyParams::new_vartime(modulus);
231        let m = to_biguint(modulus.as_ref().as_limbs());
232
233        // An unbounded REDC reference checks the complete result, including the
234        // final carry, without reproducing the limb-wise additions under test.
235        let mut expected = to_biguint(&x) * to_biguint(&y);
236        let mask = BigUint::from(Word::MAX);
237        let inverse = BigUint::from(params.mod_neg_inv().0);
238        for _ in 0..N {
239            let u = (&expected * &inverse) & &mask;
240            expected = (expected + u * &m) >> Word::BITS;
241        }
242
243        let mut out = [Limb::ZERO; N];
244        let carry = montgomery_multiply_inner(
245            &x,
246            &y,
247            &mut out,
248            modulus.as_ref().as_limbs(),
249            params.mod_neg_inv(),
250        );
251        let actual = to_biguint(&out) + (BigUint::from(carry.0) << (N * Limb::BITS as usize));
252        assert_eq!(actual, expected, "x={x:?}, y={y:?}, modulus={modulus:?}");
253    }
254
255    #[test]
256    fn multiply_carry_boundaries() {
257        fn check<const N: usize>() {
258            let values = [
259                [Limb::ZERO; N],
260                Uint::<N>::ONE.to_limbs(),
261                [Limb::MAX; N],
262                [Limb(Word::MAX - 1); N],
263                [Limb(1 << (Word::BITS - 1)); N],
264                core::array::from_fn(|i| if i & 1 == 0 { Limb::MAX } else { Limb::ZERO }),
265            ];
266            for mut modulus in values {
267                modulus[0].0 |= 1;
268                for x in values {
269                    for y in values {
270                        check_product(x, y, modulus);
271                    }
272                }
273            }
274        }
275        check::<1>();
276        check::<2>();
277        check::<3>();
278        check::<4>();
279        check::<16>();
280    }
281
282    #[test]
283    fn multiply_matches_biguint() {
284        use chacha20::ChaCha8Rng;
285        use rand_core::{Rng, SeedableRng};
286
287        let mut rng = ChaCha8Rng::from_seed([42; 32]);
288        let mut limb = || {
289            let mut bytes = [0u8; Limb::BYTES];
290            rng.fill_bytes(&mut bytes);
291            Limb(Word::from_le_bytes(bytes))
292        };
293        for _ in 0..1024 {
294            let x = core::array::from_fn(|_| limb());
295            let y = core::array::from_fn(|_| limb());
296            let mut modulus: [Limb; 4] = core::array::from_fn(|_| limb());
297            modulus[0].0 |= 1;
298            check_product(x, y, modulus);
299        }
300    }
301}
crypto_bigint/modular/mul.rs

crypto_bigint/modular/
mul.rs
