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