Sieve of Eratosthenes

eratosthenes_sieve.cpp

// Sieve of Eratosthenes
// (c) 2021 Jin X ([email protected])

#include <vector>
#include <cmath>
#include <stdexcept>

// Sieve of Eratosthenes class
class EratosthenesSieve
{
public:
  // Initialize sieve
  EratosthenesSieve(size_t limit) : m_limit(limit) {
    size_t add = limit % 30 == 0 ? 0 : 1;
    m_sieve.resize((limit / 30 + add) * 8, true);  // init sieve with "prime" flags
    m_sieve[0] = false;  // 1 is not a prime number
    for (size_t i = 7, idx_i = 1, lim = static_cast<size_t>(std::sqrt(limit)); i <= lim; i += m_delta[idx_i], idx_i = (idx_i+1) & 7) {  // I use 'sqrt' instead of 'i*i' for overflow protection (and speed up)
      for (size_t n = i*i, idx_n = m_ofs[i % 30 >> 1]; n <= limit && n >= i; n += m_delta[idx_n] * i, idx_n = (idx_n+1) & 7) {  // n is always odd; 'n > i' is overflow protection (can be omitted if you're sure that limit can't approach SIZE_MAX)
        size_t ofs = m_ofs[n % 30 >> 1];
        if (ofs != 8) {  // non-multiple of 2, 3 and 5
          m_sieve[n / 30 * 8 + ofs] = false;
        } // if (ofs != 8)
      } // for (size_t n...)
    } // for (size_t i...)
  } // EratosthenesSieve (ctor)

  // Check value with exception for out-of-limit values
  bool is_prime(size_t n) {
    if (n > m_limit) { throw std::out_of_range("value exceeds the limit"); }
    return is_prime_safe(n);
  } // is_prime

  // Check value with no exceptions (returns false for out-of-limit value)
  bool is_prime_safe(size_t n) noexcept {
    if (n <= m_limit) {  // will be omitted by a good compiler with optimization options when 'is_prime_safe' is called from 'is_prime'
      if (n & 1) {  // odd
        size_t ofs = m_ofs[n % 30 >> 1];
        if (ofs != 8) {  // non-multiple of 3 and 5
          return m_sieve[n / 30 * 8 + ofs];  // result from sieve
        } // if (ofs != 8)
      } // if (n & 1)
      return (n <= 5 && n & 3);  // special cases (2, 3, 5: less of equal to 5 and non-multiple of 4, value 1 is in sieve)
    }
    return false;  // out-of-limit value
  } // is_prime_safe

private:
  static const uint8_t m_delta[8];  // increment array
  static const uint8_t m_ofs[15];   // offsets of odd values inside block of 30 values
  size_t m_limit;                   // value limit
  std::vector<bool> m_sieve;        // sieve (8 values of each block of 30)
}; // EratosthenesSieve

const uint8_t EratosthenesSieve::m_delta[8] = { 6, 4, 2, 4, 2, 4, 6, 2 };
const uint8_t EratosthenesSieve::m_ofs[15] = { 0, 8, 8, 1, 8, 2, 3, 8, 4, 5, 8, 6, 8, 8, 7 };

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <chrono>
#include <random>

using std::cout;
using std::endl;

int main()
{
  const size_t PRIME_LIMIT = 1'000'000'000;

  cout << "Initializing Sieve of Eratosthenes..." << endl;
  using Clock = typename std::conditional<std::chrono::high_resolution_clock::is_steady, std::chrono::high_resolution_clock, std::chrono::steady_clock>::type;
  auto start_time = Clock::now();
  EratosthenesSieve es(PRIME_LIMIT);  // init sieve
  auto stop_time = Clock::now();
  double duration = std::chrono::duration_cast<std::chrono::duration<double, std::ratio<1>>>(stop_time - start_time).count();
  cout << "Elapsed time: " << std::fixed << duration << " sec" << endl << endl;

  cout << "Primes from 0 to 1000..." << endl;
  for (size_t n = 0; n <= 1000; ++n) {
    if (es.is_prime(n)) { cout << n << ' '; }
  } // for
  cout << endl << endl;

  cout << "Checking 1000 random primes up to 1 billion..." << endl;
  std::random_device rnd;
  std::mt19937 rng(rnd());
  std::uniform_int_distribution<size_t> random;
  using pt = std::uniform_int_distribution<size_t>::param_type;
  for (int i = 0; i < 1000; ++i) {
    size_t n = random(rng, pt{0, random(rng, pt{0, random(rng, pt{0, random(rng, pt{0, PRIME_LIMIT}) }) }) });  // shift random numbers to smaller values
    if (es.is_prime(n)) { cout << n << ' '; }
  } // for
  cout << endl;
} // main