From cf98424536e2eca1c54455f042dcb0112420430a Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Tue, 8 Sep 2020 20:58:36 +1000 Subject: [PATCH 1/2] Make interface of linear_sieve based on output iterator. I broke something in the logic, so it segfaults in the benchmark. --- .../math/special_functions/prime_sieve.hpp | 75 +++++++++++++------ .../performance/prime_sieve_performance.cpp | 14 +--- 2 files changed, 55 insertions(+), 34 deletions(-) diff --git a/include/boost/math/special_functions/prime_sieve.hpp b/include/boost/math/special_functions/prime_sieve.hpp index 5e50c32fcb..0f5895974f 100644 --- a/include/boost/math/special_functions/prime_sieve.hpp +++ b/include/boost/math/special_functions/prime_sieve.hpp @@ -21,30 +21,59 @@ #include #include -namespace boost::math { namespace detail + +namespace boost::math { + +template +Integer upper_bound_prime_count(Integer x) +{ + using std::floor; + using std::log; + constexpr + auto c = 30 * log(113) / 113; // Magic numbers from Wikipedia. + return floor(c * x / log(x)); +} + +namespace detail { // https://mathworld.wolfram.com/SieveofEratosthenes.html // https://www.cs.utexas.edu/users/misra/scannedPdf.dir/linearSieve.pdf -template -void linear_sieve(Integer upper_bound, Container &resultant_primes) +template +OutputIterator linear_sieve(Integer upper_bound, OutputIterator resultant_primes) { - std::size_t least_divisors_size{static_cast(upper_bound + 1)}; + auto const first = resultant_primes; + auto const least_divisors_size = upper_bound + 1; std::unique_ptr least_divisors{new Integer[least_divisors_size]{0}}; - for (std::size_t i{2}; i < upper_bound; ++i) + for (Integer i{2}; i < upper_bound; ++i) { if (least_divisors[i] == 0) { least_divisors[i] = i; - resultant_primes.emplace_back(i); + *resultant_primes++ = i; } - for (std::size_t j{}; j < resultant_primes.size() && i * resultant_primes[j] <= upper_bound && - resultant_primes[j] <= least_divisors[i] && j < least_divisors_size; ++j) + for (Integer j = 0; j < resultant_primes - first + && i * resultant_primes[j] <= upper_bound + && resultant_primes[j] <= least_divisors[i] + && j < least_divisors_size; ++j) { - least_divisors[i * static_cast(resultant_primes[j])] = resultant_primes[j]; + least_divisors[i * resultant_primes[j]] = resultant_primes[j]; } } + + return resultant_primes; +} + +// This wrapper function could possibly drop the _container suffix with the +// judicious use of SFINAE. +template +void linear_sieve_container(Integer upper_bound, Container &resultant_primes) +{ + resultant_primes.resize(upper_bound_prime_count(upper_bound)); + auto const first = std::begin(resultant_primes); + auto const last = linear_sieve(upper_bound, first); + resultant_primes.resize(last - first); } // 4096 is where benchmarked performance of linear_sieve begins to diverge @@ -95,10 +124,9 @@ template void mask_sieve(Integer lower_bound, Integer upper_bound, Container &resultant_primes) { auto limit{std::floor(std::sqrt(static_cast(upper_bound))) + 1}; - std::vector primes {}; - primes.reserve(limit / std::log(limit)); + std::vector primes; - boost::math::detail::linear_sieve(limit, primes); + boost::math::detail::linear_sieve_container(limit, primes); boost::math::detail::mask_sieve(lower_bound, upper_bound, primes, resultant_primes); } @@ -187,12 +215,12 @@ void segmented_sieve(Integer lower_bound, Integer upper_bound, Container &result // Prepare for max value so you do not have to calculate this again if(limit < linear_sieve_limit) { - boost::math::detail::linear_sieve(static_cast(limit), primes); + boost::math::detail::linear_sieve_container(static_cast(limit), primes); } else { - boost::math::detail::linear_sieve(linear_sieve_limit, primes); + boost::math::detail::linear_sieve_container(linear_sieve_limit, primes); boost::math::detail::segmented_sieve(linear_sieve_limit, limit, primes, primes); } @@ -251,22 +279,21 @@ void prime_sieve(ExecutionPolicy&& policy, Integer upper_bound, Container &prime if(upper_bound <= linear_sieve_limit) { - boost::math::detail::linear_sieve(static_cast(upper_bound), primes); + boost::math::detail::linear_sieve_container(static_cast(upper_bound), primes); } else if(typeid(policy) == typeid(std::execution::seq)) { - boost::math::detail::linear_sieve(linear_sieve_limit, primes); + boost::math::detail::linear_sieve_container(linear_sieve_limit, primes); boost::math::detail::sequential_segmented_sieve(linear_sieve_limit, upper_bound, primes); } else { - std::vector small_primes {}; - small_primes.reserve(1028); + std::vector small_primes; std::thread t1([&small_primes] { - boost::math::detail::linear_sieve(static_cast(linear_sieve_limit * 2), small_primes); + boost::math::detail::linear_sieve_container(static_cast(linear_sieve_limit * 2), small_primes); }); std::thread t2([upper_bound, &primes] { boost::math::detail::segmented_sieve(static_cast(linear_sieve_limit * 2), upper_bound, primes); @@ -298,14 +325,14 @@ void prime_range(ExecutionPolicy&& policy, Integer lower_bound, Integer upper_bo if(upper_bound <= linear_sieve_limit) { - boost::math::detail::linear_sieve(static_cast(upper_bound), primes); + boost::math::detail::linear_sieve_container(static_cast(upper_bound), primes); } else if(typeid(policy) == typeid(std::execution::seq)) { if(limit <= linear_sieve_limit) { - boost::math::detail::linear_sieve(limit, primes); + boost::math::detail::linear_sieve_container(limit, primes); if(lower_bound <= limit) { @@ -320,7 +347,7 @@ void prime_range(ExecutionPolicy&& policy, Integer lower_bound, Integer upper_bo else { - boost::math::detail::linear_sieve(linear_sieve_limit, primes); + boost::math::detail::linear_sieve_container(linear_sieve_limit, primes); boost::math::detail::sequential_segmented_sieve(linear_sieve_limit, limit, primes); boost::math::detail::sequential_segmented_sieve(lower_bound, upper_bound, primes); } @@ -335,7 +362,7 @@ void prime_range(ExecutionPolicy&& policy, Integer lower_bound, Integer upper_bo small_primes.reserve(1028); std::thread t1([limit, &small_primes] { - boost::math::detail::linear_sieve(limit, small_primes); + boost::math::detail::linear_sieve_container(limit, small_primes); }); std::thread t2([lower_bound, limit, upper_bound, &primes] { @@ -361,7 +388,7 @@ void prime_range(ExecutionPolicy&& policy, Integer lower_bound, Integer upper_bo boost::math::prime_reserve(limit, small_primes); std::thread t1([&small_primes] { - boost::math::detail::linear_sieve(static_cast(linear_sieve_limit * 2), small_primes); + boost::math::detail::linear_sieve_container(static_cast(linear_sieve_limit * 2), small_primes); }); std::thread t2([limit, &primes] { diff --git a/reporting/performance/prime_sieve_performance.cpp b/reporting/performance/prime_sieve_performance.cpp index a806e6efb9..4ba5e218d2 100644 --- a/reporting/performance/prime_sieve_performance.cpp +++ b/reporting/performance/prime_sieve_performance.cpp @@ -74,22 +74,16 @@ void interval_sieve(benchmark::State& state) state.SetComplexityN(state.range(0)); } -// Complete Implementations -template -inline auto prime_sieve_helper(ExecuitionPolicy policy, Integer upper, Container primes) -{ - boost::math::prime_sieve(policy, upper, primes); - return primes; -} - template void prime_sieve(benchmark::State& state) { Integer upper = static_cast(state.range(0)); + std::vector primes; + for(auto _ : state) { - std::vector primes; - benchmark::DoNotOptimize(prime_sieve_helper(std::execution::par, upper, primes)); + primes.clear(); + boost::math::prime_sieve(std::execution::par, upper, primes); } state.SetComplexityN(state.range(0)); } From 29cf8f4769a2293c75f78531a6725418963a1d01 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Wed, 9 Sep 2020 10:41:28 +1000 Subject: [PATCH 2/2] Simplify logic of linear_sieve slightly. Just replace the use of resultant_primes - first with i - 1. --- include/boost/math/special_functions/prime_sieve.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/prime_sieve.hpp b/include/boost/math/special_functions/prime_sieve.hpp index 0f5895974f..e20622b53d 100644 --- a/include/boost/math/special_functions/prime_sieve.hpp +++ b/include/boost/math/special_functions/prime_sieve.hpp @@ -41,7 +41,6 @@ namespace detail template OutputIterator linear_sieve(Integer upper_bound, OutputIterator resultant_primes) { - auto const first = resultant_primes; auto const least_divisors_size = upper_bound + 1; std::unique_ptr least_divisors{new Integer[least_divisors_size]{0}}; @@ -53,7 +52,7 @@ OutputIterator linear_sieve(Integer upper_bound, OutputIterator resultant_primes *resultant_primes++ = i; } - for (Integer j = 0; j < resultant_primes - first + for (Integer j = 0; j < i - 1 && i * resultant_primes[j] <= upper_bound && resultant_primes[j] <= least_divisors[i] && j < least_divisors_size; ++j)