HiveBrain v1.2.0
Get Started
← Back to all entries
patternrustMinor

Counting primes less than n

Submitted by: @import:stackexchange-codereview··
0
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 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.

  • 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 in vec![]



  • 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! (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 ]) -> 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.

  • 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.