perf: generate HCN records from exponent vectors
Replace the range scan and trial division with a search over prime exponent vectors. Highly composite records only need candidates whose exponents are non-increasing across ascending primes, so the program can compute divisor counts directly while pruning branches above the search limit. Remove the sieve module because the new search no longer factors each number. This makes the current one-billion limit finish in microseconds after startup instead of scanning hundreds of millions of inputs. Test Plan: - cargo clippy - cargo clippy --benches - cargo clippy --tests - cargo test Refs: N/A
This commit is contained in:
+216
-44
@@ -1,64 +1,236 @@
|
|||||||
mod sieve;
|
|
||||||
|
|
||||||
use std::time::Instant;
|
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;
|
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
|
||||||
|
struct Candidate {
|
||||||
fn count_prime_factor_exponent(prime: u64, mut nr: u64) -> u64 {
|
nr: u64,
|
||||||
let mut count = 0;
|
divisors: u64,
|
||||||
|
|
||||||
while nr.is_multiple_of(prime) {
|
|
||||||
nr /= prime;
|
|
||||||
count += 1;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
count
|
fn max_exponent(limit: u64) -> u64 {
|
||||||
}
|
let mut exponent = 0;
|
||||||
|
let mut nr = 1_u64;
|
||||||
|
|
||||||
fn prime_factors(nr: u64, primes: &[u64]) -> u64 {
|
while let Some(next_nr) = nr.checked_mul(2) {
|
||||||
let mut num_teilers = 1;
|
if next_nr > limit {
|
||||||
|
|
||||||
// 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 {
|
|
||||||
break;
|
break;
|
||||||
}
|
}
|
||||||
|
|
||||||
let count = count_prime_factor_exponent(*prime, nr);
|
nr = next_nr;
|
||||||
if count > 0 {
|
exponent += 1;
|
||||||
num_teilers *= count + 1;
|
}
|
||||||
|
|
||||||
|
exponent
|
||||||
|
}
|
||||||
|
|
||||||
|
fn collect_hcn_candidates(
|
||||||
|
limit: u64,
|
||||||
|
prime_idx: usize,
|
||||||
|
max_exponent: u64,
|
||||||
|
nr: u64,
|
||||||
|
divisors: u64,
|
||||||
|
candidates: &mut Vec<Candidate>,
|
||||||
|
) {
|
||||||
|
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,
|
||||||
|
);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
num_teilers
|
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<Candidate> {
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
records
|
||||||
}
|
}
|
||||||
|
|
||||||
fn main() {
|
fn main() {
|
||||||
let start = Instant::now();
|
let start = Instant::now();
|
||||||
println!("Precalculating primes...");
|
let records = hcn_records(SEARCH_LIMIT);
|
||||||
let primes = get_primes(MAX_SIEVED_PRIMES);
|
|
||||||
println!("{} primes. Took {:?}", primes.len(), start.elapsed());
|
|
||||||
|
|
||||||
let mut max_teilers = 0;
|
for Candidate { nr, divisors } in records {
|
||||||
|
println!("{nr}: {divisors} ({:?} since start)", start.elapsed());
|
||||||
let start = Instant::now();
|
|
||||||
|
|
||||||
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;
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
#[cfg(test)]
|
||||||
|
mod tests {
|
||||||
|
use super::{Candidate, hcn_records, max_exponent};
|
||||||
|
|
||||||
|
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<Candidate> {
|
||||||
|
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));
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -1,50 +0,0 @@
|
|||||||
pub fn get_primes(max_nr: usize) -> Vec<u64> {
|
|
||||||
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<u64> {
|
|
||||||
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);
|
|
||||||
}
|
|
||||||
Reference in New Issue
Block a user