Program Listing for File euler_sieve_mult.hpp

↰ Return to documentation for file (include/primes/euler_sieve_mult.hpp)

#pragma once

#include <algorithm>
#include <concepts>
#include <forward_list>
#include <type_traits>
#include <vector>

#include "../../external/gcem/include/gcem.hpp"
#include "primes_interface.hpp"

using std::forward_list;
using std::integral;
using std::signed_integral;
using std::unsigned_integral;
using std::vector;

namespace madlib {

// Member declarations

class euler_sieve_mult : public primes_interface<euler_sieve_mult> {

  // Used to tell in constant time whether a given number is prime or not.
  vector<bool> odd_sieve;

  // Used to access calculated primes in constant time.
  vector<prime_type> ordered_primes;

public:
  PRIMES_CONSTEXPR euler_sieve_mult() = default;

  explicit PRIMES_CONSTEXPR
  euler_sieve_mult(const unsigned_integral auto &initially_through);

  PRIMES_CONSTEXPR euler_sieve_mult &
  calculate_primes_through_impl(const unsigned_integral auto &through_value);

  PRIMES_CONSTEXPR bool is_prime_impl(const unsigned_integral auto &number);

  PRIMES_CONSTEXPR prime_type
  get_nth_prime_impl(const unsigned_integral auto &n);

  PRIMES_CONSTEXPR size_t number_of_mapped_primes_impl();

  PRIMES_CONSTEXPR prime_type highest_odd_number_mapped_impl();

  PRIMES_CONSTEXPR void auto_expand_mapped_primes_impl();

  PRIMES_CONSTEXPR void map_next_prime_impl();
};

// Some day when containers receive static trait checking, add that here.

// member definitions

inline PRIMES_CONSTEXPR euler_sieve_mult::euler_sieve_mult(
    const unsigned_integral auto &initially_through) {
  if (initially_through == 0) {
    return;
  }

  calculate_primes_through_impl(initially_through);
}

inline PRIMES_CONSTEXPR euler_sieve_mult &
euler_sieve_mult::calculate_primes_through_impl(
    const unsigned_integral auto &through_value) {
  if (through_value < odd_sieve.size()) {
    return *this;
  }

  // Call this close enough for reserving space.  Let the OS take on the virtual
  // memory optimizations.
  ordered_primes.reserve(static_cast<size_t>(
      static_cast<double>(through_value) / gcem::log2(through_value)));

  // prime datastructures
  if (ordered_primes.empty()) {
    ordered_primes.push_back(2);
    ordered_primes.push_back(3);
    odd_sieve.push_back(false);
    odd_sieve.push_back(true);
    if (through_value < odd_sieve.size()) {
      return *this;
    }
  }

  const size_t ordered_primes_initial_size = ordered_primes.size();
  { // Keep scope of euler_sieve only to where it is used.

    // Memory optimization
    // if 20, then checks through 19.  if 19, then checks through 19. 19 maps to
    // 10 stored odds, with a max index of 9.
    const size_t forward_list_size =
        static_cast<size_t>(
            gcem::ceil(static_cast<double>(through_value) / 2.0)) -
        odd_sieve.size();
    forward_list<prime_type> euler_sieve(forward_list_size);

    // Euler sieve
    generate(euler_sieve.begin(), euler_sieve.end(),
             [current_odd = odd_sieve.size() * 2 - 1]() mutable {
               current_odd += 2;
               return current_odd;
             });

    for (size_t i = 1; i * i <= through_value; ++i) {
      const size_t a_prime = ordered_primes.at(i);
      erase_if(euler_sieve,
               [a_prime](size_t num) { return num % a_prime == 0; });
      while (!euler_sieve.empty() && euler_sieve.front() < a_prime * a_prime) {
        ordered_primes.push_back(euler_sieve.front());
        euler_sieve.pop_front();
      }
      if (euler_sieve.empty()) {
        break;
      }
    }
  }

  // Adapt to bool vector for O(1) prime checking of precalculated values.
  odd_sieve.reserve(through_value * 2 + 1);

  for (auto new_prime_itr = ordered_primes.begin() +
                            static_cast<long>(ordered_primes_initial_size);
       new_prime_itr != ordered_primes.end(); ++new_prime_itr) {
    while (odd_sieve.size() * 2 + 1 < *new_prime_itr) {
      odd_sieve.push_back(false);
    }
    odd_sieve.push_back(true);
  }

  while (odd_sieve.size() * 2 + 1 <= through_value) {
    odd_sieve.push_back(false);
  }

  return *this;
}

inline PRIMES_CONSTEXPR bool
euler_sieve_mult::is_prime_impl(const unsigned_integral auto &number) {
  if (number == 2) {
    return true;
  }
  if (number % 2 == 0) {
    return false;
  }

  const size_t idx = (number - 1) / 2;
  if (idx >= odd_sieve.size()) {
    calculate_primes_through_impl(number);
  }
  return odd_sieve.at(idx);
}

inline PRIMES_CONSTEXPR prime_type
euler_sieve_mult::get_nth_prime_impl(const unsigned_integral auto &n) {
  while (n >= ordered_primes.size()) {
    auto_expand_mapped_primes_impl();
  }
  return ordered_primes.at(n);
}

inline PRIMES_CONSTEXPR void
euler_sieve_mult::auto_expand_mapped_primes_impl() {
  // Tuning factor for expanding primes known.  Probably overkill, but primes
  // exhibit complex, not-smooth behavior so what looks like over-engineering
  // is required.
  const size_t expansion_constant = 10;
  const size_t number_of_primes = number_of_mapped_primes_impl();
  for (size_t expansion_factor =
           (highest_odd_number_mapped_impl() + 1) * expansion_constant;
       number_of_primes == number_of_mapped_primes_impl();
       expansion_factor *= expansion_constant) {
    calculate_primes_through_impl(odd_sieve.size() * expansion_factor);
  }
}

inline PRIMES_CONSTEXPR void euler_sieve_mult::map_next_prime_impl() {
  const size_t number_of_primes = number_of_mapped_primes_impl();
  for (size_t map_through_number = highest_odd_number_mapped_impl() + 2;
       number_of_primes == ordered_primes.size(); map_through_number += 2) {
    calculate_primes_through_impl(map_through_number);
  }
}

inline PRIMES_CONSTEXPR size_t
euler_sieve_mult::number_of_mapped_primes_impl() {
  return ordered_primes.size() + 1; // 2 is hardcoded
}

inline PRIMES_CONSTEXPR prime_type
euler_sieve_mult::highest_odd_number_mapped_impl() {
  return odd_sieve.size() * 2 + 1;
}

} // namespace madlib