factor.rs

1/// Returns the greatest common divisor of two numbers by the Euclidean algorithm, zero for two zeroes.2pub fn gcd(a: usize, b: usize) -> usize {3    let (mut a, mut b) = (a, b);4    while b != 0 {5        (a, b) = (b, a % b);6    }7    a8}910/// Returns the least common multiple of two numbers, zero when either side is zero.11///12/// Wraps silently once the multiple passes the pointer width.13pub fn lcm(a: usize, b: usize) -> usize {14    if a == 0 || b == 0 {15        return 0;16    }17    a / gcd(a, b) * b18}1920/// Returns whether two numbers share no divisor above one.21pub fn coprime(a: usize, b: usize) -> bool {22    gcd(a, b) == 123}2425fn peel(prime: u64, rest: &mut u64, out: &mut Vec<(u64, u32)>) {26    let mut power = 0;27    while rest.is_multiple_of(prime) {28        *rest /= prime;29        power += 1;30    }31    if power > 0 {32        out.push((prime, power));33    }34}3536/// Returns the prime and exponent pairs of a wide number in ascending primes, by trial division on the six-step wheel.37///38/// ```39/// assert_eq!(mrlynum::factor::factorize_wide(999_999_999_989), vec![(999_999_999_989, 1)]);40/// ```41pub fn factorize_wide(number: u64) -> Vec<(u64, u32)> {42    let mut out = Vec::new();43    let mut rest = number;44    if rest < 2 {45        return out;46    }47    peel(2, &mut rest, &mut out);48    peel(3, &mut rest, &mut out);49    let mut step = 5;50    while step <= rest / step {51        peel(step, &mut rest, &mut out);52        peel(step + 2, &mut rest, &mut out);53        step += 6;54    }55    if rest > 1 {56        out.push((rest, 1));57    }58    out59}6061/// Returns the prime and exponent pairs of the number in ascending primes, by trial division on the six-step wheel.62///63/// ```64/// assert_eq!(mrlynum::factor::factorize(240), vec![(2, 4), (3, 1), (5, 1)]);65/// ```66pub fn factorize(number: usize) -> Vec<(usize, u32)> {67    factorize_wide(number as u64)68        .into_iter()69        .map(|(prime, power)| (prime as usize, power))70        .collect()71}7273/// Builds every divisor of the number from its factorization, ascending, empty for zero.74pub fn divisors(number: usize) -> Vec<usize> {75    if number == 0 {76        return Vec::new();77    }78    let mut out = vec![1];79    for (prime, power) in factorize(number) {80        let mut next = Vec::with_capacity(out.len() * (power as usize + 1));81        for &divisor in &out {82            let mut value = divisor;83            next.push(value);84            for _ in 0..power {85                value *= prime;86                next.push(value);87            }88        }89        out = next;90    }91    out.sort_unstable();92    out93}9495/// Returns the sum of every divisor of the number raised to the power, so power zero counts them.96///97/// Stays exact for a power of two or below, above which the sum can wrap past a hundred and twenty-eight bits.98pub fn sigma(number: usize, power: u32) -> u128 {99    divisors(number)100        .iter()101        .map(|&divisor| (divisor as u128).pow(power))102        .sum()103}104105/// Returns the sum of the proper divisors of the number, its divisor sum less itself, zero for zero and for one.106pub fn aliquot(number: usize) -> usize {107    (sigma(number, 1) - number as u128) as usize108}109110/// Returns the divisor sum with a periodic rhythm painted on each divisor, zero for zero and for an empty rhythm.111///112/// ```113/// assert_eq!(mrlynum::factor::twisted(5, &[0, 1, 0, -1]), 2);114/// ```115pub fn twisted(number: usize, rhythm: &[i8]) -> i64 {116    if rhythm.is_empty() {117        return 0;118    }119    divisors(number)120        .iter()121        .map(|&divisor| i64::from(rhythm[divisor % rhythm.len()]))122        .sum()123}124125/// Returns the Mobius value of the number: zero for zero or a squared factor, else minus one to the count of primes.126pub fn mobius(number: usize) -> i8 {127    if number == 0 {128        return 0;129    }130    let mut out = 1;131    for (_, power) in factorize(number) {132        if power > 1 {133            return 0;134        }135        out = -out;136    }137    out138}139140/// Sieves the Mobius values of zero through the limit in one pass.141pub fn mobius_sieve(limit: usize) -> Vec<i8> {142    let mut out = vec![1i8; limit + 1];143    let mut composite = vec![false; limit + 1];144    out[0] = 0;145    for p in 2..=limit {146        if composite[p] {147            continue;148        }149        for m in (p..=limit).step_by(p) {150            composite[m] = m != p;151            out[m] = -out[m];152        }153        if p <= limit / p {154            for m in (p * p..=limit).step_by(p * p) {155                out[m] = 0;156            }157        }158    }159    out160}161162/// Returns the Euler totient of the number from its factorization, zero for zero and one for one.163pub fn totient(number: usize) -> usize {164    if number == 0 {165        return 0;166    }167    let mut out = number;168    for (prime, _) in factorize(number) {169        out -= out / prime;170    }171    out172}173174/// Returns the radical of the number, the product of its distinct primes, zero for zero and one for one.175pub fn radical(number: usize) -> usize {176    if number == 0 {177        return 0;178    }179    factorize(number).iter().map(|&(prime, _)| prime).product()180}181182/// Returns whether no prime squares into the number, true for one and false for zero.183pub fn squarefree(number: usize) -> bool {184    number != 0 && factorize(number).iter().all(|&(_, power)| power < 2)185}186187#[cfg(test)]188mod tests {189    use super::*;190    use crate::lattice::totients;191    use crate::prime::rectangles;192    use crate::series::{chi4, chi8};193194    fn represented(m: usize, weight: usize) -> usize {195        let mut count = 0;196        for a in -20i64..=20 {197            for b in -20i64..=20 {198                if a * a + weight as i64 * b * b == m as i64 {199                    count += 1;200                }201            }202        }203        count204    }205206    #[test]207    fn factorize_pins_the_hand_checked_shapes() {208        assert_eq!(factorize(240), vec![(2, 4), (3, 1), (5, 1)]);209        assert!(factorize(1).is_empty());210        assert!(factorize(0).is_empty());211        assert_eq!(factorize(97), vec![(97, 1)]);212        assert_eq!(factorize(243), vec![(3, 5)]);213        assert_eq!(factorize(1_000_003), vec![(1_000_003, 1)]);214    }215216    #[test]217    fn a_factorization_multiplies_back_to_its_number() {218        for number in 2..=5_000usize {219            let parts = factorize(number);220            let rebuilt: usize = parts.iter().map(|&(p, e)| p.pow(e)).product();221            assert_eq!(rebuilt, number, "{number}");222            for pair in parts.windows(2) {223                assert!(pair[0].0 < pair[1].0, "{number}");224            }225        }226    }227228    #[test]229    fn the_divisor_count_is_sigma_zero_and_twice_the_rectangles_less_a_square() {230        assert!(divisors(0).is_empty());231        assert_eq!(divisors(1), vec![1]);232        assert_eq!(divisors(28), vec![1, 2, 4, 7, 14, 28]);233        for number in 1..=2_000usize {234            let count = divisors(number).len() as u128;235            assert_eq!(count, sigma(number, 0), "{number}");236            let root = number.isqrt();237            let square = u128::from(root * root == number);238            assert_eq!(239                count,240                2 * rectangles(number).len() as u128 - square,241                "{number}"242            );243        }244    }245246    #[test]247    fn sigma_pins_the_hand_checked_sums() {248        assert_eq!(sigma(6, 1), 12);249        assert_eq!(sigma(240, 0), 20);250        assert_eq!(sigma(1, 3), 1);251        assert_eq!(sigma(0, 1), 0);252    }253254    #[test]255    fn the_aliquot_sum_is_sigma_one_less_the_number() {256        assert_eq!(aliquot(0), 0);257        assert_eq!(aliquot(1), 0);258        assert_eq!(aliquot(6), 6);259        assert_eq!(aliquot(28), 28);260        assert_eq!(aliquot(12), 16);261        for number in 1..=2_000usize {262            assert_eq!(263                sigma(number, 1),264                aliquot(number) as u128 + number as u128,265                "{number}"266            );267        }268    }269270    #[test]271    fn the_mod_four_twist_counts_the_sums_of_two_squares() {272        let rhythm: Vec<i8> = (0..4).map(chi4).collect();273        for m in 1..=300usize {274            assert_eq!(represented(m, 1) as i64, 4 * twisted(m, &rhythm), "{m}");275        }276    }277278    #[test]279    fn the_mod_eight_twist_counts_a_square_plus_two_squares() {280        let rhythm: Vec<i8> = (0..8).map(chi8).collect();281        for m in 1..=300usize {282            assert_eq!(represented(m, 2) as i64, 2 * twisted(m, &rhythm), "{m}");283        }284    }285286    #[test]287    fn the_flat_twist_is_the_divisor_count() {288        assert_eq!(twisted(0, &[1]), 0);289        assert_eq!(twisted(12, &[]), 0);290        for number in 1..=1_000usize {291            assert_eq!(twisted(number, &[1]) as u128, sigma(number, 0), "{number}");292        }293    }294295    #[test]296    fn mobius_starts_the_known_sequence() {297        let start: Vec<i8> = (1..=10).map(mobius).collect();298        assert_eq!(start, vec![1, -1, -1, 0, -1, 1, -1, 0, 0, 1]);299        assert_eq!(mobius(0), 0);300        assert_eq!(mobius(240), 0);301        assert_eq!(mobius(30), -1);302    }303304    #[test]305    fn the_mobius_sieve_agrees_with_the_single_value() {306        let sieved = mobius_sieve(10_000);307        for (number, &value) in sieved.iter().enumerate() {308            assert_eq!(value, mobius(number), "{number}");309        }310    }311312    #[test]313    fn totient_agrees_with_the_lattice_sieve() {314        let sieved = totients(5_000);315        for (number, &value) in sieved.iter().enumerate() {316            assert_eq!(totient(number) as u64, value, "{number}");317        }318    }319320    #[test]321    fn the_radical_keeps_one_of_each_prime_and_equals_a_squarefree_number() {322        assert_eq!(radical(240), 30);323        assert_eq!(radical(0), 0);324        assert_eq!(radical(1), 1);325        assert!(squarefree(30));326        assert!(!squarefree(240));327        assert!(squarefree(1));328        assert!(!squarefree(0));329        for number in 1..=2_000usize {330            assert_eq!(squarefree(number), radical(number) == number, "{number}");331        }332    }333334    #[test]335    fn gcd_times_lcm_is_the_product_and_a_gcd_of_one_means_coprime() {336        assert_eq!(gcd(0, 0), 0);337        assert_eq!(lcm(0, 7), 0);338        for a in 1..=60usize {339            for b in 1..=60usize {340                assert_eq!(gcd(a, b) * lcm(a, b), a * b, "{a} {b}");341                assert_eq!(coprime(a, b), gcd(a, b) == 1, "{a} {b}");342            }343        }344    }345}