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)
