Skip to content

Commit

Permalink
Complete revamp of algorithm. Hide implementation behind detail names…
Browse files Browse the repository at this point in the history 
…pace.
  • Loading branch information
@mborland
mborland committed Jul 12, 2020
1 parent a684dbd commit 349ca6f
Show file tree
Hide file tree
Showing 3 changed files with 117 additions and 67 deletions.
140 changes: 107 additions & 33 deletions include/boost/math/special_functions/prime_sieve.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,88 +14,162 @@
#include <vector>
#include <iterator>
#include <iostream>
#include <cmath>
#include <thread>

namespace boost { namespace math
namespace boost { namespace math { namespace detail
{

// 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)
template<class Z, class Container>
void linear_sieve(Z upper_bound, Container &c)
{
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;
Z least_divisors_size{upper_bound + 1};
Z *least_divisors{new Z[least_divisors_size]{0}};

for (Z i{2}; i <= upper_bound; ++i)
{
if (least_divisors[i] == 0)
{
least_divisors[i] = i;
primes.emplace_back(i);
c.emplace_back(i);
}

for (size_t j{}; j < least_divisors.size(); ++j)
for (size_t j{}; j < least_divisors_size; ++j)
{
if (j >= primes.size())
if (j >= c.size())
{
break;
}

else if (primes[j] > least_divisors[i])
else if (c[j] > least_divisors[i])
{
break;
}

else if (i * primes[j] > upper_bound)
else if (i * c[j] > upper_bound)
{
break;
}

else
{
least_divisors[i * primes[j]] = primes[j];
least_divisors[i * c[j]] = c[j];
}
}
}

auto it{primes.begin()};
while (*it < lower_bound && it != primes.end())
delete[] least_divisors;
}

template<class Z, class Container>
void prime_table(Z upper_bound, Container &c)
{
Z i{2};
unsigned counter{};

while (i <= upper_bound && counter < 9999) // 10k elements are in the lookup table
{
primes.pop_front();
++it;
c.emplace_back(i);
++counter;
i = static_cast<Z>(boost::math::prime(counter));
}

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 Z, class Container>
void mask_sieve(Z lower_bound, Z upper_bound, Container &c)
{
if (upper_bound <= 104729)
Z limit{static_cast<Z>(std::floor(std::sqrt(upper_bound))) + 1};
std::vector<Z> primes;
primes.reserve(limit / std::log(limit));

boost::math::detail::linear_sieve(limit, primes);

const Z n{upper_bound - lower_bound + 1};
bool *mask{new bool[n + 1]{false}};

for (size_t i{}; i < primes.size(); ++i)
{
Z i{2};
unsigned counter {};
std::deque<Z> primes;
while (i <= upper_bound)
Z lower_limit = std::floor(lower_bound / primes[i]) * primes[i];

if (lower_limit < lower_bound)
{
lower_limit += primes[i];
}

if (lower_limit == primes[i])
{
lower_limit += primes[i];
}

for (Z j{lower_limit}; j <= upper_bound; j += primes[i])
{
mask[j - lower_bound] = true;
}
}

// Numbers which are not masked in range, are prime
for (Z i{lower_bound}; i <= upper_bound; i++)
{
if (!mask[i - lower_bound])
{
if (i >= lower_bound)
{
primes.emplace_back(i);
c.emplace_back(i);
}

++counter;
i = static_cast<Z>(boost::math::prime(counter));
}
}

delete[] mask;
}
} // End namespace detail

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");

return std::move(primes.begin(), primes.end(), output);
std::vector<Z> primes;
primes.reserve(upper_bound / std::log(upper_bound));

if (upper_bound <= 104729)
{
boost::math::detail::prime_table(upper_bound, primes);
}

else
{
return prime_sieve(lower_bound, upper_bound, output);
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);});

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)
{
std::vector<Z> primes;
primes.reserve(upper_bound / std::log(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);
}

template<class Z, class OutputIterator>
Expand Down
27 changes: 2 additions & 25 deletions reporting/performance/prime_sieve_performance.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,19 +15,7 @@ void prime_sieve(benchmark::State& state)
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)));
benchmark::DoNotOptimize(boost::math::prime_sieve(upper, std::back_inserter(primes)));
}
state.SetComplexityN(state.range(0));
}
Expand All @@ -40,7 +28,7 @@ void prime_sieve_partial_range(benchmark::State& state)
for(auto _ : state)
{
std::vector<Z> primes;
benchmark::DoNotOptimize(boost::math::prime_sieve(lower, upper, std::back_inserter(primes)));
benchmark::DoNotOptimize(boost::math::prime_range(lower, upper, std::back_inserter(primes)));
}
state.SetComplexityN(state.range(0));
}
Expand All @@ -51,16 +39,5 @@ BENCHMARK_TEMPLATE(prime_sieve, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 <
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();
17 changes: 8 additions & 9 deletions test/test_prime_sieve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,38 +19,38 @@ void test_prime_sieve()
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::math::prime_sieve(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::math::prime_sieve(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::math::prime_sieve(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::math::prime_sieve(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::math::prime_sieve(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::math::prime_sieve(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::math::prime_sieve(1000, std::back_inserter(d_primes));
BOOST_TEST_EQ(d_primes.size(), ref);
}

Expand Down Expand Up @@ -109,6 +109,7 @@ void test_prime_sieve_overflow()

int main()
{

test_prime_sieve<int>();
test_prime_sieve<int32_t>();
test_prime_sieve<int64_t>();
Expand All @@ -119,8 +120,6 @@ int main()
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();
Expand Down

0 comments on commit 349ca6f

Please sign in to comment.