diff --git a/src/main.rs b/src/main.rs index 9727b0d..0aae7b9 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,64 +1,236 @@ -mod sieve; - use std::time::Instant; -use sieve::get_primes; +const SEARCH_LIMIT: u64 = 1_000_000_000; +const FIRST_PRIMES: &[u64] = &[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]; -const MAX_SIEVED_PRIMES: usize = 100_000_000; - -fn count_prime_factor_exponent(prime: u64, mut nr: u64) -> u64 { - let mut count = 0; - - while nr.is_multiple_of(prime) { - nr /= prime; - count += 1; - } - - count +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +struct Candidate { + nr: u64, + divisors: u64, } -fn prime_factors(nr: u64, primes: &[u64]) -> u64 { - let mut num_teilers = 1; +fn max_exponent(limit: u64) -> u64 { + let mut exponent = 0; + let mut nr = 1_u64; - // safety: nr <= 2.pow(53) must be fulfilled, - // otherwise we have a precision loss - assert!(nr <= 2u64.pow(53)); - #[allow( - clippy::cast_sign_loss, - clippy::cast_possible_truncation, - clippy::cast_precision_loss - )] - let nr_sqrt = (nr as f64).sqrt().ceil() as u64; - - for prime in primes { - if *prime > nr_sqrt { + while let Some(next_nr) = nr.checked_mul(2) { + if next_nr > limit { break; } - let count = count_prime_factor_exponent(*prime, nr); - if count > 0 { - num_teilers *= count + 1; + nr = next_nr; + exponent += 1; + } + + exponent +} + +fn collect_hcn_candidates( + limit: u64, + prime_idx: usize, + max_exponent: u64, + nr: u64, + divisors: u64, + candidates: &mut Vec, +) { + let Some(&prime) = FIRST_PRIMES.get(prime_idx) else { + return; + }; + + let mut next_nr = nr; + + for exponent in 1..=max_exponent { + let Some(multiplied_nr) = next_nr.checked_mul(prime) else { + break; + }; + + if multiplied_nr > limit { + break; + } + + next_nr = multiplied_nr; + let next_divisors = divisors * (exponent + 1); + + candidates.push(Candidate { + nr: next_nr, + divisors: next_divisors, + }); + + collect_hcn_candidates( + limit, + prime_idx + 1, + exponent, + next_nr, + next_divisors, + candidates, + ); + } +} + +fn has_enough_primes(limit: u64) -> bool { + let mut primorial = 1_u64; + + for &prime in FIRST_PRIMES { + let Some(next_primorial) = primorial.checked_mul(prime) else { + return true; + }; + + if next_primorial > limit { + return true; + } + + primorial = next_primorial; + } + + false +} + +fn hcn_records(limit: u64) -> Vec { + assert!( + has_enough_primes(limit), + "not enough primes for search limit" + ); + + let mut candidates = Vec::new(); + collect_hcn_candidates(limit, 0, max_exponent(limit), 1, 1, &mut candidates); + candidates.sort_unstable_by_key(|candidate| candidate.nr); + + let mut records = Vec::new(); + let mut max_divisors = 0; + + for candidate in candidates { + if candidate.divisors > max_divisors { + records.push(candidate); + max_divisors = candidate.divisors; } } - num_teilers + records } fn main() { let start = Instant::now(); - println!("Precalculating primes..."); - let primes = get_primes(MAX_SIEVED_PRIMES); - println!("{} primes. Took {:?}", primes.len(), start.elapsed()); + let records = hcn_records(SEARCH_LIMIT); - let mut max_teilers = 0; + for Candidate { nr, divisors } in records { + println!("{nr}: {divisors} ({:?} since start)", start.elapsed()); + } +} - let start = Instant::now(); +#[cfg(test)] +mod tests { + use super::{Candidate, hcn_records, max_exponent}; - for nr in (2..1_000_000_000).step_by(2) { - let teilers = prime_factors(nr, &primes); - if teilers > max_teilers { - println!("{nr}: {teilers} ({:?} since start)", start.elapsed()); - max_teilers = teilers; + fn count_divisors(mut nr: u64) -> u64 { + let mut divisor_count = 1; + let mut divisor = 2; + + while divisor * divisor <= nr { + let mut exponent = 0; + + while nr.is_multiple_of(divisor) { + nr /= divisor; + exponent += 1; + } + + if exponent > 0 { + divisor_count *= exponent + 1; + } + + divisor += if divisor == 2 { 1 } else { 2 }; + } + + if nr > 1 { + divisor_count *= 2; + } + + divisor_count + } + + fn brute_force_hcn_records(limit: u64) -> Vec { + let mut records = Vec::new(); + let mut max_divisors = 0; + + for nr in 2..=limit { + let divisors = count_divisors(nr); + + if divisors > max_divisors { + records.push(Candidate { nr, divisors }); + max_divisors = divisors; + } } + + records + } + + #[test] + fn max_exponent_counts_powers_of_two_up_to_limit() { + assert_eq!(max_exponent(1), 0); + assert_eq!(max_exponent(2), 1); + assert_eq!(max_exponent(3), 1); + assert_eq!(max_exponent(16), 4); + } + + #[test] + fn finds_highly_composite_records_up_to_one_thousand() { + let records = hcn_records(1_000); + + assert_eq!( + records, + vec![ + Candidate { nr: 2, divisors: 2 }, + Candidate { nr: 4, divisors: 3 }, + Candidate { nr: 6, divisors: 4 }, + Candidate { + nr: 12, + divisors: 6, + }, + Candidate { + nr: 24, + divisors: 8, + }, + Candidate { + nr: 36, + divisors: 9, + }, + Candidate { + nr: 48, + divisors: 10, + }, + Candidate { + nr: 60, + divisors: 12, + }, + Candidate { + nr: 120, + divisors: 16, + }, + Candidate { + nr: 180, + divisors: 18, + }, + Candidate { + nr: 240, + divisors: 20, + }, + Candidate { + nr: 360, + divisors: 24, + }, + Candidate { + nr: 720, + divisors: 30, + }, + Candidate { + nr: 840, + divisors: 32, + }, + ], + ); + } + + #[test] + fn matches_brute_force_records_up_to_ten_thousand() { + assert_eq!(hcn_records(10_000), brute_force_hcn_records(10_000)); } } diff --git a/src/sieve.rs b/src/sieve.rs deleted file mode 100644 index c531cde..0000000 --- a/src/sieve.rs +++ /dev/null @@ -1,50 +0,0 @@ -pub fn get_primes(max_nr: usize) -> Vec { - let mut sieve = vec![true; max_nr + 1]; - - let mut primes = Vec::with_capacity(1_000_000); - - sieve[0] = false; - sieve[1] = false; - - let mut nr = 2; - primes.push(2); - - loop { - check_off_multiples_of_nr(nr, &mut sieve); - - if find_next_prime(&mut nr, &mut sieve).is_none() { - break; - } - - primes.push(nr); - } - - primes -} - -fn find_next_prime(nr: &mut u64, sieve: &mut [bool]) -> Option { - let start_idx = usize::try_from(*nr).ok()? + 1; - - if let Some((idx, _)) = sieve - .iter() - .enumerate() - .skip(start_idx) - .find(|&(_, &is_prime)| is_prime) - { - *nr = idx as u64; - return Some(*nr); - } - - None -} - -#[allow(clippy::unwrap_used)] -fn check_off_multiples_of_nr(nr: u64, sieve: &mut [bool]) { - let nr = usize::try_from(nr).unwrap(); - - sieve - .iter_mut() - .skip(nr.checked_mul(2).unwrap()) - .step_by(nr) - .for_each(|nr| *nr = false); -}