Prime Numbers Generator after http://www.rsdn.ru/Forum/?mid=2969898

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#include <memory>
#include <math.h>
#include <time.h>
#include "../generator.h"
typedef unsigned word_t;
const word_t bits_in_word = 8 * sizeof(word_t);
template <typename Primes>
struct gen_map
{
typedef typename Primes::nat_type nat_type;
gen_map(nat_type p0, nat_type p1)
: lo_(p0 * p0)
, hi_(p1 * p1)
, bound_((hi_ - lo_) / 2 - 1)
, arr_(new word_t [bound_ / bits_in_word + 1])
, i_(~nat_type(0))
, n_(lo_)
{
for (word_t i = 0
; i != bound_ / bits_in_word + 1
; ++i)
arr_[i] = word_t(0);
Primes ps(3);
for (ps.next(); ps.value() != p1; ps.next())
{
const nat_type i = ps.value();
for (nat_type j = correct(i - lo_ % i, i) / 2 - 1
; j < bound_
; j += i)
arr_[j / bits_in_word] |= 1 << (j % bits_in_word);
}
next();
}
~gen_map()
{ delete[] arr_; }
bool avail() const
{ return i_ != bound_; }
nat_type value() const
{ return n_; }
void next()
{
do { ++i_; n_ += 2; }
while (
avail() &&
(arr_[i_ / bits_in_word] & (1 << (i_ % bits_in_word))));
}
private:
static nat_type correct(nat_type n, nat_type i)
{ return n & 1 ? n + i : n; }
const nat_type lo_, hi_, bound_;
word_t *arr_;
nat_type i_, n_;
};
template <typename Natural>
inline Natural sqrt_int(Natural x)
{ return static_cast<Natural>(sqrt(static_cast<double>(x))); }
template <typename Natural>
generator(primes, Natural, primes<Natural>)
{
public:
typedef Natural nat_type;
primes(nat_type n = 2)
: n_(n)
, sq_n_(sqrt_int(n))
{ }
emit(nat_type)
{
if (n_ <= 2) yield(2);
if (n_ <= 3) yield(3);
if (n_ <= 5) yield(5);
if (n_ <= 7) yield(7);
ps_.reset(new primes(3));
ps_->next();
p1_ = ps_->value();
for (;;)
{
p0_ = p1_;
ps_->next();
p1_ = ps_->value();
if (p1_ <= sq_n_)
continue;
gmap_.reset(new gen_map<primes>(p0_, p1_));
for (; gmap_->avail(); gmap_->next())
if (gmap_->value() >= n_)
yield(gmap_->value());
}
}
stop_emit()
private:
const nat_type n_, sq_n_;
std::auto_ptr<gen_map<primes> > gmap_;
nat_type p0_, p1_;
std::auto_ptr<primes> ps_;
};
int main()
{
primes<word_t> ps(100000000);
word_t sum = 0;
long tim0 = clock();
for (ps.next(); ps.value() <= 120000000; ps.next(), ++sum);
printf("sum = %u, time: %0.3f sec.\n", sum, (clock() - tim0) / 1000.0);
}