diff --git a/include/boost/math/special_functions/prime_sieve.hpp b/include/boost/math/special_functions/prime_sieve.hpp index 8ec48ca2e..5d415b34c 100644 --- a/include/boost/math/special_functions/prime_sieve.hpp +++ b/include/boost/math/special_functions/prime_sieve.hpp @@ -191,6 +191,39 @@ void segmented_sieve(Integer lower_bound, Integer upper_bound, Container &result boost::math::detail::segmented_sieve(lower_bound, upper_bound, primes, resultant_primes); } + +template +void sequential_segmented_sieve(Integer lower_bound, Integer upper_bound, Container &resultant_primes) +{ + const Integer L1_SIZE {32768}; + const Integer interval {L1_SIZE * 8}; + Integer current_lower_bound{lower_bound}; + Integer current_upper_bound{current_lower_bound + interval}; + + if(current_upper_bound > upper_bound) + { + current_upper_bound = upper_bound; + } + + std::size_t ranges {static_cast((upper_bound - lower_bound) / interval)}; + + boost::math::detail::IntervalSieve sieve(current_lower_bound, current_upper_bound, resultant_primes, resultant_primes); + if(ranges == 0) + { + return; + } + + for(std::size_t i {}; i < ranges; ++i) + { + current_lower_bound = current_upper_bound; + current_upper_bound += interval; + if(current_upper_bound > upper_bound) + { + current_upper_bound = upper_bound; + } + sieve.NewRange(current_lower_bound, current_upper_bound, resultant_primes); + } +} } // End namespace detail template @@ -215,9 +248,10 @@ void prime_sieve(ExecutionPolicy&& policy, Integer upper_bound, Container &prime boost::math::detail::linear_sieve(upper_bound, primes); } - else if (std::is_same_v) + else if(typeid(policy) == typeid(std::execution::seq)) { - boost::math::detail::segmented_sieve(static_cast(2), upper_bound, primes); + boost::math::detail::linear_sieve(static_cast(4096), primes); + boost::math::detail::sequential_segmented_sieve(static_cast(4096), upper_bound, primes); } else diff --git a/test/test_prime_sieve.cpp b/test/test_prime_sieve.cpp index bf47d4eff..8f0a5ee7d 100644 --- a/test/test_prime_sieve.cpp +++ b/test/test_prime_sieve.cpp @@ -13,6 +13,7 @@ #include #include #include +#include template void test_prime_sieve() @@ -24,11 +25,6 @@ void test_prime_sieve() boost::math::prime_sieve(1000, primes); BOOST_TEST_EQ(primes.size(), ref); - // Does the sequential policy work - primes.clear(); - boost::math::prime_sieve(std::execution::seq, 1000, primes); - BOOST_TEST_EQ(primes.size(), ref); - // Tests for correctness // 2 primes.clear(); @@ -68,6 +64,27 @@ void test_prime_sieve() BOOST_TEST_EQ(d_primes.size(), ref); } +template +void test_sequential_prime_sieve() +{ + std::vector primes; + + // 10'000 + primes.clear(); + boost::math::prime_sieve(std::execution::seq, 10000, primes); + BOOST_TEST_EQ(primes.size(), 1229); + + // 100'000 + primes.clear(); + boost::math::prime_sieve(std::execution::seq, 100000, primes); + BOOST_TEST_EQ(primes.size(), 9592); + + // 1'000'000 + primes.clear(); + boost::math::prime_sieve(std::execution::seq, 1000000, primes); + BOOST_TEST_EQ(primes.size(), 78498); +} + template void test_prime_range() { @@ -121,7 +138,8 @@ void test_par_prime_sieve_large() std::vector primes; Integer ref {54400028}; // Calculated with wolfram-alpha - // Force the sieve into the multi-threading section + // Force the sieve into the multi-threading section and test reserve functionality + boost::math::prime_reserve(static_cast(1073741824), primes); boost::math::prime_sieve(static_cast(1073741824), primes); BOOST_TEST_EQ(primes.size(), ref); } @@ -149,15 +167,15 @@ void test_linear_sieve() { std::vector primes; - boost::math::detail::linear_sieve(1'000, primes); + boost::math::detail::linear_sieve(static_cast(1'000), primes); BOOST_TEST_EQ(primes.size(), 168); primes.clear(); - boost::math::detail::linear_sieve(10'000, primes); + boost::math::detail::linear_sieve(static_cast(10'000), primes); BOOST_TEST_EQ(primes.size(), 1229); primes.clear(); - boost::math::detail::linear_sieve(100'000, primes); + boost::math::detail::linear_sieve(static_cast(100'000), primes); BOOST_TEST_EQ(primes.size(), 9592); } @@ -181,6 +199,12 @@ int main() //test_prime_sieve(); //test_prime_sieve(); + //test_sequential_prime_sieve(); + //test_sequential_prime_sieve(); + test_sequential_prime_sieve(); + //test_sequential_prime_sieve(); + //test_sequential_prime_sieve(); + //test_prime_range(); //test_prime_range(); //test_prime_range();