Synopsis

#include <boost/math/special_functions/prime_sieve.hpp>

namespace boost { namespace math {

namespace execution {
   struct cuda_policy {};
   inline constexpr cuda_policy cuda {};
}

enum class prime_range_strategy { automatic, full_sieve, test_survivors };

struct prime_sieve_options
{
   std::size_t l1d_bytes {0};          // 0 = detect
   std::size_t l2_bytes {0};           // 0 = detect (per core)
   std::size_t sieve_bytes {0};        // 0 = derive from sqrt(upper) and the caches; clamped to [16 KiB, 8 MiB]
   unsigned max_threads {0};           // 0 = std::thread::hardware_concurrency()
   std::size_t chunk_primes {1u << 19};  // primes per parallel chunk when storing output
   prime_range_strategy range_strategy {prime_range_strategy::automatic};
   bool probable_prime_only {false};   // beyond 2^64: Baillie-PSW only, skip the pseudosquares test
};

// primes in [2, upper_bound)
template <class ExecutionPolicy, class Integer, class OutputIterator>
OutputIterator prime_sieve(ExecutionPolicy&& policy, Integer upper_bound, OutputIterator out, const prime_sieve_options& options = {});
template <class Integer, class OutputIterator>
OutputIterator prime_sieve(Integer upper_bound, OutputIterator out, const prime_sieve_options& options = {});
template <class ExecutionPolicy, class Integer, class T, class Alloc>
void prime_sieve(ExecutionPolicy&& policy, Integer upper_bound, std::vector<T, Alloc>& out, const prime_sieve_options& options = {});
template <class Integer, class T, class Alloc>
void prime_sieve(Integer upper_bound, std::vector<T, Alloc>& out, const prime_sieve_options& options = {});

// primes in [lower_bound, upper_bound); the bounds may have different integer types
template <class ExecutionPolicy, class Lower, class Upper, class OutputIterator>
OutputIterator prime_range(ExecutionPolicy&& policy, Lower lower_bound, Upper upper_bound, OutputIterator out, const prime_sieve_options& options = {});
template <class Lower, class Upper, class OutputIterator>
OutputIterator prime_range(Lower lower_bound, Upper upper_bound, OutputIterator out, const prime_sieve_options& options = {});
template <class ExecutionPolicy, class Lower, class Upper, class T, class Alloc>
void prime_range(ExecutionPolicy&& policy, Lower lower_bound, Upper upper_bound, std::vector<T, Alloc>& out, const prime_sieve_options& options = {});
template <class Lower, class Upper, class T, class Alloc>
void prime_range(Lower lower_bound, Upper upper_bound, std::vector<T, Alloc>& out, const prime_sieve_options& options = {});

// number of primes in [2, upper_bound) and in [lower_bound, upper_bound)
template <class ExecutionPolicy, class Integer>
std::uint64_t prime_count(ExecutionPolicy&& policy, Integer upper_bound, const prime_sieve_options& options = {});
template <class ExecutionPolicy, class Lower, class Upper>
std::uint64_t prime_count(ExecutionPolicy&& policy, Lower lower_bound, Upper upper_bound, const prime_sieve_options& options = {});
template <class Integer>
std::uint64_t prime_count(Integer upper_bound, const prime_sieve_options& options = {});
template <class Lower, class Upper>
std::uint64_t prime_count(Lower lower_bound, Upper upper_bound, const prime_sieve_options& options = {});

// reserves room for all primes below upper_bound
template <class Integer, class T, class Alloc>
void prime_reserve(Integer upper_bound, std::vector<T, Alloc>& primes);

}} // namespaces

Description

These functions require C++17.

prime_sieve produces the primes in the half-open range [2, upper_bound), prime_range the primes in [lower_bound, upper_bound), and prime_count returns how many primes such a range contains without storing them. Output is always in ascending order. The overloads taking a std::vector append to it and are the fastest way to store the primes; the output iterator overloads return the advanced iterator. A lower_bound below 2 is treated as 2, negative bounds as 0, and an empty range produces nothing. A negative bound is clamped in its own type, so mixing a signed and an unsigned bound is safe.

Integer may be any builtin integer type or a Boost.Multiprecision integer such as cpp_int or mpz_int. Bounds that fit in 64 bits are sieved with a segmented sieve of Eratosthenes in the style of primesieve: one byte represents 30 integers (the eight residues coprime to 30), multiples of the primes up to 163 are removed by pre-computed periodic patterns, and the sieving primes are split into three classes (many multiples per segment, a few, or at most one) that use a mod 30 wheel with unrolled loops, lists sorted by wheel position, and Oliveira e Silva's bucket sieve (fixed 8 KiB buckets from a memory pool) respectively. Segments are sized to the L1 and L2 data caches, which are detected on Linux and macOS and can be overridden through prime_sieve_options; when big primes are present the largest segment (8 MiB) is used. The wheel tables and the 16 pre-sieve patterns (123 KB) are computed at compile time with constexpr on GCC and Clang; MSVC and nvcc, whose constexpr step limits are lower, generate the patterns on first use (define BOOST_MATH_PRIME_SIEVE_RUNTIME_PRESIEVE to force that everywhere).

Short intervals at large magnitudes (roughly when the interval is narrower than sqrt(upper_bound) / ln(upper_bound)) switch automatically to a window method: the interval is sieved by the primes up to a modest depth and the survivors are tested with a Baillie-PSW test in Montgomery arithmetic, which is deterministic below 2⁶⁴. prime_range_strategy forces either method.

Bounds beyond 2⁶⁴ use the same window method with the Lukes-Patterson-Williams pseudosquares test (Sorenson, The pseudosquares prime sieve, 2006), which is deterministic while upper_bound divided by the sieving depth stays below the largest tabulated pseudosquare (about 2.95 * 10²³), and Baillie-PSW beyond that, for which no counterexample is known. Setting probable_prime_only uses Baillie-PSW everywhere and is roughly ten times faster. This path is intended for windows of up to about 10⁸ integers; the head of a range that straddles 2⁶⁴ is still processed by the 64-bit engine.

Execution policies

The policy overloads accept the std::execution policies and boost::math::execution::cuda:

The std::execution overloads are available when the standard library provides the policy types (libstdc++ 9 and later, MSVC, and libc++ with -fexperimental-library); no parallel algorithms and therefore no TBB are required. Threads are never used when Boost.Math is configured without thread support.

Examples

#include <boost/math/special_functions/prime_sieve.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <execution>
#include <vector>
#include <list>
#include <iterator>
#include <cstdint>

// All primes below one billion, using every core
std::vector<std::uint64_t> primes;
boost::math::prime_reserve(1'000'000'000, primes);
boost::math::prime_sieve(std::execution::par, 1'000'000'000, primes);

// Count the primes in [10^18, 10^18 + 10^9) on one thread
std::uint64_t n = boost::math::prime_count(1'000'000'000'000'000'000ull, 1'000'000'001'000'000'000ull);

// Primes just above 2^64 with an arbitrary precision type
using boost::multiprecision::cpp_int;
std::vector<cpp_int> big;
boost::math::prime_range(cpp_int(1) << 64, (cpp_int(1) << 64) + 1'000'000, big);

// A range through any output iterator
std::list<int> small;
boost::math::prime_range(100, 1000, std::back_inserter(small));

Performance

Measured with the Google Benchmark harness in reporting/performance/prime_sieve_performance.cpp on an Apple M4 Max (12 performance and 4 efficiency cores, 128 KiB L1D), Apple clang 17, -O3 -march=native, against primesieve 12.16 built from source on the same machine. Times are wall clock in milliseconds.

The same comparison on an Intel Core i9-11900K (8 cores, 48 KiB L1D, 512 KiB L2 per core) with GCC 13, -O3 -march=native, primesieve 12.16 built from source there:

Single threaded this library is 10 to 15 percent ahead on Apple silicon and within 2 to 3 percent of primesieve on x86 for whole ranges, while primesieve keeps a lead of 20 to 30 percent on x86 for windows at 10¹⁵ and beyond (the bucket sieve is memory bound there and the Intel machine has a third of the memory bandwidth). With all cores this library is ahead on both machines. When compiled with AVX-512 VBMI2 support (-march=native on Ice Lake or Zen 4 and later) primes are extracted from the sieve with byte compression instructions; define BOOST_MATH_PRIME_SIEVE_NO_SIMD to keep the portable loop.

Both sieves are O(n log log n) in time and use O(sqrt(n)) memory for the sieving primes (about four bytes per prime here) plus the segment; the multithreaded path additionally buffers at most 2 * max_threads * chunk_primes output values.

The CUDA backend was measured on an NVIDIA GeForce RTX 3060 (CUDA 13.1, nvcc -O3 -arch=sm_86) against the same code running with std::execution::par on that machine's Intel Core i9-11900K (8 cores, 16 threads). It stores 16 bytes per sieving prime on the device, so ranges up to 2⁶⁴ need about 3.3 GB of device memory.

Counting is where the GPU pays off; storing primes is dominated by the transfer and compaction of the output (3.6 GB for the primes below 10¹⁰) and is faster on the CPU. Windows near 2⁶⁴ are dominated by generating the 203 million sieving primes on the host.

References