diff --git a/content/maths/number-theory/fast-sieve.cpp b/content/maths/number-theory/fast-sieve.cpp new file mode 100644 index 0000000..b5b51c3 --- /dev/null +++ b/content/maths/number-theory/fast-sieve.cpp @@ -0,0 +1,94 @@ +/* + *Author:* min\_25 + *Description:* Fast prime sieve up to $N$. + Returns a sorted vector of all primes $<= N$. + Uses a segmented sieve with wheel factorization (mod 30). + *Time:* ~0.5s for $N = 10^9$. + *Status:* Tested on Yosupo (enumerate\_primes) +*/ +vector sieve(int N, int Q = 17, int L = 1 << 15) { + static const int rs[] = {1, 7, 11, 13, 17, 19, 23, 29}; + struct P { + P(int p) : p(p) {} + int p, pos[8]; + }; + + auto prime_count = [](int N) -> int { + return N > 60184 ? N / (log(N) - 1.1) : max(1., N / (log(N) - 1.11)) + 1; + }; + + int v = sqrt(N), vv = sqrt(v); + vector isp(v + 1, true); + for (int i = 2; i <= vv; i++) + if (isp[i]) + for (int j = i * i; j <= v; j += i) + isp[j] = false; + + int rsize = prime_count(N + 30); + vector primes = {2, 3, 5}; + int psize = 3; + primes.resize(rsize); + + vector

sprimes; + size_t pbeg = 0; + int prod = 1; + + for (int p = 7; p <= v; p++) { + if (!isp[p]) + continue; + if (p <= Q) + prod *= p, pbeg++, primes[psize++] = p; + auto pp = P(p); + for (int t = 0; t < 8; t++) { + int j = (p <= Q) ? p : p * p; + while (j % 30 != rs[t]) + j += p << 1; + pp.pos[t] = j / 30; + } + sprimes.push_back(pp); + } + + vector pre(prod, 0xFF); + for (size_t pi = 0; pi < pbeg; pi++) { + auto pp = sprimes[pi]; + int p = pp.p; + for (int t = 0; t < 8; t++) { + unsigned char m = ~(1 << t); + for (int i = pp.pos[t]; i < prod; i += p) + pre[i] &= m; + } + } + + int block_size = (L + prod - 1) / prod * prod; + vector block(block_size); + unsigned char *pblock = block.data(); + int M = (N + 29) / 30; + + for (int beg = 0; beg < M; beg += block_size, pblock -= block_size) { + int end = min(M, beg + block_size); + for (int i = beg; i < end; i += prod) + copy(pre.begin(), pre.end(), pblock + i); + if (beg == 0) + pblock[0] &= 0xFE; + for (size_t pi = pbeg; pi < sprimes.size(); pi++) { + auto &pp = sprimes[pi]; + int p = pp.p; + for (int t = 0; t < 8; t++) { + int i = pp.pos[t]; + unsigned char m = ~(1 << t); + for (; i < end; i += p) + pblock[i] &= m; + pp.pos[t] = i; + } + } + for (int i = beg; i < end; i++) + for (int m = pblock[i]; m > 0; m &= m - 1) + primes[psize++] = i * 30 + rs[__builtin_ctz(m)]; + } + + assert(psize <= rsize); + while (psize > 0 && primes[psize - 1] > N) + psize--; + primes.resize(psize); + return primes; +} diff --git a/tracker-wf.yaml b/tracker-wf.yaml index 4b3ae4a..f7c354d 100644 --- a/tracker-wf.yaml +++ b/tracker-wf.yaml @@ -99,6 +99,7 @@ content: - counting-divisors.cpp - discrete-log.cpp - erathostenes-sieve.cpp + - fast-sieve.cpp - euler-phi.cpp - euler-phi-sieve.cpp - mobius-sieve.cpp diff --git a/tracker.yaml b/tracker.yaml index 0a10bcb..513385a 100644 --- a/tracker.yaml +++ b/tracker.yaml @@ -76,6 +76,7 @@ content: - divisors.cpp - discrete-log.cpp - erathostenes-sieve.cpp + - fast-sieve.cpp - euler-phi.cpp - euler-phi-sieve.cpp - mobius-sieve.cpp