# Sieve of Eratosthenes - segmented to increase speed and range

A lot of the sieve code floating around here and on the 'net is sloppy or broken; in fact, I've yet to see a single one that is correct, clean (in the sense of avoiding superfluous iterations etc.) and that can actually sieve up to the end of its operating range. Some state at least that they cannot sieve the full range, some state limits and crash anyway even for ranges below that, the rest simply crash or return wrong results. Even the 'canonical' examples - Achim Flammenkamp's code, Tomás Oliveira e Silva's and Kim Walisch's @ primesieve.org - range over the whole spectrum in that regard.

That's why I thought it would be a good idea to have a simple, robust example that is correct and flexible enough to serve as a basis for people to experiment with their own ideas. The code I'm proposing here (full code in zrbj_sx_sieve64_v1.cpp) demonstrates two different uses of segmentation with the Sieve of Eratosthenes.

The first use of segmentation is to avoid unnecessary work when the start of the range to be sieved is higher than the square root of the upper end of the range, like in SPOJ problem #2 PRIME1. Sieving all the way up to the upper end of the range (10^9 for PRIME1) can easily exceed the time limit of 5 seconds; by contrast, sieving only the potential factors up to the square root and then the actual target range is several thousand times as fast. This is pretty much the only way of sieving close to 2^64 without taking a long holiday... The second use of segmentation is for speeding up the sieving of bigger ranges but that is explained below, with the code for extending the factor sieve.

The code here uses a packed, odds-only bitmap where set bits signify 'composite', and it implements a 32-bit capable version that limits sieving ranges - and index variables - to at most 2^32 numbers in one go. This limits the sieved range to 2^31 bits per call but that is already more than enough since smaller segments work better even at the upper end of the range (L2/L3 cache).

unsigned char odd_composites32[UINT32_MAX / (2 * CHAR_BIT) + 1];  // the small factor sieve
uint32_t sieved_bits = 0;                                         // how far it's been initialised

void extend_factor_sieve_to_cover (uint32_t max_factor_bit);      // bit, not number!

void sieve32 (unsigned char *target_segment, uint64_t bit_offset, uint32_t bit_count)
{
assert( bit_count > 0 && bit_count <= UINT32_MAX / 2 + 1 );

uint32_t max_bit = bit_count - 1;
uint64_t max_num = 2 * (bit_offset + max_bit) + 1;
uint32_t max_factor_bit = (max_factor32(max_num) - 1) / 2;

if (target_segment != odd_composites32)
{
extend_factor_sieve_to_cover(max_factor_bit);
}

std::memset(target_segment, 0, std::size_t((max_bit + CHAR_BIT) / CHAR_BIT));

for (uint32_t i = 3u >> 1; i <= max_factor_bit; ++i)
{
if (bit(odd_composites32, i))  continue;

uint32_t n = (i << 1) + 1;   // the actual prime represented by bit i (< 2^32)

uint32_t stride = n;         // == (n * 2) / 2
uint64_t start = (uint64_t(n) * n) >> 1;
uint32_t k;

if (start >= bit_offset)
{
k = uint32_t(start - bit_offset);
}
else // start < offset
{
k = stride - (bit_offset - start - 1) % stride - 1;
}

while (k <= max_bit)
{
set_bit(target_segment, k);

if ((k += stride) < stride)  // k can wrap since strides go up to almost 2^32
{
break;
}
}
}
}


The second use of segmentation is to achieve locality of access, because sieving in small segments that fit into the fast L1 cache of the CPU can speed things up significantly. At least until ranges get so high that the iteration over all potential odd factors - 203,280,220 of them, when working close to 2^64 - starts dominating. When initialising the small factors sieve this is definitely not the case, since it has only 6541 small factor primes to deal with. The speed can be more than doubled by presieving and remembering offsets, but that would make the code a hundred lines longer. Two versions of these optimisations are demoed in zrbj_sx_sieve32_v2.cpp and zrbj_sx_sieve32_v4.cpp.

This function is equivalent to the single call sieve32(odd_composites32, 0, max_factor_bit) except for the speed. For sieving no higher than 2^32 (SPOJ) the single call is more than good enough, since the maximum factor for a 32-bit number cannot exceed 2^16 and the call takes less than a millisecond.

void extend_factor_sieve_to_cover (uint32_t max_factor_bit)
{
uint32_t const SIEVE_BITS = sizeof(odd_composites32) * CHAR_BIT;

assert( max_factor_bit < SIEVE_BITS );

if (max_factor_bit >= sieved_bits)
{
assert( sieved_bits % CHAR_BIT == 0 );

uint32_t bits_to_sieve = max_factor_bit - sieved_bits + 1;
uint32_t segment_size = 1u << 18;  // == 32K * CHAR_BIT == L1 cache size
uint32_t partial_bits = bits_to_sieve % segment_size;
uint32_t segments     = bits_to_sieve / segment_size + (partial_bits != 0);

if ((SIEVE_BITS - sieved_bits) % segment_size == 0)
{
partial_bits = 0;   // always sieve full segments as long as there is enough space
}

for ( ; segments--; )
{
uint32_t bits_this_round = segments == 0 && partial_bits ? partial_bits : segment_size;

sieve32(odd_composites32 + sieved_bits / CHAR_BIT, sieved_bits, bits_this_round);

sieved_bits += bits_this_round;
}
}
}


Idiosyncratic implementation choice #1: using byte arrays instead of std::vector or std::bitset<> makes the code simpler, cleaner, faster and more robust. Not to mention more flexible. A more sophisticated implementation might allocate bitmaps as uint32_t[] or uint64_t[] and let functions take void* pointers as parameters, so that functions can cast to uint8_t or whatever (up to the allocated type size) as is most suitable.

Idiosyncratic coding choice #1: as you can see, I calculate many values into separate variables with relatively long (but hopefully accurate) names. This gives the calculation of these values several times as much source code space as elsewhere, where the stuff might be crammed into big expressions and loop headers. But elsewhere these calculations are often not very accurate. Besides, I'm using this as some sort of 'live' (executable) documentation for myself of how the stuff works. Most, if not all, Eratosthenes coding errors actually occur with indexing-related math, so I thought this might be appropriate here as well.

Here are two helper function templates to make pretty much anything usable as a bitmap, the moral equivalent to bit-twiddling macros in C:

///////////////////////////////////////////////////////////////////////////////////////////////////
// Bit manipulation functions are sometimes faster using the machine word size, sometimes using
// unsigned char. Using templated functions makes it possible to mix and match. It also makes
// it easy to experiment with intrinsics, by the simple expedient of defining specialisations
// for certain types and implementing them via intrinsics.

template<typename word_t>
void set_bit (word_t *p, uintxx_t index)
{
enum {  BITS_PER_WORD = sizeof(word_t) * CHAR_BIT  };

// we can trust the compiler to use masking and shifting instead of division; we cannot do that
// ourselves without having the log2 which cannot easily be computed as a constexpr

p[index / BITS_PER_WORD] |= word_t(1) << (index % BITS_PER_WORD);
}

//-------------------------------------------------------------------------------------------------
// Returning the result as unsigned char makes this compatible with the intrinsics for equivalent
// CPU instructions; forcing the intrinsics to bool instead would cause potential inefficiencies.

template<typename word_t>
unsigned char bit (word_t const *p, uintxx_t index)
{
enum {  BITS_PER_WORD = sizeof(word_t) * CHAR_BIT  };

return (unsigned char)((p[index / BITS_PER_WORD] >> (index % BITS_PER_WORD)) & 1u);
}


Here's part of the output for the test program, compiled with gcc 4.8.1 (MinGW64):

sizeof(uintxx_t) == 4
sizeof(void*)    == 8

initialising the small factor sieve... 4.255 s (962.6 M/s)
...
sieving some bigger ranges (10M)...
1 ..            100000000    115 ms 1658.6 M/s   5761455  OK
793877173 ..           1000000000    473 ms  831.2 M/s  10000000  OK
9769857367 ..          10000000000    634 ms  692.4 M/s  10000000  OK
99746761237 ..         100000000000    822 ms  587.6 M/s  10000000  OK
999723733787 ..        1000000000000   1037 ms  508.1 M/s  10000000  OK
9999700629011 ..       10000000000000   1254 ms  455.3 M/s  10000000  OK
99999677617663 ..      100000000000000   1504 ms  408.8 M/s  10000000  OK
999999654617569 ..     1000000000000000   1790 ms  368.0 M/s  10000000  OK
9999999631636541 ..    10000000000000000   2190 ms  320.8 M/s  10000000  OK
99999999608465399 ..   100000000000000000   2884 ms  258.9 M/s  10000000  OK
999999999585415333 ..  1000000000000000000   4260 ms  185.6 M/s  10000000  OK
9999999999562573471 .. 10000000000000000000   7213 ms  115.7 M/s  10000000  OK
18446744030316425227 .. 18446744030759878665   8734 ms   96.8 M/s  10000000  OK
18446744073265777349 .. 18446744073709551615   8793 ms   96.3 M/s  10000000  OK
sieving with maximum window size...
0 ..           4294967294  14405 ms  568.7 M/s 203280221  OK


Observe the speed difference between sieving the full 32-bit range in one go (last line) and sieving the same range in small segments when initialising the factor sieve.

The digests for verifying the sieve output have been computed from reliable sources, as indicated in the topic Checksumming large swathes of prime numbers? (for verification). Sources:

• up to 1,000,000,000,000 available as files from sites like primos.mat.br
• up to 2^64-10*2^32 in super-fast bulk via the primesieve.org console program (pipe)
• up to 2^64-1 - and beyond - via the gp/PARI program (pipe, about 1 million primes/minute)

Generally speaking, I have to admit that I don't understand everything the code is doing and that I could use a little more time in order to understand it. However, I hope I can still give you some hopefully useful pieces of advice:

• You could use a little more C++ in your code, which could be considered a bit too C-ish. What comes to my mind is BITS_PER_WORD which can actually be obtained with std::numeric_limits:

enum { BITS_PER_WORD = std::numeric_limits<word_t>::digits };


The same goes for several other constants such as UINT32_MAX which could be std::numeric_limits<uint32_t>::max(). I have to admit that it makes the code longer too, but it may be easier to change the integer type with a basic find-and-replace if you ever need to.

• You could also use the C++ algorithms instead of the C library. For example, replace std::memset by std::fill_n, which is likely to be optimized back into an std::memset anyway while being clearer:

std::fill_n(target_segment, std::size_t((max_bit + CHAR_BIT) / CHAR_BIT), 0);

• It is likely that it has already been optimized by the compiler but since you compute both index / BITS_PER_WORD and index % BITS_PER_WORD, you may as well make even more sure that these two values are computed together by using std::div. It also applies to bits_to_sieve % segment_size and bits_to_sieve / segment_size somewhere else in the program.

• If you have access to a C++11 compiler, don't hesitate to declare some variables constexpr, like SIEVE_BITS which can totally be computed at compile time.

Sorry that I'm four years late, but I'm very interested in this topic. I hope you are still interested in it as well.

WHEEL FACTORIZATION

That would greatly speed up your code. The canonical examples you reference use wheel factorization, so I'm quite sure you are aware of it. Perhaps you feel it's too large or complex to add here and might confuse people unfamiliar with it.

• Greg, I suggest you try it and come back with results. Good caching will massively improve speed, and obviously he is using a wheel { 2 } which gives most of the potential savings already. – gnasher729 Dec 17 '19 at 13:48
• Gnasher has hit the nail squarley on the head. A packed odds-only bitmap (a.k.a. mod-2 wheel) realises most of the potential of wheeled representation already, and most of the remaining potential can be realised via presieving and/or wheelish striding when reading out the numbers from the packed bitmap, all at virtually no cost. By contrast, the famous mod-30 wheel offers a further space reduction to a bit more than half but it makes the code hugely more complicated. If you were to drive it with lookup tables then you would lose speed, which means you have to use massively unrolled code at – DarthGizka Dec 24 '19 at 13:33
• some chosen phase for blasting through the in-phase whole-phase part of a segment, with setup and finisher code to get the operation into phase before the main loop and then to finish off odds and ends after the main loop. Huuugely complicated. – DarthGizka Dec 24 '19 at 13:35