Skip to content

Commit a2aedbb

Browse files
mborlandmborland
committed
Complete revamp of algorithm. Hide implementation behind detail namespace.
1 parent a684dbd commit a2aedbb

File tree

3 files changed

+117
-67
lines changed

3 files changed

+117
-67
lines changed

‎include/boost/math/special_functions/prime_sieve.hpp

+107-33
Original file line numberDiff line numberDiff line change
@@ -14,88 +14,162 @@
1414
#include <vector>
1515
#include <iterator>
1616
#include <iostream>
17+
#include <cmath>
18+
#include <thread>
1719

18-
namespace boost { namespace math
20+
namespace boost { namespace math { namespace detail
1921
{
20-
2122
// https://mathworld.wolfram.com/SieveofEratosthenes.html
2223
// https://www.cs.utexas.edu/users/misra/scannedPdf.dir/linearSieve.pdf
23-
template<class Z, class OutputIterator>
24-
auto prime_sieve(Z lower_bound, Z upper_bound, OutputIterator output) -> decltype(output)
24+
template<class Z, class Container>
25+
void linear_sieve(Z upper_bound, Container &c)
2526
{
26-
static_assert(std::is_integral<Z>::value, "No primes for floating point types");
27-
BOOST_ASSERT_MSG(upper_bound + 1 < std::numeric_limits<Z>::max(), "Type Overflow");
28-
std::vector<Z> least_divisors(upper_bound + 1, 0);
29-
std::deque<Z> primes;
27+
Z least_divisors_size{upper_bound + 1};
28+
Z *least_divisors{new Z[least_divisors_size]{0}};
3029

3130
for (Z i{2}; i <= upper_bound; ++i)
3231
{
3332
if (least_divisors[i] == 0)
3433
{
3534
least_divisors[i] = i;
36-
primes.emplace_back(i);
35+
c.emplace_back(i);
3736
}
3837

39-
for (size_t j{}; j < least_divisors.size(); ++j)
38+
for (size_t j{}; j < least_divisors_size; ++j)
4039
{
41-
if (j >= primes.size())
40+
if (j >= c.size())
4241
{
4342
break;
4443
}
4544

46-
else if (primes[j] > least_divisors[i])
45+
else if (c[j] > least_divisors[i])
4746
{
4847
break;
4948
}
5049

51-
else if (i * primes[j] > upper_bound)
50+
else if (i * c[j] > upper_bound)
5251
{
5352
break;
5453
}
5554

5655
else
5756
{
58-
least_divisors[i * primes[j]] = primes[j];
57+
least_divisors[i * c[j]] = c[j];
5958
}
6059
}
6160
}
6261

63-
auto it{primes.begin()};
64-
while (*it < lower_bound && it != primes.end())
62+
delete[] least_divisors;
63+
}
64+
65+
template<class Z, class Container>
66+
void prime_table(Z upper_bound, Container &c)
67+
{
68+
Z i{2};
69+
unsigned counter{};
70+
71+
while (i <= upper_bound && counter < 9999) // 10k elements are in the lookup table
6572
{
66-
primes.pop_front();
67-
++it;
73+
c.emplace_back(i);
74+
++counter;
75+
i = static_cast<Z>(boost::math::prime(counter));
6876
}
69-
70-
return std::move(primes.begin(), primes.end(), output);
7177
}
7278

73-
template<class Z, class OutputIterator>
74-
auto prime_range(Z lower_bound, Z upper_bound, OutputIterator output) -> decltype(output)
79+
template<class Z, class Container>
80+
void mask_sieve(Z lower_bound, Z upper_bound, Container &c)
7581
{
76-
if (upper_bound <= 104729)
82+
Z limit{static_cast<Z>(std::floor(std::sqrt(upper_bound))) + 1};
83+
std::vector<Z> primes;
84+
primes.reserve(limit / std::log(limit));
85+
86+
boost::math::detail::linear_sieve(limit, primes);
87+
88+
const Z n{upper_bound - lower_bound + 1};
89+
bool *mask{new bool[n + 1]{false}};
90+
91+
for (size_t i{}; i < primes.size(); ++i)
7792
{
78-
Z i{2};
79-
unsigned counter {};
80-
std::deque<Z> primes;
81-
while (i <= upper_bound)
93+
Z lower_limit = std::floor(lower_bound / primes[i]) * primes[i];
94+
95+
if (lower_limit < lower_bound)
96+
{
97+
lower_limit += primes[i];
98+
}
99+
100+
if (lower_limit == primes[i])
101+
{
102+
lower_limit += primes[i];
103+
}
104+
105+
for (Z j{lower_limit}; j <= upper_bound; j += primes[i])
106+
{
107+
mask[j - lower_bound] = true;
108+
}
109+
}
110+
111+
// Numbers which are not masked in range, are prime
112+
for (Z i{lower_bound}; i <= upper_bound; i++)
113+
{
114+
if (!mask[i - lower_bound])
82115
{
83116
if (i >= lower_bound)
84117
{
85-
primes.emplace_back(i);
118+
c.emplace_back(i);
86119
}
87-
88-
++counter;
89-
i = static_cast<Z>(boost::math::prime(counter));
90120
}
121+
}
122+
123+
delete[] mask;
124+
}
125+
} // End namespace detail
126+
127+
template<typename Z, class OutputIterator>
128+
auto prime_sieve(Z upper_bound, OutputIterator output) -> decltype(output)
129+
{
130+
static_assert(std::is_integral<Z>::value, "No primes for floating point types");
131+
BOOST_ASSERT_MSG(upper_bound + 1 < std::numeric_limits<Z>::max(), "Type Overflow");
91132

92-
return std::move(primes.begin(), primes.end(), output);
133+
std::vector<Z> primes;
134+
primes.reserve(upper_bound / std::log(upper_bound));
135+
136+
if (upper_bound <= 104729)
137+
{
138+
boost::math::detail::prime_table(upper_bound, primes);
93139
}
94140

95141
else
96142
{
97-
return prime_sieve(lower_bound, upper_bound, output);
143+
std::vector<Z> small_primes;
144+
small_primes.reserve(1000);
145+
146+
// Spilt into two vectors and merge after joined to avoid data races
147+
std::thread t1([upper_bound, &small_primes]{boost::math::detail::prime_table(static_cast<Z>(104729), small_primes);});
148+
std::thread t2([upper_bound, &primes]{boost::math::detail::mask_sieve(static_cast<Z>(104729), upper_bound, primes);});
149+
150+
t1.join();
151+
t2.join();
152+
primes.insert(primes.begin(), small_primes.begin(), small_primes.end());
153+
}
154+
155+
return std::move(primes.begin(), primes.end(), output);
156+
}
157+
158+
template<class Z, class OutputIterator>
159+
auto prime_range(Z lower_bound, Z upper_bound, OutputIterator output) -> decltype(output)
160+
{
161+
std::vector<Z> primes;
162+
primes.reserve(upper_bound / std::log(upper_bound));
163+
164+
boost::math::prime_sieve(upper_bound, std::back_inserter(primes));
165+
166+
auto it{primes.begin()};
167+
while(*it < lower_bound && it != primes.end())
168+
{
169+
++it;
98170
}
171+
172+
return std::move(it, primes.end(), output);
99173
}
100174

101175
template<class Z, class OutputIterator>

‎reporting/performance/prime_sieve_performance.cpp

+2-25
Original file line numberDiff line numberDiff line change
@@ -15,19 +15,7 @@ void prime_sieve(benchmark::State& state)
1515
for(auto _ : state)
1616
{
1717
std::vector<Z> primes;
18-
benchmark::DoNotOptimize(boost::math::prime_sieve(static_cast<Z>(2), upper, std::back_inserter(primes)));
19-
}
20-
state.SetComplexityN(state.range(0));
21-
}
22-
23-
template <typename Z>
24-
void prime_range(benchmark::State& state)
25-
{
26-
Z upper = static_cast<Z>(state.range(0));
27-
for(auto _ : state)
28-
{
29-
std::vector<Z> primes;
30-
benchmark::DoNotOptimize(boost::math::prime_range(static_cast<Z>(2), upper, std::back_inserter(primes)));
18+
benchmark::DoNotOptimize(boost::math::prime_sieve(upper, std::back_inserter(primes)));
3119
}
3220
state.SetComplexityN(state.range(0));
3321
}
@@ -40,7 +28,7 @@ void prime_sieve_partial_range(benchmark::State& state)
4028
for(auto _ : state)
4129
{
4230
std::vector<Z> primes;
43-
benchmark::DoNotOptimize(boost::math::prime_sieve(lower, upper, std::back_inserter(primes)));
31+
benchmark::DoNotOptimize(boost::math::prime_range(lower, upper, std::back_inserter(primes)));
4432
}
4533
state.SetComplexityN(state.range(0));
4634
}
@@ -51,16 +39,5 @@ BENCHMARK_TEMPLATE(prime_sieve, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 <
5139
BENCHMARK_TEMPLATE(prime_sieve_partial_range, int32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
5240
BENCHMARK_TEMPLATE(prime_sieve_partial_range, int64_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
5341
BENCHMARK_TEMPLATE(prime_sieve_partial_range, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
54-
BENCHMARK_TEMPLATE(prime_range, int32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
55-
BENCHMARK_TEMPLATE(prime_range, int64_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
56-
BENCHMARK_TEMPLATE(prime_range, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 22)->Complexity();
57-
58-
// Direct comparison of lookup vs sieve using only range of lookup
59-
BENCHMARK_TEMPLATE(prime_sieve, int32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
60-
BENCHMARK_TEMPLATE(prime_range, int32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
61-
BENCHMARK_TEMPLATE(prime_sieve, int64_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
62-
BENCHMARK_TEMPLATE(prime_range, int64_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
63-
BENCHMARK_TEMPLATE(prime_sieve, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
64-
BENCHMARK_TEMPLATE(prime_range, uint32_t)->RangeMultiplier(2)->Range(1 << 1, 1 << 16)->Complexity();
6542

6643
BENCHMARK_MAIN();

‎test/test_prime_sieve.cpp

+8-9
Original file line numberDiff line numberDiff line change
@@ -19,38 +19,38 @@ void test_prime_sieve()
1919
Z ref {168}; // Calculated with wolfram-alpha
2020

2121
// Does the function work with a vector
22-
boost::math::prime_sieve(2, 1000, std::back_inserter(primes));
22+
boost::math::prime_sieve(1000, std::back_inserter(primes));
2323
BOOST_TEST_EQ(primes.size(), ref);
2424

2525
// Tests for correctness
2626
// 100
2727
primes.clear();
28-
boost::math::prime_sieve(2, 100, std::back_inserter(primes));
28+
boost::math::prime_sieve(100, std::back_inserter(primes));
2929
BOOST_TEST_EQ(primes.size(), 25);
3030

3131
// 10'000
3232
primes.clear();
33-
boost::math::prime_sieve(2, 10000, std::back_inserter(primes));
33+
boost::math::prime_sieve(10000, std::back_inserter(primes));
3434
BOOST_TEST_EQ(primes.size(), 1229);
3535

3636
// 100'000
3737
primes.clear();
38-
boost::math::prime_sieve(2, 100000, std::back_inserter(primes));
38+
boost::math::prime_sieve(100000, std::back_inserter(primes));
3939
BOOST_TEST_EQ(primes.size(), 9592);
4040

4141
// 1'000'000
4242
primes.clear();
43-
boost::math::prime_sieve(2, 1000000, std::back_inserter(primes));
43+
boost::math::prime_sieve(1000000, std::back_inserter(primes));
4444
BOOST_TEST_EQ(primes.size(), 78498);
4545

4646
// Does the function work with a list?
4747
std::list<Z> l_primes;
48-
boost::math::prime_sieve(2, 1000, std::back_inserter(l_primes));
48+
boost::math::prime_sieve(1000, std::back_inserter(l_primes));
4949
BOOST_TEST_EQ(l_primes.size(), ref);
5050

5151
// Does the function work with a deque?
5252
std::deque<Z> d_primes;
53-
boost::math::prime_sieve(2, 1000, std::back_inserter(d_primes));
53+
boost::math::prime_sieve(1000, std::back_inserter(d_primes));
5454
BOOST_TEST_EQ(d_primes.size(), ref);
5555
}
5656

@@ -109,6 +109,7 @@ void test_prime_sieve_overflow()
109109

110110
int main()
111111
{
112+
112113
test_prime_sieve<int>();
113114
test_prime_sieve<int32_t>();
114115
test_prime_sieve<int64_t>();
@@ -119,8 +120,6 @@ int main()
119120
test_prime_range<int64_t>();
120121
test_prime_range<uint32_t>();
121122

122-
test_prime_sieve<boost::multiprecision::cpp_int>();
123-
124123
//test_prime_sieve_overflow<int16_t>();
125124

126125
boost::report_errors();

0 commit comments

Comments
 (0)