Added execution policies. Increased performance for dynamically linke…

…d libraries. Fixed -Wextra errors. [CI SKIP]

include/boost/math/special_functions/prime_sieve.hpp

@@ -14,7 +14,11 @@
 #include <vector>
 #include <iterator>
 #include <cmath>
+
+#if __has_include(<execution>)
 #include <thread>
+#include <execution>
+#endif
 
 namespace boost { namespace math { namespace detail
 {
@@ -23,7 +27,7 @@ namespace boost { namespace math { namespace detail
 template<class Z, class Container>
 void linear_sieve(Z upper_bound, Container &c)
 {
-    Z least_divisors_size{upper_bound + 1};
+    size_t least_divisors_size{static_cast<size_t>(upper_bound + 1)};
     Z *least_divisors{new Z[least_divisors_size]{0}};
 
     for (Z i{2}; i <= upper_bound; ++i)
@@ -123,44 +127,57 @@ void mask_sieve(Z lower_bound, Z upper_bound, Container &c)
 }
 } // End namespace detail
 
-template<typename Z, class OutputIterator>
-auto prime_sieve(Z upper_bound, OutputIterator output) -> decltype(output)
+#if __has_include(<execution>)
+template<class ExecutionPolicy, typename Z, class OutputIterator>
+auto prime_sieve(ExecutionPolicy&& policy, Z upper_bound, OutputIterator output) -> decltype(output)
 {
     static_assert(std::is_integral<Z>::value, "No primes for floating point types");
     BOOST_ASSERT_MSG(upper_bound + 1 < std::numeric_limits<Z>::max(), "Type Overflow");
 
     std::vector<Z> primes;
     primes.reserve(upper_bound / std::log(upper_bound));
 
-    if (upper_bound <= 104729)
+    if (upper_bound <= 8192)
     {
-        boost::math::detail::prime_table(upper_bound, primes);
+        boost::math::detail::linear_sieve(upper_bound, primes);
     }
 
     else
     {
-        std::vector<Z> small_primes;
-        small_primes.reserve(1000);
-
-        // Spilt into two vectors and merge after joined to avoid data races
-        std::thread t1([upper_bound, &small_primes]{boost::math::detail::prime_table(static_cast<Z>(104729), small_primes);});
-        std::thread t2([upper_bound, &primes]{boost::math::detail::mask_sieve(static_cast<Z>(104729), upper_bound, primes);});
+        if constexpr (std::is_same_v<decltype(policy), std::execution::sequenced_policy>)
+        {
+            boost::math::detail::mask_sieve(static_cast<Z>(2), upper_bound, primes);
+        }
 
-        t1.join();
-        t2.join();
-        primes.insert(primes.begin(), small_primes.begin(), small_primes.end());
+        else
+        {
+            std::vector<Z> small_primes;
+            small_primes.reserve(1000);
+
+            // Split into two vectors and merge after joined to avoid data races
+            std::thread t1([upper_bound, &small_primes] {
+                boost::math::detail::prime_table(static_cast<Z>(8192), small_primes);
+            });
+            std::thread t2([upper_bound, &primes] {
+                boost::math::detail::mask_sieve(static_cast<Z>(8192), upper_bound, primes);
+            });
+
+            t1.join();
+            t2.join();
+            primes.insert(primes.begin(), small_primes.begin(), small_primes.end());
+        }
     }
 
     return std::move(primes.begin(), primes.end(), output);
 }
 
-template<class Z, class OutputIterator>
-auto prime_range(Z lower_bound, Z upper_bound, OutputIterator output) -> decltype(output)
+template<class ExecutionPolicy, class Z, class OutputIterator>
+auto prime_range(ExecutionPolicy&& policy, Z lower_bound, Z upper_bound, OutputIterator output) -> decltype(output)
 {
     std::vector<Z> primes;
     primes.reserve(upper_bound / std::log(static_cast<double>(upper_bound)));
 
-    boost::math::prime_sieve(upper_bound, std::back_inserter(primes));
+    boost::math::prime_sieve(policy, upper_bound, std::back_inserter(primes));
 
     auto it{primes.begin()};
     while(*it < lower_bound && it != primes.end())
@@ -170,11 +187,45 @@ auto prime_range(Z lower_bound, Z upper_bound, OutputIterator output) -> decltyp
 
     return std::move(it, primes.end(), output);
 }
+#endif //__has_include(<execution>)
+
+template<typename Z, class OutputIterator>
+auto prime_sieve(Z upper_bound, OutputIterator output) -> decltype(output)
+{
+    static_assert(std::is_integral<Z>::value, "No primes for floating point types");
+    BOOST_ASSERT_MSG(upper_bound + 1 < std::numeric_limits<Z>::max(), "Type Overflow");
+
+    std::vector<Z> primes;
+    primes.reserve(upper_bound / std::log(upper_bound));
+
+    if (upper_bound <= 8192)
+    {
+        boost::math::detail::linear_sieve(upper_bound, primes);
+    }
+
+    else
+    {
+        boost::math::detail::mask_sieve(static_cast<Z>(2), upper_bound, primes);
+    }
+
+    return std::move(primes.begin(), primes.end(), output);
+}
 
 template<class Z, class OutputIterator>
-inline auto prime_range(Z upper_bound, OutputIterator output) -> decltype(output)
+auto prime_range(Z lower_bound, Z upper_bound, OutputIterator output) -> decltype(output)
 {
-    return prime_range(static_cast<Z>(2), upper_bound, output);
+    std::vector<Z> primes;
+    primes.reserve(upper_bound / std::log(static_cast<double>(upper_bound)));
+
+    boost::math::prime_sieve(upper_bound, std::back_inserter(primes));
+
+    auto it{primes.begin()};
+    while(*it < lower_bound && it != primes.end())
+    {
+        ++it;
+    }
+
+    return std::move(it, primes.end(), output);
 }
 }}
 

reporting/performance/prime_sieve_performance.cpp

@@ -15,7 +15,7 @@ void prime_sieve(benchmark::State& state)
     for(auto _ : state)
     {
         std::vector<Z> primes;
-        benchmark::DoNotOptimize(boost::math::prime_sieve(upper, std::back_inserter(primes)));
+        benchmark::DoNotOptimize(boost::math::prime_sieve(std::execution::par, upper, std::back_inserter(primes)));
     }
     state.SetComplexityN(state.range(0));
 }
@@ -28,7 +28,7 @@ void prime_sieve_partial_range(benchmark::State& state)
     for(auto _ : state)
     {
         std::vector<Z> primes;
-        benchmark::DoNotOptimize(boost::math::prime_range(lower, upper, std::back_inserter(primes)));
+        benchmark::DoNotOptimize(boost::math::prime_range(std::execution::par, lower, upper, std::back_inserter(primes)));
     }
     state.SetComplexityN(state.range(0));
 }

test/test_prime_sieve.cpp

@@ -66,17 +66,17 @@ void test_prime_range()
 
     // Does the upper bound call work
     primes.clear();
-    boost::math::prime_range(1000, std::back_inserter(primes));
+    boost::math::prime_range(2, 1000, std::back_inserter(primes));
     BOOST_TEST_EQ(primes.size(), ref);
 
     // Does it work with a deque?
     std::deque<Z> d_primes;
-    boost::math::prime_range(1000, std::back_inserter(d_primes));
+    boost::math::prime_range(2, 1000, std::back_inserter(d_primes));
     BOOST_TEST_EQ(d_primes.size(), ref);
 
     // Does it work with a list?
     std::list<Z> l_primes;
-    boost::math::prime_range(1000, std::front_inserter(l_primes));
+    boost::math::prime_range(2, 1000, std::front_inserter(l_primes));
     BOOST_TEST_EQ(l_primes.size(), ref);
 
     // Does the lower bound change the results?
@@ -93,7 +93,7 @@ void test_prime_range()
     // Will it call the sieve for large input
     ref = 78498; // Calculated with wolfram-alpha
     primes.clear();
-    boost::math::prime_range(1000000, std::back_inserter(primes));
+    boost::math::prime_range(2, 1000000, std::back_inserter(primes));
     BOOST_TEST_EQ(primes.size(), ref);
 }
 
@@ -107,6 +107,48 @@ void test_prime_sieve_overflow()
                              std::back_inserter(primes));
 }
 
+template<typename Z>
+void test_par_prime_sieve()
+{
+    std::vector<Z> primes;
+    Z ref {168}; // Calculated with wolfram-alpha
+
+    // Does the function work with a vector
+    boost::math::prime_sieve(std::execution::par, 1000, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), ref);
+
+    // Tests for correctness
+    // 100
+    primes.clear();
+    boost::math::prime_sieve(std::execution::par, 100, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), 25);
+
+    // 10'000
+    primes.clear();
+    boost::math::prime_sieve(std::execution::par, 10000, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), 1229);
+
+    // 100'000
+    primes.clear();
+    boost::math::prime_sieve(std::execution::par, 100000, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), 9592);
+
+    // 1'000'000
+    primes.clear();
+    boost::math::prime_sieve(std::execution::par, 1000000, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), 78498);
+
+    // Does the function work with a list?
+    std::list<Z> l_primes;
+    boost::math::prime_sieve(std::execution::par, 1000, std::back_inserter(l_primes));
+    BOOST_TEST_EQ(l_primes.size(), ref);
+
+    // Does the function work with a deque?
+    std::deque<Z> d_primes;
+    boost::math::prime_sieve(std::execution::par, 1000, std::back_inserter(d_primes));
+    BOOST_TEST_EQ(d_primes.size(), ref);
+}
+
 int main()
 {
 
@@ -124,5 +166,10 @@ int main()
 
     test_prime_sieve<boost::multiprecision::cpp_int>();
 
+    test_par_prime_sieve<int>();
+    test_par_prime_sieve<int32_t>();
+    test_par_prime_sieve<int64_t>();
+    test_par_prime_sieve<uint32_t>();
+
     boost::report_errors();
 }