crypto_bigint/modular/bingcd/
xgcd.rs

1//! The Binary Extended GCD algorithm.
2use super::gcd::bingcd_step;
3use crate::modular::bingcd::matrix::{DividedIntMatrix, DividedPatternMatrix, PatternMatrix, Unit};
4use crate::primitives::u32_max;
5use crate::{Choice, Int, Limb, NonZeroUint, Odd, OddUint, U64, U128, Uint};
6
7/// Binary XGCD update step.
8///
9/// This is a condensed, constant time execution of the following algorithm:
10/// ```text
11/// if a mod 2 == 1
12///    if a < b
13///        (a, b) ← (b, a)
14///        (f0, g0, f1, g1) ← (f1, g1, f0, g0)
15///    a ← a - b
16///    (f0, g0) ← (f0 - f1, g0 - g1)
17/// if a > 0
18///     a ← a/2
19///     (f1, g1) ← (2f1, 2g1)
20/// ```
21/// where `matrix` represents
22/// ```text
23///  (f0 g0)
24///  (f1 g1).
25/// ```
26///
27/// Note: this algorithm assumes `b` to be an odd integer. The algorithm will likely not yield
28/// the correct result when this is not the case.
29///
30/// Ref: Pornin, Algorithm 2, L8-17, <https://eprint.iacr.org/2020/972.pdf>.
31#[inline(always)]
32const fn binxgcd_step<const LIMBS: usize, const MATRIX_LIMBS: usize>(
33    a: &mut Uint<LIMBS>,
34    b: &mut Uint<LIMBS>,
35    matrix: &mut DividedPatternMatrix<MATRIX_LIMBS>,
36    halt_at_zero: Choice,
37) {
38    let (a_odd, swap, _) = bingcd_step(a, b);
39
40    // swap if a odd and a < b
41    matrix.conditional_swap_rows(swap);
42
43    // subtract b from a when a is odd
44    matrix.conditional_subtract_bottom_row_from_top(a_odd);
45
46    // Double the bottom row of the matrix when a was ≠ 0 and when not halting.
47    matrix.conditional_double_bottom_row(a.is_nonzero().or(halt_at_zero.not()));
48}
49
50/// Container for the raw output of the Binary XGCD algorithm.
51pub(crate) struct RawXgcdOutput<const LIMBS: usize, MATRIX> {
52    gcd: OddUint<LIMBS>,
53    matrix: MATRIX,
54}
55
56pub(crate) type DividedPatternXgcdOutput<const LIMBS: usize> =
57    RawXgcdOutput<LIMBS, DividedPatternMatrix<LIMBS>>;
58
59impl<const LIMBS: usize> DividedPatternXgcdOutput<LIMBS> {
60    /// Divide `self.matrix.inner` by `2^self.matrix.k`, allowing us to simplify `inner` from a
61    /// [`DividedPatternMatrix`] to a [`PatternMatrix`].
62    ///
63    /// The performed divisions are modulo `lhs/gcd` and `rhs/gcd` to maintain the correctness of
64    /// the XGCD state.
65    ///
66    /// This operation is 'fast' since it only applies the division to the top row of the matrix.
67    /// This is allowed since it is assumed that `self.matrix * (lhs, rhs) = (gcd, 0)`; dividing
68    /// the bottom row of the matrix by a constant has no impact since its inner-product with the
69    /// input vector is zero.
70    ///
71    /// Executes in variable time w.r.t. `k_upper_bound`.
72    pub(crate) const fn divide(self) -> PatternXgcdOutput<LIMBS> {
73        let DividedPatternMatrix {
74            inner: mut matrix,
75            k,
76            k_upper_bound,
77            ..
78        } = self.matrix;
79
80        let PatternMatrix {
81            m00: x,
82            m01: y,
83            m10: rhs_div_gcd,
84            m11: lhs_div_gcd,
85            ..
86        } = &mut matrix;
87
88        if k_upper_bound > 0 {
89            *x = x.bounded_div2k_mod_q(
90                k,
91                k_upper_bound,
92                &rhs_div_gcd.to_odd().expect_copied("odd by construction"),
93            );
94            *y = y.bounded_div2k_mod_q(
95                k,
96                k_upper_bound,
97                &lhs_div_gcd.to_odd().expect_copied("odd by construction"),
98            );
99        }
100
101        PatternXgcdOutput {
102            gcd: self.gcd,
103            matrix,
104        }
105    }
106}
107
108pub(crate) type PatternXgcdOutput<const LIMBS: usize> = RawXgcdOutput<LIMBS, PatternMatrix<LIMBS>>;
109
110impl<const LIMBS: usize> PatternXgcdOutput<LIMBS> {
111    /// Obtain the `gcd`.
112    pub(crate) const fn gcd(&self) -> OddUint<LIMBS> {
113        self.gcd
114    }
115
116    /// Obtain the bezout coefficients `(x, y)` such that `lhs * x + rhs * y = gcd`.
117    pub(crate) const fn bezout_coefficients(&self) -> (Int<LIMBS>, Int<LIMBS>) {
118        let PatternMatrix {
119            m00,
120            m01,
121            m10,
122            m11,
123            pattern,
124            ..
125        } = self.matrix;
126
127        // TODO: can we simplify this?
128        let m10_sub_m00 = m10.wrapping_sub(&m00);
129        let m11_sub_m01 = m11.wrapping_sub(&m01);
130        let apply = Uint::lte(&m10_sub_m00, &m00).and(Uint::lte(&m11_sub_m01, &m01));
131
132        let m00 = *Uint::select(&m00, &m10_sub_m00, apply)
133            .wrapping_neg_if(apply.xor(pattern.not()))
134            .as_int();
135        let m01 = *Uint::select(&m01, &m11_sub_m01, apply)
136            .wrapping_neg_if(apply.xor(pattern))
137            .as_int();
138        (m00, m01)
139    }
140
141    /// Obtain the quotients `lhs/gcd` and `rhs/gcd` from `matrix`.
142    pub(crate) const fn quotients(&self) -> (Uint<LIMBS>, Uint<LIMBS>) {
143        let PatternMatrix {
144            m10: rhs_div_gcd,
145            m11: lhs_div_gcd,
146            ..
147        } = self.matrix;
148        (lhs_div_gcd, rhs_div_gcd)
149    }
150}
151
152/// Number of bits used by [`OddUint::optimized_binxgcd`] to represent a "compact" [`Uint`].
153const SUMMARY_BITS: u32 = U64::BITS - 1;
154
155/// Number of limbs used to represent [`SUMMARY_BITS`].
156const SUMMARY_LIMBS: usize = U64::LIMBS;
157
158/// Twice the number of limbs used to represent [`SUMMARY_BITS`], i.e., two times [`SUMMARY_LIMBS`].
159const DOUBLE_SUMMARY_LIMBS: usize = U128::LIMBS;
160
161impl<const LIMBS: usize> OddUint<LIMBS> {
162    /// The minimal number of binary GCD iterations required to guarantee successful completion.
163    const MIN_BINXGCD_ITERATIONS: u32 = 2 * Self::BITS - 1;
164
165    /// Given `(self, rhs)`, computes `(g, x, y)` s.t. `self * x + rhs * y = g = gcd(self, rhs)`,
166    /// leveraging the Binary Extended GCD algorithm.
167    pub(crate) const fn binxgcd_nz(&self, rhs: &NonZeroUint<LIMBS>) -> PatternXgcdOutput<LIMBS> {
168        let (lhs_, rhs_) = (self.as_ref(), rhs.as_ref());
169
170        // The `xgcd` subroutine requires `rhs` to be odd.
171        // We leverage the equality gcd(lhs, rhs) = gcd(lhs, |lhs-rhs|) to deal with the case that
172        // `rhs` is even.
173        let rhs_is_even = rhs_.is_odd().not();
174        let (abs_diff, rhs_gt_lhs) = lhs_.abs_diff(rhs_);
175        let odd_rhs = Odd::new_unchecked(Uint::select(rhs_, &abs_diff, rhs_is_even));
176
177        let mut output = self.binxgcd_odd(&odd_rhs);
178        let matrix = &mut output.matrix;
179
180        // Modify the output to negate the transformation applied to the input.
181        let case_one = rhs_is_even.and(rhs_gt_lhs);
182        matrix.conditional_subtract_right_column_from_left(case_one);
183
184        let case_two = rhs_is_even.and(rhs_gt_lhs.not());
185        matrix.conditional_add_right_column_to_left(case_two);
186        matrix.conditional_negate(case_two);
187
188        output
189    }
190
191    /// Execute the classic Extended GCD algorithm.
192    ///
193    /// Given `(self, rhs)`, computes `(g, x, y)` s.t. `self * x + rhs * y = g = gcd(self, rhs)`.
194    #[inline]
195    pub(crate) const fn binxgcd_odd(&self, rhs: &Self) -> PatternXgcdOutput<LIMBS> {
196        if LIMBS < 4 {
197            self.classic_binxgcd(rhs).divide()
198        } else {
199            self.optimized_binxgcd(rhs).divide()
200        }
201    }
202
203    /// Execute the classic Binary Extended GCD algorithm.
204    ///
205    /// Given `(self, rhs)`, computes `(g, x, y)` s.t. `self * x + rhs * y = g = gcd(self, rhs)`.
206    ///
207    /// Ref: Pornin, Optimized Binary GCD for Modular Inversion, Algorithm 1.
208    /// <https://eprint.iacr.org/2020/972.pdf>.
209    pub(crate) const fn classic_binxgcd(&self, rhs: &Self) -> DividedPatternXgcdOutput<LIMBS> {
210        let (gcd, _, matrix) =
211            self.partial_binxgcd::<LIMBS>(rhs.as_ref(), Self::MIN_BINXGCD_ITERATIONS, Choice::TRUE);
212        DividedPatternXgcdOutput { gcd, matrix }
213    }
214
215    /// Given `(self, rhs)`, computes `(g, x, y)` s.t. `self * x + rhs * y = g = gcd(self, rhs)`,
216    /// leveraging the Binary Extended GCD algorithm.
217    ///
218    /// **Warning**: `self` and `rhs` must be contained in an [U128] or larger.
219    ///
220    /// Note: this algorithm becomes more efficient than the classical algorithm for [Uint]s with
221    /// relatively many `LIMBS`. A best-effort threshold is presented in [`Self::binxgcd`_].
222    ///
223    /// Note: the full algorithm has an additional parameter; this function selects the best-effort
224    /// value for this parameter. You might be able to further tune your performance by calling the
225    /// [`Self::optimized_bingcd`_] function directly.
226    ///
227    /// Ref: Pornin, Optimized Binary GCD for Modular Inversion, Algorithm 2.
228    /// <https://eprint.iacr.org/2020/972.pdf>.
229    pub(crate) const fn optimized_binxgcd(&self, rhs: &Self) -> DividedPatternXgcdOutput<LIMBS> {
230        assert!(Self::BITS >= U128::BITS);
231        self.optimized_binxgcd_::<SUMMARY_BITS, SUMMARY_LIMBS, DOUBLE_SUMMARY_LIMBS>(rhs)
232    }
233
234    /// Given `(self, rhs)`, computes `(g, x, y)`, s.t. `self * x + rhs * y = g = gcd(self, rhs)`,
235    /// leveraging the optimized Binary Extended GCD algorithm.
236    ///
237    /// Ref: Pornin, Optimized Binary GCD for Modular Inversion, Algorithm 2.
238    /// <https://eprint.iacr.org/2020/972.pdf>
239    ///
240    /// In summary, the optimized algorithm does not operate on `self` and `rhs` directly, but
241    /// instead of condensed summaries that fit in few registers. Based on these summaries, an
242    /// update matrix is constructed by which `self` and `rhs` are updated in larger steps.
243    ///
244    /// This function is generic over the following three values:
245    /// - `K`: the number of bits used when summarizing `self` and `rhs` for the inner loop. The
246    ///   `K+1` top bits and `K-1` least significant bits are selected. It is recommended to keep
247    ///   `K` close to a (multiple of) the number of bits that fit in a single register.
248    /// - `LIMBS_K`: should be chosen as the minimum number s.t. `Uint::<LIMBS>::BITS ≥ K`,
249    /// - `LIMBS_2K`: should be chosen as the minimum number s.t. `Uint::<LIMBS>::BITS ≥ 2K`.
250    pub(crate) const fn optimized_binxgcd_<
251        const K: u32,
252        const LIMBS_K: usize,
253        const LIMBS_2K: usize,
254    >(
255        &self,
256        rhs: &Self,
257    ) -> DividedPatternXgcdOutput<LIMBS> {
258        let (mut a, mut b) = (*self.as_ref(), *rhs.as_ref());
259        let mut state = DividedIntMatrix::UNIT;
260
261        let (mut a_is_negative, mut b_is_negative);
262        let mut i = 0;
263        while i < Self::MIN_BINXGCD_ITERATIONS.div_ceil(K - 1) {
264            // Loop invariants:
265            //  i) each iteration of this loop, `a.bits() + b.bits()` shrinks by at least K-1,
266            //     until `b = 0`.
267            // ii) `a` is odd.
268            i += 1;
269
270            // Construct compact_a and compact_b as the summary of a and b, respectively.
271            let b_bits = b.bits();
272            let n = u32_max(2 * K, u32_max(a.bits(), b_bits));
273            let compact_a = a.compact::<K, LIMBS_2K>(n);
274            let compact_b = b.compact::<K, LIMBS_2K>(n);
275            let b_eq_compact_b =
276                Choice::from_u32_le(b_bits, K - 1).or(Choice::from_u32_eq(n, 2 * K));
277
278            // Compute the K-1 iteration update matrix from a_ and b_
279            let (.., update_matrix) = compact_a
280                .to_odd()
281                .expect_copied("a is always odd")
282                .partial_binxgcd::<LIMBS_K>(&compact_b, K - 1, b_eq_compact_b);
283
284            // Update `a` and `b` using the update matrix
285            let (updated_a, updated_b) = update_matrix.extended_apply_to::<LIMBS, K>((a, b));
286            (a, a_is_negative) = updated_a.dropped_abs_sign();
287            (b, b_is_negative) = updated_b.dropped_abs_sign();
288
289            state = update_matrix.mul_int_matrix(&state);
290
291            // Correct for the sign change in a/b after dropping the extension.
292            state.conditional_negate_top_row(a_is_negative);
293            state.conditional_negate_bottom_row(b_is_negative);
294        }
295
296        let gcd = a
297            .to_odd()
298            .expect_copied("gcd of an odd value with something else is always odd");
299
300        let matrix = state.to_divided_pattern_matrix();
301        DividedPatternXgcdOutput { gcd, matrix }
302    }
303
304    /// Executes the optimized Binary GCD inner loop.
305    ///
306    /// Ref: Pornin, Optimized Binary GCD for Modular Inversion, Algorithm 2.
307    /// <https://eprint.iacr.org/2020/972.pdf>.
308    ///
309    /// The function outputs the reduced values `(a, b)` for the input values `(self, rhs)` as well
310    /// as the matrix that yields the former two when multiplied with the latter two.
311    ///
312    /// Note: this implementation deviates slightly from the paper, in that it can be instructed to
313    /// "run in place" (i.e., execute iterations that do nothing) once `a` becomes zero.
314    /// This is done by passing a truthy `halt_at_zero`.
315    ///
316    /// The function executes in time variable in `iterations`.
317    #[inline(always)]
318    pub(super) const fn partial_binxgcd<const UPDATE_LIMBS: usize>(
319        &self,
320        rhs: &Uint<LIMBS>,
321        iterations: u32,
322        halt_at_zero: Choice,
323    ) -> (Self, Uint<LIMBS>, DividedPatternMatrix<UPDATE_LIMBS>) {
324        let (mut a, mut b) = (*self.as_ref(), *rhs);
325        // This matrix corresponds with (f0, g0, f1, g1) in the paper.
326        let mut matrix = DividedPatternMatrix::UNIT;
327
328        // Compute the update matrix.
329        // Note: to be consistent with the paper, the `binxgcd_step` algorithm requires the second
330        // argument to be odd. Here, we have `a` odd, so we have to swap a and b before and after
331        // calling the subroutine. The columns of the matrix have to be swapped accordingly.
332        Uint::swap(&mut a, &mut b);
333        matrix.swap_rows();
334
335        let mut i = 0;
336
337        while i < iterations {
338            binxgcd_step::<LIMBS, UPDATE_LIMBS>(&mut a, &mut b, &mut matrix, halt_at_zero);
339            i += 1;
340        }
341
342        // Undo swap
343        Uint::swap(&mut a, &mut b);
344        matrix.swap_rows();
345
346        let a = a.to_odd().expect_copied("a is always odd");
347        (a, b, matrix)
348    }
349}
350
351impl<const LIMBS: usize> Uint<LIMBS> {
352    /// Compute the absolute difference between `self` and `rhs`.
353    /// In addition to the result, also returns whether `rhs > self`.
354    const fn abs_diff(&self, rhs: &Self) -> (Self, Choice) {
355        let (diff, borrow) = self.borrowing_sub(rhs, Limb::ZERO);
356        let rhs_gt_self = borrow.is_nonzero();
357        let abs_diff = diff.wrapping_neg_if(rhs_gt_self);
358        (abs_diff, rhs_gt_self)
359    }
360}
361
362#[cfg(all(test, not(miri)))]
363mod tests {
364    use crate::modular::bingcd::xgcd::PatternXgcdOutput;
365    use crate::{Concat, Uint};
366    use core::ops::Div;
367
368    mod test_extract_quotients {
369        use crate::modular::bingcd::matrix::DividedPatternMatrix;
370        use crate::modular::bingcd::xgcd::{DividedPatternXgcdOutput, RawXgcdOutput};
371        use crate::{Choice, U64, Uint};
372
373        fn raw_binxgcdoutput_setup<const LIMBS: usize>(
374            matrix: DividedPatternMatrix<LIMBS>,
375        ) -> DividedPatternXgcdOutput<LIMBS> {
376            RawXgcdOutput {
377                gcd: Uint::<LIMBS>::ONE.to_odd().unwrap(),
378                matrix,
379            }
380        }
381
382        #[test]
383        fn test_extract_quotients_unit() {
384            let output =
385                raw_binxgcdoutput_setup(DividedPatternMatrix::<{ U64::LIMBS }>::UNIT).divide();
386            let (lhs_on_gcd, rhs_on_gcd) = output.quotients();
387            assert_eq!(lhs_on_gcd, Uint::ONE);
388            assert_eq!(rhs_on_gcd, Uint::ZERO);
389        }
390
391        #[test]
392        fn test_extract_quotients_basic() {
393            let output = raw_binxgcdoutput_setup(DividedPatternMatrix::<{ U64::LIMBS }>::new_u64(
394                (0, 0, 5, 7),
395                Choice::FALSE,
396                0,
397                0,
398            ))
399            .divide();
400            let (lhs_on_gcd, rhs_on_gcd) = output.quotients();
401            assert_eq!(lhs_on_gcd, Uint::from(7u32));
402            assert_eq!(rhs_on_gcd, Uint::from(5u32));
403
404            let output = raw_binxgcdoutput_setup(DividedPatternMatrix::<{ U64::LIMBS }>::new_u64(
405                (0, 0, 7u64, 5u64),
406                Choice::TRUE,
407                0,
408                0,
409            ))
410            .divide();
411            let (lhs_on_gcd, rhs_on_gcd) = output.quotients();
412            assert_eq!(lhs_on_gcd, Uint::from(5u32));
413            assert_eq!(rhs_on_gcd, Uint::from(7u32));
414        }
415    }
416
417    mod test_derive_bezout_coefficients {
418        use crate::modular::bingcd::matrix::DividedPatternMatrix;
419        use crate::modular::bingcd::xgcd::RawXgcdOutput;
420        use crate::{Choice, Int, U64, Uint};
421
422        #[test]
423        fn test_derive_bezout_coefficients_unit() {
424            let output = RawXgcdOutput {
425                gcd: Uint::ONE.to_odd().unwrap(),
426                matrix: DividedPatternMatrix::<{ U64::LIMBS }>::UNIT,
427            }
428            .divide();
429            let (x, y) = output.bezout_coefficients();
430            assert_eq!(x, Int::ONE);
431            assert_eq!(y, Int::ZERO);
432        }
433
434        #[test]
435        fn test_derive_bezout_coefficients_basic() {
436            let output = RawXgcdOutput {
437                gcd: U64::ONE.to_odd().unwrap(),
438                matrix: DividedPatternMatrix::new_u64((2u64, 3u64, 5u64, 5u64), Choice::TRUE, 0, 0),
439            }
440            .divide();
441            let (x, y) = output.bezout_coefficients();
442            assert_eq!(x, Int::from(2i32));
443            assert_eq!(y, Int::from(-3i32));
444
445            let output = RawXgcdOutput {
446                gcd: U64::ONE.to_odd().unwrap(),
447                matrix: DividedPatternMatrix::new_u64(
448                    (2u64, 3u64, 3u64, 5u64),
449                    Choice::FALSE,
450                    0,
451                    1,
452                ),
453            }
454            .divide();
455            let (x, y) = output.bezout_coefficients();
456            assert_eq!(x, Int::from(1i32));
457            assert_eq!(y, Int::from(-2i32));
458        }
459
460        #[test]
461        fn test_derive_bezout_coefficients_removes_doublings_easy() {
462            let output = RawXgcdOutput {
463                gcd: U64::ONE.to_odd().unwrap(),
464                matrix: DividedPatternMatrix::new_u64((2u64, 6u64, 3u64, 5u64), Choice::TRUE, 1, 1),
465            }
466            .divide();
467            let (x, y) = output.bezout_coefficients();
468            assert_eq!(x, Int::ONE);
469            assert_eq!(y, Int::from(-3i32));
470
471            let output = RawXgcdOutput {
472                gcd: U64::ONE.to_odd().unwrap(),
473                matrix: DividedPatternMatrix::new_u64(
474                    (120u64, 64u64, 7u64, 5u64),
475                    Choice::FALSE,
476                    5,
477                    6,
478                ),
479            }
480            .divide();
481            let (x, y) = output.bezout_coefficients();
482            assert_eq!(x, Int::from(-9i32));
483            assert_eq!(y, Int::from(2i32));
484        }
485
486        #[test]
487        fn test_derive_bezout_coefficients_removes_doublings_for_odd_numbers() {
488            let output = RawXgcdOutput {
489                gcd: U64::ONE.to_odd().unwrap(),
490                matrix: DividedPatternMatrix::new_u64(
491                    (2u64, 6u64, 7u64, 5u64),
492                    Choice::FALSE,
493                    3,
494                    7,
495                ),
496            }
497            .divide();
498            let (x, y) = output.bezout_coefficients();
499            assert_eq!(x, Int::from(-2i32));
500            assert_eq!(y, Int::from(2i32));
501        }
502    }
503
504    mod test_partial_binxgcd {
505        use crate::modular::bingcd::matrix::DividedPatternMatrix;
506        use crate::{Choice, Gcd, Odd, U64};
507
508        const A: Odd<U64> = Odd::from_be_hex("CA048AFA63CD6A1F");
509        const B: U64 = U64::from_be_hex("AE693BF7BE8E5566");
510
511        #[test]
512        fn test_partial_binxgcd() {
513            let (.., matrix) = A.partial_binxgcd::<{ U64::LIMBS }>(&B, 5, Choice::TRUE);
514            assert_eq!(matrix.k, 5);
515            assert_eq!(
516                matrix,
517                DividedPatternMatrix::new_u64((8u64, 4u64, 2u64, 5u64), Choice::TRUE, 5, 5)
518            );
519        }
520
521        #[test]
522        fn test_partial_binxgcd_constructs_correct_matrix() {
523            let target_a = U64::from_be_hex("1CB3FB3FA1218FDB").to_odd().unwrap();
524            let target_b = U64::from_be_hex("0EA028AF0F8966B6");
525
526            let (new_a, new_b, matrix) = A.partial_binxgcd::<{ U64::LIMBS }>(&B, 5, Choice::TRUE);
527
528            assert_eq!(new_a, target_a);
529            assert_eq!(new_b, target_b);
530
531            let (computed_a, computed_b) =
532                matrix.extended_apply_to::<{ U64::LIMBS }, 6>((A.get(), B));
533            let computed_a = computed_a.dropped_abs_sign().0;
534            let computed_b = computed_b.dropped_abs_sign().0;
535
536            assert_eq!(computed_a, target_a);
537            assert_eq!(computed_b, target_b);
538        }
539
540        const SMALL_A: Odd<U64> = Odd::from_be_hex("0000000003CD6A1F");
541        const SMALL_B: U64 = U64::from_be_hex("000000000E8E5566");
542
543        #[test]
544        fn test_partial_binxgcd_halts() {
545            let (gcd, _, matrix) =
546                SMALL_A.partial_binxgcd::<{ U64::LIMBS }>(&SMALL_B, 60, Choice::TRUE);
547            assert_eq!(matrix.k, 35);
548            assert_eq!(matrix.k_upper_bound, 60);
549            assert_eq!(gcd.get(), SMALL_A.gcd(&SMALL_B));
550        }
551
552        #[test]
553        fn test_partial_binxgcd_does_not_halt() {
554            let (gcd, .., matrix) =
555                SMALL_A.partial_binxgcd::<{ U64::LIMBS }>(&SMALL_B, 60, Choice::FALSE);
556            assert_eq!(matrix.k, 60);
557            assert_eq!(matrix.k_upper_bound, 60);
558            assert_eq!(gcd.get(), SMALL_A.gcd(&SMALL_B));
559        }
560    }
561
562    /// Helper function to effectively test xgcd.
563    fn test_xgcd<const LIMBS: usize, const DOUBLE: usize>(
564        lhs: Uint<LIMBS>,
565        rhs: Uint<LIMBS>,
566        output: PatternXgcdOutput<LIMBS>,
567    ) where
568        Uint<LIMBS>: Concat<LIMBS, Output = Uint<DOUBLE>>,
569    {
570        // Test the gcd
571        assert_eq!(lhs.gcd(&rhs), output.gcd, "{lhs} {rhs}");
572
573        // Test the quotients
574        let (lhs_on_gcd, rhs_on_gcd) = output.quotients();
575        assert_eq!(lhs_on_gcd, lhs.div(output.gcd.as_nz_ref()));
576        assert_eq!(rhs_on_gcd, rhs.div(output.gcd.as_nz_ref()));
577
578        // Test the Bezout coefficients for correctness
579        let (x, y) = output.bezout_coefficients();
580        assert_eq!(
581            x.concatenating_mul_unsigned(&lhs) + y.concatenating_mul_unsigned(&rhs),
582            *output.gcd.resize().as_int(),
583        );
584
585        // Test the Bezout coefficients for minimality
586        assert!(x.abs() <= rhs.div(output.gcd.as_nz_ref()));
587        assert!(y.abs() <= lhs.div(output.gcd.as_nz_ref()));
588        if lhs != rhs {
589            assert!(x.abs() <= rhs_on_gcd.shr(1) || rhs_on_gcd.is_zero().to_bool());
590            assert!(y.abs() <= lhs_on_gcd.shr(1) || lhs_on_gcd.is_zero().to_bool());
591        }
592    }
593
594    mod test_binxgcd_nz {
595        use crate::modular::bingcd::xgcd::tests::test_xgcd;
596        use crate::{
597            Concat, Int, U64, U128, U192, U256, U384, U512, U768, U1024, U2048, U4096, U8192, Uint,
598        };
599
600        fn binxgcd_nz_test<const LIMBS: usize, const DOUBLE: usize>(
601            lhs: Uint<LIMBS>,
602            rhs: Uint<LIMBS>,
603        ) where
604            Uint<LIMBS>: Concat<LIMBS, Output = Uint<DOUBLE>>,
605        {
606            let output = lhs.to_odd().unwrap().binxgcd_nz(&rhs.to_nz().unwrap());
607            test_xgcd(lhs, rhs, output);
608        }
609
610        fn binxgcd_nz_tests<const LIMBS: usize, const DOUBLE: usize>()
611        where
612            Uint<LIMBS>: Concat<LIMBS, Output = Uint<DOUBLE>>,
613        {
614            let max_int = *Int::MAX.as_uint();
615            let int_abs_min = Int::MIN.abs();
616
617            binxgcd_nz_test(Uint::ONE, Uint::ONE);
618            binxgcd_nz_test(Uint::ONE, max_int);
619            binxgcd_nz_test(Uint::ONE, int_abs_min);
620            binxgcd_nz_test(Uint::ONE, Uint::MAX);
621            binxgcd_nz_test(max_int, Uint::ONE);
622            binxgcd_nz_test(max_int, max_int);
623            binxgcd_nz_test(max_int, int_abs_min);
624            binxgcd_nz_test(max_int, Uint::MAX);
625            binxgcd_nz_test(Uint::MAX, Uint::ONE);
626            binxgcd_nz_test(Uint::MAX, max_int);
627            binxgcd_nz_test(Uint::MAX, int_abs_min);
628            binxgcd_nz_test(Uint::MAX, Uint::MAX);
629        }
630
631        #[test]
632        fn test_binxgcd_nz() {
633            binxgcd_nz_tests::<{ U64::LIMBS }, { U128::LIMBS }>();
634            binxgcd_nz_tests::<{ U128::LIMBS }, { U256::LIMBS }>();
635            binxgcd_nz_tests::<{ U192::LIMBS }, { U384::LIMBS }>();
636            binxgcd_nz_tests::<{ U256::LIMBS }, { U512::LIMBS }>();
637            binxgcd_nz_tests::<{ U384::LIMBS }, { U768::LIMBS }>();
638            binxgcd_nz_tests::<{ U512::LIMBS }, { U1024::LIMBS }>();
639            binxgcd_nz_tests::<{ U1024::LIMBS }, { U2048::LIMBS }>();
640            binxgcd_nz_tests::<{ U2048::LIMBS }, { U4096::LIMBS }>();
641            binxgcd_nz_tests::<{ U4096::LIMBS }, { U8192::LIMBS }>();
642        }
643    }
644
645    mod test_classic_binxgcd {
646        use crate::modular::bingcd::xgcd::tests::test_xgcd;
647        use crate::{
648            Concat, Int, U64, U128, U192, U256, U384, U512, U768, U1024, U2048, U4096, U8192, Uint,
649        };
650
651        fn classic_binxgcd_test<const LIMBS: usize, const DOUBLE: usize>(
652            lhs: Uint<LIMBS>,
653            rhs: Uint<LIMBS>,
654        ) where
655            Uint<LIMBS>: Concat<LIMBS, Output = Uint<DOUBLE>>,
656        {
657            let output = lhs
658                .to_odd()
659                .unwrap()
660                .classic_binxgcd(&rhs.to_odd().unwrap())
661                .divide();
662            test_xgcd(lhs, rhs, output);
663        }
664
665        fn classic_binxgcd_tests<const LIMBS: usize, const DOUBLE: usize>()
666        where
667            Uint<LIMBS>: Concat<LIMBS, Output = Uint<DOUBLE>>,
668        {
669            let max_int = *Int::MAX.as_uint();
670
671            classic_binxgcd_test(Uint::ONE, Uint::ONE);
672            classic_binxgcd_test(Uint::ONE, max_int);
673            classic_binxgcd_test(Uint::ONE, Uint::MAX);
674            classic_binxgcd_test(max_int, Uint::ONE);
675            classic_binxgcd_test(max_int, max_int);
676            classic_binxgcd_test(max_int, Uint::MAX);
677            classic_binxgcd_test(Uint::MAX, Uint::ONE);
678            classic_binxgcd_test(Uint::MAX, max_int);
679            classic_binxgcd_test(Uint::MAX, Uint::MAX);
680        }
681
682        #[test]
683        fn test_classic_binxgcd() {
684            classic_binxgcd_tests::<{ U64::LIMBS }, { U128::LIMBS }>();
685            classic_binxgcd_tests::<{ U128::LIMBS }, { U256::LIMBS }>();
686            classic_binxgcd_tests::<{ U192::LIMBS }, { U384::LIMBS }>();
687            classic_binxgcd_tests::<{ U256::LIMBS }, { U512::LIMBS }>();
688            classic_binxgcd_tests::<{ U384::LIMBS }, { U768::LIMBS }>();
689            classic_binxgcd_tests::<{ U512::LIMBS }, { U1024::LIMBS }>();
690            classic_binxgcd_tests::<{ U1024::LIMBS }, { U2048::LIMBS }>();
691            classic_binxgcd_tests::<{ U2048::LIMBS }, { U4096::LIMBS }>();
692            classic_binxgcd_tests::<{ U4096::LIMBS }, { U8192::LIMBS }>();
693        }
694    }
695
696    mod test_optimized_binxgcd {
697        use crate::modular::bingcd::xgcd::tests::test_xgcd;
698        use crate::modular::bingcd::xgcd::{DOUBLE_SUMMARY_LIMBS, SUMMARY_BITS, SUMMARY_LIMBS};
699        use crate::{
700            Concat, Int, U64, U128, U192, U256, U384, U512, U768, U1024, U2048, U4096, U8192, Uint,
701        };
702
703        fn test<const LIMBS: usize, const DOUBLE: usize>(lhs: Uint<LIMBS>, rhs: Uint<LIMBS>)
704        where
705            Uint<LIMBS>: Concat<LIMBS, Output = Uint<DOUBLE>>,
706        {
707            let output = lhs
708                .to_odd()
709                .unwrap()
710                .optimized_binxgcd(&rhs.to_odd().unwrap())
711                .divide();
712            test_xgcd(lhs, rhs, output);
713        }
714
715        fn run_tests<const LIMBS: usize, const DOUBLE: usize>()
716        where
717            Uint<LIMBS>: Concat<LIMBS, Output = Uint<DOUBLE>>,
718        {
719            let upper_bound = *Int::MAX.as_uint();
720            test(Uint::ONE, Uint::ONE);
721            test(Uint::ONE, upper_bound);
722            test(Uint::ONE, Uint::MAX);
723            test(upper_bound, Uint::ONE);
724            test(upper_bound, upper_bound);
725            test(upper_bound, Uint::MAX);
726            test(Uint::MAX, Uint::ONE);
727            test(Uint::MAX, upper_bound);
728            test(Uint::MAX, Uint::MAX);
729        }
730
731        #[test]
732        fn test_optimized_binxgcd_edge_cases() {
733            // If one of these tests fails, you have probably tweaked the SUMMARY_BITS,
734            // SUMMARY_LIMBS or DOUBLE_SUMMARY_LIMBS settings. Please make sure to update these
735            // tests accordingly.
736            assert_eq!(SUMMARY_BITS, 63);
737            assert_eq!(SUMMARY_LIMBS, U64::LIMBS);
738            assert_eq!(DOUBLE_SUMMARY_LIMBS, U128::LIMBS);
739
740            // Case #1: a > b but a.compact() < b.compact()
741            let a = U256::from_be_hex(
742                "1234567890ABCDEF80000000000000000000000000000000BEDCBA0987654321",
743            );
744            let b = U256::from_be_hex(
745                "1234567890ABCDEF800000000000000000000000000000007EDCBA0987654321",
746            );
747            assert!(a > b);
748            assert!(
749                a.compact::<SUMMARY_BITS, DOUBLE_SUMMARY_LIMBS>(U256::BITS)
750                    < b.compact::<SUMMARY_BITS, DOUBLE_SUMMARY_LIMBS>(U256::BITS)
751            );
752            test(a, b);
753
754            // Case #2: a < b but a.compact() > b.compact()
755            test(b, a);
756
757            // Case #3: a > b but a.compact() = b.compact()
758            let a = U256::from_be_hex(
759                "1234567890ABCDEF80000000000000000000000000000000FEDCBA0987654321",
760            );
761            let b = U256::from_be_hex(
762                "1234567890ABCDEF800000000000000000000000000000007EDCBA0987654321",
763            );
764            assert!(a > b);
765            assert_eq!(
766                a.compact::<SUMMARY_BITS, DOUBLE_SUMMARY_LIMBS>(U256::BITS),
767                b.compact::<SUMMARY_BITS, DOUBLE_SUMMARY_LIMBS>(U256::BITS)
768            );
769            test(a, b);
770
771            // Case #4: a < b but a.compact() = b.compact()
772            test(b, a);
773        }
774
775        #[test]
776        fn optimized_binxgcd() {
777            run_tests::<{ U128::LIMBS }, { U256::LIMBS }>();
778            run_tests::<{ U192::LIMBS }, { U384::LIMBS }>();
779            run_tests::<{ U256::LIMBS }, { U512::LIMBS }>();
780            run_tests::<{ U384::LIMBS }, { U768::LIMBS }>();
781            run_tests::<{ U512::LIMBS }, { U1024::LIMBS }>();
782            run_tests::<{ U1024::LIMBS }, { U2048::LIMBS }>();
783            run_tests::<{ U2048::LIMBS }, { U4096::LIMBS }>();
784            run_tests::<{ U4096::LIMBS }, { U8192::LIMBS }>();
785        }
786    }
787}
crypto_bigint/modular/bingcd/xgcd.rs

crypto_bigint/modular/bingcd/
xgcd.rs
