forked from preda/gpuowl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathPm1Plan.cpp
314 lines (255 loc) · 7.88 KB
/
Pm1Plan.cpp
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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
#include "Pm1Plan.h"
#include <tuple>
#include <array>
#include <cassert>
#include <numeric>
vector<bool> Pm1Plan::sieve(u32 B1, u32 B2) {
assert(B1 && B1 < B2);
vector<bool> bits(B2 + 1);
bits[0] = bits[1] = true;
for (u32 p = 0; p <= B2; ++p) {
while (p <= B2 && bits[p]) { ++p; }
if (p > B2) { break; }
if (p <= B1) { bits[p] = true; }
if (p < (1u << 16) && p * p <= B2) { for (u32 i = p * p; i <= B2; i += p) { bits[i] = true; }}
}
bits.flip();
return bits;
}
namespace {
template<u32 D> constexpr bool isRelPrime(u32 j);
template<> constexpr bool isRelPrime<210>(u32 j) { return j % 2 && j % 3 && j % 5 && j % 7; }
template<> constexpr bool isRelPrime<330>(u32 j) { return j % 2 && j % 3 && j % 5 && j % 11; }
template<> constexpr bool isRelPrime<462>(u32 j) { return j % 2 && j % 3 && j % 7 && j % 11; }
template<> constexpr bool isRelPrime<770>(u32 j) { return j % 2 && j % 5 && j % 7 && j % 11; }
template<> constexpr bool isRelPrime<2310>(u32 j) { return j % 2 && j % 3 && j % 5 && j % 7 && j % 11; }
template<typename T>
u32 sum(const T &v) { return accumulate(v.begin(), v.end(), 0); }
// Repeatedly divides "pos" by "F" while it can, keeping it above B1.
template<u32 F>
u32 reduce(u32 B1, u32 pos) {
while (pos > B1 * F && pos % F == 0) { pos /= F; }
return pos;
}
constexpr u32 firstMissingFactor(u32 D) {
switch (D) {
case 210:
case 420:
return 11;
case 330:
case 660:
return 7;
case 462:
case 2*462:
return 5;
case 770:
case 2*770:
return 3;
case 2310:
return 13;
}
assert(false);
}
// Repeatedly divides "pos" by the smallest missing factor of "D".
u32 reduce(u32 D, u32 B1, u32 pos) {
switch (D) {
case 210:
case 420:
return reduce<11>(B1, pos);
case 330:
case 660:
return reduce<7>(B1, pos);
case 462:
case 2*462:
return reduce<5>(B1, pos);
case 770:
case 2*770:
return reduce<3>(B1, pos);
case 2310:
return reduce<13>(B1, pos);
}
assert(false);
}
// Returns whether GCD(D, j)==1.
bool isRelPrime(u32 D, u32 j) {
switch (D) {
case 210:
case 2*210:
return isRelPrime<210>(j);
case 330:
case 2*330:
return isRelPrime<330>(j);
case 462:
case 2*462:
return isRelPrime<462>(j);
case 770:
case 2*770:
return isRelPrime<770>(j);
case 2310:
return isRelPrime<2310>(j);
}
assert(false);
}
}
// Returns the minimal number of buffers J needed for D.
u32 Pm1Plan::minBufsFor(u32 D) {
switch (D) {
case 210: return 24;
case 2*210: return 2*24;
case 330: return 40;
case 2*330: return 2*40;
case 462: return 60;
case 2*462: return 2*60;
case 770: return 120;
case 2*770: return 2*120;
case 2310: return 240;
}
assert(false);
}
u32 Pm1Plan::reduce(u32 pos) const { return ::reduce(D, B1, pos); }
template<u32 F> u32 Pm1Plan::reduce(u32 pos) const { return ::reduce<F>(B1, pos); }
Pm1Plan::Pm1Plan(u32 argsD, u32 nBuf, u32 B1, u32 B2, vector<bool>&& primeBits)
: nBuf{min(nBuf, MAX_BUFS)}, primeBits{std::move(primeBits)}, D{getD(argsD, nBuf)}, B1{B1}, B2{B2}, jset{makeJset()} {
assert(nBuf >= 24);
assert(nBuf >= minBufsFor(D));
assert(B1 < B2);
// Extend the primeBits vector with guard zero bits, so that makePlan() does not read past-the-end.
for (u32 i = 0; i < 2 * jset.back() + 1; ++i) { this->primeBits.push_back(false); }
}
Pm1Plan::Pm1Plan(u32 D, u32 nBuf, u32 B1, u32 B2) : Pm1Plan{D, nBuf, B1, B2, sieve(B1, B2)} {
}
vector<u32> Pm1Plan::makeJset() {
vector<u32> jset;
assert(nBuf >= minBufsFor(D));
for (u32 j = 1; jset.size() < nBuf; j += 2) {
if (isRelPrime(D, j)) { jset.push_back(j); }
}
return jset;
}
u32 Pm1Plan::primeAfter(u32 b) const {
while (!primeBits[++b]);
return b;
}
u32 Pm1Plan::primeBefore(u32 b) const {
while (!primeBits[--b]);
return b;
}
// Returns the prime hit by "a", or 0.
u32 Pm1Plan::hit(const vector<bool>& primes, u32 a) {
u32 r = 0;
return primes[r=a]
|| primes[r=reduce(a)]
|| primes[r=reduce<13>(a)]
|| primes[r=reduce<17>(a)]
|| primes[r=reduce<19>(a)]
|| primes[r=reduce<23>(a)]
|| primes[r=reduce<29>(a)]
#if 0
|| primes[r=reduce<31>(a)]
|| primes[r=reduce<37>(a)]
|| primes[r=reduce<41>(a)]
|| primes[r=reduce<43>(a)]
#endif
? r : 0;
}
// The largest block that covers "b" (the initial block).
u32 Pm1Plan::lowerBlock(u32 b) const { return (b + jset.back()) / D; }
// The smallest block that cover "b" (the final block).
u32 Pm1Plan::upperBlock(u32 b) const { return (b - jset.back()) / D + 1; }
template<typename Fun>
void Pm1Plan::scan(const vector<bool>& primes, u32 beginBlock, vector<BitBlock>& selected, Fun fun) {
for (u32 block = selected.size() - 1; block >= beginBlock; --block) {
BitBlock& blockBits = selected[block];
const u32 base = block * D;
for (u32 pos = 0, end = jset.size(); pos < end; ++pos) {
u32 j = jset[pos];
u32 a = base - j;
u32 b = base + j;
u32 p1 = hit(primes, a);
u32 p2 = hit(primes, b);
if (fun(p1, p2)) {
assert(!blockBits[pos]);
blockBits[pos] = true;
}
}
}
}
pair<u32, vector<Pm1Plan::BitBlock>> Pm1Plan::makePlan() {
// In the unlikely case that either B1 or B2 is prime:
// B1 was included in P1, so excluded in P2.
// B2 is included in P2.
u32 lastPrime = primeBefore(B2 + 1);
u32 lastBlock = upperBlock(lastPrime);
u32 lastCovered = lastBlock * D + jset.back();
// log("last %u %u %u %u %u\n", lastPrime, lastBlock, lastCovered, u32(primeBits.size()), jset.back());
assert(lastCovered < primeBits.size());
// All primes <= cutValue can be transposed by mulitplying with firstMissingFactor.
u32 cutValue = lastCovered / firstMissingFactor(D);
u32 startValue = max(primeAfter(B1), cutValue + 1);
u32 beginBlock = lowerBlock(startValue);
u32 endBlock = lastBlock + 1;
assert(beginBlock < endBlock);
auto primes = primeBits; // use a copy in which we'll clear the primes as we cover them.
const u32 nPrimes = sum(primes);
vector<BitBlock> selected(endBlock);
u32 nPair = 0, nSingle = 0;
for (int rep = 0; rep < 4; ++rep) {
vector<bool> oneHit(primes.size()), twoHit(primes.size());
scan(primes, beginBlock, selected, [&oneHit, &twoHit](u32 p1, u32 p2) {
if (p1 && p2) {
assert(p1 != p2);
if (oneHit[p1]) {
twoHit.at(p1) = true;
} else {
oneHit.at(p1) = true;
}
if (oneHit[p2]) {
twoHit.at(p2) = true;
} else {
oneHit.at(p2) = true;
}
}
return false;
});
scan(primes, beginBlock, selected, [&primes, &twoHit, &nPair](u32 p1, u32 p2) {
if (p1 && p2 && (!twoHit[p1] || !twoHit[p2])) {
++nPair;
primes[p1] = primes[p2] = false;
return true;
} else {
return false;
}
});
}
scan(primes, beginBlock, selected, [&primes, &nPair](u32 p1, u32 p2) {
if (p1 && p2) {
++nPair;
primes[p1] = primes[p2] = false;
return true;
} else {
return false;
}
});
scan(primes, beginBlock, selected, [&primes, &nSingle](u32 p1, u32 p2) {
if (p1 || p2) {
assert(!(p1 && p2));
++nSingle;
primes[p1] = primes[p2] = false;
return true;
} else {
return false;
}
});
assert(sum(primes) == 0); // all primes are covered.
assert(nPair * 2 + nSingle == nPrimes);
u32 nBlocks = endBlock - beginBlock;
// The block transition cost is approximated as 2 MULs.
float cost = nPair + nSingle + 2 * nBlocks;
float percentPaired = 100 * (1 - nSingle / float(nPrimes));
u32 firstPrime = primeAfter(B1);
log("D=%u: %u primes in [%u, %u]: cost %.2fM (pair: %u, single: %u, (%.0f%% paired), blocks: %u)\n",
D, nPrimes, firstPrime, lastPrime,
cost * (1.0f / 1'000'000),
nPair, nSingle, percentPaired, nBlocks);
return {beginBlock, selected};
}