crypto_bigint/modular/bingcd/
gcd.rs

1use crate::{
2    Choice, Limb, Odd, U64, U128, Uint, Word,
3    primitives::{u32_max, u32_min},
4    word,
5};
6
7/// Binary GCD 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///    a ← a - b
15/// a ← a/2
16/// ```
17///
18/// Note: assumes `b` to be odd. Might yield an incorrect result if this is not the case.
19///
20/// Ref: Pornin, Algorithm 1, L3-9, <https://eprint.iacr.org/2020/972.pdf>.
21#[inline(always)]
22pub(super) const fn bingcd_step<const LIMBS: usize>(
23    a: &mut Uint<LIMBS>,
24    b: &mut Uint<LIMBS>,
25) -> (Choice, Choice, Word) {
26    let a_b_mod_4 = (a.limbs[0].0 & b.limbs[0].0) & 3;
27
28    let a_odd = a.is_odd();
29    let (a_sub_b, borrow) = a.borrowing_sub(&Uint::select(&Uint::ZERO, b, a_odd), Limb::ZERO);
30    let swap = borrow.is_nonzero();
31    *b = Uint::select(b, a, swap);
32    *a = a_sub_b.wrapping_neg_if(swap).shr1();
33
34    // (b|a) = -(a|b) iff a = b = 3 mod 4 (quadratic reciprocity)
35    let mut jacobi_neg = word::select(0, a_b_mod_4 & (a_b_mod_4 >> 1) & 1, swap);
36
37    // (2a|b) = -(a|b) iff b = ±3 mod 8
38    // b is always odd, so we ignore the lower bit and check that bits 2 and 3 are '01' or '10'
39    let b_lo = b.limbs[0].0;
40    jacobi_neg ^= ((b_lo >> 1) ^ (b_lo >> 2)) & 1;
41
42    (a_odd, swap, jacobi_neg)
43}
44
45impl<const LIMBS: usize> Odd<Uint<LIMBS>> {
46    /// The minimal number of binary GCD iterations required to guarantee successful completion.
47    pub(crate) const MINIMAL_BINGCD_ITERATIONS: u32 = 2 * Self::BITS - 1;
48
49    /// Computes `gcd(self, rhs)`, leveraging (a constant time implementation of) the classic
50    /// Binary GCD algorithm.
51    ///
52    /// Note: this algorithm is efficient for [Uint]s with relatively few `LIMBS`.
53    ///
54    /// Ref: Pornin, Optimized Binary GCD for Modular Inversion, Algorithm 1.
55    /// <https://eprint.iacr.org/2020/972.pdf>
56    #[inline]
57    pub(crate) const fn classic_bingcd(&self, rhs: &Uint<LIMBS>) -> Self {
58        self.classic_bingcd_(rhs).0
59    }
60
61    /// Computes `gcd(self, rhs)`, leveraging (a constant time implementation of) the classic
62    /// Binary GCD algorithm.
63    ///
64    /// Note: this algorithm is efficient for [Uint]s with relatively few `LIMBS`.
65    ///
66    /// Ref: Pornin, Optimized Binary GCD for Modular Inversion, Algorithm 1.
67    /// <https://eprint.iacr.org/2020/972.pdf>
68    ///
69    /// This method returns a pair consisting of the GCD and the sign of the Jacobi symbol,
70    /// 0 for positive and 1 for negative.
71    #[inline(always)]
72    pub(crate) const fn classic_bingcd_(&self, rhs: &Uint<LIMBS>) -> (Self, Word) {
73        // (self, rhs) corresponds to (m, y) in the Algorithm 1 notation.
74        let (mut a, mut b) = (*rhs, *self.as_ref());
75        let mut i = 0;
76        let mut jacobi_neg = 0;
77
78        while i < Self::MINIMAL_BINGCD_ITERATIONS {
79            jacobi_neg ^= bingcd_step(&mut a, &mut b).2;
80            i += 1;
81        }
82
83        let gcd = b
84            .to_odd()
85            .expect_copied("gcd of an odd value with something else is always odd");
86
87        (gcd, jacobi_neg)
88    }
89
90    /// Variable time equivalent of [`Self::classic_bingcd`].
91    #[inline]
92    pub(crate) const fn classic_bingcd_vartime(&self, rhs: &Uint<LIMBS>) -> Self {
93        self.classic_bingcd_vartime_(rhs).0
94    }
95
96    /// Variable time equivalent of [`Self::classic_bingcd_`].
97    #[inline(always)]
98    pub(crate) const fn classic_bingcd_vartime_(&self, rhs: &Uint<LIMBS>) -> (Self, Word) {
99        // (self, rhs) corresponds to (m, y) in the Algorithm 1 notation.
100        let (mut a, mut b) = (*rhs, *self.as_ref());
101        let mut jacobi_neg = 0;
102        while !a.is_zero_vartime() {
103            jacobi_neg ^= bingcd_step(&mut a, &mut b).2;
104        }
105
106        let gcd = b
107            .to_odd()
108            .expect_copied("gcd of an odd value with something else is always odd");
109
110        (gcd, jacobi_neg)
111    }
112
113    /// Computes `gcd(self, rhs)`, leveraging the optimized Binary GCD algorithm.
114    ///
115    /// Note: this algorithm becomes more efficient than the classical algorithm for [Uint]s with
116    /// relatively many `LIMBS`. A best-effort threshold is presented in [`Self::bingcd`].
117    ///
118    /// Note: the full algorithm has an additional parameter; this function selects the best-effort
119    /// value for this parameter. You might be able to further tune your performance by calling the
120    /// [`Self::optimized_bingcd`_] function directly.
121    ///
122    /// Ref: Pornin, Optimized Binary GCD for Modular Inversion, Algorithm 2.
123    /// <https://eprint.iacr.org/2020/972.pdf>
124    #[inline(always)]
125    pub(crate) const fn optimized_bingcd(&self, rhs: &Uint<LIMBS>) -> Self {
126        self.optimized_bingcd_::<{ U64::BITS }, { U64::LIMBS }, { U128::LIMBS }>(rhs, U64::BITS - 1)
127            .0
128    }
129
130    /// Computes `gcd(self, rhs)`, leveraging the optimized Binary GCD algorithm.
131    ///
132    /// Ref: Pornin, Optimized Binary GCD for Modular Inversion, Algorithm 2.
133    /// <https://eprint.iacr.org/2020/972.pdf>
134    ///
135    /// In summary, the optimized algorithm does not operate on `self` and `rhs` directly, but
136    /// instead of condensed summaries that fit in few registers. Based on these summaries, an
137    /// update matrix is constructed by which `self` and `rhs` are updated in larger steps.
138    ///
139    /// This function is generic over the following three values:
140    /// - `K`: the number of bits used when summarizing `self` and `rhs` for the inner loop. The
141    ///   `K` top bits and `K` least significant bits are selected. It is recommended to keep
142    ///   `K` close to a (multiple of) the number of bits that fit in a single register.
143    /// - `LIMBS_K`: should be chosen as the minimum number s.t. [`Uint::<LIMBS>::BITS`] `≥ K`,
144    /// - `LIMBS_2K`: should be chosen as the minimum number s.t. [`Uint::<LIMBS>::BITS`] `≥ 2K`.
145    ///
146    /// This method returns a pair consisting of the GCD and the sign of the Jacobi symbol,
147    /// 0 for positive and 1 for negative.
148    #[inline(always)]
149    pub(crate) const fn optimized_bingcd_<
150        const K: u32,
151        const LIMBS_K: usize,
152        const LIMBS_2K: usize,
153    >(
154        &self,
155        rhs: &Uint<LIMBS>,
156        batch_max: u32,
157    ) -> (Self, Word) {
158        let (mut a, mut b) = (*self.as_ref(), *rhs);
159        let mut jacobi_neg = 0;
160
161        let mut iterations = Self::MINIMAL_BINGCD_ITERATIONS;
162        while iterations > 0 {
163            let batch = u32_min(iterations, batch_max);
164            iterations -= batch;
165
166            // Construct a_ and b_ as the summary of a and b, respectively.
167            let n = u32_max(2 * K, a.bitor(&b).bits());
168            let a_ = a.compact::<K, LIMBS_2K>(n);
169            let b_ = b.compact::<K, LIMBS_2K>(n);
170
171            // Compute the batch update matrix from a_ and b_.
172            let (.., matrix, j_neg) = a_
173                .to_odd()
174                .expect_copied("a_ is always odd")
175                .partial_binxgcd::<LIMBS_K>(&b_, batch, Choice::FALSE);
176            jacobi_neg ^= j_neg;
177
178            // Update `a` and `b` using the update matrix.
179            // Safe to use vartime: the number of doublings is the same as the batch size.
180            let (updated_a, updated_b) = matrix.extended_apply_to_vartime::<LIMBS>((a, b));
181            (a, _) = updated_a.dropped_abs_sign();
182            (b, _) = updated_b.dropped_abs_sign();
183        }
184
185        (
186            a.to_odd()
187                .expect_copied("gcd of an odd value is always odd"),
188            jacobi_neg,
189        )
190    }
191
192    /// Variable time equivalent of [`Self::optimized_bingcd`].
193    #[inline(always)]
194    pub(crate) const fn optimized_bingcd_vartime(&self, rhs: &Uint<LIMBS>) -> Self {
195        self.optimized_bingcd_vartime_::<{ U64::BITS }, { U64::LIMBS }, { U128::LIMBS }>(
196            rhs,
197            U64::BITS - 1,
198        )
199        .0
200    }
201
202    /// Variable time equivalent of [`Self::optimized_bingcd_`].
203    #[inline(always)]
204    pub(crate) const fn optimized_bingcd_vartime_<
205        const K: u32,
206        const LIMBS_K: usize,
207        const LIMBS_2K: usize,
208    >(
209        &self,
210        rhs: &Uint<LIMBS>,
211        batch_max: u32,
212    ) -> (Self, Word) {
213        let (mut a, mut b) = (*self.as_ref(), *rhs);
214        let mut jacobi_neg = 0;
215
216        while !b.is_zero_vartime() {
217            // Construct a_ and b_ as the summary of a and b, respectively.
218            let n = u32_max(2 * K, u32_max(a.bits_vartime(), b.bits_vartime()));
219            let mut a_ = a.compact_vartime::<K, LIMBS_2K>(n);
220            let mut b_ = b.compact_vartime::<K, LIMBS_2K>(n);
221
222            if n <= Uint::<LIMBS_2K>::BITS {
223                loop {
224                    jacobi_neg ^= bingcd_step(&mut b_, &mut a_).2;
225                    if b_.is_zero_vartime() {
226                        break;
227                    }
228                }
229                a = a_.resize();
230                break;
231            }
232
233            // Compute the batch update matrix from a_ and b_.
234            let (.., matrix, j_neg) = a_
235                .to_odd()
236                .expect_copied("a_ is always odd")
237                .partial_binxgcd::<LIMBS_K>(&b_, batch_max, Choice::FALSE);
238            jacobi_neg ^= j_neg;
239
240            // Update `a` and `b` using the update matrix.
241            let (updated_a, updated_b) = matrix.extended_apply_to_vartime((a, b));
242            (a, _) = updated_a.dropped_abs_sign();
243            (b, _) = updated_b.dropped_abs_sign();
244        }
245
246        (
247            a.to_odd()
248                .expect_copied("gcd of an odd value with something else is always odd"),
249            jacobi_neg,
250        )
251    }
252}
253
254#[cfg(all(test, not(miri)))]
255mod tests {
256
257    mod test_classic_bingcd {
258        use crate::{Int, U64, U128, U192, U256, U384, U512, U1024, U2048, U4096, Uint};
259
260        fn classic_bingcd_test<const LIMBS: usize>(
261            lhs: Uint<LIMBS>,
262            rhs: Uint<LIMBS>,
263            target: Uint<LIMBS>,
264        ) {
265            assert_eq!(
266                lhs.to_odd().unwrap().classic_bingcd(&rhs),
267                target.to_odd().unwrap()
268            );
269            assert_eq!(
270                lhs.to_odd().unwrap().classic_bingcd_vartime(&rhs),
271                target.to_odd().unwrap()
272            );
273        }
274
275        fn classic_bingcd_tests<const LIMBS: usize>() {
276            classic_bingcd_test(Uint::<LIMBS>::ONE, Uint::ZERO, Uint::ONE);
277            classic_bingcd_test(Uint::<LIMBS>::ONE, Uint::ONE, Uint::ONE);
278            classic_bingcd_test(Uint::<LIMBS>::ONE, Int::MAX.abs(), Uint::ONE);
279            classic_bingcd_test(Uint::<LIMBS>::ONE, Int::MIN.abs(), Uint::ONE);
280            classic_bingcd_test(Uint::<LIMBS>::ONE, Uint::MAX, Uint::ONE);
281            classic_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::ZERO, Int::MAX.abs());
282            classic_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::ONE, Uint::ONE);
283            classic_bingcd_test(Int::<LIMBS>::MAX.abs(), Int::MAX.abs(), Int::MAX.abs());
284            classic_bingcd_test(Int::<LIMBS>::MAX.abs(), Int::MIN.abs(), Uint::ONE);
285            classic_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::MAX, Uint::ONE);
286            classic_bingcd_test(Uint::<LIMBS>::MAX, Uint::ZERO, Uint::MAX);
287            classic_bingcd_test(Uint::<LIMBS>::MAX, Uint::ONE, Uint::ONE);
288            classic_bingcd_test(Uint::<LIMBS>::MAX, Int::MAX.abs(), Uint::ONE);
289            classic_bingcd_test(Uint::<LIMBS>::MAX, Int::MIN.abs(), Uint::ONE);
290            classic_bingcd_test(Uint::<LIMBS>::MAX, Uint::MAX, Uint::MAX);
291        }
292
293        #[test]
294        fn test_classic_bingcd() {
295            classic_bingcd_tests::<{ U64::LIMBS }>();
296            classic_bingcd_tests::<{ U128::LIMBS }>();
297            classic_bingcd_tests::<{ U192::LIMBS }>();
298            classic_bingcd_tests::<{ U256::LIMBS }>();
299            classic_bingcd_tests::<{ U384::LIMBS }>();
300            classic_bingcd_tests::<{ U512::LIMBS }>();
301            classic_bingcd_tests::<{ U1024::LIMBS }>();
302            classic_bingcd_tests::<{ U2048::LIMBS }>();
303            classic_bingcd_tests::<{ U4096::LIMBS }>();
304        }
305    }
306
307    mod test_optimized_bingcd {
308        use crate::{Int, U128, U192, U256, U384, U512, U1024, U2048, U4096, Uint};
309
310        fn optimized_bingcd_test<const LIMBS: usize>(
311            lhs: Uint<LIMBS>,
312            rhs: Uint<LIMBS>,
313            target: Uint<LIMBS>,
314        ) {
315            assert_eq!(
316                lhs.to_odd().unwrap().optimized_bingcd(&rhs),
317                target.to_odd().unwrap()
318            );
319            assert_eq!(
320                lhs.to_odd().unwrap().optimized_bingcd_vartime(&rhs),
321                target.to_odd().unwrap()
322            );
323        }
324
325        fn optimized_bingcd_tests<const LIMBS: usize>() {
326            optimized_bingcd_test(Uint::<LIMBS>::ONE, Uint::ZERO, Uint::ONE);
327            optimized_bingcd_test(Uint::<LIMBS>::ONE, Uint::ONE, Uint::ONE);
328            optimized_bingcd_test(Uint::<LIMBS>::ONE, Int::MAX.abs(), Uint::ONE);
329            optimized_bingcd_test(Uint::<LIMBS>::ONE, Int::MIN.abs(), Uint::ONE);
330            optimized_bingcd_test(Uint::<LIMBS>::ONE, Uint::MAX, Uint::ONE);
331            optimized_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::ZERO, Int::MAX.abs());
332            optimized_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::ONE, Uint::ONE);
333            optimized_bingcd_test(Int::<LIMBS>::MAX.abs(), Int::MAX.abs(), Int::MAX.abs());
334            optimized_bingcd_test(Int::<LIMBS>::MAX.abs(), Int::MIN.abs(), Uint::ONE);
335            optimized_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::MAX, Uint::ONE);
336            optimized_bingcd_test(Uint::<LIMBS>::MAX, Uint::ZERO, Uint::MAX);
337            optimized_bingcd_test(Uint::<LIMBS>::MAX, Uint::ONE, Uint::ONE);
338            optimized_bingcd_test(Uint::<LIMBS>::MAX, Int::MAX.abs(), Uint::ONE);
339            optimized_bingcd_test(Uint::<LIMBS>::MAX, Int::MIN.abs(), Uint::ONE);
340            optimized_bingcd_test(Uint::<LIMBS>::MAX, Uint::MAX, Uint::MAX);
341        }
342
343        #[test]
344        fn test_optimized_bingcd() {
345            // Not applicable for U64
346            optimized_bingcd_tests::<{ U128::LIMBS }>();
347            optimized_bingcd_tests::<{ U192::LIMBS }>();
348            optimized_bingcd_tests::<{ U256::LIMBS }>();
349            optimized_bingcd_tests::<{ U384::LIMBS }>();
350            optimized_bingcd_tests::<{ U512::LIMBS }>();
351            optimized_bingcd_tests::<{ U1024::LIMBS }>();
352            optimized_bingcd_tests::<{ U2048::LIMBS }>();
353            optimized_bingcd_tests::<{ U4096::LIMBS }>();
354        }
355    }
356}
crypto_bigint/modular/bingcd/gcd.rs

crypto_bigint/modular/bingcd/
gcd.rs
