1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
/// compute g := \gcd(modulus, n),
/// and modular inverse of n/g in Z_{modulus/g}.
/// we convert parameters to i64 internally.
/// so be careful not to pass modulus > 2^63 because it overflows.
/// it's `trivial` that inverse of 0 is undefined, so if n = 0, it panics.
pub fn mod_gcd_inv(
modulus: u64,
n: u64,
) -> (u64, u64) {
assert!(0 < n && n < modulus);
let (mut a, mut b) = (n as i64, modulus as i64);
let (mut x00, mut x01) = (1, 0);
while b != 0 {
(x00, x01) = (x01, x00 - a / b * x01);
(a, b) = (b, a % b);
}
let gcd = a as u64;
let u = (modulus / gcd) as i64;
if x00 < 0 {
x00 += u;
}
debug_assert!(0 <= x00 && x00 < u);
(gcd, x00 as u64)
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_mod_gcd_inv() {
// euclidean_mod_gcd_inv(10, 0); // runtime error.
assert_eq!(mod_gcd_inv(5, 2), (1, 3));
assert_eq!(mod_gcd_inv(18, 12), (6, 2));
assert_eq!(mod_gcd_inv(111, 30), (3, 26));
// gcd(111, 30) = 3
// 111 / 3 = 37, 30 / 3 = 10, 10^{-1} \equiv 26 \mod 37
}
}