From ab19c870e132431f2bb86639e388bdf97848294e Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Tue, 15 Oct 2019 23:25:33 +0100 Subject: [PATCH 01/36] Commit new crate --- Cargo.toml | 1 + primal-count/Cargo.toml | 9 +++++++++ primal-count/src/main.rs | 3 +++ 3 files changed, 13 insertions(+) create mode 100644 primal-count/Cargo.toml create mode 100644 primal-count/src/main.rs diff --git a/Cargo.toml b/Cargo.toml index 0dd1697b..686d3885 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -42,5 +42,6 @@ members = [ "primal-estimate", "primal-sieve", "primal-slowsieve", + "primal-count", "generators", ] diff --git a/primal-count/Cargo.toml b/primal-count/Cargo.toml new file mode 100644 index 00000000..91108834 --- /dev/null +++ b/primal-count/Cargo.toml @@ -0,0 +1,9 @@ +[package] +name = "primal-count" +version = "0.1.0" +authors = ["Kaimyn Chapman "] +edition = "2018" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] diff --git a/primal-count/src/main.rs b/primal-count/src/main.rs new file mode 100644 index 00000000..e7a11a96 --- /dev/null +++ b/primal-count/src/main.rs @@ -0,0 +1,3 @@ +fn main() { + println!("Hello, world!"); +} From d846377dbaa0f3925c38c2cbcebfe9d1b4e02324 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Tue, 15 Oct 2019 23:31:07 +0100 Subject: [PATCH 02/36] Add raw script --- primal-count/src/raw_script.rs | 307 +++++++++++++++++++++++++++++++++ 1 file changed, 307 insertions(+) create mode 100644 primal-count/src/raw_script.rs diff --git a/primal-count/src/raw_script.rs b/primal-count/src/raw_script.rs new file mode 100644 index 00000000..5617ef68 --- /dev/null +++ b/primal-count/src/raw_script.rs @@ -0,0 +1,307 @@ +use std::io; +use std::collections::HashMap; +use std::cmp::Ordering::*; + +fn integer_square_root(value: usize) -> usize { + // Returns the largest integer at least sqrt(n) using Newton's method + let mut x = value; + let mut y = (x + 1) >> 1; + while y < x { + x = y; + y = (x + value / x) >> 1; + } + return x; +} + +fn integer_cubic_root(value: usize) -> usize { + // Returns the largest integer at least cbrt(n) using Newton's method + let mut x = value; + let mut y = (2 * value + 1) / 3; // Plus 1 only to protect against division by 0 later - not because it's right + while y < x { + x = y; + y = (2 * x + value / x / x) / 3; // Weird divide to protect against overflow + } + return x; +} + +fn integer_quartic_root(value: usize) -> usize { + // Returns the largest integer at least n^(1/4) using our fast integer sqrt + return integer_square_root(integer_square_root(value)); +} + +fn create_prime_array(bound: usize) -> Vec { + // Returns an array of primes up to and including bound + let mut vec = Vec::new(); + let max_value_to_check = integer_square_root(bound); + let mut sieve = vec![true; bound + 1]; + sieve[0] = false; + sieve[1] = false; + for value in 2..bound+1 { + if sieve[value] { + // Add the prime to our vector + vec.push(value); + + // Filter primes more if we're below the bound + if value <= max_value_to_check { + // Anything smaller than this already filtered + let mut idx = value * value; + while idx <= bound { + sieve[idx] = false; + idx += value; + } + } + } + } + return vec; +} + +fn meissel_fn_small(m: usize, n: usize, prime_array: &Vec) -> usize { + // The number of numbers <= m that are coprime to the first n prime numbers. + if n == 0 { + return m; + } else { + return meissel_fn_small(m, n-1, &prime_array) - meissel_fn_small(m / prime_array[n-1], n-1, &prime_array); + } +} + +// def _meissel_function_large(self, m, n): +// """The number of numbers <= m that are coprime to the first n prime numbers. +// Run for larger values where repeating isn't going to happen often. Use for n > 1000 or so""" +// m = int(m) +// # if m <= 10000: return meis104[n][m] +// if n == 0: +// return m +// result = self.meissel_memoize[n].get(m, None) +// if result is None: +// result = 0 +// primes_gen = (p for p in self.prime_array[:n:-1]) +// stacks = defaultdict(lambda: defaultdict(int)) +// stacks[n][m] = 1 +// for N in range(n, 0, -1): +// prime_dividing = next(primes_gen) +// for M in stacks[N]: +// # Marginal speed improvement? +// # if M <= 10000 and n<1229: +// # res += meis104[N][M]*stacks[N][M] +// # continue +// stacks[N-1][M] += stacks[N][M] +// stacks[N-1][M//prime_dividing] -= stacks[N][M] +// del stacks[N] +// for M in stacks[0]: +// result += M*stacks[0][M] +// return result +fn meissel_fn_large(m: usize, n: usize, prime_array: &Vec) -> usize { + meissel_fn_small(m, n, &prime_array) +} + +fn meissel_memoized(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { + if n == 0 { + return m; + } + match meissel_cache.get(&(m, n)).map(|entry| entry.clone()){ + Some(result) => result, + None => { + let value = meissel_memoized(m, n-1, &prime_array, meissel_cache) + - meissel_memoized(m / prime_array[n-1], n-1, &prime_array, meissel_cache); + meissel_cache.insert((m, n), value); + return value; + } + } +} + +fn num_primes_less_than(bound: usize) -> isize { + // Initialise prime array: + if bound < 1000 { + return create_prime_array(bound).len() as isize; + } + let sqrt_bound = integer_square_root(bound); + let primes = create_prime_array(sqrt_bound); + + let a = num_primes_less_than(integer_quartic_root(bound)); + let c = num_primes_less_than(integer_cubic_root(bound)); + let b = num_primes_less_than(sqrt_bound); + println!("a,b,c={},{},{}", a,b,c); + let mut result = meissel_fn_small(bound, a as usize, &primes) as isize + ((b + a - 2) * (b - a + 1)) / 2; + println!("result={}", result); + + for i in a..b { + let ith_prime = primes[i as usize]; + result -= num_primes_less_than(bound / ith_prime); + if i < c { + let bi = num_primes_less_than(integer_square_root(bound / ith_prime)); + for j in i..bi { + let jth_prime = primes[j as usize]; + result += j - num_primes_less_than(bound / ith_prime / jth_prime); + } + } + } + + return result +} + +fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: &mut HashMap, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { + // Initialise prime array: + match prime_cache.get(&bound).map(|entry| entry.clone()){ + Some(value) => value, + None => { // The meat of the function + if bound < 2 { + return 0; + } else if bound <= primes[primes.len()-1] { + let result = match primes.binary_search(&bound) + { + Ok(idx) => idx+1, + Err(idx) => idx, + }; + let result = result; + prime_cache.insert(bound, result); + return result; + } + + let sqrt_bound = integer_square_root(bound); + + let a = num_primes_less_than_memoized(integer_quartic_root(bound), primes, prime_cache, meissel_cache); + let c = num_primes_less_than_memoized(integer_cubic_root(bound), primes, prime_cache, meissel_cache); + let b = num_primes_less_than_memoized(sqrt_bound, primes, prime_cache, meissel_cache); + + // Issues with underflow here if b + a < 2 + let mut result = meissel_memoized(bound, a, &primes, meissel_cache) + ((b + a - 2) * (b - a + 1)) / 2; + + for i in a..b { + let ith_prime = primes[i]; + result -= num_primes_less_than_memoized(bound / ith_prime, primes, prime_cache, meissel_cache); + if i < c { + let bi = num_primes_less_than_memoized(integer_square_root(bound / ith_prime), primes, prime_cache, meissel_cache); + for j in i..bi { + let jth_prime = primes[j]; + result -= (num_primes_less_than_memoized(bound / ith_prime / jth_prime, primes, prime_cache, meissel_cache) - j); + } + } + } + + // Caching and shit + prime_cache.insert(bound, result); + return result; + } + } +} + +// Top level function +fn num_primes_less_than_main(bound: usize) -> usize { + let primes = create_prime_array(integer_square_root(bound)); + let mut value_cache = HashMap::new(); + + // Insert some basic primes - this is mainly to deal with underflow issues later + value_cache.insert(0, 0); + value_cache.insert(1, 0); + value_cache.insert(2, 1); + value_cache.insert(3, 2); + value_cache.insert(4, 2); + value_cache.insert(5, 3); + value_cache.insert(6, 3); + value_cache.insert(7, 4); + value_cache.insert(8, 4); + let mut meissel_cache = HashMap::new(); + println!("Primes initialised..."); + let value = num_primes_less_than_memoized(bound, &primes, &mut value_cache, &mut meissel_cache); + println!("Value = {}", value_cache[&bound]); + return value; +} + + + +fn main() { + // println!("How many primes below: "); + // let mut input = String::new(); + + // io::stdin().read_line(&mut input) + // .expect("Failed to read line"); + // let input: usize = input.trim().parse() + // .expect("Couldn't turn value into an integer..."); + + // println!("There are {} primes below {}", create_prime_array(input).len(), input); + // println!("There are {} primes below {}", num_primes_less_than_main(input), input); + println!("There are {} primes", num_primes_less_than_main(10_000_000_000)); + + +} + +#[test] +fn test_int_sqrt() { + assert_eq!(integer_square_root(0), 0); + assert_eq!(integer_square_root(1), 1); + assert_eq!(integer_square_root(16), 4); + assert_eq!(integer_square_root(17), 4); + assert_eq!(integer_square_root(24), 4); + assert_eq!(integer_square_root(587 * 587 - 1), 586); +} + +#[test] +fn test_int_cbrt() { + assert_eq!(integer_cubic_root(0), 0); + assert_eq!(integer_cubic_root(1), 1); + assert_eq!(integer_cubic_root(26), 2); + assert_eq!(integer_cubic_root(27), 3); + assert_eq!(integer_cubic_root(28), 3); + assert_eq!(integer_cubic_root(587 * 587 * 587 - 1), 586); +} + +#[test] +fn test_int_quartic_root() { + assert_eq!(integer_quartic_root(0), 0); + assert_eq!(integer_quartic_root(1), 1); + assert_eq!(integer_quartic_root(15), 1); + assert_eq!(integer_quartic_root(16), 2); + assert_eq!(integer_quartic_root(17), 2); + assert_eq!(integer_quartic_root(587 * 587 * 587 * 587 - 1), 586); + assert_eq!(integer_quartic_root(587 * 587 * 587 * 587 + 1), 587); +} + +#[test] +fn test_prime_array() { + let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; + assert_eq!(create_prime_array(19), prime_array); + assert_eq!(create_prime_array(20), prime_array); +} + +#[test] +fn test_meissel_fn_small() { + let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; + assert_eq!(meissel_fn_small(30, 8, &prime_array), 3); + assert_eq!(meissel_fn_small(100, 1, &prime_array), 50); +} + +#[test] +fn test_meissel_fn_large() { + let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; + assert_eq!(meissel_fn_large(30, 8, &prime_array), 3); + assert_eq!(meissel_fn_large(100, 1, &prime_array), 50); +} + + +#[test] +fn test_meissel_fn_memoized() { + let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; + let mut cache = HashMap::new(); + assert_eq!(meissel_memoized(30, 8, &prime_array, &mut cache), 3); + assert_eq!(meissel_memoized(100, 1, &prime_array, &mut cache), 50); +} + +#[test] +fn test_num_primes_less_than() { + assert_eq!(num_primes_less_than(7), 4); + assert_eq!(num_primes_less_than(100), 25); + assert_eq!(num_primes_less_than(2143), 324); + assert_eq!(num_primes_less_than(1_000_000), 78_498); + // assert_eq!(num_primes_less_than(1_000_000_000), 50_847_534); +} + +#[test] +fn test_num_primes_less_than_main() { + assert_eq!(num_primes_less_than_main(7), 4); + assert_eq!(num_primes_less_than_main(100), 25); + assert_eq!(num_primes_less_than_main(2143), 324); + assert_eq!(num_primes_less_than_main(1_000_000), 78_498); + assert_eq!(num_primes_less_than_main(1_000_000_000), 50_847_534); + // assert_eq!(num_primes_less_than_main(1_000_000_000_000), 37_607_912_018); + // assert_eq!(num_primes_less_than_main(1_000_000_000_000_000), 29_844_570_422_669); +} \ No newline at end of file From 83a79c7e35788104a3dfdea26f17cb5c4ccee4cf Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Tue, 15 Oct 2019 23:48:28 +0100 Subject: [PATCH 03/36] Better naming convention --- primal-count/src/lib.rs | 7 +++++ primal-count/src/main.rs | 3 -- .../src/{raw_script.rs => prime_count.rs} | 30 ++----------------- primal-count/src/util.rs | 26 ++++++++++++++++ 4 files changed, 36 insertions(+), 30 deletions(-) create mode 100644 primal-count/src/lib.rs delete mode 100644 primal-count/src/main.rs rename primal-count/src/{raw_script.rs => prime_count.rs} (92%) create mode 100644 primal-count/src/util.rs diff --git a/primal-count/src/lib.rs b/primal-count/src/lib.rs new file mode 100644 index 00000000..31e1bb20 --- /dev/null +++ b/primal-count/src/lib.rs @@ -0,0 +1,7 @@ +#[cfg(test)] +mod tests { + #[test] + fn it_works() { + assert_eq!(2 + 2, 4); + } +} diff --git a/primal-count/src/main.rs b/primal-count/src/main.rs deleted file mode 100644 index e7a11a96..00000000 --- a/primal-count/src/main.rs +++ /dev/null @@ -1,3 +0,0 @@ -fn main() { - println!("Hello, world!"); -} diff --git a/primal-count/src/raw_script.rs b/primal-count/src/prime_count.rs similarity index 92% rename from primal-count/src/raw_script.rs rename to primal-count/src/prime_count.rs index 5617ef68..34c52be5 100644 --- a/primal-count/src/raw_script.rs +++ b/primal-count/src/prime_count.rs @@ -1,34 +1,10 @@ +use crate::util::integer_square_root; +use crate::util::integer_cubic_root; +use crate::util::integer_quartic_root; use std::io; use std::collections::HashMap; use std::cmp::Ordering::*; -fn integer_square_root(value: usize) -> usize { - // Returns the largest integer at least sqrt(n) using Newton's method - let mut x = value; - let mut y = (x + 1) >> 1; - while y < x { - x = y; - y = (x + value / x) >> 1; - } - return x; -} - -fn integer_cubic_root(value: usize) -> usize { - // Returns the largest integer at least cbrt(n) using Newton's method - let mut x = value; - let mut y = (2 * value + 1) / 3; // Plus 1 only to protect against division by 0 later - not because it's right - while y < x { - x = y; - y = (2 * x + value / x / x) / 3; // Weird divide to protect against overflow - } - return x; -} - -fn integer_quartic_root(value: usize) -> usize { - // Returns the largest integer at least n^(1/4) using our fast integer sqrt - return integer_square_root(integer_square_root(value)); -} - fn create_prime_array(bound: usize) -> Vec { // Returns an array of primes up to and including bound let mut vec = Vec::new(); diff --git a/primal-count/src/util.rs b/primal-count/src/util.rs new file mode 100644 index 00000000..2cf31616 --- /dev/null +++ b/primal-count/src/util.rs @@ -0,0 +1,26 @@ +fn integer_square_root(value: usize) -> usize { + // Returns the largest integer at least sqrt(n) using Newton's method + let mut x = value; + let mut y = (x + 1) >> 1; + while y < x { + x = y; + y = (x + value / x) >> 1; + } + return x; +} + +fn integer_cubic_root(value: usize) -> usize { + // Returns the largest integer at least cbrt(n) using Newton's method + let mut x = value; + let mut y = (2 * value + 1) / 3; // Plus 1 only to protect against division by 0 later - not because it's right + while y < x { + x = y; + y = (2 * x + value / x / x) / 3; // Weird divide to protect against overflow + } + return x; +} + +fn integer_quartic_root(value: usize) -> usize { + // Returns the largest integer at least n^(1/4) using our fast integer sqrt + return integer_square_root(integer_square_root(value)); +} \ No newline at end of file From 3bf2d94119bc2834904fbf669071637d49e33f90 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Tue, 15 Oct 2019 23:56:20 +0100 Subject: [PATCH 04/36] Get code building --- primal-count/src/lib.rs | 86 ++++++++++++++++- primal-count/src/prime_count.rs | 157 ++------------------------------ primal-count/src/util.rs | 6 +- 3 files changed, 94 insertions(+), 155 deletions(-) diff --git a/primal-count/src/lib.rs b/primal-count/src/lib.rs index 31e1bb20..12ae7fde 100644 --- a/primal-count/src/lib.rs +++ b/primal-count/src/lib.rs @@ -1,7 +1,87 @@ #[cfg(test)] -mod tests { +mod tests { #[test] - fn it_works() { - assert_eq!(2 + 2, 4); + fn test_int_sqrt() { + assert_eq!(integer_square_root(0), 0); + assert_eq!(integer_square_root(1), 1); + assert_eq!(integer_square_root(16), 4); + assert_eq!(integer_square_root(17), 4); + assert_eq!(integer_square_root(24), 4); + assert_eq!(integer_square_root(587 * 587 - 1), 586); + } + + #[test] + fn test_int_cbrt() { + assert_eq!(integer_cubic_root(0), 0); + assert_eq!(integer_cubic_root(1), 1); + assert_eq!(integer_cubic_root(26), 2); + assert_eq!(integer_cubic_root(27), 3); + assert_eq!(integer_cubic_root(28), 3); + assert_eq!(integer_cubic_root(587 * 587 * 587 - 1), 586); + } + + #[test] + fn test_int_quartic_root() { + assert_eq!(integer_quartic_root(0), 0); + assert_eq!(integer_quartic_root(1), 1); + assert_eq!(integer_quartic_root(15), 1); + assert_eq!(integer_quartic_root(16), 2); + assert_eq!(integer_quartic_root(17), 2); + assert_eq!(integer_quartic_root(587 * 587 * 587 * 587 - 1), 586); + assert_eq!(integer_quartic_root(587 * 587 * 587 * 587 + 1), 587); + } + + #[test] + fn test_prime_array() { + let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; + assert_eq!(create_prime_array(19), prime_array); + assert_eq!(create_prime_array(20), prime_array); + } + + #[test] + fn test_meissel_fn_small() { + let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; + assert_eq!(meissel_fn_small(30, 8, &prime_array), 3); + assert_eq!(meissel_fn_small(100, 1, &prime_array), 50); + } + + #[test] + fn test_meissel_fn_large() { + let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; + assert_eq!(meissel_fn_large(30, 8, &prime_array), 3); + assert_eq!(meissel_fn_large(100, 1, &prime_array), 50); + } + + + #[test] + fn test_meissel_fn_memoized() { + let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; + let mut cache = HashMap::new(); + assert_eq!(meissel_memoized(30, 8, &prime_array, &mut cache), 3); + assert_eq!(meissel_memoized(100, 1, &prime_array, &mut cache), 50); + } + + #[test] + fn test_num_primes_less_than() { + assert_eq!(num_primes_less_than(7), 4); + assert_eq!(num_primes_less_than(100), 25); + assert_eq!(num_primes_less_than(2143), 324); + assert_eq!(num_primes_less_than(1_000_000), 78_498); + // assert_eq!(num_primes_less_than(1_000_000_000), 50_847_534); + } + + #[test] + fn test_num_primes_less_than_main() { + assert_eq!(num_primes_less_than_main(7), 4); + assert_eq!(num_primes_less_than_main(100), 25); + assert_eq!(num_primes_less_than_main(2143), 324); + assert_eq!(num_primes_less_than_main(1_000_000), 78_498); + assert_eq!(num_primes_less_than_main(1_000_000_000), 50_847_534); + // assert_eq!(num_primes_less_than_main(1_000_000_000_000), 37_607_912_018); + // assert_eq!(num_primes_less_than_main(1_000_000_000_000_000), 29_844_570_422_669); } } + +mod prime_count; +mod util; +pub use prime_count::num_primes_less_than_main; diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 34c52be5..d40ab2f7 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -1,7 +1,4 @@ -use crate::util::integer_square_root; -use crate::util::integer_cubic_root; -use crate::util::integer_quartic_root; -use std::io; +use crate::util::{integer_square_root, integer_cubic_root, integer_quartic_root}; use std::collections::HashMap; use std::cmp::Ordering::*; @@ -31,15 +28,6 @@ fn create_prime_array(bound: usize) -> Vec { return vec; } -fn meissel_fn_small(m: usize, n: usize, prime_array: &Vec) -> usize { - // The number of numbers <= m that are coprime to the first n prime numbers. - if n == 0 { - return m; - } else { - return meissel_fn_small(m, n-1, &prime_array) - meissel_fn_small(m / prime_array[n-1], n-1, &prime_array); - } -} - // def _meissel_function_large(self, m, n): // """The number of numbers <= m that are coprime to the first n prime numbers. // Run for larger values where repeating isn't going to happen often. Use for n > 1000 or so""" @@ -66,55 +54,25 @@ fn meissel_fn_small(m: usize, n: usize, prime_array: &Vec) -> usize { // for M in stacks[0]: // result += M*stacks[0][M] // return result -fn meissel_fn_large(m: usize, n: usize, prime_array: &Vec) -> usize { - meissel_fn_small(m, n, &prime_array) +fn meissel_fn_large(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { + meissel_fn_small(m, n, &prime_array, meissel_cache) } -fn meissel_memoized(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { +fn meissel_fn_small(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { if n == 0 { return m; } match meissel_cache.get(&(m, n)).map(|entry| entry.clone()){ Some(result) => result, None => { - let value = meissel_memoized(m, n-1, &prime_array, meissel_cache) - - meissel_memoized(m / prime_array[n-1], n-1, &prime_array, meissel_cache); + let value = meissel_fn_small(m, n-1, &prime_array, meissel_cache) + - meissel_fn_small(m / prime_array[n-1], n-1, &prime_array, meissel_cache); meissel_cache.insert((m, n), value); return value; } } } -fn num_primes_less_than(bound: usize) -> isize { - // Initialise prime array: - if bound < 1000 { - return create_prime_array(bound).len() as isize; - } - let sqrt_bound = integer_square_root(bound); - let primes = create_prime_array(sqrt_bound); - - let a = num_primes_less_than(integer_quartic_root(bound)); - let c = num_primes_less_than(integer_cubic_root(bound)); - let b = num_primes_less_than(sqrt_bound); - println!("a,b,c={},{},{}", a,b,c); - let mut result = meissel_fn_small(bound, a as usize, &primes) as isize + ((b + a - 2) * (b - a + 1)) / 2; - println!("result={}", result); - - for i in a..b { - let ith_prime = primes[i as usize]; - result -= num_primes_less_than(bound / ith_prime); - if i < c { - let bi = num_primes_less_than(integer_square_root(bound / ith_prime)); - for j in i..bi { - let jth_prime = primes[j as usize]; - result += j - num_primes_less_than(bound / ith_prime / jth_prime); - } - } - } - - return result -} - fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: &mut HashMap, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { // Initialise prime array: match prime_cache.get(&bound).map(|entry| entry.clone()){ @@ -140,7 +98,7 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: let b = num_primes_less_than_memoized(sqrt_bound, primes, prime_cache, meissel_cache); // Issues with underflow here if b + a < 2 - let mut result = meissel_memoized(bound, a, &primes, meissel_cache) + ((b + a - 2) * (b - a + 1)) / 2; + let mut result = meissel_fn_small(bound, a, &primes, meissel_cache) + ((b + a - 2) * (b - a + 1)) / 2; for i in a..b { let ith_prime = primes[i]; @@ -162,7 +120,7 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: } // Top level function -fn num_primes_less_than_main(bound: usize) -> usize { +pub fn num_primes_less_than_main(bound: usize) -> usize { let primes = create_prime_array(integer_square_root(bound)); let mut value_cache = HashMap::new(); @@ -182,102 +140,3 @@ fn num_primes_less_than_main(bound: usize) -> usize { println!("Value = {}", value_cache[&bound]); return value; } - - - -fn main() { - // println!("How many primes below: "); - // let mut input = String::new(); - - // io::stdin().read_line(&mut input) - // .expect("Failed to read line"); - // let input: usize = input.trim().parse() - // .expect("Couldn't turn value into an integer..."); - - // println!("There are {} primes below {}", create_prime_array(input).len(), input); - // println!("There are {} primes below {}", num_primes_less_than_main(input), input); - println!("There are {} primes", num_primes_less_than_main(10_000_000_000)); - - -} - -#[test] -fn test_int_sqrt() { - assert_eq!(integer_square_root(0), 0); - assert_eq!(integer_square_root(1), 1); - assert_eq!(integer_square_root(16), 4); - assert_eq!(integer_square_root(17), 4); - assert_eq!(integer_square_root(24), 4); - assert_eq!(integer_square_root(587 * 587 - 1), 586); -} - -#[test] -fn test_int_cbrt() { - assert_eq!(integer_cubic_root(0), 0); - assert_eq!(integer_cubic_root(1), 1); - assert_eq!(integer_cubic_root(26), 2); - assert_eq!(integer_cubic_root(27), 3); - assert_eq!(integer_cubic_root(28), 3); - assert_eq!(integer_cubic_root(587 * 587 * 587 - 1), 586); -} - -#[test] -fn test_int_quartic_root() { - assert_eq!(integer_quartic_root(0), 0); - assert_eq!(integer_quartic_root(1), 1); - assert_eq!(integer_quartic_root(15), 1); - assert_eq!(integer_quartic_root(16), 2); - assert_eq!(integer_quartic_root(17), 2); - assert_eq!(integer_quartic_root(587 * 587 * 587 * 587 - 1), 586); - assert_eq!(integer_quartic_root(587 * 587 * 587 * 587 + 1), 587); -} - -#[test] -fn test_prime_array() { - let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; - assert_eq!(create_prime_array(19), prime_array); - assert_eq!(create_prime_array(20), prime_array); -} - -#[test] -fn test_meissel_fn_small() { - let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; - assert_eq!(meissel_fn_small(30, 8, &prime_array), 3); - assert_eq!(meissel_fn_small(100, 1, &prime_array), 50); -} - -#[test] -fn test_meissel_fn_large() { - let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; - assert_eq!(meissel_fn_large(30, 8, &prime_array), 3); - assert_eq!(meissel_fn_large(100, 1, &prime_array), 50); -} - - -#[test] -fn test_meissel_fn_memoized() { - let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; - let mut cache = HashMap::new(); - assert_eq!(meissel_memoized(30, 8, &prime_array, &mut cache), 3); - assert_eq!(meissel_memoized(100, 1, &prime_array, &mut cache), 50); -} - -#[test] -fn test_num_primes_less_than() { - assert_eq!(num_primes_less_than(7), 4); - assert_eq!(num_primes_less_than(100), 25); - assert_eq!(num_primes_less_than(2143), 324); - assert_eq!(num_primes_less_than(1_000_000), 78_498); - // assert_eq!(num_primes_less_than(1_000_000_000), 50_847_534); -} - -#[test] -fn test_num_primes_less_than_main() { - assert_eq!(num_primes_less_than_main(7), 4); - assert_eq!(num_primes_less_than_main(100), 25); - assert_eq!(num_primes_less_than_main(2143), 324); - assert_eq!(num_primes_less_than_main(1_000_000), 78_498); - assert_eq!(num_primes_less_than_main(1_000_000_000), 50_847_534); - // assert_eq!(num_primes_less_than_main(1_000_000_000_000), 37_607_912_018); - // assert_eq!(num_primes_less_than_main(1_000_000_000_000_000), 29_844_570_422_669); -} \ No newline at end of file diff --git a/primal-count/src/util.rs b/primal-count/src/util.rs index 2cf31616..c42cc81d 100644 --- a/primal-count/src/util.rs +++ b/primal-count/src/util.rs @@ -1,4 +1,4 @@ -fn integer_square_root(value: usize) -> usize { +pub fn integer_square_root(value: usize) -> usize { // Returns the largest integer at least sqrt(n) using Newton's method let mut x = value; let mut y = (x + 1) >> 1; @@ -9,7 +9,7 @@ fn integer_square_root(value: usize) -> usize { return x; } -fn integer_cubic_root(value: usize) -> usize { +pub fn integer_cubic_root(value: usize) -> usize { // Returns the largest integer at least cbrt(n) using Newton's method let mut x = value; let mut y = (2 * value + 1) / 3; // Plus 1 only to protect against division by 0 later - not because it's right @@ -20,7 +20,7 @@ fn integer_cubic_root(value: usize) -> usize { return x; } -fn integer_quartic_root(value: usize) -> usize { +pub fn integer_quartic_root(value: usize) -> usize { // Returns the largest integer at least n^(1/4) using our fast integer sqrt return integer_square_root(integer_square_root(value)); } \ No newline at end of file From 998c7d8ae687e2a6f9b10a677df3be28cc1a39fc Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 16 Oct 2019 00:05:21 +0100 Subject: [PATCH 05/36] Clean up code + warnings --- primal-count/src/lib.rs | 45 ++++++++------------------------- primal-count/src/prime_count.rs | 9 ++++--- 2 files changed, 16 insertions(+), 38 deletions(-) diff --git a/primal-count/src/lib.rs b/primal-count/src/lib.rs index 12ae7fde..e2442fb5 100644 --- a/primal-count/src/lib.rs +++ b/primal-count/src/lib.rs @@ -2,6 +2,7 @@ mod tests { #[test] fn test_int_sqrt() { + use crate::util::integer_square_root; assert_eq!(integer_square_root(0), 0); assert_eq!(integer_square_root(1), 1); assert_eq!(integer_square_root(16), 4); @@ -12,6 +13,7 @@ mod tests { #[test] fn test_int_cbrt() { + use crate::util::integer_cubic_root; assert_eq!(integer_cubic_root(0), 0); assert_eq!(integer_cubic_root(1), 1); assert_eq!(integer_cubic_root(26), 2); @@ -22,6 +24,7 @@ mod tests { #[test] fn test_int_quartic_root() { + use crate::util::integer_quartic_root; assert_eq!(integer_quartic_root(0), 0); assert_eq!(integer_quartic_root(1), 1); assert_eq!(integer_quartic_root(15), 1); @@ -32,46 +35,18 @@ mod tests { } #[test] - fn test_prime_array() { + fn test_meissel_fn() { + use crate::prime_count::meissel_fn; + use std::collections::HashMap; let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; - assert_eq!(create_prime_array(19), prime_array); - assert_eq!(create_prime_array(20), prime_array); - } - - #[test] - fn test_meissel_fn_small() { - let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; - assert_eq!(meissel_fn_small(30, 8, &prime_array), 3); - assert_eq!(meissel_fn_small(100, 1, &prime_array), 50); - } - - #[test] - fn test_meissel_fn_large() { - let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; - assert_eq!(meissel_fn_large(30, 8, &prime_array), 3); - assert_eq!(meissel_fn_large(100, 1, &prime_array), 50); - } - - - #[test] - fn test_meissel_fn_memoized() { - let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; - let mut cache = HashMap::new(); - assert_eq!(meissel_memoized(30, 8, &prime_array, &mut cache), 3); - assert_eq!(meissel_memoized(100, 1, &prime_array, &mut cache), 50); - } - - #[test] - fn test_num_primes_less_than() { - assert_eq!(num_primes_less_than(7), 4); - assert_eq!(num_primes_less_than(100), 25); - assert_eq!(num_primes_less_than(2143), 324); - assert_eq!(num_primes_less_than(1_000_000), 78_498); - // assert_eq!(num_primes_less_than(1_000_000_000), 50_847_534); + let mut meissel_cache = HashMap::new(); + assert_eq!(meissel_fn(30, 8, &prime_array, &mut meissel_cache), 3); + assert_eq!(meissel_fn(100, 1, &prime_array, &mut meissel_cache), 50); } #[test] fn test_num_primes_less_than_main() { + use crate::prime_count::num_primes_less_than_main; assert_eq!(num_primes_less_than_main(7), 4); assert_eq!(num_primes_less_than_main(100), 25); assert_eq!(num_primes_less_than_main(2143), 324); diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index d40ab2f7..e69d0f73 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -1,6 +1,5 @@ use crate::util::{integer_square_root, integer_cubic_root, integer_quartic_root}; use std::collections::HashMap; -use std::cmp::Ordering::*; fn create_prime_array(bound: usize) -> Vec { // Returns an array of primes up to and including bound @@ -73,6 +72,10 @@ fn meissel_fn_small(m: usize, n: usize, prime_array: &Vec, meissel_cache: } } +pub fn meissel_fn(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { + meissel_fn_small(m, n, &prime_array, meissel_cache) +} + fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: &mut HashMap, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { // Initialise prime array: match prime_cache.get(&bound).map(|entry| entry.clone()){ @@ -98,7 +101,7 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: let b = num_primes_less_than_memoized(sqrt_bound, primes, prime_cache, meissel_cache); // Issues with underflow here if b + a < 2 - let mut result = meissel_fn_small(bound, a, &primes, meissel_cache) + ((b + a - 2) * (b - a + 1)) / 2; + let mut result = meissel_fn(bound, a, &primes, meissel_cache) + ((b + a - 2) * (b - a + 1)) / 2; for i in a..b { let ith_prime = primes[i]; @@ -107,7 +110,7 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: let bi = num_primes_less_than_memoized(integer_square_root(bound / ith_prime), primes, prime_cache, meissel_cache); for j in i..bi { let jth_prime = primes[j]; - result -= (num_primes_less_than_memoized(bound / ith_prime / jth_prime, primes, prime_cache, meissel_cache) - j); + result -= num_primes_less_than_memoized(bound / ith_prime / jth_prime, primes, prime_cache, meissel_cache) - j; } } } From 44fce0876f5e631974a6ab7e252ffa558a956d2d Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 16 Oct 2019 00:18:49 +0100 Subject: [PATCH 06/36] Tactical sed --- primal-sieve/src/sieve.rs | 12 ++++++------ primal-sieve/src/streaming/mod.rs | 8 ++++---- primal-sieve/src/wheel/wheel210.rs | 16 ++++++++-------- primal-sieve/src/wheel/wheel30.rs | 16 ++++++++-------- 4 files changed, 26 insertions(+), 26 deletions(-) diff --git a/primal-sieve/src/sieve.rs b/primal-sieve/src/sieve.rs index 8be885ef..39fb24f3 100644 --- a/primal-sieve/src/sieve.rs +++ b/primal-sieve/src/sieve.rs @@ -183,11 +183,11 @@ impl Sieve { pub fn prime_pi(&self, n: usize) -> usize { assert!(n <= self.upper_bound()); match n { - 0...1 => 0, + 0..=1 => 0, 2 => 1, - 3...4 => 2, - 5...6 => 3, - 7...10 => 4, + 3..=4 => 2, + 5..=6 => 3, + 7..=10 => 4, _ => { let (includes, base, tweak) = self.index_for(n); let mut count = match wheel::BYTE_MODULO { @@ -354,9 +354,9 @@ impl Sieve { pub fn primes_from<'a>(&'a self, n: usize) -> SievePrimes<'a> { assert!(n <= self.upper_bound()); let early = match n { - 0...2 => Early::Two, + 0..=2 => Early::Two, 3 => Early::Three, - 4...5 => Early::Five, + 4..=5 => Early::Five, _ => Early::Done }; let (_, base, tweak) = self.index_for(n); diff --git a/primal-sieve/src/streaming/mod.rs b/primal-sieve/src/streaming/mod.rs index 25bd0d31..ae3dd0f7 100644 --- a/primal-sieve/src/streaming/mod.rs +++ b/primal-sieve/src/streaming/mod.rs @@ -108,11 +108,11 @@ impl StreamingSieve { /// ``` pub fn prime_pi(n: usize) -> usize { match n { - 0...1 => 0, + 0..=1 => 0, 2 => 1, - 3...4 => 2, - 5...6 => 3, - 7...10 => 4, + 3..=4 => 2, + 5..=6 => 3, + 7..=10 => 4, _ => { let mut sieve = StreamingSieve::new(n); let (includes, base, tweak) = sieve.index_for(n); diff --git a/primal-sieve/src/wheel/wheel210.rs b/primal-sieve/src/wheel/wheel210.rs index a07a1f63..d68499d9 100755 --- a/primal-sieve/src/wheel/wheel210.rs +++ b/primal-sieve/src/wheel/wheel210.rs @@ -689,7 +689,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize 'outer: loop { match wi { - 0...47 => { // 30 * x + 1 + 0..=47 => { // 30 * x + 1 loop { 'label47: loop { 'label46: loop { @@ -1180,7 +1180,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 0 } } - 48...95 => { // 30 * x + 11 + 48..=95 => { // 30 * x + 11 loop { 'label95: loop { 'label94: loop { @@ -1671,7 +1671,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 48 } } - 96...143 => { // 30 * x + 13 + 96..=143 => { // 30 * x + 13 loop { 'label143: loop { 'label142: loop { @@ -2162,7 +2162,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 96 } } - 144...191 => { // 30 * x + 17 + 144..=191 => { // 30 * x + 17 loop { 'label191: loop { 'label190: loop { @@ -2653,7 +2653,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 144 } } - 192...239 => { // 30 * x + 19 + 192..=239 => { // 30 * x + 19 loop { 'label239: loop { 'label238: loop { @@ -3144,7 +3144,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 192 } } - 240...287 => { // 30 * x + 23 + 240..=287 => { // 30 * x + 23 loop { 'label287: loop { 'label286: loop { @@ -3635,7 +3635,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 240 } } - 288...335 => { // 30 * x + 29 + 288..=335 => { // 30 * x + 29 loop { 'label335: loop { 'label334: loop { @@ -4126,7 +4126,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 288 } } - 336...383 => { // 30 * x + 31 + 336..=383 => { // 30 * x + 31 loop { 'label383: loop { 'label382: loop { diff --git a/primal-sieve/src/wheel/wheel30.rs b/primal-sieve/src/wheel/wheel30.rs index bd85c76a..f544f838 100755 --- a/primal-sieve/src/wheel/wheel30.rs +++ b/primal-sieve/src/wheel/wheel30.rs @@ -195,7 +195,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize 'outer: loop { match wi { - 0...7 => { // 30 * x + 1 + 0..=7 => { // 30 * x + 1 loop { 'label7: loop { 'label6: loop { @@ -286,7 +286,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 0 } } - 8...15 => { // 30 * x + 7 + 8..=15 => { // 30 * x + 7 loop { 'label15: loop { 'label14: loop { @@ -377,7 +377,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 8 } } - 16...23 => { // 30 * x + 11 + 16..=23 => { // 30 * x + 11 loop { 'label23: loop { 'label22: loop { @@ -468,7 +468,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 16 } } - 24...31 => { // 30 * x + 13 + 24..=31 => { // 30 * x + 13 loop { 'label31: loop { 'label30: loop { @@ -559,7 +559,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 24 } } - 32...39 => { // 30 * x + 17 + 32..=39 => { // 30 * x + 17 loop { 'label39: loop { 'label38: loop { @@ -650,7 +650,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 32 } } - 40...47 => { // 30 * x + 19 + 40..=47 => { // 30 * x + 19 loop { 'label47: loop { 'label46: loop { @@ -741,7 +741,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 40 } } - 48...55 => { // 30 * x + 23 + 48..=55 => { // 30 * x + 23 loop { 'label55: loop { 'label54: loop { @@ -832,7 +832,7 @@ pub unsafe fn hardcoded_sieve(bytes: &mut [u8], si_: &mut usize, wi_: &mut usize wi = 48 } } - 56...63 => { // 30 * x + 29 + 56..=63 => { // 30 * x + 29 loop { 'label63: loop { 'label62: loop { From 7d7ff8bd095422749ae4588657da34d4a8787088 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 16 Oct 2019 00:23:27 +0100 Subject: [PATCH 07/36] Commit minor textual change to function I'm removing anyway --- primal-count/src/prime_count.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index e69d0f73..9552e8af 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -8,7 +8,7 @@ fn create_prime_array(bound: usize) -> Vec { let mut sieve = vec![true; bound + 1]; sieve[0] = false; sieve[1] = false; - for value in 2..bound+1 { + for value in 2..=bound { if sieve[value] { // Add the prime to our vector vec.push(value); From c675fe3f28f930f8710222f0e4049f688049e68c Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 16 Oct 2019 00:24:40 +0100 Subject: [PATCH 08/36] Minor textual change --- primal-count/src/lib.rs | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/primal-count/src/lib.rs b/primal-count/src/lib.rs index e2442fb5..4e6ee7b2 100644 --- a/primal-count/src/lib.rs +++ b/primal-count/src/lib.rs @@ -45,18 +45,18 @@ mod tests { } #[test] - fn test_num_primes_less_than_main() { - use crate::prime_count::num_primes_less_than_main; - assert_eq!(num_primes_less_than_main(7), 4); - assert_eq!(num_primes_less_than_main(100), 25); - assert_eq!(num_primes_less_than_main(2143), 324); - assert_eq!(num_primes_less_than_main(1_000_000), 78_498); - assert_eq!(num_primes_less_than_main(1_000_000_000), 50_847_534); - // assert_eq!(num_primes_less_than_main(1_000_000_000_000), 37_607_912_018); - // assert_eq!(num_primes_less_than_main(1_000_000_000_000_000), 29_844_570_422_669); + fn test_primes_below() { + use crate::prime_count::primes_below; + assert_eq!(primes_below(7), 4); + assert_eq!(primes_below(100), 25); + assert_eq!(primes_below(2143), 324); + assert_eq!(primes_below(1_000_000), 78_498); + assert_eq!(primes_below(1_000_000_000), 50_847_534); + // assert_eq!(primes_below(1_000_000_000_000), 37_607_912_018); + // assert_eq!(primes_below(1_000_000_000_000_000), 29_844_570_422_669); } } mod prime_count; mod util; -pub use prime_count::num_primes_less_than_main; +pub use prime_count::primes_below; From fc2792f6c02ce844616100b3d624fd0d215d6daf Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 16 Oct 2019 00:25:18 +0100 Subject: [PATCH 09/36] *sigh* --- primal-count/src/prime_count.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 9552e8af..194cab26 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -123,7 +123,7 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: } // Top level function -pub fn num_primes_less_than_main(bound: usize) -> usize { +pub fn primes_below(bound: usize) -> usize { let primes = create_prime_array(integer_square_root(bound)); let mut value_cache = HashMap::new(); From 4c01281ffa22bbe2e7265da004c42a40c340a205 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 16 Oct 2019 22:45:31 +0100 Subject: [PATCH 10/36] Improve readability --- primal-count/src/prime_count.rs | 15 ++++++++------- primal-count/src/util.rs | 1 + 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 194cab26..c8a64394 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -96,17 +96,18 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: let sqrt_bound = integer_square_root(bound); - let a = num_primes_less_than_memoized(integer_quartic_root(bound), primes, prime_cache, meissel_cache); - let c = num_primes_less_than_memoized(integer_cubic_root(bound), primes, prime_cache, meissel_cache); - let b = num_primes_less_than_memoized(sqrt_bound, primes, prime_cache, meissel_cache); + let nprimes_below_4thr = num_primes_less_than_memoized(integer_quartic_root(bound), primes, prime_cache, meissel_cache); + let nprimes_below_3rdr = num_primes_less_than_memoized(integer_cubic_root(bound), primes, prime_cache, meissel_cache); + let nprimes_below_2ndr = num_primes_less_than_memoized(sqrt_bound, primes, prime_cache, meissel_cache); - // Issues with underflow here if b + a < 2 - let mut result = meissel_fn(bound, a, &primes, meissel_cache) + ((b + a - 2) * (b - a + 1)) / 2; + // Issues with underflow here if nprimes_below_2ndr + nprimes_below_4thr < 2 + let mut result = ((nprimes_below_2ndr + nprimes_below_4thr - 2) * (nprimes_below_2ndr - nprimes_below_4thr + 1)) / 2; + result += meissel_fn(bound, nprimes_below_4thr, &primes, meissel_cache); - for i in a..b { + for i in nprimes_below_4thr..nprimes_below_2ndr { let ith_prime = primes[i]; result -= num_primes_less_than_memoized(bound / ith_prime, primes, prime_cache, meissel_cache); - if i < c { + if i < nprimes_below_3rdr { let bi = num_primes_less_than_memoized(integer_square_root(bound / ith_prime), primes, prime_cache, meissel_cache); for j in i..bi { let jth_prime = primes[j]; diff --git a/primal-count/src/util.rs b/primal-count/src/util.rs index c42cc81d..7d23ce07 100644 --- a/primal-count/src/util.rs +++ b/primal-count/src/util.rs @@ -22,5 +22,6 @@ pub fn integer_cubic_root(value: usize) -> usize { pub fn integer_quartic_root(value: usize) -> usize { // Returns the largest integer at least n^(1/4) using our fast integer sqrt + // N.b. it's faster to use two sqrts than naively apply Newton here return integer_square_root(integer_square_root(value)); } \ No newline at end of file From 3ff61d8c2a6724eae9bac1bdab246b764546d994 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 16 Oct 2019 23:46:13 +0100 Subject: [PATCH 11/36] Fix merge conflict --- primal-count/Cargo.toml | 12 ++++++++++++ primal-count/benches/bench.rs | 27 +++++++++++++++++++++++++++ primal-count/src/prime_count.rs | 1 - 3 files changed, 39 insertions(+), 1 deletion(-) create mode 100644 primal-count/benches/bench.rs diff --git a/primal-count/Cargo.toml b/primal-count/Cargo.toml index 91108834..57767e83 100644 --- a/primal-count/Cargo.toml +++ b/primal-count/Cargo.toml @@ -7,3 +7,15 @@ edition = "2018" # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html [dependencies] +<<<<<<< Updated upstream +======= +primal-sieve = { path = "../primal-sieve", version = "0.2.9" } + +[dev-dependencies] +primal = { path = "..", version = "0.2" } +bencher = "0.1.5" + +[[bench]] +name = "bench" +harness = false +>>>>>>> Stashed changes diff --git a/primal-count/benches/bench.rs b/primal-count/benches/bench.rs new file mode 100644 index 00000000..4c1d813f --- /dev/null +++ b/primal-count/benches/bench.rs @@ -0,0 +1,27 @@ +#[macro_use] +extern crate bencher; + +use bencher::Bencher; + +use primal_count::primes_below; + +fn a(bench: &mut Bencher) { + bench.iter(|| { + primes_below(7) + }) +} + +fn b(bench: &mut Bencher) { + bench.iter(|| { + primes_below(1_000) + }); +} + +fn c(bench: &mut Bencher) { + bench.iter(|| { + primes_below(10_000_000_000) + }); +} + +benchmark_group!(benches, a, b, c); +benchmark_main!(benches); \ No newline at end of file diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index c8a64394..d5bc833f 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -141,6 +141,5 @@ pub fn primes_below(bound: usize) -> usize { let mut meissel_cache = HashMap::new(); println!("Primes initialised..."); let value = num_primes_less_than_memoized(bound, &primes, &mut value_cache, &mut meissel_cache); - println!("Value = {}", value_cache[&bound]); return value; } From 992c3fc8b82fc1c96496a781109598fef19fd435 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 16 Oct 2019 23:47:16 +0100 Subject: [PATCH 12/36] Remove logging line --- primal-count/src/prime_count.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index d5bc833f..83e3bccc 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -139,7 +139,6 @@ pub fn primes_below(bound: usize) -> usize { value_cache.insert(7, 4); value_cache.insert(8, 4); let mut meissel_cache = HashMap::new(); - println!("Primes initialised..."); let value = num_primes_less_than_memoized(bound, &primes, &mut value_cache, &mut meissel_cache); return value; } From 1922499e6ca48d896632e712041544e7147eb9d8 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 16 Oct 2019 23:56:52 +0100 Subject: [PATCH 13/36] Whoops, fix merge --- primal-count/Cargo.toml | 3 --- 1 file changed, 3 deletions(-) diff --git a/primal-count/Cargo.toml b/primal-count/Cargo.toml index 57767e83..bc832350 100644 --- a/primal-count/Cargo.toml +++ b/primal-count/Cargo.toml @@ -7,8 +7,6 @@ edition = "2018" # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html [dependencies] -<<<<<<< Updated upstream -======= primal-sieve = { path = "../primal-sieve", version = "0.2.9" } [dev-dependencies] @@ -18,4 +16,3 @@ bencher = "0.1.5" [[bench]] name = "bench" harness = false ->>>>>>> Stashed changes From 1f7368a5d3a66d55bdcc2469b4959a51a9ad87d7 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Thu, 17 Oct 2019 22:56:48 +0100 Subject: [PATCH 14/36] Cleanup and some comments --- primal-count/src/prime_count.rs | 55 +++++++++------------------------ 1 file changed, 14 insertions(+), 41 deletions(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 83e3bccc..c5622f87 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -6,8 +6,6 @@ fn create_prime_array(bound: usize) -> Vec { let mut vec = Vec::new(); let max_value_to_check = integer_square_root(bound); let mut sieve = vec![true; bound + 1]; - sieve[0] = false; - sieve[1] = false; for value in 2..=bound { if sieve[value] { // Add the prime to our vector @@ -27,57 +25,26 @@ fn create_prime_array(bound: usize) -> Vec { return vec; } -// def _meissel_function_large(self, m, n): -// """The number of numbers <= m that are coprime to the first n prime numbers. -// Run for larger values where repeating isn't going to happen often. Use for n > 1000 or so""" -// m = int(m) -// # if m <= 10000: return meis104[n][m] -// if n == 0: -// return m -// result = self.meissel_memoize[n].get(m, None) -// if result is None: -// result = 0 -// primes_gen = (p for p in self.prime_array[:n:-1]) -// stacks = defaultdict(lambda: defaultdict(int)) -// stacks[n][m] = 1 -// for N in range(n, 0, -1): -// prime_dividing = next(primes_gen) -// for M in stacks[N]: -// # Marginal speed improvement? -// # if M <= 10000 and n<1229: -// # res += meis104[N][M]*stacks[N][M] -// # continue -// stacks[N-1][M] += stacks[N][M] -// stacks[N-1][M//prime_dividing] -= stacks[N][M] -// del stacks[N] -// for M in stacks[0]: -// result += M*stacks[0][M] -// return result -fn meissel_fn_large(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { - meissel_fn_small(m, n, &prime_array, meissel_cache) -} - -fn meissel_fn_small(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { +pub fn meissel_fn(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { + // The number of numbers less than m that are coprime to the first n prime numbers + // Get recurrence m_fn(m, n) = m_fn(m, n - 1) - m_fn(m/p_n, n-1) by thinking about p_n, the nth prime if n == 0 { return m; } match meissel_cache.get(&(m, n)).map(|entry| entry.clone()){ Some(result) => result, None => { - let value = meissel_fn_small(m, n-1, &prime_array, meissel_cache) - - meissel_fn_small(m / prime_array[n-1], n-1, &prime_array, meissel_cache); + let value = meissel_fn(m, n-1, &prime_array, meissel_cache) + - meissel_fn(m / prime_array[n-1], n-1, &prime_array, meissel_cache); meissel_cache.insert((m, n), value); return value; } } } -pub fn meissel_fn(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { - meissel_fn_small(m, n, &prime_array, meissel_cache) -} - fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: &mut HashMap, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { - // Initialise prime array: + // Memoized combinatorial prime counting function + // Basic idea here: https://en.wikipedia.org/wiki/Meissel%E2%80%93Lehmer_algorithm match prime_cache.get(&bound).map(|entry| entry.clone()){ Some(value) => value, None => { // The meat of the function @@ -101,6 +68,7 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: let nprimes_below_2ndr = num_primes_less_than_memoized(sqrt_bound, primes, prime_cache, meissel_cache); // Issues with underflow here if nprimes_below_2ndr + nprimes_below_4thr < 2 + // Dealt with by populating the offending (small) values in the cache at the top level let mut result = ((nprimes_below_2ndr + nprimes_below_4thr - 2) * (nprimes_below_2ndr - nprimes_below_4thr + 1)) / 2; result += meissel_fn(bound, nprimes_below_4thr, &primes, meissel_cache); @@ -116,7 +84,7 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: } } - // Caching and shit + // Caching prime_cache.insert(bound, result); return result; } @@ -125,6 +93,11 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: // Top level function pub fn primes_below(bound: usize) -> usize { + // Designed for one call + // Should refactor into a class so that multiple calls share the caches + + // N.b. we generate primes by a naive sieve, because it doesn't take much time and + // it's speed of access of primes that really matters, not generation let primes = create_prime_array(integer_square_root(bound)); let mut value_cache = HashMap::new(); From 7b518cfc0c21ce0c678cf863dc72f1860f2d3b1c Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Thu, 17 Oct 2019 23:04:49 +0100 Subject: [PATCH 15/36] More comments --- primal-count/src/prime_count.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index c5622f87..55de5284 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -45,6 +45,7 @@ pub fn meissel_fn(m: usize, n: usize, prime_array: &Vec, meissel_cache: & fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: &mut HashMap, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { // Memoized combinatorial prime counting function // Basic idea here: https://en.wikipedia.org/wiki/Meissel%E2%80%93Lehmer_algorithm + // What I've called the "Meissel Function" is phi on that Wikipedia page match prime_cache.get(&bound).map(|entry| entry.clone()){ Some(value) => value, None => { // The meat of the function From 92c7d76a57bf2f618754362f2c091fd2bc359a86 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sat, 19 Oct 2019 11:07:50 +0100 Subject: [PATCH 16/36] Add object for multiple calls --- primal-count/src/lib.rs | 16 +++++++++ primal-count/src/prime_count.rs | 60 +++++++++++++++++++++++++++------ 2 files changed, 66 insertions(+), 10 deletions(-) diff --git a/primal-count/src/lib.rs b/primal-count/src/lib.rs index 4e6ee7b2..873af4a3 100644 --- a/primal-count/src/lib.rs +++ b/primal-count/src/lib.rs @@ -55,8 +55,24 @@ mod tests { // assert_eq!(primes_below(1_000_000_000_000), 37_607_912_018); // assert_eq!(primes_below(1_000_000_000_000_000), 29_844_570_422_669); } + + #[test] + fn test_primes_below_struct() { + use crate::prime_count::PrimeCounter; + let mut pc = PrimeCounter::new(10_000); + assert_eq!(pc.primes_below(7), 4); + assert_eq!(pc.primes_below(100), 25); + assert_eq!(pc.primes_below(2143), 324); + + pc.update_limit(1_000_000_000); + assert_eq!(pc.primes_below(1_000_000), 78_498); + assert_eq!(pc.primes_below(1_000_000_000), 50_847_534); + // assert_eq!(primes_below(1_000_000_000_000), 37_607_912_018); + // assert_eq!(primes_below(1_000_000_000_000_000), 29_844_570_422_669); + } } mod prime_count; mod util; pub use prime_count::primes_below; +pub use prime_count::PrimeCounter; diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 55de5284..be91d6ef 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -102,17 +102,57 @@ pub fn primes_below(bound: usize) -> usize { let primes = create_prime_array(integer_square_root(bound)); let mut value_cache = HashMap::new(); - // Insert some basic primes - this is mainly to deal with underflow issues later - value_cache.insert(0, 0); - value_cache.insert(1, 0); - value_cache.insert(2, 1); - value_cache.insert(3, 2); - value_cache.insert(4, 2); - value_cache.insert(5, 3); - value_cache.insert(6, 3); - value_cache.insert(7, 4); - value_cache.insert(8, 4); + // Insert primes <= 10 - this is mainly to deal with underflow issues later + for n in 0..=10 { + let nprimes = match n { + 2 => 1, + 3..=4 => 2, + 5..=6 => 3, + 7..=10 => 4, + 0..=1 | _ => 0, // N.B. _ never hit + }; + value_cache.insert(n, nprimes); + } + let mut meissel_cache = HashMap::new(); let value = num_primes_less_than_memoized(bound, &primes, &mut value_cache, &mut meissel_cache); return value; } + +pub struct PrimeCounter { + limit: usize, + primes: Vec, + prime_cache: HashMap, + meissel_cache: HashMap<(usize, usize), usize> +} + +impl PrimeCounter { + pub fn new(limit: usize) -> PrimeCounter { + let primes = create_prime_array(integer_square_root(limit)); + let mut prime_cache = HashMap::new(); + + // Insert primes <= 10 - this is mainly to deal with underflow issues later + for n in 0..=10 { + let nprimes = match n { + 2 => 1, + 3..=4 => 2, + 5..=6 => 3, + 7..=10 => 4, + 0..=1 | _ => 0, // N.B. _ never hit + }; + prime_cache.insert(n, nprimes); + } + let meissel_cache = HashMap::new(); + + PrimeCounter {limit, primes, prime_cache, meissel_cache} + } + + pub fn update_limit(&mut self, limit: usize) { + self.limit = limit; + self.primes = create_prime_array(integer_square_root(limit)); + } + + pub fn primes_below(&mut self, bound: usize) -> usize { + num_primes_less_than_memoized(bound, &self.primes, &mut self.prime_cache, &mut self.meissel_cache) + } +} \ No newline at end of file From e78d1d56af4dd5b3c409d6fd6b4df2808bceb47e Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sat, 19 Oct 2019 11:08:03 +0100 Subject: [PATCH 17/36] Add some benchmarks --- primal-count/benches/bench.rs | 43 +++++++++++++++++++++++++++++++++-- 1 file changed, 41 insertions(+), 2 deletions(-) diff --git a/primal-count/benches/bench.rs b/primal-count/benches/bench.rs index 4c1d813f..c82b2772 100644 --- a/primal-count/benches/bench.rs +++ b/primal-count/benches/bench.rs @@ -4,6 +4,9 @@ extern crate bencher; use bencher::Bencher; use primal_count::primes_below; +use primal_count::PrimeCounter; + +const LARGE_VAL: usize = 10_000_000_000; fn a(bench: &mut Bencher) { bench.iter(|| { @@ -19,9 +22,45 @@ fn b(bench: &mut Bencher) { fn c(bench: &mut Bencher) { bench.iter(|| { - primes_below(10_000_000_000) + primes_below(LARGE_VAL) + }); +} + +fn d(bench: &mut Bencher) { + bench.iter(|| { + let mut pc = PrimeCounter::new(LARGE_VAL); + pc.primes_below(LARGE_VAL); + }); +} + +fn e(bench: &mut Bencher) { + bench.iter(|| { + primes_below(LARGE_VAL); + primes_below(LARGE_VAL >> 4); }); } -benchmark_group!(benches, a, b, c); +fn f(bench: &mut Bencher) { + bench.iter(|| { + let mut pc = PrimeCounter::new(LARGE_VAL); + pc.primes_below(LARGE_VAL); + pc.primes_below(LARGE_VAL >> 4); + }); +} + +fn g(bench: &mut Bencher) { + bench.iter(|| { + primes_below(LARGE_VAL * 100); + }); +} + +fn h(bench: &mut Bencher) { + bench.iter(|| { + primes_below(1_000_000_000_000_000); + }); +} + + + +benchmark_group!(benches, a, b, c, d, e, f, g, h); benchmark_main!(benches); \ No newline at end of file From 0de5956935c33138f5dd75553e19dfa3ad8dc2e3 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sat, 19 Oct 2019 12:39:01 +0100 Subject: [PATCH 18/36] Update benches and .toml --- primal-count/Cargo.toml | 12 +++- primal-count/benches/bench.rs | 121 ++++++++++++++++------------------ 2 files changed, 68 insertions(+), 65 deletions(-) diff --git a/primal-count/Cargo.toml b/primal-count/Cargo.toml index bc832350..786bb014 100644 --- a/primal-count/Cargo.toml +++ b/primal-count/Cargo.toml @@ -2,9 +2,16 @@ name = "primal-count" version = "0.1.0" authors = ["Kaimyn Chapman "] -edition = "2018" -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html +homepage = "https://github.com/huonw/primal" +repository = "https://github.com/huonw/primal" +documentation = "http://docs.rs/primal-sieve/" +license = "MIT/Apache-2.0" +keywords = ["math", "mathematics", "primes", "number-theory"] + +description = """ +A combinatorial prime counting library +""" [dependencies] primal-sieve = { path = "../primal-sieve", version = "0.2.9" } @@ -12,6 +19,7 @@ primal-sieve = { path = "../primal-sieve", version = "0.2.9" } [dev-dependencies] primal = { path = "..", version = "0.2" } bencher = "0.1.5" +criterion = "0.2" [[bench]] name = "bench" diff --git a/primal-count/benches/bench.rs b/primal-count/benches/bench.rs index c82b2772..e169ee5f 100644 --- a/primal-count/benches/bench.rs +++ b/primal-count/benches/bench.rs @@ -1,66 +1,61 @@ #[macro_use] -extern crate bencher; - -use bencher::Bencher; - +extern crate criterion; +extern crate primal_count; use primal_count::primes_below; use primal_count::PrimeCounter; - -const LARGE_VAL: usize = 10_000_000_000; - -fn a(bench: &mut Bencher) { - bench.iter(|| { - primes_below(7) - }) -} - -fn b(bench: &mut Bencher) { - bench.iter(|| { - primes_below(1_000) - }); -} - -fn c(bench: &mut Bencher) { - bench.iter(|| { - primes_below(LARGE_VAL) - }); -} - -fn d(bench: &mut Bencher) { - bench.iter(|| { - let mut pc = PrimeCounter::new(LARGE_VAL); - pc.primes_below(LARGE_VAL); - }); -} - -fn e(bench: &mut Bencher) { - bench.iter(|| { - primes_below(LARGE_VAL); - primes_below(LARGE_VAL >> 4); - }); -} - -fn f(bench: &mut Bencher) { - bench.iter(|| { - let mut pc = PrimeCounter::new(LARGE_VAL); - pc.primes_below(LARGE_VAL); - pc.primes_below(LARGE_VAL >> 4); - }); -} - -fn g(bench: &mut Bencher) { - bench.iter(|| { - primes_below(LARGE_VAL * 100); - }); -} - -fn h(bench: &mut Bencher) { - bench.iter(|| { - primes_below(1_000_000_000_000_000); - }); -} - - - -benchmark_group!(benches, a, b, c, d, e, f, g, h); -benchmark_main!(benches); \ No newline at end of file +use criterion::{Criterion, ParameterizedBenchmark}; + +const SIZES: [usize; 6] = [100, 10_000, 100_000, 1_000_000, 10_000_000, 10_000_000_000]; + +macro_rules! create_benchmarks { + ($( + fn $group_id: ident($input: expr) { + $first_name: expr => $first_func: expr, + $($rest_name: expr => $rest_func: expr,)* + } + )*) => { + $( + fn $group_id(c: &mut Criterion) { + let input = $input; + + let bench = ParameterizedBenchmark::new( + $first_name, $first_func, input.into_iter().cloned()) + $( .with_function($rest_name, $rest_func) )*; + c.bench(stringify!($group_id), bench); + } + )* + } +} + +create_benchmarks! { + fn new(SIZES) { + "PrimeCounter" => |b, upto: &usize| b.iter(|| PrimeCounter::new(*upto)), + } + + fn prime_pi(SIZES) { + "PrimeCounter" => |b, upto: &usize| { + let mut s = PrimeCounter::new(*upto + 1); + b.iter(|| s.primes_below(*upto)); + }, + + "PrimeCounter with init" => |b, upto: &usize| { + b.iter(|| { + let mut s = PrimeCounter::new(*upto + 1); + s.primes_below(*upto) + }); + }, + + "PrimesBelow" => |b, upto: &usize| { + b.iter(|| primes_below(*upto)); + }, + // "Sieve with init" => |b, upto: &usize| { + // b.iter(|| { + // let s = Sieve::new(*upto + 1); + // s.prime_pi(*upto) + // }); + // }, + } +} + +criterion_group!(benches, new, prime_pi); +criterion_main!(benches); From 9a9d44d36e570e1b41a0578f2d2029bfef557de9 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sat, 19 Oct 2019 17:08:00 +0100 Subject: [PATCH 19/36] Cleanup and start using primal library --- primal-count/src/prime_count.rs | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index be91d6ef..d7f659a2 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -1,3 +1,6 @@ +extern crate primal_sieve; +use std::iter::FromIterator; + use crate::util::{integer_square_root, integer_cubic_root, integer_quartic_root}; use std::collections::HashMap; @@ -128,7 +131,13 @@ pub struct PrimeCounter { impl PrimeCounter { pub fn new(limit: usize) -> PrimeCounter { - let primes = create_prime_array(integer_square_root(limit)); + // let primes = create_prime_array(integer_square_root(limit)); + let int_sqrt = integer_square_root(limit); + let sieve = primal_sieve::Sieve::new(int_sqrt); + let sieve_iter = sieve.primes_from(2).take_while(|x| *x <= int_sqrt); + let primes = Vec::from_iter(sieve_iter); + // let primes = primal_sieve::Primes::all().take_while(|x| *x <= integer_square_root(limit)); + let mut prime_cache = HashMap::new(); // Insert primes <= 10 - this is mainly to deal with underflow issues later @@ -153,6 +162,6 @@ impl PrimeCounter { } pub fn primes_below(&mut self, bound: usize) -> usize { - num_primes_less_than_memoized(bound, &self.primes, &mut self.prime_cache, &mut self.meissel_cache) + num_primes_less_than_memoized2(bound, &self.primes, &mut self.prime_cache, &mut self.meissel_cache) } } \ No newline at end of file From f76d2445106ee00e002f548f5e4012ae167761b4 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sat, 19 Oct 2019 18:30:55 +0100 Subject: [PATCH 20/36] Refactor --- primal-count/benches/bench.rs | 10 --- primal-count/src/lib.rs | 75 ---------------- primal-count/src/prime_count.rs | 151 ++++++++++++++++---------------- primal-count/src/util.rs | 60 ++++++++++--- 4 files changed, 124 insertions(+), 172 deletions(-) diff --git a/primal-count/benches/bench.rs b/primal-count/benches/bench.rs index e169ee5f..12171f58 100644 --- a/primal-count/benches/bench.rs +++ b/primal-count/benches/bench.rs @@ -44,16 +44,6 @@ create_benchmarks! { s.primes_below(*upto) }); }, - - "PrimesBelow" => |b, upto: &usize| { - b.iter(|| primes_below(*upto)); - }, - // "Sieve with init" => |b, upto: &usize| { - // b.iter(|| { - // let s = Sieve::new(*upto + 1); - // s.prime_pi(*upto) - // }); - // }, } } diff --git a/primal-count/src/lib.rs b/primal-count/src/lib.rs index 873af4a3..147db5ac 100644 --- a/primal-count/src/lib.rs +++ b/primal-count/src/lib.rs @@ -1,78 +1,3 @@ -#[cfg(test)] -mod tests { - #[test] - fn test_int_sqrt() { - use crate::util::integer_square_root; - assert_eq!(integer_square_root(0), 0); - assert_eq!(integer_square_root(1), 1); - assert_eq!(integer_square_root(16), 4); - assert_eq!(integer_square_root(17), 4); - assert_eq!(integer_square_root(24), 4); - assert_eq!(integer_square_root(587 * 587 - 1), 586); - } - - #[test] - fn test_int_cbrt() { - use crate::util::integer_cubic_root; - assert_eq!(integer_cubic_root(0), 0); - assert_eq!(integer_cubic_root(1), 1); - assert_eq!(integer_cubic_root(26), 2); - assert_eq!(integer_cubic_root(27), 3); - assert_eq!(integer_cubic_root(28), 3); - assert_eq!(integer_cubic_root(587 * 587 * 587 - 1), 586); - } - - #[test] - fn test_int_quartic_root() { - use crate::util::integer_quartic_root; - assert_eq!(integer_quartic_root(0), 0); - assert_eq!(integer_quartic_root(1), 1); - assert_eq!(integer_quartic_root(15), 1); - assert_eq!(integer_quartic_root(16), 2); - assert_eq!(integer_quartic_root(17), 2); - assert_eq!(integer_quartic_root(587 * 587 * 587 * 587 - 1), 586); - assert_eq!(integer_quartic_root(587 * 587 * 587 * 587 + 1), 587); - } - - #[test] - fn test_meissel_fn() { - use crate::prime_count::meissel_fn; - use std::collections::HashMap; - let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; - let mut meissel_cache = HashMap::new(); - assert_eq!(meissel_fn(30, 8, &prime_array, &mut meissel_cache), 3); - assert_eq!(meissel_fn(100, 1, &prime_array, &mut meissel_cache), 50); - } - - #[test] - fn test_primes_below() { - use crate::prime_count::primes_below; - assert_eq!(primes_below(7), 4); - assert_eq!(primes_below(100), 25); - assert_eq!(primes_below(2143), 324); - assert_eq!(primes_below(1_000_000), 78_498); - assert_eq!(primes_below(1_000_000_000), 50_847_534); - // assert_eq!(primes_below(1_000_000_000_000), 37_607_912_018); - // assert_eq!(primes_below(1_000_000_000_000_000), 29_844_570_422_669); - } - - #[test] - fn test_primes_below_struct() { - use crate::prime_count::PrimeCounter; - let mut pc = PrimeCounter::new(10_000); - assert_eq!(pc.primes_below(7), 4); - assert_eq!(pc.primes_below(100), 25); - assert_eq!(pc.primes_below(2143), 324); - - pc.update_limit(1_000_000_000); - assert_eq!(pc.primes_below(1_000_000), 78_498); - assert_eq!(pc.primes_below(1_000_000_000), 50_847_534); - // assert_eq!(primes_below(1_000_000_000_000), 37_607_912_018); - // assert_eq!(primes_below(1_000_000_000_000_000), 29_844_570_422_669); - } -} - mod prime_count; mod util; -pub use prime_count::primes_below; pub use prime_count::PrimeCounter; diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index d7f659a2..8de857a9 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -1,36 +1,25 @@ extern crate primal_sieve; use std::iter::FromIterator; -use crate::util::{integer_square_root, integer_cubic_root, integer_quartic_root}; +use crate::util::{int_square_root, int_cubic_root, int_quartic_root}; use std::collections::HashMap; -fn create_prime_array(bound: usize) -> Vec { - // Returns an array of primes up to and including bound - let mut vec = Vec::new(); - let max_value_to_check = integer_square_root(bound); - let mut sieve = vec![true; bound + 1]; - for value in 2..=bound { - if sieve[value] { - // Add the prime to our vector - vec.push(value); - - // Filter primes more if we're below the bound - if value <= max_value_to_check { - // Anything smaller than this already filtered - let mut idx = value * value; - while idx <= bound { - sieve[idx] = false; - idx += value; - } - } - } - } - return vec; +/// Generate a vec of primes from 2 up to and including limit +/// Leverages the fast sieve in primal to do so +fn generate_primes(limit: usize) -> Vec { + let sieve = primal_sieve::Sieve::new(limit); + let sieve_iter = sieve.primes_from(2).take_while(|x| *x <= limit); + // Note that we put the primes into a vec here because later we want to have both + // 1) Very quick access to the nth prime + // 2) Quick counting of number of primes below a value, achieved with a binary search + // Experiments replacing 1) or 2) with the methods in sieve seem to significantly + // slow things down for larger numbers + return Vec::from_iter(sieve_iter); } -pub fn meissel_fn(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { - // The number of numbers less than m that are coprime to the first n prime numbers - // Get recurrence m_fn(m, n) = m_fn(m, n - 1) - m_fn(m/p_n, n-1) by thinking about p_n, the nth prime +/// The number of numbers less than m that are coprime to the first n prime numbers +/// Get recurrence m_fn(m, n) = m_fn(m, n - 1) - m_fn(m/p_n, n-1) by thinking about p_n, the nth prime +fn meissel_fn(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { if n == 0 { return m; } @@ -45,13 +34,14 @@ pub fn meissel_fn(m: usize, n: usize, prime_array: &Vec, meissel_cache: & } } -fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: &mut HashMap, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { - // Memoized combinatorial prime counting function - // Basic idea here: https://en.wikipedia.org/wiki/Meissel%E2%80%93Lehmer_algorithm - // What I've called the "Meissel Function" is phi on that Wikipedia page +/// Find the number of primes less than bound using the Meissel-Lehmer method +/// Leverages caching to speed up the recursive calls +fn primes_less_than(bound: usize, primes: &Vec, prime_cache: &mut HashMap, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { + // First check if it's in the cache already match prime_cache.get(&bound).map(|entry| entry.clone()){ Some(value) => value, - None => { // The meat of the function + None => { + // The meat of the function if bound < 2 { return 0; } else if bound <= primes[primes.len()-1] { @@ -65,11 +55,11 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: return result; } - let sqrt_bound = integer_square_root(bound); + let sqrt_bound = int_square_root(bound); - let nprimes_below_4thr = num_primes_less_than_memoized(integer_quartic_root(bound), primes, prime_cache, meissel_cache); - let nprimes_below_3rdr = num_primes_less_than_memoized(integer_cubic_root(bound), primes, prime_cache, meissel_cache); - let nprimes_below_2ndr = num_primes_less_than_memoized(sqrt_bound, primes, prime_cache, meissel_cache); + let nprimes_below_4thr = primes_less_than(int_quartic_root(bound), primes, prime_cache, meissel_cache); + let nprimes_below_3rdr = primes_less_than(int_cubic_root(bound), primes, prime_cache, meissel_cache); + let nprimes_below_2ndr = primes_less_than(sqrt_bound, primes, prime_cache, meissel_cache); // Issues with underflow here if nprimes_below_2ndr + nprimes_below_4thr < 2 // Dealt with by populating the offending (small) values in the cache at the top level @@ -78,12 +68,12 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: for i in nprimes_below_4thr..nprimes_below_2ndr { let ith_prime = primes[i]; - result -= num_primes_less_than_memoized(bound / ith_prime, primes, prime_cache, meissel_cache); + result -= primes_less_than(bound / ith_prime, primes, prime_cache, meissel_cache); if i < nprimes_below_3rdr { - let bi = num_primes_less_than_memoized(integer_square_root(bound / ith_prime), primes, prime_cache, meissel_cache); + let bi = primes_less_than(int_square_root(bound / ith_prime), primes, prime_cache, meissel_cache); for j in i..bi { let jth_prime = primes[j]; - result -= num_primes_less_than_memoized(bound / ith_prime / jth_prime, primes, prime_cache, meissel_cache) - j; + result -= primes_less_than(bound / ith_prime / jth_prime, primes, prime_cache, meissel_cache) - j; } } } @@ -95,33 +85,12 @@ fn num_primes_less_than_memoized(bound: usize, primes: &Vec, prime_cache: } } -// Top level function -pub fn primes_below(bound: usize) -> usize { - // Designed for one call - // Should refactor into a class so that multiple calls share the caches - - // N.b. we generate primes by a naive sieve, because it doesn't take much time and - // it's speed of access of primes that really matters, not generation - let primes = create_prime_array(integer_square_root(bound)); - let mut value_cache = HashMap::new(); - - // Insert primes <= 10 - this is mainly to deal with underflow issues later - for n in 0..=10 { - let nprimes = match n { - 2 => 1, - 3..=4 => 2, - 5..=6 => 3, - 7..=10 => 4, - 0..=1 | _ => 0, // N.B. _ never hit - }; - value_cache.insert(n, nprimes); - } - - let mut meissel_cache = HashMap::new(); - let value = num_primes_less_than_memoized(bound, &primes, &mut value_cache, &mut meissel_cache); - return value; -} - +/// Memoized combinatorial prime counting function +/// Basic idea here: https://en.wikipedia.org/wiki/Meissel%E2%80%93Lehmer_algorithm +/// The "Meissel Function" here is phi on that Wikipedia page +/// +/// # Examples +/// pub struct PrimeCounter { limit: usize, primes: Vec, @@ -130,15 +99,10 @@ pub struct PrimeCounter { } impl PrimeCounter { + /// Create a new PrimeCounter instance, which generates all the primes up to sqrt(limit) pub fn new(limit: usize) -> PrimeCounter { - // let primes = create_prime_array(integer_square_root(limit)); - let int_sqrt = integer_square_root(limit); - let sieve = primal_sieve::Sieve::new(int_sqrt); - let sieve_iter = sieve.primes_from(2).take_while(|x| *x <= int_sqrt); - let primes = Vec::from_iter(sieve_iter); - // let primes = primal_sieve::Primes::all().take_while(|x| *x <= integer_square_root(limit)); - let mut prime_cache = HashMap::new(); + let primes = generate_primes(limit); // Insert primes <= 10 - this is mainly to deal with underflow issues later for n in 0..=10 { @@ -152,16 +116,51 @@ impl PrimeCounter { prime_cache.insert(n, nprimes); } let meissel_cache = HashMap::new(); - + PrimeCounter {limit, primes, prime_cache, meissel_cache} } + /// Bla pub fn update_limit(&mut self, limit: usize) { self.limit = limit; - self.primes = create_prime_array(integer_square_root(limit)); + self.primes = generate_primes(int_square_root(limit)); } - + + /// Bla pub fn primes_below(&mut self, bound: usize) -> usize { - num_primes_less_than_memoized2(bound, &self.primes, &mut self.prime_cache, &mut self.meissel_cache) + primes_less_than(bound, &self.primes, &mut self.prime_cache, &mut self.meissel_cache) + } + + /// Bla + pub fn meissel_fn(&mut self, m: usize, n: usize) -> usize { + meissel_fn(m, n, &self.primes, &mut self.meissel_cache) + } +} + +#[cfg(test)] +mod tests { + #[test] + fn test_meissel_fn() { + use prime_count::meissel_fn; + use std::collections::HashMap; + let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; + let mut meissel_cache = HashMap::new(); + assert_eq!(meissel_fn(30, 8, &prime_array, &mut meissel_cache), 3); + assert_eq!(meissel_fn(100, 1, &prime_array, &mut meissel_cache), 50); + } + + #[test] + fn test_primes_below() { + use crate::prime_count::PrimeCounter; + let mut pc = PrimeCounter::new(10_000); + assert_eq!(pc.primes_below(7), 4); + assert_eq!(pc.primes_below(100), 25); + assert_eq!(pc.primes_below(2143), 324); + + pc.update_limit(1_000_000_000); + assert_eq!(pc.primes_below(1_000_000), 78_498); + assert_eq!(pc.primes_below(1_000_000_000), 50_847_534); + // assert_eq!(primes_below(1_000_000_000_000), 37_607_912_018); + // assert_eq!(primes_below(1_000_000_000_000_000), 29_844_570_422_669); } } \ No newline at end of file diff --git a/primal-count/src/util.rs b/primal-count/src/util.rs index 7d23ce07..a2682fa8 100644 --- a/primal-count/src/util.rs +++ b/primal-count/src/util.rs @@ -1,5 +1,5 @@ -pub fn integer_square_root(value: usize) -> usize { - // Returns the largest integer at least sqrt(n) using Newton's method +/// Returns the largest integer at least sqrt(n) using Newton's method +pub fn int_square_root(value: usize) -> usize { let mut x = value; let mut y = (x + 1) >> 1; while y < x { @@ -9,19 +9,57 @@ pub fn integer_square_root(value: usize) -> usize { return x; } -pub fn integer_cubic_root(value: usize) -> usize { - // Returns the largest integer at least cbrt(n) using Newton's method +/// Returns the largest integer at least cbrt(n) using Newton's method +pub fn int_cubic_root(value: usize) -> usize { + let mut x = value; - let mut y = (2 * value + 1) / 3; // Plus 1 only to protect against division by 0 later - not because it's right + let mut y = (2 * value + 1) / 3; // Plus 1 only to protect against division by 0 later while y < x { x = y; - y = (2 * x + value / x / x) / 3; // Weird divide to protect against overflow + y = (2 * x + value / x / x) / 3; // Divide split to protect against overflow } return x; } -pub fn integer_quartic_root(value: usize) -> usize { - // Returns the largest integer at least n^(1/4) using our fast integer sqrt - // N.b. it's faster to use two sqrts than naively apply Newton here - return integer_square_root(integer_square_root(value)); -} \ No newline at end of file +// Returns the largest integer at least n^(1/4) using our fast integer sqrt +// N.b. it's faster to use two sqrts than naively apply Newton here +pub fn int_quartic_root(value: usize) -> usize { + return int_square_root(int_square_root(value)); +} + +#[cfg(test)] +mod tests { + #[test] + fn test_int_sqrt() { + use util::int_square_root; + assert_eq!(int_square_root(0), 0); + assert_eq!(int_square_root(1), 1); + assert_eq!(int_square_root(16), 4); + assert_eq!(int_square_root(17), 4); + assert_eq!(int_square_root(24), 4); + assert_eq!(int_square_root(587 * 587 - 1), 586); + } + + #[test] + fn test_int_cbrt() { + use util::int_cubic_root; + assert_eq!(int_cubic_root(0), 0); + assert_eq!(int_cubic_root(1), 1); + assert_eq!(int_cubic_root(26), 2); + assert_eq!(int_cubic_root(27), 3); + assert_eq!(int_cubic_root(28), 3); + assert_eq!(int_cubic_root(587 * 587 * 587 - 1), 586); + } + + #[test] + fn test_int_quartic_root() { + use util::int_quartic_root; + assert_eq!(int_quartic_root(0), 0); + assert_eq!(int_quartic_root(1), 1); + assert_eq!(int_quartic_root(15), 1); + assert_eq!(int_quartic_root(16), 2); + assert_eq!(int_quartic_root(17), 2); + assert_eq!(int_quartic_root(587 * 587 * 587 * 587 - 1), 586); + assert_eq!(int_quartic_root(587 * 587 * 587 * 587 + 1), 587); + } +} From d928a07897a13296f3496188e2cb5b3d036acf81 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sat, 19 Oct 2019 18:31:37 +0100 Subject: [PATCH 21/36] Update Cargo.toml --- primal-count/Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/primal-count/Cargo.toml b/primal-count/Cargo.toml index 786bb014..0e53f6b5 100644 --- a/primal-count/Cargo.toml +++ b/primal-count/Cargo.toml @@ -1,7 +1,7 @@ [package] name = "primal-count" version = "0.1.0" -authors = ["Kaimyn Chapman "] +authors = ["Kaimyn Chapman "] homepage = "https://github.com/huonw/primal" repository = "https://github.com/huonw/primal" From 11bf5c8547cb5e6a88d37e48533c3c905b3fa6e7 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sat, 19 Oct 2019 18:37:48 +0100 Subject: [PATCH 22/36] Clean up tests --- primal-count/src/prime_count.rs | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 8de857a9..d7abe88d 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -141,12 +141,11 @@ impl PrimeCounter { mod tests { #[test] fn test_meissel_fn() { - use prime_count::meissel_fn; - use std::collections::HashMap; - let prime_array = vec![2, 3, 5, 7, 11, 13, 17, 19]; - let mut meissel_cache = HashMap::new(); - assert_eq!(meissel_fn(30, 8, &prime_array, &mut meissel_cache), 3); - assert_eq!(meissel_fn(100, 1, &prime_array, &mut meissel_cache), 50); + use crate::prime_count::PrimeCounter; + let mut pc = PrimeCounter::new(100); + assert_eq!(pc.meissel_fn(30, 8), 3); + assert_eq!(pc.meissel_fn(100, 1), 50); + assert_eq!(pc.meissel_fn(100, 25), 1); } #[test] From 54c2ea647a5f3164394213dc283d1383be3d2614 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sat, 19 Oct 2019 18:45:36 +0100 Subject: [PATCH 23/36] Clean up tests --- primal-count/src/prime_count.rs | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index d7abe88d..c11008ad 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -120,7 +120,10 @@ impl PrimeCounter { PrimeCounter {limit, primes, prime_cache, meissel_cache} } - /// Bla + /// Updates the limit - to be used when you want to make the prime cache larger + /// May be slow for large values of limit - it's recommended that you don't call + /// this and instead ensure that your first call to construct the PrimeCounter + /// object is large enough. pub fn update_limit(&mut self, limit: usize) { self.limit = limit; self.primes = generate_primes(int_square_root(limit)); @@ -131,7 +134,16 @@ impl PrimeCounter { primes_less_than(bound, &self.primes, &mut self.prime_cache, &mut self.meissel_cache) } - /// Bla + /// The number of numbers less than m that are coprime to the first n prime numbers + /// + /// # Examples + /// + /// ```rust + /// # extern crate primal; + /// let pc = primal::Sieve::PrimeCounter(10_000); + /// assert_eq!(pc.meissel_fn(100, 10), 30) + /// assert_eq!(pc.meissel_fn(1234, 5), 300) + /// ``` pub fn meissel_fn(&mut self, m: usize, n: usize) -> usize { meissel_fn(m, n, &self.primes, &mut self.meissel_cache) } From 605a225c205198cd21b88b1ebcff41c03c16e019 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sat, 19 Oct 2019 19:31:46 +0100 Subject: [PATCH 24/36] Fix bug and update documentation --- primal-count/benches/bench.rs | 21 +++++++++++++++++--- primal-count/src/lib.rs | 1 + primal-count/src/prime_count.rs | 35 +++++++++++++++++++++++---------- 3 files changed, 44 insertions(+), 13 deletions(-) diff --git a/primal-count/benches/bench.rs b/primal-count/benches/bench.rs index 12171f58..192ab6e7 100644 --- a/primal-count/benches/bench.rs +++ b/primal-count/benches/bench.rs @@ -1,8 +1,7 @@ #[macro_use] extern crate criterion; extern crate primal_count; -use primal_count::primes_below; -use primal_count::PrimeCounter; +use primal_count::{PrimeCounter, int_quartic_root}; use criterion::{Criterion, ParameterizedBenchmark}; const SIZES: [usize; 6] = [100, 10_000, 100_000, 1_000_000, 10_000_000, 10_000_000_000]; @@ -45,7 +44,23 @@ create_benchmarks! { }); }, } + + fn meissel_fn(SIZES) { + "PrimeCounter with init, n=10" => |b, upto: &usize| { + b.iter(|| { + let mut s = PrimeCounter::new(*upto + 1); + s.meissel_fn(*upto, 10) + }); + }, + + "PrimeCounter with init, n=4th root" => |b, upto: &usize| { + b.iter(|| { + let mut s = PrimeCounter::new(*upto + 1); + s.meissel_fn(*upto, int_quartic_root(*upto)) + }); + }, + } } -criterion_group!(benches, new, prime_pi); +criterion_group!(benches, new, prime_pi, meissel_fn); criterion_main!(benches); diff --git a/primal-count/src/lib.rs b/primal-count/src/lib.rs index 147db5ac..2bb864ad 100644 --- a/primal-count/src/lib.rs +++ b/primal-count/src/lib.rs @@ -1,3 +1,4 @@ mod prime_count; mod util; pub use prime_count::PrimeCounter; +pub use util::{int_square_root, int_cubic_root, int_quartic_root}; \ No newline at end of file diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index c11008ad..f212ed41 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -88,9 +88,6 @@ fn primes_less_than(bound: usize, primes: &Vec, prime_cache: &mut HashMap /// Memoized combinatorial prime counting function /// Basic idea here: https://en.wikipedia.org/wiki/Meissel%E2%80%93Lehmer_algorithm /// The "Meissel Function" here is phi on that Wikipedia page -/// -/// # Examples -/// pub struct PrimeCounter { limit: usize, primes: Vec, @@ -102,7 +99,7 @@ impl PrimeCounter { /// Create a new PrimeCounter instance, which generates all the primes up to sqrt(limit) pub fn new(limit: usize) -> PrimeCounter { let mut prime_cache = HashMap::new(); - let primes = generate_primes(limit); + let primes = generate_primes(int_square_root(limit)); // Insert primes <= 10 - this is mainly to deal with underflow issues later for n in 0..=10 { @@ -129,20 +126,38 @@ impl PrimeCounter { self.primes = generate_primes(int_square_root(limit)); } - /// Bla + /// The number of primes that are at least `bound` + /// + /// # Panics + /// + /// If the limit in the constructor is smaller than the input of primes_below + /// (Unless it's later been updated with update_limit) + /// + /// # Examples + /// + /// ```rust + /// # extern crate primal; + /// let mut pc = primal::PrimeCounter::new(10_000); + /// assert_eq!(pc.primes_below(100), 25); + /// assert_eq!(pc.primes_below(8166), 1024); + /// ``` pub fn primes_below(&mut self, bound: usize) -> usize { primes_less_than(bound, &self.primes, &mut self.prime_cache, &mut self.meissel_cache) } - /// The number of numbers less than m that are coprime to the first n prime numbers + /// The number of numbers less than `m` that are coprime to the first `n` prime numbers + /// + /// # Panics + /// + /// If the `n`th prime is larger than `limit` /// /// # Examples /// /// ```rust /// # extern crate primal; - /// let pc = primal::Sieve::PrimeCounter(10_000); - /// assert_eq!(pc.meissel_fn(100, 10), 30) - /// assert_eq!(pc.meissel_fn(1234, 5), 300) + /// let mut pc = primal::PrimeCounter::new(10_000); + /// assert_eq!(pc.meissel_fn(100, 10), 16); + /// assert_eq!(pc.meissel_fn(1234, 5), 256); /// ``` pub fn meissel_fn(&mut self, m: usize, n: usize) -> usize { meissel_fn(m, n, &self.primes, &mut self.meissel_cache) @@ -154,7 +169,7 @@ mod tests { #[test] fn test_meissel_fn() { use crate::prime_count::PrimeCounter; - let mut pc = PrimeCounter::new(100); + let mut pc = PrimeCounter::new(10_000); assert_eq!(pc.meissel_fn(30, 8), 3); assert_eq!(pc.meissel_fn(100, 1), 50); assert_eq!(pc.meissel_fn(100, 25), 1); From 7ec5a825f3d9b8ba5d8c51069722e8fd9fe2ade9 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sun, 20 Oct 2019 14:56:10 +0100 Subject: [PATCH 25/36] Update global Cargo.toml and lib.rs --- Cargo.toml | 1 + src/lib.rs | 2 ++ 2 files changed, 3 insertions(+) diff --git a/Cargo.toml b/Cargo.toml index 686d3885..31a791e2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -22,6 +22,7 @@ prime). primal-check = { path = "primal-check", version = "0.2" } primal-estimate = { path = "primal-estimate", version = "0.2" } primal-sieve = { path = "primal-sieve", version = "0.2" } +primal-count = { path = "primal-count", version = "0.1" } [dev-dependencies] primal-slowsieve = { path = "primal-slowsieve", version = "0.2" } diff --git a/src/lib.rs b/src/lib.rs index c902837b..5e654f14 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -254,6 +254,7 @@ extern crate primal_estimate; extern crate primal_check; extern crate primal_sieve; +extern crate primal_count; pub use primal_estimate::prime_pi as estimate_prime_pi; pub use primal_estimate::nth_prime as estimate_nth_prime; @@ -261,3 +262,4 @@ pub use primal_check::miller_rabin as is_prime; pub use primal_check::{as_perfect_power, as_prime_power}; pub use primal_sieve::{StreamingSieve, Sieve, SievePrimes, Primes}; +pub use primal_count::PrimeCounter; From a1a678d883a3520e1527647722b0a596a7464d47 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sun, 20 Oct 2019 14:56:28 +0100 Subject: [PATCH 26/36] Update benchmarks --- benches/bench.rs | 59 +++++++++++++++++++++++++++++++++-- primal-count/benches/bench.rs | 18 +---------- 2 files changed, 57 insertions(+), 20 deletions(-) diff --git a/benches/bench.rs b/benches/bench.rs index 23acda6e..ba4b02a9 100644 --- a/benches/bench.rs +++ b/benches/bench.rs @@ -1,10 +1,31 @@ #[macro_use] extern crate criterion; extern crate primal; -use criterion::{Criterion, Fun}; +use criterion::{Criterion, Fun, ParameterizedBenchmark}; const N: usize = 1_000_000; const STEP: usize = 101; +const SIZES: [usize; 8] = [100, 10_000, 100_000, 1_000_000, 10_000_000, 100_000_000, 1_000_000_000, 10_000_000_000]; + +macro_rules! create_benchmarks { + ($( + fn $group_id: ident($input: expr) { + $first_name: expr => $first_func: expr, + $($rest_name: expr => $rest_func: expr,)* + } + )*) => { + $( + fn $group_id(c: &mut Criterion) { + let input = $input; + + let bench = ParameterizedBenchmark::new( + $first_name, $first_func, input.into_iter().cloned()) + $( .with_function($rest_name, $rest_func) )*; + c.bench(stringify!($group_id), bench); + } + )* + } +} fn is_prime(c: &mut Criterion) { let miller_rabin = Fun::new( @@ -35,6 +56,38 @@ fn is_prime(c: &mut Criterion) { "is_prime", vec![miller_rabin, sieve, sieve_with_init], ()); } -criterion_group!(benches, is_prime); -criterion_main!(benches); +create_benchmarks!{ + fn prime_pi(SIZES) { + "PrimeCounter" => |b, upto: &usize| { + let mut s = primal::PrimeCounter::new(*upto + 1); + b.iter(|| s.primes_below(*upto)); + }, + "Sieve" => |b, upto: &usize| { + let s = primal::Sieve::new(*upto + 1); + b.iter(|| s.prime_pi(*upto)); + }, + + "PrimeCounter with init" => |b, upto: &usize| { + b.iter(|| { + let mut s = primal::PrimeCounter::new(*upto + 1); + s.primes_below(*upto) + }); + }, + "Sieve with init" => |b, upto: &usize| { + b.iter(|| { + let s = primal::Sieve::new(*upto + 1); + s.prime_pi(*upto) + }); + }, + "StreamingSieve" => |b, upto: &usize| { + b.iter(|| primal::StreamingSieve::prime_pi(*upto)) + }, + "Primes" => |b, upto: &usize| { + b.iter(|| primal::Primes::all().take_while(|x| *x <= *upto).count()) + }, + } +} + +criterion_group!(benches, is_prime, prime_pi); +criterion_main!(benches); diff --git a/primal-count/benches/bench.rs b/primal-count/benches/bench.rs index 192ab6e7..4251e934 100644 --- a/primal-count/benches/bench.rs +++ b/primal-count/benches/bench.rs @@ -44,23 +44,7 @@ create_benchmarks! { }); }, } - - fn meissel_fn(SIZES) { - "PrimeCounter with init, n=10" => |b, upto: &usize| { - b.iter(|| { - let mut s = PrimeCounter::new(*upto + 1); - s.meissel_fn(*upto, 10) - }); - }, - - "PrimeCounter with init, n=4th root" => |b, upto: &usize| { - b.iter(|| { - let mut s = PrimeCounter::new(*upto + 1); - s.meissel_fn(*upto, int_quartic_root(*upto)) - }); - }, - } } -criterion_group!(benches, new, prime_pi, meissel_fn); +criterion_group!(benches, new, prime_pi); criterion_main!(benches); From 7febca939a0bc9e98979c001e4f1b45520c24b99 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sun, 20 Oct 2019 21:23:29 +0100 Subject: [PATCH 27/36] Significantly speed up computation of meissel for large values --- primal-count/benches/bench.rs | 2 +- primal-count/src/prime_count.rs | 19 ++++++++++++++++--- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/primal-count/benches/bench.rs b/primal-count/benches/bench.rs index 4251e934..7c0c3964 100644 --- a/primal-count/benches/bench.rs +++ b/primal-count/benches/bench.rs @@ -1,7 +1,7 @@ #[macro_use] extern crate criterion; extern crate primal_count; -use primal_count::{PrimeCounter, int_quartic_root}; +use primal_count::PrimeCounter; use criterion::{Criterion, ParameterizedBenchmark}; const SIZES: [usize; 6] = [100, 10_000, 100_000, 1_000_000, 10_000_000, 10_000_000_000]; diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index f212ed41..8ef40ea1 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -4,6 +4,10 @@ use std::iter::FromIterator; use crate::util::{int_square_root, int_cubic_root, int_quartic_root}; use std::collections::HashMap; +const MEISSEL_LOOKUP_SIZE: usize = 8; // Number of primes we do the reduce trick for +const SMALL_PRIME_PRODUCTS: [usize; MEISSEL_LOOKUP_SIZE + 1] = [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690]; +const SMALL_TOTIENT_VALUES: [usize; MEISSEL_LOOKUP_SIZE + 1] = [0, 1, 2, 8, 48, 480, 5760, 92160, 1658880]; + /// Generate a vec of primes from 2 up to and including limit /// Leverages the fast sieve in primal to do so fn generate_primes(limit: usize) -> Vec { @@ -20,14 +24,23 @@ fn generate_primes(limit: usize) -> Vec { /// The number of numbers less than m that are coprime to the first n prime numbers /// Get recurrence m_fn(m, n) = m_fn(m, n - 1) - m_fn(m/p_n, n-1) by thinking about p_n, the nth prime fn meissel_fn(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { - if n == 0 { + if n == 0 || m < 2 { return m; } match meissel_cache.get(&(m, n)).map(|entry| entry.clone()){ Some(result) => result, None => { - let value = meissel_fn(m, n-1, &prime_array, meissel_cache) - - meissel_fn(m / prime_array[n-1], n-1, &prime_array, meissel_cache); + // For small values of n, we can decrease the size of m by noting that + // the meissel function is almost periodic with period p_1 * .. * p_n + let mut value = 0; + let mut m_shrunk = m; + if n <= MEISSEL_LOOKUP_SIZE { + value = (m / SMALL_PRIME_PRODUCTS[n]) * SMALL_TOTIENT_VALUES[n]; + m_shrunk = m_shrunk % SMALL_PRIME_PRODUCTS[n]; + } + + // After shrinkage, just apply the recursion + value += meissel_fn(m_shrunk, n-1, &prime_array, meissel_cache) - meissel_fn(m_shrunk / prime_array[n-1], n-1, &prime_array, meissel_cache); meissel_cache.insert((m, n), value); return value; } From 5f18d7c5fd854cca27c54ab7c6130c6323706e2f Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sun, 20 Oct 2019 21:45:41 +0100 Subject: [PATCH 28/36] Paranoia bugfix --- primal-count/src/prime_count.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 8ef40ea1..b14f7ea4 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -6,7 +6,7 @@ use std::collections::HashMap; const MEISSEL_LOOKUP_SIZE: usize = 8; // Number of primes we do the reduce trick for const SMALL_PRIME_PRODUCTS: [usize; MEISSEL_LOOKUP_SIZE + 1] = [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690]; -const SMALL_TOTIENT_VALUES: [usize; MEISSEL_LOOKUP_SIZE + 1] = [0, 1, 2, 8, 48, 480, 5760, 92160, 1658880]; +const SMALL_TOTIENT_VALUES: [usize; MEISSEL_LOOKUP_SIZE + 1] = [1, 1, 2, 8, 48, 480, 5760, 92160, 1658880]; /// Generate a vec of primes from 2 up to and including limit /// Leverages the fast sieve in primal to do so From bd169f4947a1f458ff19ae3274c8733f18e43030 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sun, 20 Oct 2019 22:17:44 +0100 Subject: [PATCH 29/36] Minor improvement --- primal-count/src/prime_count.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index b14f7ea4..df1182dd 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -27,6 +27,9 @@ fn meissel_fn(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut if n == 0 || m < 2 { return m; } + if prime_array[n-1] >= m { + return 1; + } match meissel_cache.get(&(m, n)).map(|entry| entry.clone()){ Some(result) => result, None => { From 3dc98307ed76b680b23ca135c3e7cbff1d1e5fe7 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sun, 20 Oct 2019 22:59:06 +0100 Subject: [PATCH 30/36] Unify api --- benches/bench.rs | 4 ++-- primal-count/benches/bench.rs | 4 ++-- primal-count/src/prime_count.rs | 24 ++++++++++++------------ 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/benches/bench.rs b/benches/bench.rs index ba4b02a9..07a132fd 100644 --- a/benches/bench.rs +++ b/benches/bench.rs @@ -60,7 +60,7 @@ create_benchmarks!{ fn prime_pi(SIZES) { "PrimeCounter" => |b, upto: &usize| { let mut s = primal::PrimeCounter::new(*upto + 1); - b.iter(|| s.primes_below(*upto)); + b.iter(|| s.prime_pi(*upto)); }, "Sieve" => |b, upto: &usize| { let s = primal::Sieve::new(*upto + 1); @@ -70,7 +70,7 @@ create_benchmarks!{ "PrimeCounter with init" => |b, upto: &usize| { b.iter(|| { let mut s = primal::PrimeCounter::new(*upto + 1); - s.primes_below(*upto) + s.prime_pi(*upto) }); }, "Sieve with init" => |b, upto: &usize| { diff --git a/primal-count/benches/bench.rs b/primal-count/benches/bench.rs index 7c0c3964..c8ca6305 100644 --- a/primal-count/benches/bench.rs +++ b/primal-count/benches/bench.rs @@ -34,13 +34,13 @@ create_benchmarks! { fn prime_pi(SIZES) { "PrimeCounter" => |b, upto: &usize| { let mut s = PrimeCounter::new(*upto + 1); - b.iter(|| s.primes_below(*upto)); + b.iter(|| s.prime_pi(*upto)); }, "PrimeCounter with init" => |b, upto: &usize| { b.iter(|| { let mut s = PrimeCounter::new(*upto + 1); - s.primes_below(*upto) + s.prime_pi(*upto) }); }, } diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index df1182dd..8371d5fc 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -146,7 +146,7 @@ impl PrimeCounter { /// /// # Panics /// - /// If the limit in the constructor is smaller than the input of primes_below + /// If the limit in the constructor is smaller than the input of prime_pi /// (Unless it's later been updated with update_limit) /// /// # Examples @@ -154,10 +154,10 @@ impl PrimeCounter { /// ```rust /// # extern crate primal; /// let mut pc = primal::PrimeCounter::new(10_000); - /// assert_eq!(pc.primes_below(100), 25); - /// assert_eq!(pc.primes_below(8166), 1024); + /// assert_eq!(pc.prime_pi(100), 25); + /// assert_eq!(pc.prime_pi(8166), 1024); /// ``` - pub fn primes_below(&mut self, bound: usize) -> usize { + pub fn prime_pi(&mut self, bound: usize) -> usize { primes_less_than(bound, &self.primes, &mut self.prime_cache, &mut self.meissel_cache) } @@ -192,17 +192,17 @@ mod tests { } #[test] - fn test_primes_below() { + fn test_prime_pi() { use crate::prime_count::PrimeCounter; let mut pc = PrimeCounter::new(10_000); - assert_eq!(pc.primes_below(7), 4); - assert_eq!(pc.primes_below(100), 25); - assert_eq!(pc.primes_below(2143), 324); + assert_eq!(pc.prime_pi(7), 4); + assert_eq!(pc.prime_pi(100), 25); + assert_eq!(pc.prime_pi(2143), 324); pc.update_limit(1_000_000_000); - assert_eq!(pc.primes_below(1_000_000), 78_498); - assert_eq!(pc.primes_below(1_000_000_000), 50_847_534); - // assert_eq!(primes_below(1_000_000_000_000), 37_607_912_018); - // assert_eq!(primes_below(1_000_000_000_000_000), 29_844_570_422_669); + assert_eq!(pc.prime_pi(1_000_000), 78_498); + assert_eq!(pc.prime_pi(1_000_000_000), 50_847_534); + // assert_eq!(prime_pi(1_000_000_000_000), 37_607_912_018); + // assert_eq!(prime_pi(1_000_000_000_000_000), 29_844_570_422_669); } } \ No newline at end of file From 2a85ef25b0a99f7b9cd9827a83b58a127ffe6c6a Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Sun, 20 Oct 2019 23:02:16 +0100 Subject: [PATCH 31/36] Get rid of unneeded line --- primal-count/src/prime_count.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 8371d5fc..733f1f5f 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -66,7 +66,6 @@ fn primes_less_than(bound: usize, primes: &Vec, prime_cache: &mut HashMap Ok(idx) => idx+1, Err(idx) => idx, }; - let result = result; prime_cache.insert(bound, result); return result; } From 8079235691cd6af0ae5c3956cbd9aaa937268592 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 30 Oct 2019 00:28:07 +0000 Subject: [PATCH 32/36] Speed things up by better evaluation of meissel and cache a larger sieve --- primal-count/src/prime_count.rs | 218 +++++++++++++++++++------------- 1 file changed, 127 insertions(+), 91 deletions(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 733f1f5f..2cd63080 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -1,8 +1,9 @@ extern crate primal_sieve; use std::iter::FromIterator; -use crate::util::{int_square_root, int_cubic_root, int_quartic_root}; +use super::util; use std::collections::HashMap; +use std::cmp; const MEISSEL_LOOKUP_SIZE: usize = 8; // Number of primes we do the reduce trick for const SMALL_PRIME_PRODUCTS: [usize; MEISSEL_LOOKUP_SIZE + 1] = [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690]; @@ -21,100 +22,21 @@ fn generate_primes(limit: usize) -> Vec { return Vec::from_iter(sieve_iter); } -/// The number of numbers less than m that are coprime to the first n prime numbers -/// Get recurrence m_fn(m, n) = m_fn(m, n - 1) - m_fn(m/p_n, n-1) by thinking about p_n, the nth prime -fn meissel_fn(m: usize, n: usize, prime_array: &Vec, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { - if n == 0 || m < 2 { - return m; - } - if prime_array[n-1] >= m { - return 1; - } - match meissel_cache.get(&(m, n)).map(|entry| entry.clone()){ - Some(result) => result, - None => { - // For small values of n, we can decrease the size of m by noting that - // the meissel function is almost periodic with period p_1 * .. * p_n - let mut value = 0; - let mut m_shrunk = m; - if n <= MEISSEL_LOOKUP_SIZE { - value = (m / SMALL_PRIME_PRODUCTS[n]) * SMALL_TOTIENT_VALUES[n]; - m_shrunk = m_shrunk % SMALL_PRIME_PRODUCTS[n]; - } - - // After shrinkage, just apply the recursion - value += meissel_fn(m_shrunk, n-1, &prime_array, meissel_cache) - meissel_fn(m_shrunk / prime_array[n-1], n-1, &prime_array, meissel_cache); - meissel_cache.insert((m, n), value); - return value; - } - } -} - -/// Find the number of primes less than bound using the Meissel-Lehmer method -/// Leverages caching to speed up the recursive calls -fn primes_less_than(bound: usize, primes: &Vec, prime_cache: &mut HashMap, meissel_cache: &mut HashMap<(usize, usize), usize>) -> usize { - // First check if it's in the cache already - match prime_cache.get(&bound).map(|entry| entry.clone()){ - Some(value) => value, - None => { - // The meat of the function - if bound < 2 { - return 0; - } else if bound <= primes[primes.len()-1] { - let result = match primes.binary_search(&bound) - { - Ok(idx) => idx+1, - Err(idx) => idx, - }; - prime_cache.insert(bound, result); - return result; - } - - let sqrt_bound = int_square_root(bound); - - let nprimes_below_4thr = primes_less_than(int_quartic_root(bound), primes, prime_cache, meissel_cache); - let nprimes_below_3rdr = primes_less_than(int_cubic_root(bound), primes, prime_cache, meissel_cache); - let nprimes_below_2ndr = primes_less_than(sqrt_bound, primes, prime_cache, meissel_cache); - - // Issues with underflow here if nprimes_below_2ndr + nprimes_below_4thr < 2 - // Dealt with by populating the offending (small) values in the cache at the top level - let mut result = ((nprimes_below_2ndr + nprimes_below_4thr - 2) * (nprimes_below_2ndr - nprimes_below_4thr + 1)) / 2; - result += meissel_fn(bound, nprimes_below_4thr, &primes, meissel_cache); - - for i in nprimes_below_4thr..nprimes_below_2ndr { - let ith_prime = primes[i]; - result -= primes_less_than(bound / ith_prime, primes, prime_cache, meissel_cache); - if i < nprimes_below_3rdr { - let bi = primes_less_than(int_square_root(bound / ith_prime), primes, prime_cache, meissel_cache); - for j in i..bi { - let jth_prime = primes[j]; - result -= primes_less_than(bound / ith_prime / jth_prime, primes, prime_cache, meissel_cache) - j; - } - } - } - - // Caching - prime_cache.insert(bound, result); - return result; - } - } -} - /// Memoized combinatorial prime counting function /// Basic idea here: https://en.wikipedia.org/wiki/Meissel%E2%80%93Lehmer_algorithm /// The "Meissel Function" here is phi on that Wikipedia page pub struct PrimeCounter { limit: usize, primes: Vec, - prime_cache: HashMap, + pi_cache: HashMap, meissel_cache: HashMap<(usize, usize), usize> } impl PrimeCounter { /// Create a new PrimeCounter instance, which generates all the primes up to sqrt(limit) pub fn new(limit: usize) -> PrimeCounter { - let mut prime_cache = HashMap::new(); - let primes = generate_primes(int_square_root(limit)); + let mut pi_cache = HashMap::new(); + let primes = generate_primes(util::int_cubic_root(limit * limit)); // Insert primes <= 10 - this is mainly to deal with underflow issues later for n in 0..=10 { @@ -125,11 +47,12 @@ impl PrimeCounter { 7..=10 => 4, 0..=1 | _ => 0, // N.B. _ never hit }; - prime_cache.insert(n, nprimes); + pi_cache.insert(n, nprimes); } + // let meissel_cache = meissel::generate_meissel_lookup(6); let meissel_cache = HashMap::new(); - PrimeCounter {limit, primes, prime_cache, meissel_cache} + PrimeCounter {limit, primes, pi_cache, meissel_cache} } /// Updates the limit - to be used when you want to make the prime cache larger @@ -138,7 +61,7 @@ impl PrimeCounter { /// object is large enough. pub fn update_limit(&mut self, limit: usize) { self.limit = limit; - self.primes = generate_primes(int_square_root(limit)); + self.primes = generate_primes(util::int_square_root(limit)); } /// The number of primes that are at least `bound` @@ -157,7 +80,7 @@ impl PrimeCounter { /// assert_eq!(pc.prime_pi(8166), 1024); /// ``` pub fn prime_pi(&mut self, bound: usize) -> usize { - primes_less_than(bound, &self.primes, &mut self.prime_cache, &mut self.meissel_cache) + self.primes_less_than(bound) } /// The number of numbers less than `m` that are coprime to the first `n` prime numbers @@ -175,7 +98,118 @@ impl PrimeCounter { /// assert_eq!(pc.meissel_fn(1234, 5), 256); /// ``` pub fn meissel_fn(&mut self, m: usize, n: usize) -> usize { - meissel_fn(m, n, &self.primes, &mut self.meissel_cache) + self.meissel_fn_large(m, n) + } + + /// The number of numbers less than m that are coprime to the first n prime numbers + /// Get recurrence m_fn(m, n) = m_fn(m, n - 1) - m_fn(m/p_n, n-1) by thinking about p_n, the nth prime + fn meissel_fn_small(&mut self, m: usize, n: usize) -> usize { + if n == 0 || m < 2 { + return m; + } + if self.primes[n-1] >= m { + return 1; + } + match self.meissel_cache.get(&(m, n)).map(|entry| entry.clone()){ + Some(result) => result, + None => { + // For small values of n, we can decrease the size of m by noting that + // the meissel function is almost periodic with period p_1 * .. * p_n + let mut value = 0; + let mut m_shrunk = m; + if n <= MEISSEL_LOOKUP_SIZE { + value = (m / SMALL_PRIME_PRODUCTS[n]) * SMALL_TOTIENT_VALUES[n]; + m_shrunk = m_shrunk % SMALL_PRIME_PRODUCTS[n]; + } + + // After shrinkage, just apply the recursion + value += self.meissel_fn_small(m_shrunk, n-1) - self.meissel_fn_small(m_shrunk / self.primes[n-1], n-1); + self.meissel_cache.insert((m, n), value); + return value; + } + } + } + + // Here we make use of the fact that m_fn(m, n) can be evaluated with fewer calls by + // taking advantage of n being large + fn meissel_fn_large(&mut self, m: usize, n: usize) -> usize { + if n <= MEISSEL_LOOKUP_SIZE { + return self.meissel_fn_small(m, n); + } + if self.primes[n-1] >= m { + return 1; + } + match self.meissel_cache.get(&(m, n)).map(|entry| entry.clone()){ + Some(result) => result, + None => { + let mut result = self.meissel_fn_small(m, MEISSEL_LOOKUP_SIZE); + let sqrt_m = util::int_square_root(m); + let bound = self.primes_less_than(sqrt_m); + let largest_prime_index = cmp::min(cmp::max(bound, MEISSEL_LOOKUP_SIZE), n); + result = result + largest_prime_index - n; + + for idx in MEISSEL_LOOKUP_SIZE..largest_prime_index { + result -= self.meissel_fn_large(m / self.primes[idx], idx); + } + self.meissel_cache.insert((m, n), result); + return result; + } + } + } + + /// Find the number of primes less than bound using the Meissel-Lehmer method + /// Leverages caching to speed up the recursive calls + fn primes_less_than(&mut self, bound: usize) -> usize { + // First check if it's in the cache already + match self.pi_cache.get(&bound).map(|entry| entry.clone()){ + Some(value) => value, + None => { + // The meat of the function + if bound < 2 { + return 0; + } else if bound <= self.primes[self.primes.len()-1] { + let result = match self.primes.binary_search(&bound) + { + Ok(idx) => idx+1, + Err(idx) => idx, + }; + self.pi_cache.insert(bound, result); + return result; + } + + let sqrt_bound = util::int_square_root(bound); + let quartic_bound = util::int_quartic_root(bound); + + let nprimes_below_4thr = self.primes_less_than(quartic_bound); + let nprimes_below_3rdr = self.primes_less_than(util::int_cubic_root(bound)); + let nprimes_below_2ndr = self.primes_less_than(sqrt_bound); + + // Issues with underflow here if nprimes_below_2ndr + nprimes_below_4thr < 2 + // Dealt with by populating the offending (small) values in the cache at the top level + let mut result = ((nprimes_below_2ndr + nprimes_below_4thr - 2) * (nprimes_below_2ndr - nprimes_below_4thr + 1)) / 2; + result += self.meissel_fn_large(bound, nprimes_below_4thr); + + for i in nprimes_below_4thr..nprimes_below_2ndr { + let ith_prime = self.primes[i]; + // Need to make this faster + result -= self.primes_less_than(bound / ith_prime); + if i < nprimes_below_3rdr { + let bi = self.primes_less_than(util::int_square_root(bound / ith_prime)); + for j in i..bi { + let jth_prime = self.primes[j]; + // p_i, p_j > bound^(1/4), so bound / (p_i * p_j) < bound ^ (1/2) + // Hence, this is a cheap lookup, and thus why we don't bother + // optimising this further... + result -= self.primes_less_than(bound / ith_prime / jth_prime) - j; + } + } + } + + // Caching + self.pi_cache.insert(bound, result); + return result; + } + } } } @@ -184,10 +218,11 @@ mod tests { #[test] fn test_meissel_fn() { use crate::prime_count::PrimeCounter; - let mut pc = PrimeCounter::new(10_000); + let mut pc = PrimeCounter::new(10_000_000); assert_eq!(pc.meissel_fn(30, 8), 3); assert_eq!(pc.meissel_fn(100, 1), 50); assert_eq!(pc.meissel_fn(100, 25), 1); + assert_eq!(pc.meissel_fn(10_000_000, 50), 1025747); } #[test] @@ -201,7 +236,8 @@ mod tests { pc.update_limit(1_000_000_000); assert_eq!(pc.prime_pi(1_000_000), 78_498); assert_eq!(pc.prime_pi(1_000_000_000), 50_847_534); - // assert_eq!(prime_pi(1_000_000_000_000), 37_607_912_018); - // assert_eq!(prime_pi(1_000_000_000_000_000), 29_844_570_422_669); + // pc.update_limit(1_000_000_000_000); + // assert_eq!(pc.prime_pi(1_000_000_000_000), 37_607_912_018); + // assert_eq!(pc.prime_pi(1_000_000_000_000_000), 29_844_570_422_669); } } \ No newline at end of file From 405d160f5a6a8b186a10aefe12b0cb92d50e59c7 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Wed, 30 Oct 2019 00:30:01 +0000 Subject: [PATCH 33/36] Whitespace --- primal-count/src/prime_count.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 2cd63080..23659478 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -130,10 +130,10 @@ impl PrimeCounter { } } - // Here we make use of the fact that m_fn(m, n) can be evaluated with fewer calls by + // Here we make use of the fact that m_fn(m, n) can be evaluated with fewer calls by // taking advantage of n being large fn meissel_fn_large(&mut self, m: usize, n: usize) -> usize { - if n <= MEISSEL_LOOKUP_SIZE { + if n <= MEISSEL_LOOKUP_SIZE { return self.meissel_fn_small(m, n); } if self.primes[n-1] >= m { @@ -176,7 +176,7 @@ impl PrimeCounter { self.pi_cache.insert(bound, result); return result; } - + let sqrt_bound = util::int_square_root(bound); let quartic_bound = util::int_quartic_root(bound); @@ -198,7 +198,7 @@ impl PrimeCounter { for j in i..bi { let jth_prime = self.primes[j]; // p_i, p_j > bound^(1/4), so bound / (p_i * p_j) < bound ^ (1/2) - // Hence, this is a cheap lookup, and thus why we don't bother + // Hence, this is a cheap lookup, and thus why we don't bother // optimising this further... result -= self.primes_less_than(bound / ith_prime / jth_prime) - j; } From 91a00ff02f51654caf2735102192f581fe8be4b8 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Fri, 1 Nov 2019 22:06:18 +0000 Subject: [PATCH 34/36] Rust fmt --- primal-count/benches/bench.rs | 14 +++++-- primal-count/src/lib.rs | 2 +- primal-count/src/prime_count.rs | 65 +++++++++++++++++++-------------- primal-count/src/util.rs | 1 - 4 files changed, 49 insertions(+), 33 deletions(-) diff --git a/primal-count/benches/bench.rs b/primal-count/benches/bench.rs index c8ca6305..598299b1 100644 --- a/primal-count/benches/bench.rs +++ b/primal-count/benches/bench.rs @@ -1,10 +1,18 @@ #[macro_use] extern crate criterion; extern crate primal_count; -use primal_count::PrimeCounter; use criterion::{Criterion, ParameterizedBenchmark}; +use primal_count::PrimeCounter; -const SIZES: [usize; 6] = [100, 10_000, 100_000, 1_000_000, 10_000_000, 10_000_000_000]; +const SIZES: [usize; 7] = [ + 100, + 10_000, + 100_000, + 1_000_000, + 10_000_000, + 10_000_000_000, + 100_000_000_000, +]; macro_rules! create_benchmarks { ($( @@ -36,7 +44,7 @@ create_benchmarks! { let mut s = PrimeCounter::new(*upto + 1); b.iter(|| s.prime_pi(*upto)); }, - + "PrimeCounter with init" => |b, upto: &usize| { b.iter(|| { let mut s = PrimeCounter::new(*upto + 1); diff --git a/primal-count/src/lib.rs b/primal-count/src/lib.rs index 2bb864ad..817c341b 100644 --- a/primal-count/src/lib.rs +++ b/primal-count/src/lib.rs @@ -1,4 +1,4 @@ mod prime_count; mod util; pub use prime_count::PrimeCounter; -pub use util::{int_square_root, int_cubic_root, int_quartic_root}; \ No newline at end of file +pub use util::{int_cubic_root, int_quartic_root, int_square_root}; diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 23659478..511f8dd6 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -2,12 +2,14 @@ extern crate primal_sieve; use std::iter::FromIterator; use super::util; -use std::collections::HashMap; use std::cmp; +use std::collections::HashMap; -const MEISSEL_LOOKUP_SIZE: usize = 8; // Number of primes we do the reduce trick for -const SMALL_PRIME_PRODUCTS: [usize; MEISSEL_LOOKUP_SIZE + 1] = [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690]; -const SMALL_TOTIENT_VALUES: [usize; MEISSEL_LOOKUP_SIZE + 1] = [1, 1, 2, 8, 48, 480, 5760, 92160, 1658880]; +const MEISSEL_LOOKUP_SIZE: usize = 8; // Number of primes we do the reduce trick for +const SMALL_PRIME_PRODUCTS: [usize; MEISSEL_LOOKUP_SIZE + 1] = + [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690]; +const SMALL_TOTIENT_VALUES: [usize; MEISSEL_LOOKUP_SIZE + 1] = + [1, 1, 2, 8, 48, 480, 5760, 92160, 1658880]; /// Generate a vec of primes from 2 up to and including limit /// Leverages the fast sieve in primal to do so @@ -29,7 +31,7 @@ pub struct PrimeCounter { limit: usize, primes: Vec, pi_cache: HashMap, - meissel_cache: HashMap<(usize, usize), usize> + meissel_cache: HashMap<(usize, usize), usize>, } impl PrimeCounter { @@ -41,18 +43,23 @@ impl PrimeCounter { // Insert primes <= 10 - this is mainly to deal with underflow issues later for n in 0..=10 { let nprimes = match n { - 2 => 1, - 3..=4 => 2, - 5..=6 => 3, - 7..=10 => 4, - 0..=1 | _ => 0, // N.B. _ never hit + 2 => 1, + 3..=4 => 2, + 5..=6 => 3, + 7..=10 => 4, + 0..=1 | _ => 0, // N.B. _ never hit }; pi_cache.insert(n, nprimes); } // let meissel_cache = meissel::generate_meissel_lookup(6); let meissel_cache = HashMap::new(); - - PrimeCounter {limit, primes, pi_cache, meissel_cache} + + PrimeCounter { + limit, + primes, + pi_cache, + meissel_cache, + } } /// Updates the limit - to be used when you want to make the prime cache larger @@ -63,9 +70,9 @@ impl PrimeCounter { self.limit = limit; self.primes = generate_primes(util::int_square_root(limit)); } - + /// The number of primes that are at least `bound` - /// + /// /// # Panics /// /// If the limit in the constructor is smaller than the input of prime_pi @@ -88,7 +95,7 @@ impl PrimeCounter { /// # Panics /// /// If the `n`th prime is larger than `limit` - /// + /// /// # Examples /// /// ```rust @@ -107,10 +114,10 @@ impl PrimeCounter { if n == 0 || m < 2 { return m; } - if self.primes[n-1] >= m { + if self.primes[n - 1] >= m { return 1; } - match self.meissel_cache.get(&(m, n)).map(|entry| entry.clone()){ + match self.meissel_cache.get(&(m, n)).map(|entry| entry.clone()) { Some(result) => result, None => { // For small values of n, we can decrease the size of m by noting that @@ -123,7 +130,8 @@ impl PrimeCounter { } // After shrinkage, just apply the recursion - value += self.meissel_fn_small(m_shrunk, n-1) - self.meissel_fn_small(m_shrunk / self.primes[n-1], n-1); + value += self.meissel_fn_small(m_shrunk, n - 1) + - self.meissel_fn_small(m_shrunk / self.primes[n - 1], n - 1); self.meissel_cache.insert((m, n), value); return value; } @@ -136,10 +144,10 @@ impl PrimeCounter { if n <= MEISSEL_LOOKUP_SIZE { return self.meissel_fn_small(m, n); } - if self.primes[n-1] >= m { + if self.primes[n - 1] >= m { return 1; } - match self.meissel_cache.get(&(m, n)).map(|entry| entry.clone()){ + match self.meissel_cache.get(&(m, n)).map(|entry| entry.clone()) { Some(result) => result, None => { let mut result = self.meissel_fn_small(m, MEISSEL_LOOKUP_SIZE); @@ -161,16 +169,15 @@ impl PrimeCounter { /// Leverages caching to speed up the recursive calls fn primes_less_than(&mut self, bound: usize) -> usize { // First check if it's in the cache already - match self.pi_cache.get(&bound).map(|entry| entry.clone()){ + match self.pi_cache.get(&bound).map(|entry| entry.clone()) { Some(value) => value, None => { // The meat of the function if bound < 2 { return 0; - } else if bound <= self.primes[self.primes.len()-1] { - let result = match self.primes.binary_search(&bound) - { - Ok(idx) => idx+1, + } else if bound <= self.primes[self.primes.len() - 1] { + let result = match self.primes.binary_search(&bound) { + Ok(idx) => idx + 1, Err(idx) => idx, }; self.pi_cache.insert(bound, result); @@ -186,7 +193,9 @@ impl PrimeCounter { // Issues with underflow here if nprimes_below_2ndr + nprimes_below_4thr < 2 // Dealt with by populating the offending (small) values in the cache at the top level - let mut result = ((nprimes_below_2ndr + nprimes_below_4thr - 2) * (nprimes_below_2ndr - nprimes_below_4thr + 1)) / 2; + let mut result = ((nprimes_below_2ndr + nprimes_below_4thr - 2) + * (nprimes_below_2ndr - nprimes_below_4thr + 1)) + / 2; result += self.meissel_fn_large(bound, nprimes_below_4thr); for i in nprimes_below_4thr..nprimes_below_2ndr { @@ -214,7 +223,7 @@ impl PrimeCounter { } #[cfg(test)] -mod tests { +mod tests { #[test] fn test_meissel_fn() { use crate::prime_count::PrimeCounter; @@ -240,4 +249,4 @@ mod tests { // assert_eq!(pc.prime_pi(1_000_000_000_000), 37_607_912_018); // assert_eq!(pc.prime_pi(1_000_000_000_000_000), 29_844_570_422_669); } -} \ No newline at end of file +} diff --git a/primal-count/src/util.rs b/primal-count/src/util.rs index a2682fa8..96443bc7 100644 --- a/primal-count/src/util.rs +++ b/primal-count/src/util.rs @@ -11,7 +11,6 @@ pub fn int_square_root(value: usize) -> usize { /// Returns the largest integer at least cbrt(n) using Newton's method pub fn int_cubic_root(value: usize) -> usize { - let mut x = value; let mut y = (2 * value + 1) / 3; // Plus 1 only to protect against division by 0 later while y < x { From e9169854fd6c512f8e5837e326f5a64c202a5783 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Fri, 1 Nov 2019 22:09:18 +0000 Subject: [PATCH 35/36] More cleaning --- primal-count/src/prime_count.rs | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 511f8dd6..9617c658 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -21,7 +21,7 @@ fn generate_primes(limit: usize) -> Vec { // 2) Quick counting of number of primes below a value, achieved with a binary search // Experiments replacing 1) or 2) with the methods in sieve seem to significantly // slow things down for larger numbers - return Vec::from_iter(sieve_iter); + sieve_iter.collect() } /// Memoized combinatorial prime counting function @@ -132,8 +132,9 @@ impl PrimeCounter { // After shrinkage, just apply the recursion value += self.meissel_fn_small(m_shrunk, n - 1) - self.meissel_fn_small(m_shrunk / self.primes[n - 1], n - 1); + self.meissel_cache.insert((m, n), value); - return value; + value } } } @@ -159,8 +160,9 @@ impl PrimeCounter { for idx in MEISSEL_LOOKUP_SIZE..largest_prime_index { result -= self.meissel_fn_large(m / self.primes[idx], idx); } + self.meissel_cache.insert((m, n), result); - return result; + result } } } @@ -216,7 +218,7 @@ impl PrimeCounter { // Caching self.pi_cache.insert(bound, result); - return result; + result } } } From 6f8eebfa25a94d4acc1801c26485240c21c42c52 Mon Sep 17 00:00:00 2001 From: Kaimyn Chapman Date: Fri, 1 Nov 2019 22:10:42 +0000 Subject: [PATCH 36/36] silly me --- primal-count/src/prime_count.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/primal-count/src/prime_count.rs b/primal-count/src/prime_count.rs index 9617c658..b811d282 100644 --- a/primal-count/src/prime_count.rs +++ b/primal-count/src/prime_count.rs @@ -1,5 +1,4 @@ extern crate primal_sieve; -use std::iter::FromIterator; use super::util; use std::cmp;