Problem
I have a project in which I implemented variants of the Sieve of Eratosthenes as well as benchmarking and profiling harnesses for these (which can be ran with ./run.sh benchmark and ./run.sh profile, respectively; the results can then be accessed from the "./target/criterion/" directory). I'd like some assistance in trying to figure out why the SegmentedWheelFactorisedCountedSieve variant seems to produce primes at a slower rate relative to the other "counted" sieves (CountedSieve, WheelFactorisedCountedSieve and SegmentedCountedSieve), even though it should theoretically be faster than all of them (since it uses all of their optimisations).
For context, this page has some explanations about what wheel factorisation is, and this page has some explanations about what a segmented sieve is.
Benchmarking Results
Here are some graphs showing the average time per iteration for each of CountedSieve, WheelFactorisedCountedSieve, SegmentedCountedSieve and SegmentedWheelFactorisedCountedSieve (all having calculated 1,000,000 primes, with 3 seed primes used by the wheel factorised variants and sample size of 100 each), respectively:
And this was the terminal's output:
$ ./run.sh benchmark
Finished release [optimized] target(s) in 0.22s
Running `target/release/examples/criterion --bench`
WARNING: HTML report generation will become a non-default optional feature in Criterion.rs 0.4.0.
This feature is being moved to cargo-criterion (https://github.com/bheisler/cargo-criterion) and will be optional in a future version of Criterion.rs. To silence this warning, either switch to cargo-criterion or enable the 'html_reports' feature in your Cargo.toml.
Gnuplot not found, using plotters backend
Benchmarking Sieves/Counted/Amount = 1000000: Warming up for 3.0000 s
Warning: Unable to complete 100 samples in 5.0s. You may wish to increase target time to 2050.5s, enable flat sampling, or reduce sample count to 10.
Sieves/Counted/Amount = 1000000
time: [393.87 ms 394.04 ms 394.28 ms]
Found 16 outliers among 100 measurements (16.00%)
1 (1.00%) low severe
1 (1.00%) low mild
3 (3.00%) high mild
11 (11.00%) high severe
Benchmarking Sieves/Segmented Counted/Amount = 1000000: Warming up for 3.0000 s
Warning: Unable to complete 100 samples in 5.0s. You may wish to increase target time to 1398.4s, enable flat sampling, or reduce sample count to 10.
Sieves/Segmented Counted/Amount = 1000000
time: [276.96 ms 277.28 ms 277.70 ms]
Found 14 outliers among 100 measurements (14.00%)
6 (6.00%) low mild
1 (1.00%) high mild
7 (7.00%) high severe
Benchmarking Sieves/Wheel Factorised Counted/Amount = 1000000, Seeds = 3: Warming up for 3.0000 s
Warning: Unable to complete 100 samples in 5.0s. You may wish to increase target time to 946.7s, enable flat sampling, or reduce sample count to 10.
Sieves/Wheel Factorised Counted/Amount = 1000000, Seeds = 3
time: [187.42 ms 187.44 ms 187.47 ms]
Found 12 outliers among 100 measurements (12.00%)
3 (3.00%) low severe
3 (3.00%) low mild
4 (4.00%) high mild
2 (2.00%) high severe
Benchmarking Sieves/Segmented Wheel Factorised Counted/Amount = 1000000, Seeds = 3: Warming up for 3.0000 s
Warning: Unable to complete 100 samples in 5.0s. You may wish to increase target time to 5399.2s, enable flat sampling, or reduce sample count to 10.
Sieves/Segmented Wheel Factorised Counted/Amount = 1000000, Seeds = 3
time: [1.0644 s 1.0656 s 1.0669 s]
Found 1 outliers among 100 measurements (1.00%)
1 (1.00%) high mild
Evidently, the segmented, wheel factorised variant took the longest time per run.
Profiling Results
Here's a flame graph generated from profiling SegmentedWheelFactorisedCountedSieve for 1 minute, having it attempt to calculate 1,000,000 primes with 3 seed primes for wheel factorisation:
Evidently, most of the time spent by SegmentedWheelFactorisedCountedSieve is on calls to WheelFactorisedSegment::new, but I can't tell if that's primarily because of the number of segments or of the performance of WheelFactorisedSegment::new.
I've been reviewing the code for a couple of days but I haven't been able to figure out what I'm doing wrong with the segmented, wheel factorised variant for it to run so much slower, or how I might improve its performance to at least be on par with the other variants. Some help diagnosing this would be greatly appreciated.
Code
To review the project, it's probably easiest just to clone from the repository, and perhaps run cargo doc --no-deps to generate and peruse the documentation (which can then be accessed from "./target/doc/sieve/index.html") for the most important parts. That said, I've included the most important parts of the code here (taken from "src/lib.rs") for convenience:
use std::{
cell::RefCell,
cmp::{max, min},
mem::swap,
rc::Rc,
};
const CACHE_BYTES: usize = 8_000;
fn prime_count_overestimate(a: usize) -> usize {
let ln_a = (a as f64).ln();
((a as f64 / ln_a) * (1.0 + 1.2762 / ln_a)).ceil() as usize
}
fn prime_count_underestimate(a: usize) -> usize {
let f_a = a as f64;
(f_a / f_a.ln()).floor() as usize
}
fn inverse_prime_count(a: usize) -> usize {
let f_a = a as f64;
let ln_a = f_a.ln();
(f_a * (ln_a + ln_a.ln() - 1.0 + 1.8 * ln_a.ln() / ln_a)).ceil() as usize + 16
}
pub struct SizedSieve {
list: Vec<bool>,
index: usize,
limit: usize,
count: usize,
}
impl SizedSieve {
pub fn new(size: usize) -> Self {
Self {
list: vec![true; size],
index: 2,
limit: (size as f64).sqrt().floor() as usize,
count: 0,
}
}
}
impl Iterator for SizedSieve {
type Item = usize;
fn next(&mut self) -> Option<Self::Item> {
while self.index <= self.limit {
if self.list[self.index] {
let prime = self.index;
self.count += 1;
self.list
.iter_mut()
.skip(self.index.pow(2))
.step_by(self.index)
.for_each(|multiple| {
*multiple = false;
});
self.index += 1;
return Some(prime);
}
self.index += 1;
}
while self.index < self.list.len() {
if self.list[self.index] {
let prime = self.index;
self.count += 1;
self.index += 1;
return Some(prime);
}
self.index += 1;
}
None
}
fn size_hint(&self) -> (usize, Option<usize>) {
if !self.list.is_empty() {
let (underestimate, overestimate) = (
prime_count_underestimate(self.list.len() - 1),
prime_count_overestimate(self.list.len() - 1),
);
match (underestimate < self.count, overestimate < self.count) {
(false, false) => (underestimate - self.count, Some(overestimate - self.count)),
(false, true) => (underestimate - self.count, None),
(true, false) => (0, Some(overestimate - self.count)),
_ => (0, None),
}
} else {
(0, None)
}
}
}
pub struct CountedSieve {
list: Vec<bool>,
index: usize,
limit: usize,
count: usize,
max_count: usize,
}
impl CountedSieve {
pub fn new(max_count: usize) -> Self {
let size = inverse_prime_count(max_count);
Self {
list: vec![true; size],
index: 2,
limit: (size as f64).sqrt().floor() as usize,
count: 0,
max_count,
}
}
}
impl Iterator for CountedSieve {
type Item = usize;
fn next(&mut self) -> Option<Self::Item> {
if self.count > self.max_count {
return None;
}
while self.index <= self.limit {
if self.list[self.index] {
let prime = self.index;
self.count += 1;
self.list
.iter_mut()
.skip(self.index.pow(2))
.step_by(self.index)
.for_each(|multiple| {
*multiple = false;
});
self.index += 1;
return Some(prime);
}
self.index += 1;
}
while self.index < self.list.len() {
if self.list[self.index] {
let prime = self.index;
self.count += 1;
self.index += 1;
return Some(prime);
}
self.index += 1;
}
None
}
fn size_hint(&self) -> (usize, Option<usize>) {
let size = 1 + self.max_count - self.count;
(size, Some(size))
}
}
impl ExactSizeIterator for CountedSieve {}
fn gcd(mut u: usize, mut v: usize) -> usize {
if u == 0 {
return v;
} else if v == 0 {
return u;
}
let k = {
let i = u.trailing_zeros();
let j = v.trailing_zeros();
u >>= i;
v >>= j;
min(i, j)
};
loop {
debug_assert!(u % 2 == 1, "u = {} is even", u);
debug_assert!(v % 2 == 1, "v = {} is even", v);
if u > v {
swap(&mut u, &mut v);
}
v -= u;
if v == 0 {
return u << k;
}
v >>= v.trailing_zeros();
}
}
struct Totatives {
number: usize,
current: usize,
}
impl Totatives {
fn new(number: usize) -> Self {
Self { number, current: 1 }
}
fn is_coprime(&self, number: usize) -> bool {
gcd(self.number, number).eq(&1)
}
}
impl Iterator for Totatives {
type Item = usize;
fn next(&mut self) -> Option<Self::Item> {
while self.current < self.number {
if self.is_coprime(self.current) {
let current = self.current;
self.current += 1;
return Some(current);
}
self.current += 1;
}
None
}
fn size_hint(&self) -> (usize, Option<usize>) {
(0, Some(self.number - self.current))
}
}
struct Wheel {
seeds: Rc<Vec<usize>>,
modulus: usize,
totatives: Vec<usize>,
inverse_totatives: Vec<Option<usize>>,
}
impl Wheel {
fn new(amount: usize) -> Self {
let seeds: Vec<usize> = CountedSieve::new(amount).collect();
let modulus = seeds.iter().product();
let totatives: Vec<usize> = Totatives::new(modulus).collect();
let inverse_totatives = {
let mut out: Vec<Option<usize>> = vec![None; *totatives.last().unwrap() + 1];
totatives.iter().enumerate().for_each(|(index, &number)| {
out[number] = Some(index);
});
out
};
Self {
seeds: Rc::new(seeds),
modulus,
inverse_totatives,
totatives,
}
}
fn inverse_totative(&self, index: usize) -> Option<usize> {
*self.inverse_totatives.get(index)?
}
fn to_number(&self, index: usize) -> usize {
self.modulus * (index / self.totatives.len()) + self.totatives[index % self.totatives.len()]
}
fn to_index(&self, number: usize) -> Option<usize> {
Some(
self.totatives.len() * (number / self.modulus)
+ self.inverse_totative(number % self.modulus)?,
)
}
}
pub struct WheelFactorisedSizedSieve {
list: Vec<bool>,
index: usize,
limit: usize,
count: usize,
wheel: Wheel,
seed_index: usize,
end_number: usize,
}
impl WheelFactorisedSizedSieve {
pub fn new(size: usize, seed_amount: usize) -> Self {
let wheel = Wheel::new(seed_amount);
Self {
list: vec![
true;
(size as f64 * wheel.totatives.len() as f64 / wheel.modulus as f64).ceil()
as usize + seed_amount
+ 1
],
index: 1,
limit: (size as f64).sqrt().floor() as usize,
count: 0,
seed_index: 0,
end_number: size,
wheel,
}
}
}
impl Iterator for WheelFactorisedSizedSieve {
type Item = usize;
fn next(&mut self) -> Option<Self::Item> {
if self.seed_index < self.wheel.seeds.len() {
let prime = self.wheel.seeds[self.seed_index];
if prime < self.end_number {
self.count += 1;
self.seed_index += 1;
return Some(prime);
}
return None;
}
while self.wheel.to_number(self.index) <= self.limit {
if self.list[self.index] {
let prime = self.wheel.to_number(self.index);
if prime < self.end_number {
self.count += 1;
for factor_index in self.index..self.list.len() {
if let Some(multiple_index) = self
.wheel
.to_index(prime * self.wheel.to_number(factor_index))
{
if multiple_index < self.list.len() {
self.list[multiple_index] = false;
} else {
break;
}
}
}
self.index += 1;
return Some(prime);
}
return None;
}
self.index += 1;
}
while self.index < self.list.len() {
if self.list[self.index] {
let prime = self.wheel.to_number(self.index);
if prime < self.end_number {
self.count += 1;
self.index += 1;
return Some(prime);
}
return None;
}
self.index += 1;
}
None
}
fn size_hint(&self) -> (usize, Option<usize>) {
let (underestimate, overestimate) = (
prime_count_underestimate(self.end_number),
prime_count_overestimate(self.end_number),
);
match (underestimate < self.count, overestimate < self.count) {
(false, false) => (underestimate - self.count, Some(overestimate - self.count)),
(false, true) => (underestimate - self.count, None),
(true, false) => (0, Some(overestimate - self.count)),
_ => (0, None),
}
}
}
pub struct WheelFactorisedCountedSieve {
list: Vec<bool>,
index: usize,
limit: usize,
count: usize,
wheel: Wheel,
seed_index: usize,
max_count: usize,
}
impl WheelFactorisedCountedSieve {
pub fn new(max_count: usize, seed_amount: usize) -> Self {
let wheel = Wheel::new(seed_amount);
let size = inverse_prime_count(max_count);
Self {
list: vec![
true;
(size as f64 * wheel.totatives.len() as f64 / wheel.modulus as f64).floor()
as usize + seed_amount
+ 1
],
index: 1,
limit: (size as f64).sqrt().floor() as usize,
count: 0,
wheel,
seed_index: 0,
max_count,
}
}
}
impl Iterator for WheelFactorisedCountedSieve {
type Item = usize;
fn next(&mut self) -> Option<Self::Item> {
if self.count > self.max_count {
return None;
}
if self.seed_index < self.wheel.seeds.len() {
let prime = self.wheel.seeds[self.seed_index];
self.count += 1;
self.seed_index += 1;
return Some(prime);
}
while self.wheel.to_number(self.index) <= self.limit {
if self.list[self.index] {
let prime = self.wheel.to_number(self.index);
self.count += 1;
for factor_index in self.index..self.list.len() {
if let Some(multiple_index) = self
.wheel
.to_index(prime * self.wheel.to_number(factor_index))
{
if multiple_index < self.list.len() {
self.list[multiple_index] = false;
} else {
break;
}
}
}
self.index += 1;
return Some(prime);
}
self.index += 1;
}
while self.index < self.list.len() {
if self.list[self.index] {
let prime = self.wheel.to_number(self.index);
self.count += 1;
self.index += 1;
return Some(prime);
}
self.index += 1;
}
None
}
fn size_hint(&self) -> (usize, Option<usize>) {
let size = 1 + self.max_count - self.count;
(size, Some(size))
}
}
impl ExactSizeIterator for WheelFactorisedCountedSieve {}
fn smallest_multiple_at_least(factor: usize, minimum: usize) -> usize {
let remainder = (minimum + factor) % factor;
if remainder != 0 {
minimum + factor - remainder
} else {
minimum
}
}
struct Segment {
list: Vec<bool>,
index: usize,
start: usize,
}
impl Segment {
fn new(
size: usize,
start: usize,
primes: Rc<Vec<usize>>,
multiples: Rc<RefCell<Vec<usize>>>,
) -> Self {
let (mut list, end_sqrt) = (
vec![true; size],
((size + start) as f64).sqrt().floor() as usize,
);
primes
.iter()
.take_while(|&prime| *prime <= end_sqrt)
.enumerate()
.for_each(|(multiple_index, &prime)| {
while multiples.borrow_mut()[multiple_index] < start {
multiples.borrow_mut()[multiple_index] += prime;
}
let skip = multiples.borrow_mut()[multiple_index] - start;
list.iter_mut()
.skip(skip)
.step_by(prime)
.for_each(|multiple| {
*multiple = false;
multiples.borrow_mut()[multiple_index] += prime;
});
});
Self {
list,
index: 0,
start,
}
}
}
impl Iterator for Segment {
type Item = usize;
fn next(&mut self) -> Option<Self::Item> {
while self.index < self.list.len() {
if self.list[self.index] {
let prime = self.start + self.index;
self.index += 1;
return Some(prime);
}
self.index += 1;
}
None
}
fn size_hint(&self) -> (usize, Option<usize>) {
(
0,
Some(if !self.list.is_empty() {
let end_index = self.list.len() - 1;
if end_index >= self.index {
end_index - self.index
} else {
0
}
} else {
0
}),
)
}
}
pub struct SegmentedCountedSieve {
count: usize,
size: usize,
segment_size: usize,
segment: Segment,
segment_end_index: usize,
primes: Rc<Vec<usize>>,
prime_index: usize,
multiples: Rc<RefCell<Vec<usize>>>,
max_count: usize,
}
impl SegmentedCountedSieve {
pub fn new(max_count: usize) -> Self {
let size = inverse_prime_count(max_count);
let limit = (size as f64).sqrt().floor() as usize;
let primes: Rc<Vec<usize>> = Rc::new(SizedSieve::new(limit).collect());
let start_index = if let Some(prime) = primes.last() {
prime + 1
} else {
2
};
let segment_size = if size < CACHE_BYTES {
size
} else {
CACHE_BYTES
};
let segment_end_index = start_index + segment_size;
let multiples = Rc::new(RefCell::new(
primes
.iter()
.map(|&prime| smallest_multiple_at_least(prime, max(prime.pow(2), start_index)))
.collect::<Vec<usize>>(),
));
Self {
count: 0,
size,
segment_size,
segment: Segment::new(
segment_size,
start_index,
Rc::clone(&primes),
Rc::clone(&multiples),
),
segment_end_index,
primes,
prime_index: 0,
multiples,
max_count,
}
}
}
impl Iterator for SegmentedCountedSieve {
type Item = usize;
fn next(&mut self) -> Option<Self::Item> {
if self.count > self.max_count {
return None;
}
if self.prime_index < self.primes.len() {
let prime = self.primes[self.prime_index];
self.count += 1;
self.prime_index += 1;
return Some(prime);
}
if let Some(prime) = self.segment.next() {
self.count += 1;
return Some(prime);
}
let start_index = self.segment_end_index;
self.segment_end_index = start_index + self.segment_size;
if self.segment_end_index < self.size {
self.segment = Segment::new(
self.segment_size,
start_index,
Rc::clone(&self.primes),
Rc::clone(&self.multiples),
);
} else {
self.segment_end_index = self.size - 1;
self.segment = Segment::new(
self.segment_end_index - start_index,
start_index,
Rc::clone(&self.primes),
Rc::clone(&self.multiples),
);
}
if let Some(prime) = self.segment.next() {
self.count += 1;
return Some(prime);
}
None
}
fn size_hint(&self) -> (usize, Option<usize>) {
let size = 1 + self.max_count - self.count;
(size, Some(size))
}
}
impl ExactSizeIterator for SegmentedCountedSieve {}
struct WheelFactorisedSegment {
list: Vec<bool>,
index: usize,
start: usize,
end: usize,
wheel: Rc<Wheel>,
}
impl WheelFactorisedSegment {
fn new(
size: usize,
start: usize,
primes: Rc<Vec<usize>>,
multiples: Rc<RefCell<Vec<usize>>>,
wheel: Rc<Wheel>,
) -> Self {
let (mut list, start_number, end) = (
vec![true; size + wheel.seeds.len()],
wheel.to_number(start),
start + size,
);
let end_number = wheel.to_number(end);
let end_number_sqrt = (end_number as f64).sqrt().floor() as usize;
primes
.iter()
.take_while(|&prime| *prime <= end_number_sqrt)
.enumerate()
.for_each(|(prime_index, &prime)| {
for (totative_index, totative) in wheel.totatives.iter().enumerate() {
let multiple_index = totative_index + wheel.totatives.len() * prime_index;
if multiples.borrow_mut()[multiple_index] > end_number {
break;
}
let addend = prime * totative;
while multiples.borrow_mut()[multiple_index] < start_number {
multiples.borrow_mut()[multiple_index] += addend;
}
while multiples.borrow_mut()[multiple_index] <= end_number {
if let Some(index) = wheel.to_index(multiples.borrow_mut()[multiple_index])
{
list[index - start] = false;
}
multiples.borrow_mut()[multiple_index] += addend;
}
}
});
Self {
list,
index: 0,
start,
end,
wheel,
}
}
}
impl Iterator for WheelFactorisedSegment {
type Item = usize;
fn next(&mut self) -> Option<Self::Item> {
while self.index + self.start <= self.end {
if self.list[self.index] {
let prime = self.wheel.to_number(self.index + self.start);
self.index += 1;
return Some(prime);
}
self.index += 1;
}
None
}
fn size_hint(&self) -> (usize, Option<usize>) {
(
0,
Some(if !self.list.is_empty() {
let end_index = self.list.len() - 1;
if end_index >= self.index {
end_index - self.index
} else {
0
}
} else {
0
}),
)
}
}
pub struct SegmentedWheelFactorisedCountedSieve {
count: usize,
end_index: usize,
segment_size: usize,
segment: WheelFactorisedSegment,
segment_end_index: usize,
primes: Rc<Vec<usize>>,
prime_index: usize,
multiples: Rc<RefCell<Vec<usize>>>,
wheel: Rc<Wheel>,
max_count: usize,
}
impl SegmentedWheelFactorisedCountedSieve {
pub fn new(max_count: usize, seed_amount: usize) -> Self {
let wheel = Rc::new(Wheel::new(seed_amount));
let size = inverse_prime_count(max_count) as f64;
let limit = (size as f64).sqrt().floor() as usize;
let primes = {
let primes: Vec<usize> = WheelFactorisedSizedSieve::new(limit, seed_amount).collect();
if primes.len() > seed_amount + 1 {
Rc::new(primes)
} else {
Rc::clone(&wheel.seeds)
}
};
let start_index = {
let mut start_number = *primes.last().unwrap() + 1;
while wheel
.inverse_totative(start_number % wheel.modulus)
.is_none()
{
start_number += 1;
}
wheel.to_index(start_number).unwrap()
};
let size_factored =
(size as f64 * wheel.totatives.len() as f64 / wheel.modulus as f64).floor() as usize
+ seed_amount + 1;
let segment_size = if size_factored < CACHE_BYTES {
size_factored
} else {
CACHE_BYTES
};
let segment_end_index = start_index + segment_size;
let multiples = Rc::new(RefCell::new(
(0..primes.len() * wheel.totatives.len())
.map(|index| {
let (prime_index, totative_index) =
(index / wheel.totatives.len(), index % wheel.totatives.len());
let prime = primes[prime_index];
smallest_multiple_at_least(
prime * wheel.totatives[totative_index],
max(prime.pow(2), wheel.to_number(start_index)),
)
})
.collect::<Vec<usize>>(),
));
Self {
count: 0,
end_index: size_factored,
segment_size,
segment: WheelFactorisedSegment::new(
segment_size,
start_index,
Rc::clone(&primes),
Rc::clone(&multiples),
Rc::clone(&wheel),
),
segment_end_index,
multiples,
primes,
prime_index: 0,
wheel,
max_count,
}
}
}
impl Iterator for SegmentedWheelFactorisedCountedSieve {
type Item = usize;
fn next(&mut self) -> Option<Self::Item> {
if self.count > self.max_count {
return None;
}
if self.prime_index < self.primes.len() {
let prime = self.primes[self.prime_index];
self.count += 1;
self.prime_index += 1;
return Some(prime);
}
if let Some(prime) = self.segment.next() {
self.count += 1;
return Some(prime);
}
let start_index = self.segment_end_index + 1;
self.segment_end_index = start_index + self.segment_size;
if self.segment_end_index <= self.end_index {
self.segment = WheelFactorisedSegment::new(
self.segment_size,
start_index,
Rc::clone(&self.primes),
Rc::clone(&self.multiples),
Rc::clone(&self.wheel),
);
} else {
self.segment_end_index = self.end_index;
self.segment = WheelFactorisedSegment::new(
self.end_index - start_index,
start_index,
Rc::clone(&self.primes),
Rc::clone(&self.multiples),
Rc::clone(&self.wheel),
);
}
if let Some(prime) = self.segment.next() {
self.count += 1;
return Some(prime);
}
None
}
fn size_hint(&self) -> (usize, Option<usize>) {
let size = 1 + self.max_count - self.count;
(size, Some(size))
}
}
impl ExactSizeIterator for SegmentedWheelFactorisedCountedSieve {}
Remarks
Aside from that, I'd appreciate any other advice about diagnosing slow code or constructive feedback about improving code in general that you might be able to provide. Thanks.
P.S. This post was moved from Stack Overflow as it's apparently more suitable for this website instead. The original post is here.
SegmentedCountedSieveandSegmentedWheelFactorisedCountedSievehave improved, but the latter is still slower than all other variants, so this remains an open problem. \$\endgroup\$