qualia_core_db/solvers/number_theory/
primes.rs

1//! Primality and factorization. `is_prime` is a **deterministic** Miller–Rabin (the
2//! witness set `{2,3,5,…,37}` is proven correct for all of `u64`), and `prime_factors`
3//! uses trial division for small factors then **Pollard's rho** (Brent's variant) for
4//! the rest, so factorization is correct across the full `u64` range — not just up to
5//! `√n`.
6
7use super::modular::mod_pow;
8
9#[inline]
10fn mulmod(a: u64, b: u64, m: u64) -> u64 {
11    ((a as u128 * b as u128) % m as u128) as u64
12}
13
14/// Deterministic Miller–Rabin primality test, exact for all `u64`.
15pub fn is_prime(n: u64) -> bool {
16    if n < 2 {
17        return false;
18    }
19    for &p in &[2u64, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37] {
20        if n == p {
21            return true;
22        }
23        if n % p == 0 {
24            return false;
25        }
26    }
27    // n − 1 = d · 2^r
28    let mut d = n - 1;
29    let mut r = 0u32;
30    while d & 1 == 0 {
31        d >>= 1;
32        r += 1;
33    }
34    'witness: for &a in &[2u64, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37] {
35        let mut x = mod_pow(a, d, n);
36        if x == 1 || x == n - 1 {
37            continue;
38        }
39        for _ in 0..r - 1 {
40            x = mulmod(x, x, n);
41            if x == n - 1 {
42                continue 'witness;
43            }
44        }
45        return false; // composite
46    }
47    true
48}
49
50/// The smallest prime strictly greater than `n`.
51pub fn next_prime(n: u64) -> u64 {
52    let mut c = n.saturating_add(1);
53    if c <= 2 {
54        return 2;
55    }
56    if c % 2 == 0 {
57        c += 1;
58    }
59    loop {
60        if is_prime(c) {
61            return c;
62        }
63        c = c.saturating_add(2);
64    }
65}
66
67/// Pollard's rho (Brent) — returns a non-trivial factor of a composite `n`.
68fn pollard_rho(n: u64) -> u64 {
69    if n % 2 == 0 {
70        return 2;
71    }
72    let mut c = 1u64;
73    loop {
74        let f = |x: u64| (mulmod(x, x, n) + c) % n;
75        let mut x = 2u64;
76        let mut y = 2u64;
77        let mut d = 1u64;
78        while d == 1 {
79            x = f(x);
80            y = f(f(y));
81            d = super::modular::gcd(x.abs_diff(y), n);
82        }
83        if d != n {
84            return d;
85        }
86        c += 1; // cycle hit n itself; retry with a different constant
87    }
88}
89
90fn factor_into(n: u64, out: &mut Vec<u64>) {
91    if n == 1 {
92        return;
93    }
94    if is_prime(n) {
95        out.push(n);
96        return;
97    }
98    let d = pollard_rho(n);
99    factor_into(d, out);
100    factor_into(n / d, out);
101}
102
103/// Prime factorization as `(prime, exponent)` pairs, ascending by prime. Empty for
104/// `n < 2` (0 and 1 have no prime factorization).
105pub fn prime_factors(n: u64) -> Vec<(u64, u32)> {
106    if n < 2 {
107        return Vec::new();
108    }
109    let mut flat = Vec::new();
110    factor_into(n, &mut flat);
111    flat.sort_unstable();
112    let mut out: Vec<(u64, u32)> = Vec::new();
113    for p in flat {
114        if let Some(last) = out.last_mut() {
115            if last.0 == p {
116                last.1 += 1;
117                continue;
118            }
119        }
120        out.push((p, 1));
121    }
122    out
123}
124
125/// All positive divisors of `n`, ascending. `[1]` for `n = 1`; empty for `n = 0`.
126pub fn divisors(n: u64) -> Vec<u64> {
127    if n == 0 {
128        return Vec::new();
129    }
130    let mut divs = vec![1u64];
131    for (p, e) in prime_factors(n) {
132        let mut pk = 1u64;
133        let base = divs.clone();
134        for _ in 0..e {
135            pk *= p;
136            for &d in &base {
137                divs.push(d * pk);
138            }
139        }
140    }
141    divs.sort_unstable();
142    divs
143}
144
145#[cfg(test)]
146mod tests {
147    use super::*;
148
149    #[test]
150    fn primality_known_cases() {
151        for p in [2u64, 3, 5, 7, 13, 97, 7919, 104729, 1_000_000_007] {
152            assert!(is_prime(p), "{p} should be prime");
153        }
154        for c in [0u64, 1, 4, 9, 15, 100, 7917, 1_000_000_011] {
155            assert!(!is_prime(c), "{c} should be composite");
156        }
157        // A large Carmichael number (561 = 3·11·17) fools Fermat but not Miller–Rabin.
158        assert!(!is_prime(561));
159        // A 64-bit semiprime is composite.
160        assert!(!is_prime(10_000_000_000_000_061 * 3));
161    }
162
163    #[test]
164    fn factorization_is_correct_and_reconstructs() {
165        assert_eq!(prime_factors(360), vec![(2, 3), (3, 2), (5, 1)]); // 2³·3²·5
166        assert_eq!(prime_factors(97), vec![(97, 1)]);
167        assert_eq!(prime_factors(1), vec![]);
168        // A hard semiprime that needs Pollard's rho (beyond √n trial division).
169        let n = 1_000_000_007u64 * 1_000_000_009u64;
170        let f = prime_factors(n);
171        assert_eq!(f, vec![(1_000_000_007, 1), (1_000_000_009, 1)]);
172        // Product of (prime^exp) reconstructs n.
173        let prod: u64 = f.iter().map(|&(p, e)| p.pow(e)).product();
174        assert_eq!(prod, n);
175    }
176
177    #[test]
178    fn divisors_of_28_are_perfect() {
179        let d = divisors(28);
180        assert_eq!(d, vec![1, 2, 4, 7, 14, 28]);
181        // 28 is perfect: its proper divisors sum to itself.
182        let proper: u64 = d.iter().filter(|&&x| x != 28).sum();
183        assert_eq!(proper, 28);
184    }
185
186    #[test]
187    fn next_prime_walks_forward() {
188        assert_eq!(next_prime(13), 17);
189        assert_eq!(next_prime(0), 2);
190        assert_eq!(next_prime(1), 2);
191        assert_eq!(next_prime(89), 97);
192    }
193}
qualia_core_db/solvers/number_theory/primes.rs

qualia_core_db/solvers/number_theory/
primes.rs
