crypto_bigint/modular/safegcd/
boxed.rs

1//! Implementation of Bernstein-Yang modular inversion and GCD algorithm (a.k.a. safegcd)
2//! as described in: <https://eprint.iacr.org/2019/266>.
3//!
4//! See parent module for more information.
5
6use super::{GCD_BATCH_SIZE, Matrix, iterations, jump};
7use crate::{
8    BoxedUint, Choice, ConcatenatingMul, CtAssign, CtOption, CtSelect, I64, Int, Limb, NonZero,
9    Odd, Resize, U64, Uint,
10    primitives::{u32_max, u32_min},
11};
12use core::fmt;
13
14/// Modular multiplicative inverter based on the Bernstein-Yang method.
15///
16/// See [`super::SafeGcdInverter`] for more information.
17#[derive(Clone, Debug)]
18pub(crate) struct BoxedSafeGcdInverter {
19    /// Modulus
20    pub(crate) modulus: Odd<BoxedUint>,
21
22    /// Multiplicative inverse of the modulus modulo 2^62
23    inverse: u64,
24
25    /// Adjusting parameter (see toplevel documentation).
26    adjuster: BoxedUint,
27}
28
29impl BoxedSafeGcdInverter {
30    /// Creates the inverter for specified modulus and adjusting parameter.
31    ///
32    /// Modulus must be odd. Returns `None` if it is not.
33    #[cfg(test)]
34    pub fn new(modulus: Odd<BoxedUint>, adjuster: BoxedUint) -> Self {
35        let inverse = U64::from_u64(modulus.as_uint_ref().invert_mod_u64());
36        Self::new_with_inverse(modulus, inverse, adjuster)
37    }
38
39    /// Creates the inverter for specified modulus and adjusting parameter.
40    ///
41    /// Modulus must be odd. Returns `None` if it is not.
42    pub(crate) fn new_with_inverse(
43        modulus: Odd<BoxedUint>,
44        inverse: U64,
45        mut adjuster: BoxedUint,
46    ) -> Self {
47        adjuster = adjuster.resize(modulus.bits_precision());
48        Self {
49            modulus,
50            inverse: inverse.as_uint_ref().lowest_u64(),
51            adjuster,
52        }
53    }
54
55    /// Perform constant-time modular inversion.
56    pub(crate) fn invert(&self, value: &BoxedUint) -> CtOption<BoxedUint> {
57        invert_odd_mod_precomp::<false>(
58            value,
59            &self.modulus,
60            self.inverse,
61            Some(self.adjuster.clone()),
62        )
63    }
64
65    /// Perform variable-time modular inversion.
66    pub(crate) fn invert_vartime(&self, value: &BoxedUint) -> CtOption<BoxedUint> {
67        invert_odd_mod_precomp::<true>(
68            value,
69            &self.modulus,
70            self.inverse,
71            Some(self.adjuster.clone()),
72        )
73    }
74}
75
76#[inline]
77pub fn invert_odd_mod<const VARTIME: bool>(
78    a: &BoxedUint,
79    m: &Odd<BoxedUint>,
80) -> CtOption<BoxedUint> {
81    let mi = m.as_uint_ref().invert_mod_u64();
82    invert_odd_mod_precomp::<VARTIME>(a, m, mi, None)
83}
84
85/// Calculate the multiplicative inverse of `a` modulo `m`.
86///
87fn invert_odd_mod_precomp<const VARTIME: bool>(
88    a: &BoxedUint,
89    m: &Odd<BoxedUint>,
90    mi: u64,
91    e: Option<BoxedUint>,
92) -> CtOption<BoxedUint> {
93    let a_nonzero = a.is_nonzero();
94    let bits_precision = u32_max(a.bits_precision(), m.as_ref().bits_precision());
95    let m = m.as_ref().resize(bits_precision);
96    let (mut f, mut g) = (
97        SignedBoxedInt::from_uint(m.clone()),
98        SignedBoxedInt::from_uint_with_precision(a, bits_precision),
99    );
100    let (mut d, mut e) = (
101        SignedBoxedInt::zero_with_precision(bits_precision),
102        SignedBoxedInt::from_uint(e.map_or_else(
103            || BoxedUint::one_with_precision(bits_precision),
104            |e| e.resize(bits_precision),
105        )),
106    );
107    let mut steps = iterations(bits_precision);
108    let mut delta = 1;
109    let mut t;
110
111    while steps > 0 {
112        if VARTIME && g.is_zero_vartime() {
113            break;
114        }
115        let batch = u32_min(steps, GCD_BATCH_SIZE);
116        (delta, t) = jump::<VARTIME>(f.lowest(), g.lowest(), delta, batch);
117        (f, g) = update_fg(&f, &g, t, batch);
118        (d, e) = update_de(&d, &e, &m, mi, t, batch);
119        steps -= batch;
120    }
121
122    let d = d.norm(f.is_negative(), &m);
123    CtOption::new(d, f.magnitude().is_one() & a_nonzero)
124}
125
126/// Calculate the greatest common denominator of `f` and `g`.
127pub fn gcd<const VARTIME: bool>(f: &BoxedUint, g: &BoxedUint) -> BoxedUint {
128    if VARTIME {
129        if let Some(f_nz) = f.as_nz_vartime() {
130            gcd_nz::<VARTIME>(f_nz, g).get()
131        } else {
132            // gcd of (0, g) is g
133            g.resize(u32_max(f.bits_precision(), g.bits_precision()))
134        }
135    } else {
136        let (f_nz, f_is_nonzero) = f.to_nz_or_one();
137        let mut r = gcd_nz::<VARTIME>(&f_nz, g).get();
138        r.ct_assign(g, !f_is_nonzero);
139        r
140    }
141}
142
143/// Calculate the greatest common denominator of nonzero `f`, and `g`.
144pub fn gcd_nz<const VARTIME: bool>(f: &NonZero<BoxedUint>, g: &BoxedUint) -> NonZero<BoxedUint> {
145    if VARTIME {
146        if let Some(f_odd) = f.as_odd_vartime() {
147            return gcd_odd::<VARTIME>(f_odd, g).into_nz();
148        }
149    }
150
151    // Note the following two GCD identity rules:
152    // 1) gcd(2f, 2g) = 2•gcd(f, g), and
153    // 2) gcd(a, 2g) = gcd(f, g) if f is odd.
154    //
155    // Combined, these rules imply that
156    // 3) gcd(2^i•f, 2^j•g) = 2^k•gcd(f, g), with k = min(i, j).
157    //
158    // However, to save ourselves having to divide out 2^j, we also note that
159    // 4) 2^k•gcd(f, g) = 2^k•gcd(a, 2^j•b)
160
161    let i = f.as_ref().trailing_zeros();
162    let j = g
163        .is_nonzero()
164        .select_u32(f.bits_precision(), g.trailing_zeros());
165    let k = u32_min(i, j);
166
167    let f_odd = Odd::new_unchecked(f.as_ref().shr(i));
168    let mut r = gcd_odd::<VARTIME>(&f_odd, g).get();
169    r.shl_assign(k);
170    NonZero::new_unchecked(r)
171}
172
173/// Calculate the greatest common denominator of odd `f`, and `g`.
174pub fn gcd_odd<const VARTIME: bool>(f: &Odd<BoxedUint>, g: &BoxedUint) -> Odd<BoxedUint> {
175    let bits_precision = u32_max(f.as_ref().bits_precision(), g.bits_precision());
176    let (mut f, mut g) = (
177        SignedBoxedInt::from_uint_with_precision(f.as_ref(), bits_precision),
178        SignedBoxedInt::from_uint_with_precision(g, bits_precision),
179    );
180    let mut steps = iterations(bits_precision);
181    let mut delta = 1;
182    let mut t;
183
184    while steps > 0 {
185        if VARTIME && g.is_zero_vartime() {
186            break;
187        }
188        let batch = u32_min(steps, GCD_BATCH_SIZE);
189        (delta, t) = jump::<VARTIME>(f.lowest(), g.lowest(), delta, batch);
190        (f, g) = update_fg(&f, &g, t, batch);
191        steps -= batch;
192    }
193
194    f.magnitude()
195        .resize_unchecked(bits_precision)
196        .to_odd()
197        .expect("odd by construction")
198}
199
200#[inline]
201fn update_fg(
202    a: &SignedBoxedInt,
203    b: &SignedBoxedInt,
204    t: Matrix,
205    shift: u32,
206) -> (SignedBoxedInt, SignedBoxedInt) {
207    (
208        SignedBoxedInt::lincomb_int_reduce_shift(
209            a,
210            b,
211            &I64::from_i64(t[0][0]),
212            &I64::from_i64(t[0][1]),
213            shift,
214        ),
215        SignedBoxedInt::lincomb_int_reduce_shift(
216            a,
217            b,
218            &I64::from_i64(t[1][0]),
219            &I64::from_i64(t[1][1]),
220            shift,
221        ),
222    )
223}
224
225#[inline]
226fn update_de(
227    d: &SignedBoxedInt,
228    e: &SignedBoxedInt,
229    m: &BoxedUint,
230    mi: u64,
231    t: Matrix,
232    shift: u32,
233) -> (SignedBoxedInt, SignedBoxedInt) {
234    (
235        SignedBoxedInt::lincomb_int_reduce_shift_mod(
236            d,
237            e,
238            &Int::from_i64(t[0][0]),
239            &Int::from_i64(t[0][1]),
240            shift,
241            m,
242            U64::from_u64(mi),
243        ),
244        SignedBoxedInt::lincomb_int_reduce_shift_mod(
245            d,
246            e,
247            &Int::from_i64(t[1][0]),
248            &Int::from_i64(t[1][1]),
249            shift,
250            m,
251            U64::from_u64(mi),
252        ),
253    )
254}
255
256/// A `Uint` which carries a separate sign in order to maintain the same range.
257#[derive(Clone)]
258struct SignedBoxedInt {
259    sign: Choice,
260    magnitude: BoxedUint,
261}
262
263impl SignedBoxedInt {
264    pub fn zero_with_precision(bits_precision: u32) -> Self {
265        Self::from_uint(BoxedUint::zero_with_precision(bits_precision))
266    }
267
268    /// Construct a new `SignedInt` from a `Uint`.
269    pub const fn from_uint(uint: BoxedUint) -> Self {
270        Self {
271            sign: Choice::FALSE,
272            magnitude: uint,
273        }
274    }
275
276    /// Construct a new `SignedInt` from a `Uint`.
277    pub fn from_uint_with_precision(uint: &BoxedUint, bits_precision: u32) -> Self {
278        Self {
279            sign: Choice::FALSE,
280            magnitude: uint.resize(bits_precision),
281        }
282    }
283
284    /// Construct a new `SignedInt` from a `Uint` and a sign flag.
285    pub const fn from_uint_sign(magnitude: BoxedUint, sign: Choice) -> Self {
286        Self { sign, magnitude }
287    }
288
289    /// Obtain the magnitude of the `SignedInt`, i.e. its absolute value.
290    pub const fn magnitude(&self) -> &BoxedUint {
291        &self.magnitude
292    }
293
294    /// Determine if the `SignedInt` is non-zero.
295    pub fn is_nonzero(&self) -> Choice {
296        self.magnitude.is_nonzero()
297    }
298
299    /// Determine if the `SignedInt` is zero in variable time.
300    pub fn is_zero_vartime(&self) -> bool {
301        self.magnitude.is_zero_vartime()
302    }
303
304    /// Determine if the `SignedInt` is negative.
305    /// Note: `-0` is representable in this type, so it may be necessary
306    /// to check `self.is_nonzero()` as well.
307    pub const fn is_negative(&self) -> Choice {
308        self.sign
309    }
310
311    /// Extract the lowest 63 bits and convert to its signed representation.
312    #[allow(clippy::cast_possible_wrap)]
313    pub fn lowest(&self) -> i64 {
314        let mag = (self.magnitude.as_uint_ref().lowest_u64() & (u64::MAX >> 1)) as i64;
315        self.sign.select_i64(mag, mag.wrapping_neg())
316    }
317
318    /// Compute the linear combination `a•b + c•d`, returning a widened result.
319    #[inline]
320    pub(crate) fn lincomb_int<const RHS: usize>(
321        a: &Self,
322        b: &Self,
323        c: &Int<RHS>,
324        d: &Int<RHS>,
325    ) -> Self {
326        debug_assert!(a.magnitude.bits_precision() == b.magnitude.bits_precision());
327        let (c, c_sign) = c.abs_sign();
328        let (d, d_sign) = d.abs_sign();
329        // Each SignedBoxedInt • abs(Int) product leaves an empty upper bit.
330        let mut x = a.magnitude.concatenating_mul(c);
331        let x_neg = a.sign.xor(c_sign);
332        let mut y = b.magnitude.concatenating_mul(d);
333        let y_neg = b.sign.xor(d_sign);
334        let odd_neg = x_neg.xor(y_neg);
335
336        // Negate y if none or both of the multiplication results are negative.
337        y.conditional_wrapping_neg_assign(odd_neg.not());
338
339        let borrow;
340        (x, borrow) = x.borrowing_sub(&y, Limb::ZERO);
341        let swap = borrow.is_nonzero().and(odd_neg);
342
343        // Negate the result if we did not negate y and there was a borrow,
344        // indicating that |y| > |x|.
345        x.conditional_wrapping_neg_assign(swap);
346
347        let sign = x_neg.and(swap.not()).or(y_neg.and(swap));
348        Self::from_uint_sign(x, sign)
349    }
350
351    /// Compute the linear combination `a•b + c•d`, and shift the result
352    /// `shift` bits to the right, returning a signed value in the same range
353    /// as the `SignedInt` inputs.
354    pub(crate) fn lincomb_int_reduce_shift<const S: usize>(
355        a: &Self,
356        b: &Self,
357        c: &Int<S>,
358        d: &Int<S>,
359        shift: u32,
360    ) -> Self {
361        debug_assert!(shift < Uint::<S>::BITS);
362        let SignedBoxedInt {
363            sign,
364            mut magnitude,
365        } = Self::lincomb_int(a, b, c, d);
366        magnitude.shr_assign(shift);
367        Self::from_uint_sign(
368            magnitude.resize_unchecked(a.magnitude.bits_precision()),
369            sign,
370        )
371    }
372
373    /// Compute the linear combination `a•b + c•d`, and shift the result
374    /// `shift` bits to the right modulo `m`, returning a signed value in the
375    /// same range as the `SignedInt` inputs.
376    pub(crate) fn lincomb_int_reduce_shift_mod<const S: usize>(
377        a: &Self,
378        b: &Self,
379        c: &Int<S>,
380        d: &Int<S>,
381        shift: u32,
382        m: &BoxedUint,
383        mi: Uint<S>,
384    ) -> Self {
385        debug_assert!(shift < Uint::<S>::BITS);
386        let SignedBoxedInt {
387            sign: mut x_sign,
388            magnitude: mut x,
389        } = Self::lincomb_int(a, b, c, d);
390
391        // Compute the multiple of m that will clear the low N bits of x.
392        let mut xs = Uint::<S>::ZERO;
393        xs.limbs.copy_from_slice(&x.limbs[..S]);
394        let mut mf = xs.wrapping_mul(&mi);
395        mf = mf.bitand(&Uint::MAX.shr_vartime(Uint::<S>::BITS - shift));
396        let xa = m.concatenating_mul(mf);
397
398        // Subtract the adjustment from x potentially producing a borrow.
399        let borrow = x.borrowing_sub_assign(&xa, Limb::ZERO);
400
401        // Negate x if the subtraction borrowed.
402        let swap = borrow.is_nonzero();
403        x.conditional_wrapping_neg_assign(swap);
404        x_sign = x_sign.xor(swap);
405
406        // Shift the result, eliminating the trailing zeros.
407        x.shr_assign(shift);
408
409        // The magnitude x is now in the range [0, 2m). We conditionally subtract
410        // m in order to keep the output in (-m, m).
411        let x_hi = x.limbs[m.nlimbs()];
412        x = x.resize_unchecked(m.bits_precision());
413        x.sub_assign_mod_with_carry(x_hi, m, m);
414
415        Self::from_uint_sign(x, x_sign)
416    }
417
418    /// Normalize the value to a `BoxedUint` in the range `[0, m)`.
419    fn norm(&self, f_sign: Choice, m: &BoxedUint) -> BoxedUint {
420        let swap = f_sign.xor(self.sign) & self.is_nonzero();
421        BoxedUint::ct_select(&self.magnitude, &m.wrapping_sub(&self.magnitude), swap)
422    }
423}
424
425impl fmt::Debug for SignedBoxedInt {
426    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
427        f.write_fmt(format_args!(
428            "{}0x{}",
429            if self.sign.to_bool() { "-" } else { "+" },
430            self.magnitude
431        ))
432    }
433}
434
435#[cfg(test)]
436mod tests {
437    use super::{BoxedSafeGcdInverter, gcd_odd, invert_odd_mod};
438    use crate::{BoxedUint, Odd, Resize};
439    use num_bigint::BigUint;
440    use num_integer::Integer;
441
442    /// The smallest precision at which `45907 * bits + 30179` no longer fits in `u32`.
443    const LARGE_BITS: u32 = 93_568;
444
445    /// A deterministic, odd, full-width `BoxedUint` of `bits` bits; different seeds give
446    /// unrelated values.
447    fn large_odd(bits: u32, seed: u32) -> BoxedUint {
448        let mut bytes: alloc::vec::Vec<u8> = (0..bits >> 3)
449            .map(|i| {
450                let h = (i ^ seed.wrapping_mul(0x9E37_79B9)).wrapping_mul(0x85EB_CA6B);
451                let h = (h ^ (h >> 13)).wrapping_mul(0xC2B2_AE35);
452                (h ^ (h >> 16)).to_le_bytes()[0]
453            })
454            .collect();
455        bytes[0] |= 0x80;
456        *bytes.last_mut().unwrap() |= 1;
457        BoxedUint::from_be_slice(&bytes, bits).unwrap()
458    }
459
460    fn to_big(x: &BoxedUint) -> BigUint {
461        BigUint::from_bytes_be(&x.to_be_bytes())
462    }
463
464    #[test]
465    fn invert_odd_mod_large_precision() {
466        let m = large_odd(LARGE_BITS, 1);
467        let a = BoxedUint::from(65537u32).resize(LARGE_BITS);
468        assert_eq!(to_big(&a).gcd(&to_big(&m)), BigUint::from(1u32));
469        let m = Odd::new(m).unwrap();
470
471        for inv in [
472            invert_odd_mod::<false>(&a, &m).unwrap(),
473            invert_odd_mod::<true>(&a, &m).unwrap(),
474        ] {
475            assert_eq!(
476                (to_big(&a) * to_big(&inv)) % to_big(m.as_ref()),
477                BigUint::from(1u32)
478            );
479        }
480    }
481
482    #[test]
483    fn gcd_odd_large_precision() {
484        // Two unrelated full-width values: the reduction needs the whole iteration budget,
485        // unlike a pair with a large common factor and small cofactors.
486        let f = Odd::new(large_odd(LARGE_BITS, 2)).unwrap();
487        let g = large_odd(LARGE_BITS, 4);
488        let expected = to_big(f.as_ref()).gcd(&to_big(&g));
489
490        assert_eq!(to_big(gcd_odd::<false>(&f, &g).as_ref()), expected);
491        assert_eq!(to_big(gcd_odd::<true>(&f, &g).as_ref()), expected);
492    }
493
494    #[test]
495    fn invert() {
496        let g = BoxedUint::from_be_hex(
497            "00000000CBF9350842F498CE441FC2DC23C7BF47D3DE91C327B2157C5E4EED77",
498            256,
499        )
500        .unwrap();
501        let modulus = BoxedUint::from_be_hex(
502            "FFFFFFFF00000000FFFFFFFFFFFFFFFFBCE6FAADA7179E84F3B9CAC2FC632551",
503            256,
504        )
505        .unwrap()
506        .to_odd()
507        .unwrap();
508        let inverter = BoxedSafeGcdInverter::new(modulus, BoxedUint::one());
509        let result = inverter.invert(&g).unwrap();
510        assert_eq!(
511            BoxedUint::from_be_hex(
512                "FB668F8F509790BC549B077098918604283D42901C92981062EB48BC723F617B",
513                256
514            )
515            .unwrap(),
516            result
517        );
518    }
519}
crypto_bigint/modular/safegcd/boxed.rs

crypto_bigint/modular/safegcd/
boxed.rs
