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    }
128
129    /// Computes `gcd(self, rhs)`, leveraging the optimized Binary GCD algorithm.
130    ///
131    /// Ref: Pornin, Optimized Binary GCD for Modular Inversion, Algorithm 2.
132    /// <https://eprint.iacr.org/2020/972.pdf>
133    ///
134    /// In summary, the optimized algorithm does not operate on `self` and `rhs` directly, but
135    /// instead of condensed summaries that fit in few registers. Based on these summaries, an
136    /// update matrix is constructed by which `self` and `rhs` are updated in larger steps.
137    ///
138    /// This function is generic over the following three values:
139    /// - `K`: the number of bits used when summarizing `self` and `rhs` for the inner loop. The
140    ///   `K` top bits and `K` least significant bits are selected. It is recommended to keep
141    ///   `K` close to a (multiple of) the number of bits that fit in a single register.
142    /// - `LIMBS_K`: should be chosen as the minimum number s.t. [`Uint::<LIMBS>::BITS`] `≥ K`,
143    /// - `LIMBS_2K`: should be chosen as the minimum number s.t. [`Uint::<LIMBS>::BITS`] `≥ 2K`.
144    ///
145    /// This method returns a pair consisting of the GCD and the sign of the Jacobi symbol,
146    /// 0 for positive and 1 for negative.
147    #[inline(always)]
148    pub(crate) const fn optimized_bingcd_<
149        const K: u32,
150        const LIMBS_K: usize,
151        const LIMBS_2K: usize,
152    >(
153        &self,
154        rhs: &Uint<LIMBS>,
155        batch_max: u32,
156    ) -> Self {
157        let (mut a, mut b) = (*self.as_ref(), *rhs);
158
159        let mut iterations = Self::MINIMAL_BINGCD_ITERATIONS;
160        while iterations > 0 {
161            let batch = u32_min(iterations, batch_max);
162            iterations -= batch;
163
164            // Construct a_ and b_ as the summary of a and b, respectively.
165            let n = u32_max(2 * K, a.bitor(&b).bits());
166            let a_ = a.compact::<K, LIMBS_2K>(n);
167            let b_ = b.compact::<K, LIMBS_2K>(n);
168
169            // Compute the batch update matrix from a_ and b_.
170            let (.., matrix) = a_
171                .to_odd()
172                .expect_copied("a_ is always odd")
173                .partial_binxgcd::<LIMBS_K>(&b_, batch, Choice::FALSE);
174
175            // Update `a` and `b` using the update matrix.
176            // Safe to use vartime: the number of doublings is the same as the batch size.
177            let (updated_a, updated_b) = matrix.extended_apply_to_vartime::<LIMBS>((a, b));
178            (a, _) = updated_a.dropped_abs_sign();
179            (b, _) = updated_b.dropped_abs_sign();
180        }
181
182        a.to_odd()
183            .expect_copied("gcd of an odd value is always odd")
184    }
185
186    /// Variable time equivalent of [`Self::optimized_bingcd`].
187    #[inline(always)]
188    pub(crate) const fn optimized_bingcd_vartime(&self, rhs: &Uint<LIMBS>) -> Self {
189        self.optimized_bingcd_vartime_::<{ U64::BITS }, { U64::LIMBS }, { U128::LIMBS }>(
190            rhs,
191            U64::BITS - 1,
192        )
193    }
194
195    /// Variable time equivalent of [`Self::optimized_bingcd_`].
196    #[inline(always)]
197    pub(crate) const fn optimized_bingcd_vartime_<
198        const K: u32,
199        const LIMBS_K: usize,
200        const LIMBS_2K: usize,
201    >(
202        &self,
203        rhs: &Uint<LIMBS>,
204        batch_max: u32,
205    ) -> Self {
206        let (mut a, mut b) = (*self.as_ref(), *rhs);
207
208        while !b.is_zero_vartime() {
209            // Construct a_ and b_ as the summary of a and b, respectively.
210            let n = u32_max(2 * K, u32_max(a.bits_vartime(), b.bits_vartime()));
211            let mut a_ = a.compact_vartime::<K, LIMBS_2K>(n);
212            let mut b_ = b.compact_vartime::<K, LIMBS_2K>(n);
213
214            if n <= Uint::<LIMBS_2K>::BITS {
215                loop {
216                    bingcd_step(&mut b_, &mut a_);
217                    if b_.is_zero_vartime() {
218                        break;
219                    }
220                }
221                a = a_.resize();
222                break;
223            }
224
225            // Compute the batch update matrix from a_ and b_.
226            let (.., matrix) = a_
227                .to_odd()
228                .expect_copied("a_ is always odd")
229                .partial_binxgcd::<LIMBS_K>(&b_, batch_max, Choice::FALSE);
230
231            // Update `a` and `b` using the update matrix.
232            let (updated_a, updated_b) = matrix.extended_apply_to_vartime((a, b));
233            (a, _) = updated_a.dropped_abs_sign();
234            (b, _) = updated_b.dropped_abs_sign();
235        }
236
237        a.to_odd()
238            .expect_copied("gcd of an odd value with something else is always odd")
239    }
240}
241
242#[cfg(all(test, not(miri)))]
243mod tests {
244
245    mod test_classic_bingcd {
246        use crate::{Int, U64, U128, U192, U256, U384, U512, U1024, U2048, U4096, Uint};
247
248        fn classic_bingcd_test<const LIMBS: usize>(
249            lhs: Uint<LIMBS>,
250            rhs: Uint<LIMBS>,
251            target: Uint<LIMBS>,
252        ) {
253            assert_eq!(
254                lhs.to_odd().unwrap().classic_bingcd(&rhs),
255                target.to_odd().unwrap()
256            );
257            assert_eq!(
258                lhs.to_odd().unwrap().classic_bingcd_vartime(&rhs),
259                target.to_odd().unwrap()
260            );
261        }
262
263        fn classic_bingcd_tests<const LIMBS: usize>() {
264            classic_bingcd_test(Uint::<LIMBS>::ONE, Uint::ZERO, Uint::ONE);
265            classic_bingcd_test(Uint::<LIMBS>::ONE, Uint::ONE, Uint::ONE);
266            classic_bingcd_test(Uint::<LIMBS>::ONE, Int::MAX.abs(), Uint::ONE);
267            classic_bingcd_test(Uint::<LIMBS>::ONE, Int::MIN.abs(), Uint::ONE);
268            classic_bingcd_test(Uint::<LIMBS>::ONE, Uint::MAX, Uint::ONE);
269            classic_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::ZERO, Int::MAX.abs());
270            classic_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::ONE, Uint::ONE);
271            classic_bingcd_test(Int::<LIMBS>::MAX.abs(), Int::MAX.abs(), Int::MAX.abs());
272            classic_bingcd_test(Int::<LIMBS>::MAX.abs(), Int::MIN.abs(), Uint::ONE);
273            classic_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::MAX, Uint::ONE);
274            classic_bingcd_test(Uint::<LIMBS>::MAX, Uint::ZERO, Uint::MAX);
275            classic_bingcd_test(Uint::<LIMBS>::MAX, Uint::ONE, Uint::ONE);
276            classic_bingcd_test(Uint::<LIMBS>::MAX, Int::MAX.abs(), Uint::ONE);
277            classic_bingcd_test(Uint::<LIMBS>::MAX, Int::MIN.abs(), Uint::ONE);
278            classic_bingcd_test(Uint::<LIMBS>::MAX, Uint::MAX, Uint::MAX);
279        }
280
281        #[test]
282        fn test_classic_bingcd() {
283            classic_bingcd_tests::<{ U64::LIMBS }>();
284            classic_bingcd_tests::<{ U128::LIMBS }>();
285            classic_bingcd_tests::<{ U192::LIMBS }>();
286            classic_bingcd_tests::<{ U256::LIMBS }>();
287            classic_bingcd_tests::<{ U384::LIMBS }>();
288            classic_bingcd_tests::<{ U512::LIMBS }>();
289            classic_bingcd_tests::<{ U1024::LIMBS }>();
290            classic_bingcd_tests::<{ U2048::LIMBS }>();
291            classic_bingcd_tests::<{ U4096::LIMBS }>();
292        }
293    }
294
295    mod test_optimized_bingcd {
296        use crate::{Int, U128, U192, U256, U384, U512, U1024, U2048, U4096, Uint};
297
298        fn optimized_bingcd_test<const LIMBS: usize>(
299            lhs: Uint<LIMBS>,
300            rhs: Uint<LIMBS>,
301            target: Uint<LIMBS>,
302        ) {
303            assert_eq!(
304                lhs.to_odd().unwrap().optimized_bingcd(&rhs),
305                target.to_odd().unwrap()
306            );
307            assert_eq!(
308                lhs.to_odd().unwrap().optimized_bingcd_vartime(&rhs),
309                target.to_odd().unwrap()
310            );
311        }
312
313        fn optimized_bingcd_tests<const LIMBS: usize>() {
314            optimized_bingcd_test(Uint::<LIMBS>::ONE, Uint::ZERO, Uint::ONE);
315            optimized_bingcd_test(Uint::<LIMBS>::ONE, Uint::ONE, Uint::ONE);
316            optimized_bingcd_test(Uint::<LIMBS>::ONE, Int::MAX.abs(), Uint::ONE);
317            optimized_bingcd_test(Uint::<LIMBS>::ONE, Int::MIN.abs(), Uint::ONE);
318            optimized_bingcd_test(Uint::<LIMBS>::ONE, Uint::MAX, Uint::ONE);
319            optimized_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::ZERO, Int::MAX.abs());
320            optimized_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::ONE, Uint::ONE);
321            optimized_bingcd_test(Int::<LIMBS>::MAX.abs(), Int::MAX.abs(), Int::MAX.abs());
322            optimized_bingcd_test(Int::<LIMBS>::MAX.abs(), Int::MIN.abs(), Uint::ONE);
323            optimized_bingcd_test(Int::<LIMBS>::MAX.abs(), Uint::MAX, Uint::ONE);
324            optimized_bingcd_test(Uint::<LIMBS>::MAX, Uint::ZERO, Uint::MAX);
325            optimized_bingcd_test(Uint::<LIMBS>::MAX, Uint::ONE, Uint::ONE);
326            optimized_bingcd_test(Uint::<LIMBS>::MAX, Int::MAX.abs(), Uint::ONE);
327            optimized_bingcd_test(Uint::<LIMBS>::MAX, Int::MIN.abs(), Uint::ONE);
328            optimized_bingcd_test(Uint::<LIMBS>::MAX, Uint::MAX, Uint::MAX);
329        }
330
331        #[test]
332        fn test_optimized_bingcd() {
333            // Not applicable for U64
334            optimized_bingcd_tests::<{ U128::LIMBS }>();
335            optimized_bingcd_tests::<{ U192::LIMBS }>();
336            optimized_bingcd_tests::<{ U256::LIMBS }>();
337            optimized_bingcd_tests::<{ U384::LIMBS }>();
338            optimized_bingcd_tests::<{ U512::LIMBS }>();
339            optimized_bingcd_tests::<{ U1024::LIMBS }>();
340            optimized_bingcd_tests::<{ U2048::LIMBS }>();
341            optimized_bingcd_tests::<{ U4096::LIMBS }>();
342        }
343    }
344}
crypto_bigint/modular/bingcd/gcd.rs

crypto_bigint/modular/bingcd/
gcd.rs
