nekolib/math/
gcd_recip.rs

1//! 最大公約数と逆元。
2
3/// 最大公約数と逆元。
4///
5/// 次の条件を満たす唯一の $(g, r)$ を返す。
6/// - $g = \\gcd(a, b)$
7/// - $a\\cdot r \\equiv g \\pmod{b}$
8/// - $0\\le r \\lt b/g$
9///
10/// $a = 0$ のとき $g = b$ であり、$0\\le g \\lt b$ とはならないことに注意せよ[^1]。
11///
12/// [^1]: $g = 0$ とすると $b/g$ が定義できないため。
13///
14/// # Complexity
15/// $O(\\log(\\min\\{a, b\\}))$ time.
16///
17/// # Examples
18/// ```
19/// use nekolib::math::GcdRecip;
20///
21/// let (g, r) = 27_u32.gcd_recip(30);
22/// assert_eq!(g, 3);
23/// assert_eq!(r, 9);
24/// assert_eq!((27 * r) % 30, g);
25/// ```
26pub trait GcdRecip: Sized {
27    fn gcd_recip(self, other: Self) -> (Self, Self);
28}
29
30macro_rules! impl_uint {
31    ($t:ty) => {
32        impl GcdRecip for $t {
33            fn gcd_recip(self, other: Self) -> (Self, Self) {
34                assert!(other > 0);
35                let a = self % other;
36                if a == 0 {
37                    return (other, 0);
38                }
39
40                let mut s = other;
41                let mut t = a;
42                let mut m0 = 0;
43                let mut m1 = 1;
44                while t > 0 {
45                    let u = s / t;
46                    s -= t * u;
47                    // m0 -= m1 * u;
48                    let v = (m1 * u) % other;
49                    m0 = if m0 < v { m0 + other - v } else { m0 - v };
50                    std::mem::swap(&mut s, &mut t);
51                    std::mem::swap(&mut m0, &mut m1);
52                }
53                (s, m0 % (other / s))
54            }
55        }
56    };
57    ( $($t:ty)* ) => { $(impl_uint!($t);)* };
58}
59
60macro_rules! impl_int {
61    ($t:ty) => {
62        impl GcdRecip for $t {
63            fn gcd_recip(self, other: Self) -> (Self, Self) {
64                assert!(other > 0);
65                let a = self.rem_euclid(other);
66                if a == 0 {
67                    return (other, 0);
68                }
69
70                let mut s = other;
71                let mut t = a;
72                let mut m0 = 0;
73                let mut m1 = 1;
74                while t > 0 {
75                    let u = s / t;
76                    s -= t * u;
77                    m0 -= m1 * u;
78                    std::mem::swap(&mut s, &mut t);
79                    std::mem::swap(&mut m0, &mut m1);
80                }
81                if m0 < 0 {
82                    m0 += other / s;
83                }
84                (s, m0)
85            }
86        }
87    };
88    ( $($t:ty)* ) => { $(impl_int!($t);)* };
89}
90
91impl_uint!(u8 u16 u32 u64 u128 usize);
92impl_int!(i8 i16 i32 i64 i128 isize);
93
94#[test]
95fn test() {
96    for b in 1_i32..=1000 {
97        for a in 0..b {
98            let (g, r) = a.gcd_recip(b);
99            assert!(0 <= r && r < b / g);
100            assert_eq!(a * r % b, g % b);
101        }
102    }
103    for b in 1_u32..=1000 {
104        for a in 0..b {
105            let (g, r) = a.gcd_recip(b);
106            assert!(r < b / g);
107            assert_eq!(a * r % b, g % b);
108        }
109    }
110}
nekolib/math/gcd_recip.rs

nekolib/math/
gcd_recip.rs
