factor.rs

1// SIEVES23pub fn small_primes(limit: usize) -> Vec<u64> {4    let mut comp = vec![false; limit + 1];5    let mut primes = Vec::new();6    for i in 2..=limit {7        if !comp[i] {8            primes.push(i as u64);9            let mut j = i * i;10            while j <= limit {11                comp[j] = true;12                j += i;13            }14        }15    }16    primes17}1819pub fn mu_sieve(limit: usize) -> Vec<i8> {20    let mut mu = vec![0i8; limit + 1];21    let mut comp = vec![false; limit + 1];22    let mut primes: Vec<u32> = Vec::new();23    if limit >= 1 {24        mu[1] = 1;25    }26    for i in 2..=limit {27        if !comp[i] {28            primes.push(i as u32);29            mu[i] = -1;30        }31        for &p in &primes {32            let ip = i * p as usize;33            if ip > limit {34                break;35            }36            comp[ip] = true;37            if i % p as usize == 0 {38                mu[ip] = 0;39                break;40            }41            mu[ip] = -mu[i];42        }43    }44    mu45}4647// MODULAR ARITHMETIC4849fn mulmod(a: u64, b: u64, n: u64) -> u64 {50    ((a as u128 * b as u128) % n as u128) as u6451}5253fn powmod(mut a: u64, mut e: u64, n: u64) -> u64 {54    let mut r = 1u64 % n;55    a %= n;56    while e > 0 {57        if e & 1 == 1 {58            r = mulmod(r, a, n);59        }60        a = mulmod(a, a, n);61        e >>= 1;62    }63    r64}6566fn gcd(mut a: u64, mut b: u64) -> u64 {67    while b != 0 {68        let t = a % b;69        a = b;70        b = t;71    }72    a73}7475// PRIMALITY7677const WITNESSES: [u64; 12] = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37];7879pub fn is_prime(n: u64) -> bool {80    if n < 2 {81        return false;82    }83    for p in WITNESSES {84        if n == p {85            return true;86        }87        if n % p == 0 {88            return false;89        }90    }91    let mut d = n - 1;92    let mut s = 0;93    while d & 1 == 0 {94        d >>= 1;95        s += 1;96    }97    'witness: for a in WITNESSES {98        let mut x = powmod(a, d, n);99        if x == 1 || x == n - 1 {100            continue;101        }102        for _ in 1..s {103            x = mulmod(x, x, n);104            if x == n - 1 {105                continue 'witness;106            }107        }108        return false;109    }110    true111}112113// FACTOR SPLITTING114115pub fn isqrt(n: u64) -> u64 {116    let mut r = (n as f64).sqrt() as u64;117    while r > 0 && r as u128 * r as u128 > n as u128 {118        r -= 1;119    }120    while (r as u128 + 1) * (r as u128 + 1) <= n as u128 {121        r += 1;122    }123    r124}125126fn brent(n: u64) -> u64 {127    if n & 1 == 0 {128        return 2;129    }130    let mut c = 1u64;131    loop {132        let step = |x: u64| {133            let y = mulmod(x, x, n) + c;134            if y >= n {135                y - n136            } else {137                y138            }139        };140        let mut y = 2u64;141        let mut x = y;142        let mut ys = y;143        let mut q = 1u64;144        let mut d = 1u64;145        let mut r = 1u64;146        let m = 128u64;147        while d == 1 {148            x = y;149            for _ in 0..r {150                y = step(y);151            }152            let mut k = 0u64;153            while k < r && d == 1 {154                ys = y;155                for _ in 0..m.min(r - k) {156                    y = step(y);157                    q = mulmod(q, x.abs_diff(y), n);158                }159                d = gcd(q, n);160                k += m;161            }162            r <<= 1;163        }164        if d == n {165            loop {166                ys = step(ys);167                d = gcd(x.abs_diff(ys), n);168                if d > 1 {169                    break;170                }171            }172        }173        if d < n {174            return d;175        }176        c += 1;177    }178}179180// MOBIUS181182pub fn mobius(n: u64, primes: &[u64]) -> i8 {183    if n == 1 {184        return 1;185    }186    let mut m = n;187    let mut parity = 0u32;188    let mut settled = false;189    for &p in primes {190        if m == 1 || p * p > m {191            settled = true;192            break;193        }194        if m % p == 0 {195            m /= p;196            if m % p == 0 {197                return 0;198            }199            parity ^= 1;200        }201    }202    if m > 1 {203        if settled {204            parity ^= 1;205        } else {206            let mut stack = vec![m];207            let mut found: Vec<u64> = Vec::new();208            while let Some(x) = stack.pop() {209                if is_prime(x) {210                    found.push(x);211                    continue;212                }213                let r = isqrt(x);214                if r * r == x {215                    return 0;216                }217                let d = brent(x);218                stack.push(d);219                stack.push(x / d);220            }221            found.sort_unstable();222            for w in found.windows(2) {223                if w[0] == w[1] {224                    return 0;225                }226            }227            parity ^= found.len() as u32 & 1;228        }229    }230    if parity & 1 == 0 {231        1232    } else {233        -1234    }235}