ruint/algorithms/gcd/
mod.rs

1#![allow(clippy::module_name_repetitions)]
2
3// TODO: https://github.com/bitcoin-core/secp256k1/blob/master/doc/safegcd_implementation.md
4
5// TODO: Make these algorithms work on limb slices.
6mod matrix;
7
8pub use self::matrix::Matrix as LehmerMatrix;
9use crate::Uint;
10use core::mem::swap;
11
12/// ⚠️ Lehmer's GCD algorithms.
13#[doc = crate::algorithms::unstable_warning!()]
14/// See [`gcd_extended`] for documentation.
15#[inline]
16#[must_use]
17pub fn gcd<const BITS: usize, const LIMBS: usize>(
18    mut a: Uint<BITS, LIMBS>,
19    mut b: Uint<BITS, LIMBS>,
20) -> Uint<BITS, LIMBS> {
21    if b > a {
22        swap(&mut a, &mut b);
23    }
24    while b != Uint::ZERO {
25        debug_assert!(a >= b);
26        let m = LehmerMatrix::from(a, b);
27        if m == LehmerMatrix::IDENTITY {
28            // Lehmer step failed to find a factor, which happens when
29            // the factor is very large. We do a regular Euclidean step, which
30            // will make a lot of progress since `q` will be large.
31            a %= b;
32            swap(&mut a, &mut b);
33        } else {
34            m.apply(&mut a, &mut b);
35        }
36    }
37    a
38}
39
40/// ⚠️ Lehmer's extended GCD.
41#[doc = crate::algorithms::unstable_warning!()]
42/// Returns `(gcd, x, y, sign)` such that `gcd = a * x + b * y`.
43///
44/// # Algorithm
45///
46/// A variation of Euclids algorithm where repeated 64-bit approximations are
47/// used to make rapid progress on.
48///
49/// See Jebelean (1994) "A Double-Digit Lehmer-Euclid Algorithm for Finding the
50/// GCD of Long Integers".
51///
52/// The function `lehmer_double` takes two `U256`'s and returns a 64-bit matrix.
53///
54/// The function `lehmer_update` updates state variables using this matrix. If
55/// the matrix makes no progress (because 64 bit precision is not enough) a full
56/// precision Euclid step is done, but this happens rarely.
57///
58/// See also `mpn_gcdext_lehmer_n` in GMP.
59/// <https://gmplib.org/repo/gmp-6.1/file/tip/mpn/generic/gcdext_lehmer.c#l146>
60#[inline]
61#[must_use]
62pub fn gcd_extended<const BITS: usize, const LIMBS: usize>(
63    mut a: Uint<BITS, LIMBS>,
64    mut b: Uint<BITS, LIMBS>,
65) -> (
66    Uint<BITS, LIMBS>,
67    Uint<BITS, LIMBS>,
68    Uint<BITS, LIMBS>,
69    bool,
70) {
71    if BITS == 0 {
72        return (Uint::ZERO, Uint::ZERO, Uint::ZERO, false);
73    }
74    let swapped = a < b;
75    if swapped {
76        swap(&mut a, &mut b);
77    }
78
79    // Initialize state matrix to identity.
80    let mut s0 = Uint::ONE;
81    let mut s1 = Uint::ZERO;
82    let mut t0 = Uint::ZERO;
83    let mut t1 = Uint::ONE;
84    let mut even = true;
85    while b != Uint::ZERO {
86        debug_assert!(a >= b);
87        let m = LehmerMatrix::from(a, b);
88        if m == LehmerMatrix::IDENTITY {
89            // Lehmer step failed to find a factor, which happens when
90            // the factor is very large. We do a regular Euclidean step, which
91            // will make a lot of progress since `q` will be large.
92            let q = a / b;
93            a -= q * b;
94            swap(&mut a, &mut b);
95            s0 -= q * s1;
96            swap(&mut s0, &mut s1);
97            t0 -= q * t1;
98            swap(&mut t0, &mut t1);
99            even = !even;
100        } else {
101            m.apply(&mut a, &mut b);
102            m.apply(&mut s0, &mut s1);
103            m.apply(&mut t0, &mut t1);
104            even ^= !m.4;
105        }
106    }
107    // TODO: Compute using absolute value instead of patching sign.
108    if even {
109        // t negative
110        t0 = Uint::ZERO - t0;
111    } else {
112        // s negative
113        s0 = Uint::ZERO - s0;
114    }
115    if swapped {
116        swap(&mut s0, &mut t0);
117        even = !even;
118    }
119    (a, s0, t0, even)
120}
121
122/// ⚠️ Modular inversion using extended GCD.
123#[doc = crate::algorithms::unstable_warning!()]
124/// It uses the Bezout identity
125///
126/// ```text
127///    a * modulus + b * num = gcd(modulus, num)
128/// ````
129///
130/// where `a` and `b` are the cofactors from the extended Euclidean algorithm.
131/// A modular inverse only exists if `modulus` and `num` are coprime, in which
132/// case `gcd(modulus, num)` is one. Reducing both sides by the modulus then
133/// results in the equation `b * num = 1 (mod modulus)`. In other words, the
134/// cofactor `b` is the modular inverse of `num`.
135///
136/// It differs from `gcd_extended` in that it only computes the required
137/// cofactor, and returns `None` if the GCD is not one (i.e. when `num` does
138/// not have an inverse).
139#[inline]
140#[must_use]
141pub fn inv_mod<const BITS: usize, const LIMBS: usize>(
142    num: Uint<BITS, LIMBS>,
143    modulus: Uint<BITS, LIMBS>,
144) -> Option<Uint<BITS, LIMBS>> {
145    if BITS == 0 || modulus.is_zero() {
146        return None;
147    }
148    let mut a = modulus;
149    let mut b = num;
150    if b >= a {
151        b %= a;
152    }
153    if b.is_zero() {
154        return None;
155    }
156
157    let mut t0 = Uint::ZERO;
158    let mut t1 = Uint::ONE;
159    let mut even = true;
160    while b != Uint::ZERO {
161        debug_assert!(a >= b);
162        let m = LehmerMatrix::from(a, b);
163        if m == LehmerMatrix::IDENTITY {
164            // Lehmer step failed to find a factor, which happens when
165            // the factor is very large. We do a regular Euclidean step, which
166            // will make a lot of progress since `q` will be large.
167            let q = a / b;
168            a -= q * b;
169            swap(&mut a, &mut b);
170            t0 -= q * t1;
171            swap(&mut t0, &mut t1);
172            even = !even;
173        } else {
174            m.apply(&mut a, &mut b);
175            m.apply(&mut t0, &mut t1);
176            even ^= !m.4;
177        }
178    }
179    if a == Uint::ONE {
180        // When `even` t0 is negative and in twos-complement form
181        Some(if even { modulus + t0 } else { t0 })
182    } else {
183        None
184    }
185}
186
187#[cfg(test)]
188#[allow(clippy::cast_lossless)]
189mod tests {
190    use super::*;
191    use crate::{const_for, nlimbs};
192    use proptest::{proptest, test_runner::Config};
193
194    #[test]
195    fn test_gcd_one() {
196        use core::str::FromStr;
197        const BITS: usize = 129;
198        const LIMBS: usize = nlimbs(BITS);
199        type U = Uint<BITS, LIMBS>;
200        let a = U::from_str("0x006d7c4641f88b729a97889164dd8d07db").unwrap();
201        let b = U::from_str("0x01de6ef6f3caa963a548d7a411b05b9988").unwrap();
202        assert_eq!(gcd(a, b), gcd_ref(a, b));
203    }
204
205    // Reference implementation
206    fn gcd_ref<const BITS: usize, const LIMBS: usize>(
207        mut a: Uint<BITS, LIMBS>,
208        mut b: Uint<BITS, LIMBS>,
209    ) -> Uint<BITS, LIMBS> {
210        while b != Uint::ZERO {
211            a %= b;
212            swap(&mut a, &mut b);
213        }
214        a
215    }
216
217    #[test]
218    #[allow(clippy::absurd_extreme_comparisons)] // Generated code
219    fn test_gcd() {
220        const_for!(BITS in SIZES {
221            const LIMBS: usize = nlimbs(BITS);
222            type U = Uint<BITS, LIMBS>;
223            let config = Config { cases: 10, ..Default::default()};
224            proptest!(config, |(a: U, b: U)| {
225                assert_eq!(gcd(a, b), gcd_ref(a, b));
226            });
227        });
228    }
229
230    #[test]
231    fn test_gcd_extended() {
232        const_for!(BITS in SIZES {
233            const LIMBS: usize = nlimbs(BITS);
234            type U = Uint<BITS, LIMBS>;
235            let config = Config { cases: 5, ..Default::default() };
236            proptest!(config, |(a: U, b: U)| {
237                let (g, x, y, sign) = gcd_extended(a, b);
238                assert_eq!(g, gcd_ref(a, b));
239                if sign {
240                    assert_eq!(a * x - b * y, g);
241                } else {
242                    assert_eq!(b * y - a * x, g);
243                }
244            });
245        });
246    }
247}
ruint/algorithms/gcd/mod.rs

ruint/algorithms/gcd/
mod.rs
