Multithreaded segmented Sieve of Eratosthenes

Question

I am fairly new to Rust and thought a good way to practice would be to write a multithreaded segmented Sieve of Eratosthenes. It performs ok (searches ten-billion numbers in about 11 seconds on my system, I tried one-hundred-billion but it couldn't allocate enough memory). I was just wondering if there was any clear improvements that could be made in my code? I'm mostly hoping to get some Rust and formatting tips, but am very open to math improvements as well. Thank you in advance!

use std::{
    sync::{Arc, RwLock},
    thread, time::Instant, cmp::min
};

pub fn main() {
    let now = Instant::now();
    let primes = threaded_segmented_sieve_of_eratosthenes(10000000000);
    let finished =now.elapsed().as_secs_f64();
    
    //println!("{:?}\n", primes);
    println!("found {} primes in {}s", primes.len(), finished);
}

fn threaded_segmented_sieve_of_eratosthenes(limit:usize) -> Vec<usize> {
    let threads = num_cpus::get();
    explicit_threaded_segmented_sieve_of_eratosthenes(limit, threads)
}

fn explicit_threaded_segmented_sieve_of_eratosthenes(limit:usize, threads:usize) -> Vec<usize> {
    let sqrt_of_limit = (limit as f64).sqrt().ceil() as usize;
    let early_primes = if limit <= 230 {
        Arc::new(RwLock::new(vec![2,3,5,7,11,13,17]))
    } else {
        Arc::new(RwLock::new(threaded_segmented_sieve_of_eratosthenes(sqrt_of_limit)))
    };

    let mut thread_handles = Vec::new();

    let thread_spacing = (limit - sqrt_of_limit) / threads;
    let segment_size = min(100000, sqrt_of_limit);

    for i in 0..threads {
        let early_primes = early_primes.clone();
        let lowest_checked = sqrt_of_limit + i *thread_spacing;
        let mut highest_checked = lowest_checked + thread_spacing;

        if i == threads - 1 {
            highest_checked = limit;
        }
        
        thread_handles.push(thread::spawn(move|| {
            eratosthenes_segment_thread(early_primes, lowest_checked, highest_checked, segment_size)
        }));
    }

    let mut new_primes = Vec::new();

    for handle in thread_handles {
        new_primes.append(&mut handle.join().unwrap());
    }

    let mut early_primes = early_primes.write().unwrap();
    early_primes.append(&mut new_primes);
    return early_primes.to_owned();
}

fn eratosthenes_segment_thread(early_primes:Arc<RwLock<Vec<usize>>>, lowest_checked: usize, highest_checked: usize, segment_size: usize) -> Vec<usize> {
    let mut returned_primes = Vec::new();

    let mut lower = lowest_checked;
    let mut higher = lowest_checked + segment_size;

    let early_primes = early_primes.read().unwrap();

    while lower < highest_checked {
        if higher > highest_checked {
            higher = highest_checked;
        }
        
        let mut new_primes = vec![true; segment_size];

        for i in 0..early_primes.len() {
            let mut lolim = (lower / early_primes[i]) * early_primes[i];
            if lolim < lower {
                lolim += early_primes[i]
            }

            let mut j = lolim;
            while j < higher {
                new_primes[j - lower] = false;
                j += early_primes[i];
            }
        }
        let mut p = lower;
        while p < higher {
            if new_primes[p - lower] {
                returned_primes.push(p);
            }
            p += 1
        }
        
        lower += segment_size;
        higher += segment_size;
    }
    return returned_primes;
}

github link: https://github.com/knot427/Primes_To_N

Answer 1

Rust has great tooling. Use it!

Formatting

I'm mostly hoping to get some Rust and formatting tips

Formatting tip: cargo fmt will format your code in the idiomatic Rust style.

I don't necessarily always agree with the format... but it's good enough that I love automating that problem away.

Linting

In the same vein, cargo clippy will run a linter-on-steroids. In this case, it points out that using return x; as the last statement of a function is unnecessary, and you can just type x instead.

---

Let us start with some Rust style.

Learn your iterators

You are using indexes and arrays as if you were using C. This is not idiomatic, and may be costing you performance.

The loop while lower < highest_checked can be rewritten as:

for lower in (lowest_checked..highest_checked).step(segment_size) {
    let higher = cmp::min(lower + segment_size, highest_checked);

     // ...
}

The loop for i in 0..early_primes.len() can be rewritten as:

for early_prime in early_primes.iter().copied() {
    ...
}

Note: the copied ensures we have a local copy instead of a reference.

Similarly, the inner loop while j < higher can be rewritten as:

new_primes[(lolim - lower)..(higher - lower)]
    .iter_mut()
    .step(early_prime)
    .for_each(|is_prime| *is_prime = false);

Note: unlike your version, - lower is only applied at the beginning, rather than repeatedly, and bounds are not checked repeatedly.

And finally that while p < higher can be rewritten as:

new_primes
    .into_iter()
    .enumerate()
    .filter(|(is_prime, _)| *is_prime)
    .for_each(|(_, offset)| returned_primes.push(lower + offset as u64));

Using iterator code will avoid faffing about with indexes as much as possible:

• This avoids the risk of getting them wrong.
• This avoids the cost of bounds checks.

---

Let us focus on the core algorithm now.

Numeric Types

Your use of numeric types is problematic. You use usize when the number wouldn't fit in u32 (and 32-bits platforms are a thing), you convert usize to f64 which could lose precision, etc...

First of all, I'd argue for replacing all sieve numbers by u64, it's guaranteed to be large enough to hold the values you want regardless of the platform.

Secondly, your calculation of the square root of a u64 using a f64 is flawed for any u64 with at least 54 significant bits. You should either guard against that with an assert (53 significant bits is pretty large, already) or you should use the floating square as estimate and refine it.

Sieve on Stack

Your while lower < highest_checked loop will repeatedly create and throw away your sieve.

Because you limit the sieve to 100 K elements, if you made that 100 K a constant you could use an array of 100 K elements on the stack, without issues.

Flip Sieve

At the moment, you initialize your sieve with 1s. Initializing with 0s may be faster, so you may want trying to flip the meaning of those booleans.

(No guarantee there)

Sieve Memory Trashing

Cutting down on the sieve size will also cut down on the sieve cache usage, but it won't change the fact that you are repeatedly (in while j < higher) looping over that memory from one end to the other.

At the moment, you limited the segment size to 100 KB, which is a good start, but the L1 data cache is only 32 KB. It may be that limiting the segment size to 16 KB (half of L1) would be more cache-friendly, speeding up this inner loop execution.

And while at it, I'd advise putting that early_primes[i] in a local variable¹, just in case the compiler doesn't figure out it doesn't have to re-read it from memory every single time.

Altogether, this will keep the inner loop with minimal memory access, and hopefully it will make it easier for the compiler to unroll that loop.

¹ Using Rust Iterators properly will do that.

---

Let us look at threading next.

Locking or not locking

Your use of RwLock appears superfluous:

• You create the vector of primes.
• You then run multiple threads in parallel which only ever read the data.
• After joining, you write to the data.

There's never, actually, concurrent attempts at reading and writing, since all writes only happen when a single thread accesses the data.

You can thus, instead:

1. Use Arc<Vec<u64>>, directly.
2. Use Arc::try_unwrap(early_primes).unwrap() to get back the vector at the end.

This will eschew locking altogether, and avoid the to_owned() call to clone the vector again at the end.

Thread Creation is a costly endeavor.

Creating and Joining a thread are NOT trivial operations. Really not.

You would need to measure the cost of creating and joining compared to the cost of actually running the sieve, but I would not be surprised to learn it's a significant overhead especially at the beginning when sqrt_of_limit is low and early_primes is small.

The typical answer to this is to use a thread-pool, or some message-passing, so that rather than spinning up a new thread per "job to do", you spin up a few threads, have each of them perform all the jobs that need doing, and only then unwind them all.

I've never used thread pool libraries, so no idea which are good ones, but it may be worth investigating.