Compare commits

4 Commits

Author SHA1 Message Date
ddidderr dc82024ed0 [deps] cargo update 2025-05-29 11:33:43 +02:00
ddidderr fd36faae13 [code] remove unnecessary reference 2024-09-28 12:59:31 +02:00
ddidderr 490c53a280 [fix] but primes behind an Arc
- don't clone primes for every thread, instead put it behind an Arc
  (atomic reference counted)
- print the pure calculation time at the end of the program
- don't show thread debug prints
2024-09-28 12:57:29 +02:00
ddidderr 939e239b36 [feature] use multiple threads to calculate chunks in parallel 2024-09-12 20:35:23 +02:00
5 changed files with 167 additions and 502 deletions
Generated
+9 -170
View File
@@ -3,184 +3,23 @@
version = 4
[[package]]
name = "anstream"
version = "1.0.0"
name = "either"
version = "1.15.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "824a212faf96e9acacdbd09febd34438f8f711fb84e09a8916013cd7815ca28d"
dependencies = [
"anstyle",
"anstyle-parse",
"anstyle-query",
"anstyle-wincon",
"colorchoice",
"is_terminal_polyfill",
"utf8parse",
]
[[package]]
name = "anstyle"
version = "1.0.14"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "940b3a0ca603d1eade50a4846a2afffd5ef57a9feac2c0e2ec2e14f9ead76000"
[[package]]
name = "anstyle-parse"
version = "1.0.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "52ce7f38b242319f7cabaa6813055467063ecdc9d355bbb4ce0c68908cd8130e"
dependencies = [
"utf8parse",
]
[[package]]
name = "anstyle-query"
version = "1.1.5"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "40c48f72fd53cd289104fc64099abca73db4166ad86ea0b4341abe65af83dadc"
dependencies = [
"windows-sys",
]
[[package]]
name = "anstyle-wincon"
version = "3.0.11"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "291e6a250ff86cd4a820112fb8898808a366d8f9f58ce16d1f538353ad55747d"
dependencies = [
"anstyle",
"once_cell_polyfill",
"windows-sys",
]
[[package]]
name = "clap"
version = "4.6.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "1ddb117e43bbf7dacf0a4190fef4d345b9bad68dfc649cb349e7d17d28428e51"
dependencies = [
"clap_builder",
"clap_derive",
]
[[package]]
name = "clap_builder"
version = "4.6.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "714a53001bf66416adb0e2ef5ac857140e7dc3a0c48fb28b2f10762fc4b5069f"
dependencies = [
"anstream",
"anstyle",
"clap_lex",
"strsim",
]
[[package]]
name = "clap_derive"
version = "4.6.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "f2ce8604710f6733aa641a2b3731eaa1e8b3d9973d5e3565da11800813f997a9"
dependencies = [
"heck",
"proc-macro2",
"quote",
"syn",
]
[[package]]
name = "clap_lex"
version = "1.1.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "c8d4a3bb8b1e0c1050499d1815f5ab16d04f0959b233085fb31653fbfc9d98f9"
[[package]]
name = "colorchoice"
version = "1.0.5"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "1d07550c9036bf2ae0c684c4297d503f838287c83c53686d05370d0e139ae570"
checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719"
[[package]]
name = "hcn"
version = "1.0.0"
version = "1.0.0-multithread"
dependencies = [
"clap",
"itertools",
]
[[package]]
name = "heck"
version = "0.5.0"
name = "itertools"
version = "0.14.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea"
[[package]]
name = "is_terminal_polyfill"
version = "1.70.2"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "a6cb138bb79a146c1bd460005623e142ef0181e3d0219cb493e02f7d08a35695"
[[package]]
name = "once_cell_polyfill"
version = "1.70.2"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "384b8ab6d37215f3c5301a95a4accb5d64aa607f1fcb26a11b5303878451b4fe"
[[package]]
name = "proc-macro2"
version = "1.0.106"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "8fd00f0bb2e90d81d1044c2b32617f68fcb9fa3bb7640c23e9c748e53fb30934"
checksum = "2b192c782037fadd9cfa75548310488aabdbf3d2da73885b31bd0abd03351285"
dependencies = [
"unicode-ident",
]
[[package]]
name = "quote"
version = "1.0.45"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "41f2619966050689382d2b44f664f4bc593e129785a36d6ee376ddf37259b924"
dependencies = [
"proc-macro2",
]
[[package]]
name = "strsim"
version = "0.11.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "7da8b5736845d9f2fcb837ea5d9e2628564b3b043a70948a3f0b778838c5fb4f"
[[package]]
name = "syn"
version = "2.0.117"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "e665b8803e7b1d2a727f4023456bbbbe74da67099c585258af0ad9c5013b9b99"
dependencies = [
"proc-macro2",
"quote",
"unicode-ident",
]
[[package]]
name = "unicode-ident"
version = "1.0.24"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "e6e4313cd5fcd3dad5cafa179702e2b244f760991f45397d14d4ebf38247da75"
[[package]]
name = "utf8parse"
version = "0.2.2"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "06abde3611657adf66d383f00b093d7faecc7fa57071cce2578660c9f1010821"
[[package]]
name = "windows-link"
version = "0.2.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5"
[[package]]
name = "windows-sys"
version = "0.61.2"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "ae137229bcbd6cdf0f7b80a31df61766145077ddf49416a728b02cb3921ff3fc"
dependencies = [
"windows-link",
"either",
]
+3 -3
View File
@@ -1,10 +1,10 @@
[package]
name = "hcn"
version = "1.0.0"
edition = "2024"
version = "1.0.0-multithread"
edition = "2021"
[dependencies]
clap = { version = "4", features = ["derive"] }
itertools = "0.14"
[lints.rust]
unsafe_code = "forbid"
-72
View File
@@ -1,72 +0,0 @@
# hcn
Find highly composite numbers up to a search limit.
```bash
cargo run --release -- [SEARCH_LIMIT]
```
If no limit is given, the default is `1_000_000_000`.
Use `--help` to see the generated command-line help and `--version` to print
the package version.
## Structure
- `src/main.rs`
- Command-line interface: parses the optional `SEARCH_LIMIT` argument with
`clap` and rejects non-positive limits before the search starts.
- Search setup: defines the default limit and the fixed prime list used for
candidate generation.
- Candidate generation: recursively builds only exponent sequences that can
produce highly composite record candidates.
- Record filtering: sorts candidates and prints each new divisor-count record.
## Algorithm
A number with prime factorization
```text
n = 2^a * 3^b * 5^c * ...
```
has
```text
(a + 1) * (b + 1) * (c + 1) * ...
```
divisors.
The exact primes do not affect the divisor count; only the exponents do.
For a fixed exponent list, the smallest possible number is made by putting the
largest exponent on the smallest prime:
```text
a >= b >= c >= ...
```
If a larger exponent appears on a larger prime, swapping those two exponents
keeps the divisor count unchanged but makes the number smaller.
So a highly composite number cannot have that shape. If it did, the smaller
swapped number would already have the same number of divisors, so the original
number would not be the first record.
That means record candidates are only numbers like:
```text
2^a * 3^b * 5^c * ... where a >= b >= c >= ...
```
All other numbers are redundant for finding new records.
The program:
1. Recursively generates numbers from the first primes.
2. Only tries exponent sequences where each exponent is no larger than the
previous one.
3. Stops a branch as soon as the number exceeds the search limit.
4. Computes the divisor count directly from the exponents.
5. Sorts candidates by number and prints each new divisor-count record.
This avoids factoring every number in the range.
+95 -247
View File
@@ -1,278 +1,126 @@
use std::{num::NonZeroU64, time::Instant};
mod sieve;
use clap::Parser;
use std::{
collections::HashMap,
env,
sync::Arc,
thread::{self, available_parallelism},
time::Instant,
};
const DEFAULT_SEARCH_LIMIT: u64 = 1_000_000_000;
const FIRST_PRIMES: &[u64] = &[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31];
use itertools::Itertools as _;
use sieve::get_primes;
#[derive(Debug, Parser)]
#[command(version, about = "Find highly composite numbers up to a search limit")]
struct Cli {
#[arg(value_name = "SEARCH_LIMIT", default_value_t = NonZeroU64::new(DEFAULT_SEARCH_LIMIT).expect("default search limit is non-zero"))]
search_limit: NonZeroU64,
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 % prime == 0 {
nr /= prime;
count += 1;
}
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
struct Candidate {
nr: u64,
divisors: u64,
count
}
fn max_exponent(limit: u64) -> u64 {
let mut exponent = 0;
let mut nr = 1_u64;
fn prime_factors(nr: u64, primes: &[u64]) -> u64 {
let mut num_teilers = 1;
while let Some(next_nr) = nr.checked_mul(2) {
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;
}
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<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,
);
let count = count_prime_factor_exponent(*prime, nr);
if count > 0 {
num_teilers *= count + 1;
}
}
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;
num_teilers
}
primorial = next_primorial;
fn the_thread(mut start: u64, end: u64, primes: &[u64]) -> HashMap<u64, u64> {
let mut max_teilers = 0;
let mut max_teilers_hashmap = HashMap::new();
let mut step_size = 10;
if start == 0 {
start = 2;
step_size = 2;
}
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;
for nr in (start..=end).step_by(step_size) {
let teilers = prime_factors(nr, primes);
if teilers > max_teilers {
max_teilers = teilers;
max_teilers_hashmap.insert(nr, teilers);
}
}
records
max_teilers_hashmap
}
fn calculate_chunk_bounds(i: usize, num_threads: usize, max_nr: usize) -> (usize, usize) {
let chunk_size = max_nr / num_threads;
let start = i * chunk_size;
let end = ((i + 1) * chunk_size).min(max_nr);
(start, end)
}
fn main() {
let search_limit = Cli::parse().search_limit.get();
let max_nr = env::args()
.nth(1)
.expect("Usage: hcn <max_nr>")
.parse::<usize>()
.expect("Invalid max_nr");
let start = Instant::now();
let records = hcn_records(search_limit);
println!("Precalculating primes...");
let primes = get_primes(MAX_SIEVED_PRIMES);
println!("{} primes. Took {:?}", primes.len(), start.elapsed());
for Candidate { nr, divisors } in records {
println!("{nr}: {divisors} ({:?} since start)", start.elapsed());
let primes = Arc::new(primes);
#[allow(clippy::unwrap_used)]
let num_threads = available_parallelism().unwrap().get();
let mut threads = Vec::with_capacity(num_threads);
let now = Instant::now();
for i in 0..num_threads {
let (start, end) = calculate_chunk_bounds(i, num_threads, max_nr);
let primes = primes.clone();
threads.push(thread::spawn(move || {
the_thread(start as u64, end as u64, &primes)
}));
}
let mut results = HashMap::new();
for thread in threads {
let thread_result = thread.join().expect("Thread failed to join (panicked?)");
results.extend(thread_result);
}
let mut max_teilers = 0;
for (nr, teilers) in results.into_iter().sorted() {
if teilers > max_teilers {
println!("{nr}: {teilers}");
max_teilers = teilers;
}
}
#[cfg(test)]
mod tests {
use clap::Parser;
use super::{Candidate, Cli, DEFAULT_SEARCH_LIMIT, 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));
}
#[test]
fn parses_default_search_limit_without_arg() {
assert_eq!(
Cli::try_parse_from(["hcn"])
.expect("CLI should parse without an explicit search limit")
.search_limit
.get(),
DEFAULT_SEARCH_LIMIT
);
}
#[test]
fn parses_search_limit_arg() {
assert_eq!(
Cli::try_parse_from(["hcn", "42000"])
.expect("CLI should parse an explicit search limit")
.search_limit
.get(),
42_000
);
}
#[test]
fn rejects_invalid_search_limit_arg() {
assert!(Cli::try_parse_from(["hcn", "nope"]).is_err());
assert!(Cli::try_parse_from(["hcn", "0"]).is_err());
assert!(Cli::try_parse_from(["hcn", "100", "200"]).is_err());
}
println!("Took {:?}", now.elapsed());
}
+50
View File
@@ -0,0 +1,50 @@
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);
}