Skip to content

Prime Number Functions #400

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
mborland wants to merge 117 commits into
base: develop
Choose a base branch
from
Draft

Prime Number Functions #400

Changes from 1 commit
Commits
Show all changes
117 commits
Select commit Hold shift + click to select a range
4015fbc
Initial Commit
mborland Jul 9, 2020
3ee737b
Changed init of least_divisors
mborland Jul 10, 2020
9512bb6
Cleanup
mborland Jul 10, 2020
d762398
Cleanup
mborland Jul 10, 2020
5375a1d
Added additional tests, benchmarks, and overflow checks
mborland Jul 10, 2020
a684dbd
Fix include guard naming
mborland Jul 10, 2020
3e4db8a
Complete revamp of algorithm. Hide implementation behind detail names…
mborland Jul 12, 2020
dd8a61c
Re-added support and tests for boost::multiprecision::cpp_int [CI SKIP]
mborland Jul 12, 2020
4debd0d
Changed benchmarks to support threading [CI SKIP]
mborland Jul 13, 2020
2a7e031
Added execution policies. Increased performance for dynamically linke…
mborland Jul 14, 2020
a68910e
Added massively parallel section to prime_sieve. Increased length of …
mborland Jul 14, 2020
386cc4d
Add test for multi-threading section and add to Jamfile [CI SKIP]
mborland Jul 15, 2020
d79eddb
Added prime sieve to existing prime numbers documentation [CI SKIP]
mborland Jul 15, 2020
3fdf917
Revisions to documentation and send to CI
mborland Jul 15, 2020
6ca245b
Changed include guards to be compatible with C++11 and 14.
mborland Jul 15, 2020
7d3a520
Fixed doc, and ensured that primes are sorted
mborland Jul 15, 2020
a1ac504
Fixed documentation. Complete re-design of mask_sieve algo. Pre-gener…
mborland Jul 16, 2020
35d2aa1
All vectors now init {}. Change include guards. Replace raw pointer.
mborland Jul 17, 2020
243a299
Changed from [lower_bound, upper_bound] to [lower_bound, upper_bound)…
mborland Jul 18, 2020
2eadeff
Fixed documentation [CI SKIP]
mborland Jul 18, 2020
173ce0d
segmented_sieve now runs using std::async. [CI SKIP]
mborland Jul 18, 2020
23fba36
Removed extraneous operations [CI SKIP]
mborland Jul 19, 2020
f7b45fd
Cleanup style, and delete unused function. Enable par_unseq mask_siev…
mborland Jul 19, 2020
d687d5e
Small Cleanup and limit reduction for sieving methods.
mborland Jul 20, 2020
e51d727
Add pritchard's sub-linear algo [WIP][CI SKIP]
mborland Jul 25, 2020
1245d27
Added pritchards segmented sieve [WIP][CI SKIP]
mborland Jul 26, 2020
2052081
Implemented binary search in SetS [WIP][CI SKIP]
mborland Jul 27, 2020
9f81e0d
Tests and performance compairsons [WIP][CI SKIP]
mborland Jul 27, 2020
55ea045
More tests. Segmented sieve small cases fix
mborland Jul 27, 2020
31a2105
Imporve segmented performance [WIP][CI SKIP]
mborland Jul 29, 2020
f545928
Replaced SetS members with stl algos [WIP][CI SKIP]
mborland Jul 29, 2020
d780db5
Replace searching with tracked index. General performance improvement…
mborland Jul 31, 2020
2639bed
Minor change to SetS remove [WIP][CI SKIP]
mborland Aug 1, 2020
94fc1ac
Pritchard segmented performance improvements [WIP][CI SKIP]
mborland Aug 1, 2020
6ebb906
Fixed failed test [WIP][CI SKIP]
mborland Aug 1, 2020
0521854
Various performance improvements [WIP][CI SKIP]
mborland Aug 2, 2020
b233a80
Build primes in situ [WIP][CI SKIP]
mborland Aug 2, 2020
f95c2cf
Refactoring. Now requires C++17
mborland Aug 2, 2020
8e2e29a
Add segmented bit sieve [WIP][CI SKIP]
mborland Aug 4, 2020
1a24f16
Removed segmented bit sieve and excess headers [CI SKIP]
mborland Aug 6, 2020
c63f1f1
Added wheel class [WIP][CI SKIP]
Aug 15, 2020
b7d4256
Added fixed mod 210 wheel [WIP][CI SKIP]
Aug 23, 2020
fbc38c8
New segmented sieve algorithm [WIP][CI SKIP]
Aug 23, 2020
2e46b81
Added interval sieve to performance test [CI SKIP]
Aug 23, 2020
6759ede
Added unit tests [WIP][CI SKIP]
Aug 23, 2020
1b40403
Implemented interval sieve [CI SKIP]
Aug 23, 2020
97244be
Performance improvements and bug fixes [CI SKIP]
Aug 24, 2020
fa04133
Significant refactoring [WIP][CI SKIP]
Aug 25, 2020
c4a89c8
Seq policy actually sequential [CI SKIP]
Aug 25, 2020
91836f6
Fixes for multiprecision and policies [CI SKIP]
Aug 25, 2020
6d6b19f
cpp_int now passes all tests [CI SKIP]
Aug 27, 2020
81e4a6c
mpz_int passes unit tests [CI SKIP]
Aug 27, 2020
0b8f1d5
Documentation edits [CI SKIP]
mborland Aug 30, 2020
5ba0a1d
Replace magic number w variable template [CI SKIP]
mborland Aug 31, 2020
8e240f7
Better prime_range implementation [WIP][CI SKIP]
mborland Aug 31, 2020
7c0cbbf
Merge branch 'develop' into prime_functions
mborland Sep 1, 2020
3d9b77c
Minor changes and doc updates
mborland Sep 2, 2020
302fb5f
Merge branch 'develop' into prime_functions
mborland Sep 5, 2020
9792a23
Minor change to policy handling [CI SKIP]
mborland Sep 6, 2020
84a69f0
diffs from @jzmaddock [CI SKIP]
mborland Sep 8, 2020
5dc3523
Implemented linear sieve with iterators
mborland Sep 9, 2020
0dbe69c
Sanitize linear prime sieve and add wheel
mborland Sep 9, 2020
eee2c86
Added container method
mborland Sep 26, 2020
f5d789a
Improved Linear Algo and testing
mborland Sep 27, 2020
f2277e3
Merge branch 'prime_iterator' into mborland/prime_iterator
mborland Sep 27, 2020
1d2f03c
Linear sieve with output iterator and refactoring [CI SKIP][WIP]
mborland Sep 27, 2020
0d9d31b
Merge branch 'prime_functions' into mborland/prime_functions
mborland Sep 27, 2020
c361cde
Linear output iterator and refactoring [CI SKIP]
mborland Sep 27, 2020
66c2642
Added prime approximation function
mborland Sep 27, 2020
de1f331
prime approximation function for a range
mborland Sep 27, 2020
eaea5f9
Added interval sieve for output iterators
mborland Sep 27, 2020
e8196f3
Merge remote-tracking branch 'origin/develop' into prime_iterator
mborland Sep 27, 2020
cb5d978
Interval sieve with OI passes unit tests [CI SKIP]
mborland Sep 27, 2020
830ccc4
Merge branch 'prime_iterator' into prime_functions
mborland Sep 27, 2020
4dc3eb2
Add sequential composite sieve with OI [CI SKIP]
mborland Sep 30, 2020
90be100
Add threaded method. Currently INOP [CI SKIP]
mborland Sep 30, 2020
a772782
Threaded method completed and validated [CI SKIP]
mborland Oct 1, 2020
2d1461f
Resolves issue #439 [CI SKIP]
mborland Oct 1, 2020
e7cdb32
Updated benchmarks [CI SKIP]
mborland Oct 2, 2020
29eef88
Correct benchmark memory allocation [CI SKIP]
mborland Oct 3, 2020
b5a28b5
Add linear sieve direct from stepanov [CI SKIP]
mborland Oct 4, 2020
6c26b53
First draft of eratosthenes w/ wheel [CI SKIP]
mborland Oct 4, 2020
980bfe7
Fixes for multiprecision types [CI SKIP]
mborland Oct 4, 2020
07e6f58
wheel sieve passes standard battery [CI SKIP]
mborland Oct 4, 2020
86b9e5a
Add prime sieve wrapper [CI SKIP]
mborland Oct 5, 2020
f19149e
Implement dual interface for prime sieve [CI SKIP]
mborland Oct 7, 2020
1ad0d51
Implement dual interface for prime_range [CI SKIP]
mborland Oct 7, 2020
5f6d06a
Remove uneeded vals from wheel [CI SKIP]
mborland Nov 22, 2020
670b06d
Internal data structure changes [CI SKIP]
mborland Nov 22, 2020
7c7a491
Rename type alias [CI SKIP]
mborland Nov 23, 2020
e507ba5
Swap mod 210 wheel for mod 30 wheel [CI SKIP]
mborland Nov 23, 2020
c9cb41c
Add simple bitset file [CI SKIP]
mborland Nov 23, 2020
e8b71a0
Reduced unnecessary traversing of bitset [CI SKIP]
mborland Nov 24, 2020
d649f65
Updated bitset [CI SKIP]
mborland Nov 24, 2020
ed98892
Add resize to simple_bitset [CI SKIP]
mborland Nov 26, 2020
c3b3934
Implement simple_bitset as internal data structure
mborland Nov 26, 2020
e9aa05d
Add prime-composite ratio to wheel to size bitset
mborland Nov 27, 2020
7c5d792
minor refactoring [CI SKIP]
mborland Nov 27, 2020
2bbfad2
Fix PrimeRatio [CI SKIP]
mborland Nov 27, 2020
ddc5c8a
Improved prime range [CI SKIP]
mborland Nov 27, 2020
f517c00
Add hamming weight calc to simple bitset [CI SKIP]
mborland Nov 28, 2020
a356b47
Add forward and reverse bit scan to simple bitset
mborland Nov 28, 2020
2212e45
Fix count, enable bit scan by type, cleanup
mborland Nov 28, 2020
b12e2dc
Add dual interface. Avoid using multiprecision
mborland Nov 28, 2020
f39a4d8
Seperate approximation for multiprecision types
mborland Nov 29, 2020
eafbefc
Remove all mulitples of wheel basis from bitset
mborland Nov 29, 2020
cb36b89
IntervalSieve performance improvements [CI SKIP]
mborland Nov 29, 2020
16c2354
Add small prime optimization [CI SKIP]
mborland Dec 5, 2020
0b1a690
Add small primes to full sieve [CI SKIP]
mborland Dec 8, 2020
8c883d1
Merge remote-tracking branch 'origin/develop' into prime_functions
mborland Dec 12, 2020
f357e0f
Implement bit scanning in sieve [CI SKIP]
mborland Dec 12, 2020
a4f2d89
Fix minor errors
mborland Dec 12, 2020
d16e562
Add trial divison to verify tests [CI SKIP]
mborland Dec 12, 2020
91be2c3
Fix small primes [CI SKIP]
mborland Dec 12, 2020
d9ae8c2
Revert interval_prime_sieve.hpp to f39a4d8
mborland Dec 13, 2020
7d22010
Remove cruft [CI SKIP]
mborland Dec 13, 2020
b541987
Continue cleanup [CI SKIP]
mborland Dec 13, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
All vectors now init {}. Change include guards. Replace raw pointer.
  • Loading branch information
@mborland
mborland committed Jul 17, 2020
commit 35d2aa1ec7230f24e1d6a21efe81f15103426a5a
47 changes: 23 additions & 24 deletions include/boost/math/special_functions/prime_sieve.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,17 @@
#include <iterator>
#include <cmath>
#include <thread>
#include <memory>

#ifdef _MSVC_LANG
#if _MSVC_LANG >= 201703 // _MSVC_LANG == __cplusplus: https://devblogs.microsoft.com/cppblog/msvc-now-correctly-reports-__cplusplus/
#include <execution>
#endif
#else
#if __cplusplus >= 201703
#include <execution>
#endif
#endif

namespace boost { namespace math { namespace detail
{
Expand All @@ -28,7 +35,7 @@ template<class Z, class Container>
void linear_sieve(Z upper_bound, Container &c)
Copy link
Contributor

@jeremy-murphy jeremy-murphy Jul 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please use something more descriptive than c here. :) Surely the values in this container have some meaning, especially since it's the prime numbers result, right?

{
size_t least_divisors_size{static_cast<size_t>(upper_bound + 1)};
Z *least_divisors{new Z[least_divisors_size]{0}};
std::unique_ptr<Z[]> least_divisors{new Z[least_divisors_size]{0}};
Copy link
Contributor

@jeremy-murphy jeremy-murphy Jul 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can't this just be vector<Z>?

Copy link
Member Author

@mborland mborland Jul 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It could be, but using a C-style array instead of an std::vector yielded lower memory usage and complexity.

Copy link
Contributor

@jeremy-murphy jeremy-murphy Jul 27, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmmm, that doesn't sound right, though, does it? Like, memory usage could be marginally less with the raw array, but complexity should be identical. What kind of difference are we talking here?

Copy link
Member Author

@mborland mborland Jul 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is in micro-optimization territory. I just re-ran it and the difference is 3.72N with the C-style array vs 3.86N for the vector.

Copy link
Contributor

@jeremy-murphy jeremy-murphy Jul 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are those values consistent between different runs of the benchmark?

Copy link
Contributor

@jeremy-murphy jeremy-murphy Jul 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Although I would generally argue not to leave any performance improvement on the table, I feel like std::vector shouldn't be any worse here. At the risk of sounding demanding... have you tried another standard library implementation?

Copy link
Member Author

@mborland mborland Jul 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

They are; benchmark's variation is usually pretty small. These were about +-0.03N.

Copy link
Member Author

@mborland mborland Jul 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Compiling with Intel C++ and their implementation of the standard library yields the same complexity for the c-style array and the std::vector. I guess I would err on the side of more users compiling with gcc/clang and libstdc++ than Intel C++.


for (Z i{2}; i <= upper_bound; ++i)
{
Expand Down Expand Up @@ -61,8 +68,6 @@ void linear_sieve(Z upper_bound, Container &c)
}
}
}

delete[] least_divisors;
}

template<class Z, class Container>
Expand All @@ -85,20 +90,14 @@ void mask_sieve(Z lower_bound, Z upper_bound, const PrimeContainer& primes, Cont
{
Z limit {static_cast<Z>(std::floor(std::sqrt(static_cast<double>(upper_bound)))) + 1};

size_t primes_size;
size_t primes_size {};
auto it{primes.begin()};
while(*it < limit && it != primes.end())
while(it != primes.end() && *it < limit)
{
++primes_size;
++it;
}

// Needed for thread sanitizer. Throws FPE without this check.
if(primes_size > primes.size())
{
primes_size = primes.size();
}

// Faster and safer than std::vector<bool> which does not behave like all other vectors
const size_t n {static_cast<size_t>(upper_bound - lower_bound + 1)};
bool* is_prime {new bool[n]};
Expand Down Expand Up @@ -132,7 +131,7 @@ template<typename Z, class Container>
void mask_sieve(Z lower_bound, Z upper_bound, Container &c)
{
Z limit{static_cast<Z>(std::floor(std::sqrt(static_cast<double>(upper_bound)))) + 1};
std::vector<Z> primes;
std::vector<Z> primes {};
primes.reserve(limit / std::log(limit));

boost::math::detail::linear_sieve(limit, primes);
Expand All @@ -146,7 +145,7 @@ template <typename Z, class PrimesContainer, class Container>
void segmented_sieve(Z lower_bound, Z upper_bound, const PrimesContainer &primes, Container &c)
{
const Z L1_SIZE {32648};
const Z interval {static_cast<Z>(std::floor(L1_SIZE / sizeof(Z)))};
const Z interval {L1_SIZE * 4};
Z current_lower_bound{lower_bound};
Z current_upper_bound{current_lower_bound + interval};

Expand All @@ -169,7 +168,7 @@ template <typename Z, class Container>
void segmented_sieve(Z lower_bound, Z upper_bound, Container &c)
{
Z limit{static_cast<Z>(std::floor(std::sqrt(static_cast<double>(upper_bound)))) + 1};
std::vector<Z> primes;
std::vector<Z> primes {};
primes.reserve(limit / std::log(limit));

// Prepare for max value so you do not have to calculate this again
Expand All @@ -194,7 +193,7 @@ auto prime_sieve(ExecutionPolicy&& policy, Z upper_bound, OutputIterator 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;
std::vector<Z> primes {};
primes.reserve(upper_bound / std::log(static_cast<double>(upper_bound)));

// Range for when the linear sieve is no longer faster than threading
Expand All @@ -219,7 +218,7 @@ auto prime_sieve(ExecutionPolicy&& policy, Z upper_bound, OutputIterator output)
processor_count = 2;
}

std::vector<Z> small_primes;
std::vector<Z> small_primes {};
small_primes.reserve(1000);

// Threshold for when 2 thread performance begins to be non-linear, or when the system can only support two threads
Expand All @@ -243,7 +242,7 @@ auto prime_sieve(ExecutionPolicy&& policy, Z upper_bound, OutputIterator output)
{
//Pre-generate all of the primes so that each thread does not have to
Z limit{static_cast<Z>(std::floor(std::sqrt(static_cast<double>(upper_bound)))) + 1};
std::vector<Z> pre_generated_primes;
std::vector<Z> pre_generated_primes {};
pre_generated_primes.reserve(limit / std::log(limit));

if(limit < 8192)
Expand All @@ -265,7 +264,7 @@ auto prime_sieve(ExecutionPolicy&& policy, Z upper_bound, OutputIterator output)
pre_generated_primes.insert(pre_generated_primes.begin(), small_primes.begin(), small_primes.end());
}

std::vector<std::thread> thread_manager;
std::vector<std::thread> thread_manager {};
std::vector<std::vector<Z>> prime_vectors(processor_count - 1);
const Z range_per_thread = upper_bound / (processor_count - 1);
Z current_lower_bound {8192};
Expand All @@ -275,7 +274,7 @@ auto prime_sieve(ExecutionPolicy&& policy, Z upper_bound, OutputIterator output)

for(size_t i{1}; i < processor_count - 1; ++i)
{
std::vector<Z> temp;
std::vector<Z> temp {};
temp.reserve(primes_in_range);
prime_vectors.emplace_back(temp);

Expand All @@ -292,7 +291,7 @@ auto prime_sieve(ExecutionPolicy&& policy, Z upper_bound, OutputIterator output)
current_lower_bound / std::log(static_cast<double>(current_lower_bound)));
}

std::vector<Z> temp;
std::vector<Z> temp {};
temp.reserve(primes_in_range);
prime_vectors.emplace_back(temp);

Expand Down Expand Up @@ -320,7 +319,7 @@ auto prime_sieve(ExecutionPolicy&& policy, Z upper_bound, OutputIterator 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;
std::vector<Z> primes {};
primes.reserve(upper_bound / std::log(static_cast<double>(upper_bound)));

boost::math::prime_sieve(policy, upper_bound, std::back_inserter(primes));
Expand All @@ -341,7 +340,7 @@ 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;
std::vector<Z> primes{};
primes.reserve(upper_bound / std::log(upper_bound));

if (upper_bound <= 8192)
Expand All @@ -351,7 +350,7 @@ auto prime_sieve(Z upper_bound, OutputIterator output) -> decltype(output)

else
{
boost::math::detail::mask_sieve(static_cast<Z>(2), upper_bound, primes);
boost::math::detail::segmented_sieve(static_cast<Z>(2), upper_bound, primes);
}

return std::move(primes.begin(), primes.end(), output);
Expand All @@ -360,7 +359,7 @@ auto prime_sieve(Z upper_bound, OutputIterator output) -> decltype(output)
template<class Z, class OutputIterator>
auto prime_range(Z lower_bound, Z upper_bound, OutputIterator output) -> decltype(output)
Copy link
Contributor

@jeremy-murphy jeremy-murphy Jul 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure it's worth including this functionality if it's not computationally more efficient.
After all, this will do the same amount of work, but also allocate extra temporary memory just to omit the primes lower than the lower_bound. The user can remove the primes less than the lower_bound without using extra temporary memory.

Copy link
Member Author

@mborland mborland Jul 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My thought process was it would make life easier in that use case. If the user were to remove the primes less than lower_bound themselves the memory usage should be the same as moving them the ones they wanted out of the full range [2, upper_bound) they would have received anyway.

Copy link
Contributor

@jeremy-murphy jeremy-murphy Jul 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To be precise, the peak memory allocation will be worse this way. Let's say,
x = upper_bound / std::log(static_cast<double>(upper_bound)), and
y = (upper_bound - lower_bound) / std::log(static_cast<double>(upper_bound - lower_bound))
Then using this function, the user allocates size y for the final result, and this function allocates x for the temporary result of prime_range(upper). So the peak memory usage is x + y (ignoring other temporary allocations).

Without this function, the user allocates size x for the result, then calls prime_range(upper) directly, so the peak memory usage is x (again, ignoring other temporary allocations). The user then executes something like move(begin(result) + lower_bound, end(result), begin(result)); result.resize(upper_bound - lower_bound); result.shrink_to_fit(); to remove the unwanted values and return unused capacity to the system. It's messier and more work for the user, but that's the cost of having lower peak memory overhead.

Is it worth it? Who knows. On reflection, there is no problem with leaving this convenience wrapper in -- users concerned with peak memory allocation more than convenience will just do the second method.

{
std::vector<Z> primes;
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));
Expand Down