Loading...
Searching...
No Matches
prime_sieve.cpp
1#include <symengine/prime_sieve.h>
2#include <ciso646>
3#include <cmath>
4#include <valarray>
5#include <algorithm>
6#include <vector>
7#include <iterator>
8
9namespace SymEngine
10{
11
12static std::vector<unsigned> &sieve_primes()
13{
14 static std::vector<unsigned> primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
15 return primes;
16}
17
18bool Sieve::_clear = true;
19unsigned Sieve::_sieve_size = 32 * 1024 * 8; // 32K in bits
20
21void Sieve::set_clear(bool clear)
22{
23 _clear = clear;
24}
25
26void Sieve::clear()
27{
28 std::vector<unsigned> &_primes = sieve_primes();
29 _primes.erase(_primes.begin() + 10, _primes.end());
30}
31
32void Sieve::set_sieve_size(unsigned size)
33{
34#ifdef HAVE_SYMENGINE_PRIMESIEVE
35 primesieve::set_sieve_size(size);
36#else
37 _sieve_size = size * 1024 * 8; // size in bits
38#endif
39}
40
41void Sieve::_extend(unsigned limit)
42{
43 std::vector<unsigned> &_primes = sieve_primes();
44#ifdef HAVE_SYMENGINE_PRIMESIEVE
45 if (_primes.back() < limit)
46 primesieve::generate_primes(_primes.back() + 1, limit, &_primes);
47#else
48 const unsigned sqrt_limit
49 = static_cast<unsigned>(std::floor(std::sqrt(limit)));
50 unsigned start = _primes.back() + 1;
51 if (limit <= start)
52 return;
53 if (sqrt_limit >= start) {
54 _extend(sqrt_limit);
55 start = _primes.back() + 1;
56 }
57
58 unsigned segment = _sieve_size;
59 std::valarray<bool> is_prime(segment);
60 for (; start <= limit; start += 2 * segment) {
61 unsigned finish = std::min(start + segment * 2 + 1, limit);
62 is_prime[std::slice(0, segment, 1)] = true;
63 // considering only odd integers. An odd number n corresponds to
64 // n-start/2 in the array.
65 for (unsigned index = 1; index < _primes.size()
66 and _primes[index] * _primes[index] <= finish;
67 ++index) {
68 unsigned n = _primes[index];
69 unsigned multiple = (start / n + 1) * n;
70 if (multiple % 2 == 0)
71 multiple += n;
72 if (multiple > finish)
73 continue;
74 std::slice sl = std::slice((multiple - start) / 2,
75 1 + (finish - multiple) / (2 * n), n);
76 // starting from n*n, all the odd multiples of n are marked not
77 // prime.
78 is_prime[sl] = false;
79 }
80 for (unsigned n = start + 1; n <= finish; n += 2) {
81 if (is_prime[(n - start) / 2])
82 _primes.push_back(n);
83 }
84 }
85#endif
86}
87
88void Sieve::generate_primes(std::vector<unsigned> &primes, unsigned limit)
89{
90 _extend(limit);
91 std::vector<unsigned> &_primes = sieve_primes();
92 auto it = std::upper_bound(_primes.begin(), _primes.end(), limit);
93 // find the first position greater than limit and reserve space for the
94 // primes
95 primes.reserve(it - _primes.begin());
96 std::copy(_primes.begin(), it, std::back_inserter(primes));
97 if (_clear)
98 clear();
99}
100
101Sieve::iterator::iterator(unsigned max)
102{
103 _limit = max;
104 _index = 0;
105}
106
107Sieve::iterator::iterator()
108{
109 _limit = 0;
110 _index = 0;
111}
112
113Sieve::iterator::~iterator()
114{
115 if (_clear)
116 Sieve::clear();
117}
118
119unsigned Sieve::iterator::next_prime()
120{
121 std::vector<unsigned> &_primes = sieve_primes();
122 if (_index >= _primes.size()) {
123 unsigned extend_to = _primes[_index - 1] * 2;
124 if (_limit > 0 and _limit < extend_to) {
125 extend_to = _limit;
126 }
127 _extend(extend_to);
128 if (_index >= _primes.size()) { // the next prime is greater than _limit
129 return _limit + 1;
130 }
131 }
132 return _primes[_index++];
133}
134
135}; // namespace SymEngine
T back(T... args)
T back_inserter(T... args)
T begin(T... args)
static void generate_primes(std::vector< unsigned > &primes, unsigned limit)
Definition: prime_sieve.cpp:88
T copy(T... args)
T end(T... args)
T erase(T... args)
T floor(T... args)
T min(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
RCP< const Basic > max(const vec_basic &arg)
Canonicalize Max:
Definition: functions.cpp:3555
T push_back(T... args)
T reserve(T... args)
T size(T... args)
T sqrt(T... args)
T upper_bound(T... args)