primesieve

For low-level optimizations, primesieve::iterator provides access to the underlying 64-bit primes array and the generate_next_primes() method. This allows you to use SIMD instructions (like AVX512 or NEON) to process batches of primes.

Use it.generate_next_primes() to fill the internal buffer, then iterate over it.primes_ up to it.size_ using vector intrinsics.

#include <primesieve.hpp>
#include <immintrin.h>
#include <iostream>

int main()
{
  primesieve::iterator it;
  it.generate_next_primes();

  uint64_t limit = 10000000000;
  __m512i sums = _mm512_setzero_si512();

  while (it.primes_[it.size_ - 1] <= limit)
  {
    // Sum 64-bit primes using AVX512
    for (std::size_t i = 0; i < it.size_; i += 8) {
      __mmask8 mask = (i + 8 < it.size_) ? 0xff : 0xff >> (i + 8 - it.size_);
      __m512i primes = _mm512_maskz_loadu_epi64(mask, (__m512i*) &it.primes_[i]);
      sums = _mm512_add_epi64(sums, primes);
    }

    // Generate up to 2^10 new primes
    it.generate_next_primes();
  }

  // Sum the 8 partial sums
  uint64_t sum = _mm512_reduce_add_epi64(sums);

  // Process the remaining primes (at most 2^10)
  for (std::size_t i = 0; it.primes_[i] <= limit; i++)
    sum += it.primes_[i];

  std::cout << "Sum of the primes <= " << limit << ": " << sum << std::endl;

  return 0;
}

To achieve maximum performance when using libprimesieve, follow these best practices:

• Array vs. Iterator: If you need to iterate over the same set of primes multiple times, use primesieve_generate_primes() or primesieve_generate_n_primes() to store them in an array (if RAM permits) instead of using a primesieve_iterator.
• Direction of Iteration: primesieve_next_prime() is up to 2x faster and uses half the memory of primesieve_prev_prime(). Prefer rewriting algorithms to move forward through primes.
• Iterator Threading: The primesieve_iterator is single-threaded. To parallelize, subdivide your sieving range into chunks and process each chunk in its own thread using its own primesieve_iterator object.
• SIMD/Vectorization: You can access the underlying 64-bit primes array within a primesieve_iterator to apply SIMD instructions (like AVX512 or NEON) for low-level optimizations.
• Jump Speedups: Use primesieve_jump_to() with an optional stop_hint parameter. If the sieving distance is small (e.g., < sqrt(start)), providing a stop_hint can significantly speed up the process by limiting the buffer size.
• Avoid Repeated Initialization: Functions like primesieve_count_primes() and primesieve_nth_prime() have an $O(\sqrt{start})$ initialization overhead. Do not call them repeatedly in a loop for small sieving distances; use a primesieve_iterator instead to avoid recurring overhead.

The primesieve::iterator is single-threaded. To parallelize an algorithm, you must manually subdivide the sieving distance into chunks and process each chunk in a separate thread using its own primesieve::iterator object.

Pattern for parallelization:

1. Subdivide the total distance into equally sized chunks.
2. Assign each chunk to a thread.
3. Each thread creates a primesieve::iterator(start, stop) for its specific range.
4. Combine results (e.g., using OpenMP reduction).

#include <primesieve.hpp>
#include <iostream>
#include <omp.h>

int main()
{
  uint64_t sum = 0;
  uint64_t dist = 1e10;
  int threads = omp_get_max_threads();
  uint64_t thread_dist = (dist / threads) + 1;

  #pragma omp parallel for reduction(+: sum)
  for (int i = 0; i < threads; i++)
  {
    uint64_t start = i * thread_dist;
    uint64_t stop = std::min(start + thread_dist, dist + 1);
    primesieve::iterator it(start, stop);
    uint64_t prime = it.next_prime();

    // Sum primes inside [start, stop[
    for (; prime < stop; prime = it.next_prime())
      sum += prime;
  }

  std::cout << "Sum of the primes <= " << dist << ": " << sum << std::endl;

  return 0;
}

primesieve uses a segmented sieve of Eratosthenes combined with wheel factorization to generate prime numbers efficiently.

Core Algorithms

• Segmentation: Instead of sieving the entire interval $[2, n]$ at once, the interval is subdivided into equal-sized segments sieved consecutively. This reduces memory requirements from $O(n)$ to $O(\sqrt{n})$ and allows segment sizes to fit into CPU L1 or L2 caches for faster access.
• Wheel Factorization: This technique skips multiples of small primes. primesieve uses a modulo 210 wheel, which skips multiples of 2, 3, 5, and 7.
• Bucket Sieve: For primes $> 2^{32}$, primesieve employs Tomás Oliveira e Silva's cache-friendly bucket list algorithm. This stores sieving primes in lists of buckets associated with specific segments, allowing the use of segments smaller than $\sqrt{n}$ without losing efficiency.