Prime Number Functions

include/boost/math/special_functions/prime_sieve.hpp

@@ -0,0 +1,108 @@
+// Copyright 2020 Matt Borland
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_SPECIAL_FUNCTIONS_PRIME_SIEVE_HPP
+#define BOOST_MATH_SPECIAL_FUNCTIONS_PRIME_SIEVE_HPP
+
+#include <boost/math/special_functions/prime.hpp>
+#include <boost/assert.hpp>
+#include <deque>
+#include <vector>
+#include <iterator>
+#include <iostream>
+
+namespace boost { namespace math
+{
+
+// https://mathworld.wolfram.com/SieveofEratosthenes.html
+// https://www.cs.utexas.edu/users/misra/scannedPdf.dir/linearSieve.pdf
+template<class Z, class OutputIterator>
+auto prime_sieve(Z lower_bound, 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> least_divisors(upper_bound + 1, 0);
+    std::deque<Z> primes;
+
+    for (Z i{2}; i <= upper_bound; ++i)
+    {
+        if (least_divisors[i] == 0)
+        {
+            least_divisors[i] = i;
+            primes.emplace_back(i);
+        }
+
+        for (size_t j{}; j < least_divisors.size(); ++j)
+        {
+            if (j >= primes.size())
+            {
+                break;
+            }
+
+            else if (primes[j] > least_divisors[i])
+            {
+                break;
+            }
+
+            else if (i * primes[j] > upper_bound)
+            {
+                break;
+            }
+
+            else
+            {
+                least_divisors[i * primes[j]] = primes[j];
+            }
+        }
+    }
+
+    auto it{primes.begin()};
+    while (*it < lower_bound && it != primes.end())
+    {
+        primes.pop_front();
+        ++it;
+    }
+
+    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)
+{
+    if (upper_bound <= 104729)
+    {
+        Z i{2};
+        unsigned counter {};
+        std::deque<Z> primes;
+        while (i <= upper_bound)
+        {
+            if (i >= lower_bound)
+            {
+                primes.emplace_back(i);
+            }
+
+            ++counter;
+            i = static_cast<Z>(boost::math::prime(counter));
+        }
+
+        return std::move(primes.begin(), primes.end(), output);
+    }
+
+    else
+    {
+        return prime_sieve(lower_bound, upper_bound, output);
+    }
+}
+
+template<class Z, class OutputIterator>
+inline auto prime_range(Z upper_bound, OutputIterator output) -> decltype(output)
+{
+    return prime_range(static_cast<Z>(2), upper_bound, output);
+}
+}}
+
+#endif //BOOST_MATH_SPECIAL_FUNCTIONS_PRIME_SIEVE_HPP

reporting/performance/prime_sieve_performance.cpp

@@ -0,0 +1,66 @@
+// Copyright 2020 Matt Borland
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/math/special_functions/prime_sieve.hpp>
+#include <benchmark/benchmark.h>
+
+template <class Z>
+void prime_sieve(benchmark::State& state)
+{
+    Z upper = static_cast<Z>(state.range(0));
+    for(auto _ : state)
+    {
+        std::vector<Z> primes;
+        benchmark::DoNotOptimize(boost::math::prime_sieve(static_cast<Z>(2), upper, std::back_inserter(primes)));
+    }
+    state.SetComplexityN(state.range(0));
+}
+
+template <typename Z>
+void prime_range(benchmark::State& state)
+{
+    Z upper = static_cast<Z>(state.range(0));
+    for(auto _ : state)
+    {
+        std::vector<Z> primes;
+        benchmark::DoNotOptimize(boost::math::prime_range(static_cast<Z>(2), upper, std::back_inserter(primes)));
+    }
+    state.SetComplexityN(state.range(0));
+}
+
+template <class Z>
+void prime_sieve_partial_range(benchmark::State& state)
+{
+    Z upper = static_cast<Z>(state.range(0));
+    Z lower = static_cast<Z>(state.range(0)) > 2 ? static_cast<Z>(state.range(0)) : 2;
+    for(auto _ : state)
+    {
+        std::vector<Z> primes;
+        benchmark::DoNotOptimize(boost::math::prime_sieve(lower, upper, std::back_inserter(primes)));
+    }
+    state.SetComplexityN(state.range(0));
+}
+
+BENCHMARK_TEMPLATE(prime_sieve, int32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
+BENCHMARK_TEMPLATE(prime_sieve, int64_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
+BENCHMARK_TEMPLATE(prime_sieve, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
+BENCHMARK_TEMPLATE(prime_sieve_partial_range, int32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
+BENCHMARK_TEMPLATE(prime_sieve_partial_range, int64_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
+BENCHMARK_TEMPLATE(prime_sieve_partial_range, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
+BENCHMARK_TEMPLATE(prime_range, int32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
+BENCHMARK_TEMPLATE(prime_range, int64_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
+BENCHMARK_TEMPLATE(prime_range, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
+
+// Direct comparison of lookup vs sieve using only range of lookup
+BENCHMARK_TEMPLATE(prime_sieve, int32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
+BENCHMARK_TEMPLATE(prime_range, int32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
+BENCHMARK_TEMPLATE(prime_sieve, int64_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
+BENCHMARK_TEMPLATE(prime_range, int64_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
+BENCHMARK_TEMPLATE(prime_sieve, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
+BENCHMARK_TEMPLATE(prime_range, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
+
+BENCHMARK_MAIN();

test/test_prime_sieve.cpp

@@ -0,0 +1,127 @@
+// Copyright 2020 Matt Borland
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/math/special_functions/prime_sieve.hpp>
+#include <boost/core/lightweight_test.hpp>
+#include <boost/multiprecision/cpp_int.hpp>
+#include <list>
+#include <deque>
+#include <array>
+
+template<typename Z>
+void test_prime_sieve()
+{
+    std::vector<Z> primes;
+    Z ref {168}; // Calculated with wolfram-alpha
+
+    // Does the function work with a vector
+    boost::math::prime_sieve(2, 1000, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), ref);
+
+    // Tests for correctness
+    // 100
+    primes.clear();
+    boost::math::prime_sieve(2, 100, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), 25);
+
+    // 10'000
+    primes.clear();
+    boost::math::prime_sieve(2, 10000, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), 1229);
+
+    // 100'000
+    primes.clear();
+    boost::math::prime_sieve(2, 100000, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), 9592);
+
+    // 1'000'000
+    primes.clear();
+    boost::math::prime_sieve(2, 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(2, 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(2, 1000, std::back_inserter(d_primes));
+    BOOST_TEST_EQ(d_primes.size(), ref);
+}
+
+template<typename Z>
+void test_prime_range()
+{
+    std::vector<Z> primes;
+    Z ref {168}; // Calculated with wolfram-alpha
+
+    // Does the upper and lower bound call work
+    boost::math::prime_range(2, 1000, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), ref);
+
+    // Does the upper bound call work
+    primes.clear();
+    boost::math::prime_range(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_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_TEST_EQ(l_primes.size(), ref);
+
+    // Does the lower bound change the results?
+    ref = 143; // Calculated with wolfram-alpha
+    primes.clear();
+    boost::math::prime_range(100, 1000, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), ref);
+
+    // Does it work with 0 difference?
+    primes.clear();
+    boost::math::prime_range(2, 2, std::back_inserter(primes));
+    BOOST_TEST_EQ(primes.size(), 1);
+
+    // 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_TEST_EQ(primes.size(), ref);
+}
+
+template<typename Z>
+void test_prime_sieve_overflow()
+{
+    std::vector<Z> primes;
+
+    // Should die with call to BOOST_ASSERT
+    boost::math::prime_sieve(static_cast<Z>(2), static_cast<Z>(std::numeric_limits<Z>::max()),
+                             std::back_inserter(primes));
+}
+
+int main()
+{
+    test_prime_sieve<int>();
+    test_prime_sieve<int32_t>();
+    test_prime_sieve<int64_t>();
+    test_prime_sieve<uint32_t>();
+
+    test_prime_range<int>();
+    test_prime_range<int32_t>();
+    test_prime_range<int64_t>();
+    test_prime_range<uint32_t>();
+
+    test_prime_sieve<boost::multiprecision::cpp_int>();
+
+    //test_prime_sieve_overflow<int16_t>();
+
+    boost::report_errors();
+}