Counting primes less than n

Question

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 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<bool> = vec![true; sieve_size ];
    let mut primes: Vec<usize> = 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<= root {
                if i<=y {
                    a+=1;
                }
                num_primes_smaller_root+=1;
            }


            primes[num_primes] = i;
            num_primes +=1;
            let mut multiples: usize = i;
            while multiples < sieve_size {
                sieve[multiples] = false;
                multiples += i;
            }
        }
    }



    let mut phi_results: Vec<Vec<usize>> = 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]< x/primes[i] {
            //println!("checking {} vs {}...",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!("y={}; a={}",y,a);
    println!("P2 = {}",p_2);
    println!("phi({},{}) = {}",x, a, phi(x,a,&primes));
    */ 
    println!("pi({}) = {}", x, phi(x,a,&primes,&mut phi_results)+a-1-p_2);
}

fn phi(x:usize, b:usize, primes:&Vec<usize>, phi_results:&mut Vec<Vec<usize>>)->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

Answer 1

Here's my first pass on stylistic feedback.

1. Unneeded parenthesis around args[1]
2. Redundant type specifications, everywhere
3. Remove commented-out code
4. Add spaces around binary operators (/, +, <=, +=)
5. Don't have double semicolons at the end of the line
6. Don't add space inside [], like in vec![]
7. Split multiple statements on different lines
8. Add spaces after commas
9. Add spaces after colons
10. Basically never use &Vec<T>, use &[T] instead
11. Don't use an explicit return at the end of function
12. Typo in the println! (pi vs phi)

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 <= 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
}

---

When I run this code with 1000000 it takes ~300 milliseconds:

$ cargo build --release && time ./target/release/sieve 1000000
phi(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.

1. 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);

Answer 2

This review is on top of Shepmaster's.

Your main problem is when you do

let mut phi_results = vec![vec![0; a + 1]; x + 1];

you're conflating "returns 0" with "not cached", which means you can't cache the many returned zero values. A quick fix is to use a better sentinel

const EMPTY: usize = std::usize::MAX;
let mut phi_results = vec![vec![EMPTY; a + 1]; x + 1];

The other problem is that in the above only 20852 of 479982048 entries end up actually being used. A sparse data structure is a lot better of a bet. Using a HashMap made it possible to calculate very large values in reasonable time frames (aka. 30 seconds for 100 billion).

It also helps to conditional your insert on the size of x, given large x are unlikely to be repeated frequently enough to be worth the cost of caching them. This only gives small gains but it enables using a more direct caching layout, like vec![EMPTY; (a + 1) * MAX_X].

Then I realized I'd forgotten to add Shepmaster's iterator-based calculation of p_2, which removed most of the rest of the overhead. It takes less than 3 seconds for 100 billion for me.

use std::env;

const EMPTY: usize = std::usize::MAX;
const MAX_X: usize = 800;

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 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);

    let mut phi_results = vec![EMPTY; (a + 1) * MAX_X];
    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 [usize]) -> usize {
    if b == 0 {
        return x;
    }

    if x < MAX_X && phi_results[x + b * MAX_X] != EMPTY {
        return phi_results[x + b * MAX_X];
    }

    let value = phi(x, b - 1, primes, phi_results) - phi(x / primes[b], b - 1, primes, phi_results);
    if x < MAX_X {
        phi_results[x + b * MAX_X] = value;
    }
    value
}