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