patternrustMinor
Counting primes less than n
Viewed 0 times
thancountinglessprimes
Problem
The following code was written to be a submission to this challenge on PPCG.
It uses this algorithm, the Meissal-Lehmer method. I follow the Wikipedia entry pretty naively, except I tried to memoize the phi function.
Please see @Shepmaster's answer for cleaned up, more idiomatic Rust.
I am looking for ways to improve speed. Most improvements I believe will come from optimizing the phi function, which is used recursively. Answers about making my code more idiomatic Rust are also welcome.
I've been compiling with
I'm kind of new to the site, so please let me know how to improve this post as a whole too, if necessary
It uses this algorithm, the Meissal-Lehmer method. I follow the Wikipedia entry pretty naively, except I tried to memoize the phi function.
Please see @Shepmaster's answer for cleaned up, more idiomatic Rust.
I am looking for ways to improve speed. Most improvements I believe will come from optimizing the phi function, which is used recursively. Answers about making my code more idiomatic Rust are also welcome.
I've been compiling with
cargo build --release. The code runs in ~3 sec for input 1,000,000 and a lot longer for 10,000,000.use std::env;
fn main() {
let args: Vec = env::args().collect();
let x: usize = (args[1]).trim().parse()
.expect("expected a number");
let root: usize = (x as f64).sqrt() as usize;
let y: usize = (x as f64).powf(0.3333333333333) as usize + 1;
let sieve_size: usize = x/y+2;//y+1;
let mut sieve: Vec = vec![true; sieve_size ];
let mut primes: Vec = vec![0; sieve_size ];
sieve[0] = false; sieve[1] = false;
let mut a: usize = 0;
let mut num_primes = 1;
let mut num_primes_smaller_root: usize = 0;
// find all primes up to x/y ~ x^2/3 aka sieve_size
for i in 2..sieve_size {
if sieve[i] {
if i> = vec![ vec![0; a+1 ];x+1 ];
let mut p_2: usize = 0;
for i in a+1..num_primes_smaller_root+1 {
let mut first_term: usize = 0;
while primes[first_term], phi_results:&mut Vec>)->usize {
if b==0 {
return x;
}
if phi_results[x][b] != 0 {
return phi_results[x][b];
}
let value: usize = phi(x,b-1,primes,phi_results) - phi(x/primes[b],b-1,primes,phi_results);
phi_results[x][b] = value;
return value;
}I'm kind of new to the site, so please let me know how to improve this post as a whole too, if necessary
Solution
Here's my first pass on stylistic feedback.
I believe that you can change some of the loops to use more iterators, but I haven't figured those out yet.
When I run this code with 1000000
phi(1000000) = 78498
real 0m0.307s
user 0m0.204s
sys 0m0.102s
- Unneeded parenthesis around
args[1]
- Redundant type specifications, everywhere
- Remove commented-out code
- Add spaces around binary operators (
/,+,
- Don't have double semicolons at the end of the line
- Don't add space inside []
, like invec![]
- Split multiple statements on different lines
- Add spaces after commas
- Add spaces after colons
- Basically never use &Vec
, use&[T]instead
- Don't use an explicit return
at the end of function
- Typo in the println!
(pivsphi)
I believe that you can change some of the loops to use more iterators, but I haven't figured those out yet.
use std::env;
fn main() {
let args: Vec = env::args().collect();
let x: usize = args[1].trim().parse().expect("expected a number");
let root = (x as f64).sqrt() as usize;
let y = (x as f64).powf(0.3333333333333) as usize + 1;
let sieve_size = x / y + 2;
let mut sieve = vec![true; sieve_size];
let mut primes = vec![0; sieve_size];
sieve[0] = false;
sieve[1] = false;
let mut a = 0;
let mut num_primes = 1;
let mut num_primes_smaller_root = 0;
// find all primes up to x/y ~ x^2/3 aka sieve_size
for i in 2..sieve_size {
if sieve[i] {
if i ]) -> usize {
if b == 0 {
return x;
}
if phi_results[x][b] != 0 {
return phi_results[x][b];
}
let value = phi(x, b - 1, primes, phi_results) - phi(x / primes[b], b - 1, primes, phi_results);
phi_results[x][b] = value;
value
}When I run this code with 1000000
it takes ~300 milliseconds:
$ cargo build --release && time ./target/release/sieve 1000000phi(1000000) = 78498
real 0m0.307s
user 0m0.204s
sys 0m0.102s
Note that you need to build the code and time it in a separate pass, to not include the time spent compiling.
Performance-wise, I'm not seeing many obvious changes.
- Try to avoid direct slice access when feasible (
vector[0]). These accesses are checked at run time to avoid walking off into undefined memory. Iterators guarantee that no undefined memory will be accessed and can sidestep that small performance bit.
Iterator adapters can also allow removing mutability:
let mut p_2 = 0;
for (i, zz) in primes[a + 1..num_primes_smaller_root + 1].iter().enumerate() {
let i = i + a + 1;
let first_term = primes.iter().take_while(|&&p| p <= x / zz).count();
p_2 += (first_term - 1) - i + 1;
}
Then, move that arithmetic on i` to the end, and apply more adapters to remove more mutability:let interesting_primes = primes[a + 1..num_primes_smaller_root + 1].iter();
let p_2 =
interesting_primes
.map(|ip| primes.iter().take_while(|&&p| p <= x / ip).count())
.enumerate()
.map(|(i, v)| v - 1 - i - a)
.fold(0, |acc, v| acc + v);Code Snippets
use std::env;
fn main() {
let args: Vec<_> = env::args().collect();
let x: usize = args[1].trim().parse().expect("expected a number");
let root = (x as f64).sqrt() as usize;
let y = (x as f64).powf(0.3333333333333) as usize + 1;
let sieve_size = x / y + 2;
let mut sieve = vec![true; sieve_size];
let mut primes = vec![0; sieve_size];
sieve[0] = false;
sieve[1] = false;
let mut a = 0;
let mut num_primes = 1;
let mut num_primes_smaller_root = 0;
// find all primes up to x/y ~ x^2/3 aka sieve_size
for i in 2..sieve_size {
if sieve[i] {
if i <= root {
if i <= y {
a += 1;
}
num_primes_smaller_root += 1;
}
primes[num_primes] = i;
num_primes += 1;
let mut multiples = i;
while multiples < sieve_size {
sieve[multiples] = false;
multiples += i;
}
}
}
let mut phi_results = vec![vec![0; a + 1]; x + 1];
let mut p_2 = 0;
for i in a + 1..num_primes_smaller_root + 1 {
let mut first_term = 0;
while primes[first_term] < x / primes[i] {
first_term += 1;
}
if primes[first_term] == x / primes[i] {
first_term += 1;
}
p_2 += (first_term - 1) - i + 1;
}
println!("pi({}) = {}", x, phi(x, a, &primes, &mut phi_results) + a - 1 - p_2);
}
fn phi(x: usize, b: usize, primes: &[usize], phi_results: &mut [Vec<usize>]) -> usize {
if b == 0 {
return x;
}
if phi_results[x][b] != 0 {
return phi_results[x][b];
}
let value = phi(x, b - 1, primes, phi_results) - phi(x / primes[b], b - 1, primes, phi_results);
phi_results[x][b] = value;
value
}let mut p_2 = 0;
for (i, zz) in primes[a + 1..num_primes_smaller_root + 1].iter().enumerate() {
let i = i + a + 1;
let first_term = primes.iter().take_while(|&&p| p <= x / zz).count();
p_2 += (first_term - 1) - i + 1;
}let interesting_primes = primes[a + 1..num_primes_smaller_root + 1].iter();
let p_2 =
interesting_primes
.map(|ip| primes.iter().take_while(|&&p| p <= x / ip).count())
.enumerate()
.map(|(i, v)| v - 1 - i - a)
.fold(0, |acc, v| acc + v);Context
StackExchange Code Review Q#122120, answer score: 4
Revisions (0)
No revisions yet.