qualia_core_db/solvers/number_theory/
modular.rs

1//! Modular arithmetic and the Euclidean algorithms — the backbone the rest of the
2//! library (primality, totient, CRT) is built on. Overflow-safe via `u128`/`i128`
3//! intermediates.
4
5/// Greatest common divisor (binary/Euclid). `gcd(0, n) = n`.
6pub fn gcd(mut a: u64, mut b: u64) -> u64 {
7    while b != 0 {
8        let t = b;
9        b = a % b;
10        a = t;
11    }
12    a
13}
14
15/// Least common multiple. `0` if either argument is `0`.
16pub fn lcm(a: u64, b: u64) -> u64 {
17    if a == 0 || b == 0 {
18        return 0;
19    }
20    a / gcd(a, b) * b
21}
22
23/// Extended Euclid: returns `(g, x, y)` with `a·x + b·y = g = gcd(a, b)`.
24pub fn extended_gcd(a: i64, b: i64) -> (i64, i64, i64) {
25    if b == 0 {
26        return (a.abs(), if a < 0 { -1 } else { 1 }, 0);
27    }
28    let (g, x, y) = extended_gcd(b, a % b);
29    (g, y, x - (a / b) * y)
30}
31
32/// `(base^exp) mod modulus` by repeated squaring. `modulus = 0` → `0` (degenerate).
33pub fn mod_pow(mut base: u64, mut exp: u64, modulus: u64) -> u64 {
34    if modulus <= 1 {
35        return 0;
36    }
37    let m = modulus as u128;
38    let mut result: u128 = 1;
39    base %= modulus;
40    let mut b = base as u128;
41    while exp > 0 {
42        if exp & 1 == 1 {
43            result = result * b % m;
44        }
45        b = b * b % m;
46        exp >>= 1;
47    }
48    result as u64
49}
50
51/// Modular multiplicative inverse of `a` mod `m`: the `x` with `a·x ≡ 1 (mod m)`.
52/// `None` when `gcd(a, m) ≠ 1` (no inverse exists) — fail closed.
53pub fn mod_inverse(a: u64, m: u64) -> Option<u64> {
54    if m == 0 {
55        return None;
56    }
57    let (g, x, _) = extended_gcd((a % m) as i64, m as i64);
58    if g != 1 {
59        return None;
60    }
61    Some(((x % m as i64 + m as i64) % m as i64) as u64)
62}
63
64/// Chinese Remainder: solve `x ≡ r1 (mod m1)`, `x ≡ r2 (mod m2)` for **coprime**
65/// moduli, returning `(x, m1·m2)`. `None` if the moduli are not coprime.
66pub fn crt(r1: u64, m1: u64, r2: u64, m2: u64) -> Option<(u64, u64)> {
67    let inv = mod_inverse(m1 % m2, m2)?;
68    let m = m1.checked_mul(m2)?;
69    // x = r1 + m1 * ((r2 - r1) * inv mod m2)
70    let diff = (r2 as i128 - r1 as i128).rem_euclid(m2 as i128) as u128;
71    let t = diff * inv as u128 % m2 as u128;
72    let x = (r1 as u128 + m1 as u128 * t) % m as u128;
73    Some((x as u64, m))
74}
75
76#[cfg(test)]
77mod tests {
78    use super::*;
79
80    #[test]
81    fn gcd_lcm_basics() {
82        assert_eq!(gcd(54, 24), 6);
83        assert_eq!(gcd(17, 5), 1);
84        assert_eq!(gcd(0, 9), 9);
85        assert_eq!(lcm(4, 6), 12);
86        assert_eq!(lcm(0, 5), 0);
87    }
88
89    #[test]
90    fn extended_gcd_satisfies_bezout() {
91        let (g, x, y) = extended_gcd(240, 46);
92        assert_eq!(g, 2);
93        assert_eq!(240 * x + 46 * y, g);
94    }
95
96    #[test]
97    fn mod_pow_matches_known_values() {
98        assert_eq!(mod_pow(2, 10, 1000), 24); // 1024 mod 1000
99        assert_eq!(mod_pow(3, 0, 7), 1);
100        // Fermat: a^(p-1) ≡ 1 (mod p) for prime p ∤ a.
101        assert_eq!(mod_pow(2, 12, 13), 1);
102    }
103
104    #[test]
105    fn mod_inverse_exists_iff_coprime() {
106        // 3·4 = 12 ≡ 1 (mod 11)
107        assert_eq!(mod_inverse(3, 11), Some(4));
108        assert_eq!((3 * 4) % 11, 1);
109        // No inverse for 4 mod 8 (gcd 4).
110        assert_eq!(mod_inverse(4, 8), None);
111    }
112
113    #[test]
114    fn crt_combines_congruences() {
115        // x ≡ 2 (mod 3), x ≡ 3 (mod 5) → x = 8 (mod 15).
116        let (x, m) = crt(2, 3, 3, 5).unwrap();
117        assert_eq!(m, 15);
118        assert_eq!(x, 8);
119        assert_eq!(x % 3, 2);
120        assert_eq!(x % 5, 3);
121        // Non-coprime moduli → None.
122        assert!(crt(1, 4, 2, 6).is_none());
123    }
124}
qualia_core_db/solvers/number_theory/modular.rs

qualia_core_db/solvers/number_theory/
modular.rs
