Program Listing for File euler_sieve_sqrt.hpp

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

#pragma once

#include <algorithm>
#include <execution>
#include <type_traits>
#include <unordered_set>
#include <vector>

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

using std::integral;
using std::signed_integral;
using std::unsigned_integral;
using std::vector;
using std::execution::par_unseq;

namespace madlib {

// Member declarations

class euler_sieve_sqrt : public primes_interface<euler_sieve_sqrt> {

  // Used to tell in constant time whether a given number is prime or not.
  std::unordered_set<prime_type> odd_sieve;
  size_t highest_number_calculated_through = 3;

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

public:
  PRIMES_CONSTEXPR euler_sieve_sqrt();

  explicit PRIMES_CONSTEXPR
  euler_sieve_sqrt(const unsigned_integral auto &initially_through);

  PRIMES_CONSTEXPR euler_sieve_sqrt &
  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_sqrt::euler_sieve_sqrt() {
  odd_sieve.max_load_factor(0.5);
}

inline PRIMES_CONSTEXPR euler_sieve_sqrt::euler_sieve_sqrt(
    const unsigned_integral auto &initially_through) {
  odd_sieve.max_load_factor(0.5);
  if (initially_through == 0) {
    return;
  }

  calculate_primes_through_impl(initially_through);
}

inline PRIMES_CONSTEXPR euler_sieve_sqrt &
euler_sieve_sqrt::calculate_primes_through_impl(
    const unsigned_integral auto &through_value) {
  if (through_value < highest_number_calculated_through) {
    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);
  }

  { // 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 - highest_number_calculated_through)));
    vector<prime_type> euler_sieve;
    euler_sieve.reserve(forward_list_size);

    // Euler sieve
    for (auto current_odd = highest_number_calculated_through + 2;
         current_odd <= through_value; current_odd += 2) {
      if (std::all_of(par_unseq, ordered_primes.begin() + 1,
                      ordered_primes.end(), [current_odd](auto known_prime) {
                        return current_odd % known_prime != 0;
                      })) {
        euler_sieve.push_back(current_odd);
      }
    }
    euler_sieve.shrink_to_fit();

    while (!euler_sieve.empty()) {
      {
        auto itr = euler_sieve.begin();
        for (; itr != euler_sieve.end() &&
               *itr < ordered_primes.back() * ordered_primes.back();
             ++itr) {
          ordered_primes.push_back(*itr);
        }
        euler_sieve.erase(euler_sieve.begin(), itr);
      }

      if (euler_sieve.empty()) {
        break;
      }

      auto assign_itr = euler_sieve.begin();
      for (auto maybe_prime : euler_sieve) {
        if (std::all_of(
                par_unseq, ordered_primes.begin() + 1, ordered_primes.end(),
                [maybe_prime](auto known_prime) { return maybe_prime % known_prime != 0; })) {
          *assign_itr = maybe_prime;
          ++assign_itr;
        }
      }
      euler_sieve.erase(assign_itr, euler_sieve.end());
    }
  }

  highest_number_calculated_through = through_value;

  return *this;
}

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

  if (number >= highest_number_calculated_through) {
    calculate_primes_through_impl(number);
  }

  if(odd_sieve.size() < ordered_primes.size() && ordered_primes[odd_sieve.size()] <= number){
    odd_sieve.reserve(ordered_primes.size());

    for (auto prime_itr = ordered_primes.begin() + odd_sieve.size();
         prime_itr != ordered_primes.end(); ++prime_itr) {
      odd_sieve.insert(*prime_itr);
    }
  }
  return odd_sieve.contains(number);
}

inline PRIMES_CONSTEXPR prime_type
euler_sieve_sqrt::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_sqrt::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(highest_number_calculated_through *
                                  expansion_factor);
  }
}

inline PRIMES_CONSTEXPR void euler_sieve_sqrt::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_sqrt::number_of_mapped_primes_impl() {
  return ordered_primes.size() + 1; // 2 is hardcoded
}

inline PRIMES_CONSTEXPR prime_type
euler_sieve_sqrt::highest_odd_number_mapped_impl() {
  return highest_number_calculated_through;
}

} // namespace madlib