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);
}
|
