segmented_sieve.hpp¶
分段枚举不超过 \(n\) 的所有素数并按序回调；适合 \(n\) 较大但只需区间块内存。
\[
\displaystyle \pi(n)=\#\{p\le n: p\text{ prime}\}
\]
Complexity: Time: \(O(n \log \log n + \lceil n / \mathrm{bs}\rceil \pi (\sqrt{n}))\). Space: \(O(\sqrt{n} + \mathrm{bs})\).
AC 记录：enumerate_primes。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(n log log n + ceil(n/bs) pi(sqrt(n))).
/// Space: O(sqrt(n) + bs).
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <vector>
namespace noya {
namespace segmented_sieve_internal {
inline std::uint64_t integer_sqrt(std::uint64_t val) {
using u128 = unsigned __int128;
std::uint64_t rt = std::uint64_t(std::sqrt(static_cast<long double>(val)));
while (u128(rt + 1) * (rt + 1) <= val) {
rt++;
}
while (u128(rt) * rt > val) {
rt--;
}
return rt;
}
// Offsets count odd values. Return the next multiple's offset in the next block.
inline std::uint64_t strike(std::vector<unsigned char> &com, std::size_t cnt,
std::uint64_t p, std::uint64_t off) {
for (; off < cnt; off += p) com[off] = true;
return off - cnt;
}
} // namespace segmented_sieve_internal
/// @brief Call callback(index, prime) for every prime not exceeding n and
/// return pi(n). The callback is invoked in increasing prime order.
/// bs counts odd values per block; zero is treated as one.
template <class Callback>
std::uint64_t for_each_prime(std::uint64_t n, Callback &&cal,
std::size_t bs = 1U << 20) {
if (n < 2) {
return 0;
}
bs = std::max<std::size_t>(bs, 1);
std::uint64_t lim = segmented_sieve_internal::integer_sqrt(n);
std::vector<std::uint64_t> ps;
{
std::vector<unsigned char> cmp(lim + 1);
for (std::uint64_t val = 2; val <= lim; val++) {
if (cmp[val]) continue;
if (val > 2) ps.push_back(val);
if (val <= lim / val) {
for (std::uint64_t mul = val * val; mul <= lim; mul += val)
cmp[mul] = true;
}
}
}
std::uint64_t np = 0;
cal(np++, std::uint64_t(2));
std::vector<std::uint64_t> nxt(ps.size());
std::size_t active = 0;
std::vector<unsigned char> com(std::min<std::uint64_t>(bs, (n - 1) / 2));
for (std::uint64_t l = 3; l <= n;) {
std::uint64_t ava = (n - l) / 2 + 1;
std::size_t cnt =
std::size_t(std::min<std::uint64_t>(ava, std::uint64_t(bs)));
std::uint64_t r = l + 2 * (cnt - 1);
std::fill_n(com.begin(), cnt, 0);
while (active < ps.size() && ps[active] <= r / ps[active]) {
// A prime becomes active in the block containing its square.
nxt[active] = (ps[active] * ps[active] - l) / 2;
++active;
}
for (std::size_t j = 0; j < active; ++j)
nxt[j] = segmented_sieve_internal::strike(com, cnt, ps[j], nxt[j]);
for (std::size_t off = 0; off < cnt; off++) {
if (!com[off]) {
cal(np++, l + 2 * off);
}
}
if (r == n || r + 2 < r) {
break;
}
l = r + 2;
}
return np;
}
} // namespace noya
// NOYA_SEGMENTED_SIEVE_HPP
#ifndef NOYA_SEGMENTED_SIEVE_HPP
#define NOYA_SEGMENTED_SIEVE_HPP 1
/// @complexity Time: O(n log log n + ceil(n/bs) pi(sqrt(n))).
/// Space: O(sqrt(n) + bs).
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <vector>
namespace noya {
namespace segmented_sieve_internal {
inline std::uint64_t integer_sqrt(std::uint64_t val) {
using u128 = unsigned __int128;
std::uint64_t rt = std::uint64_t(std::sqrt(static_cast<long double>(val)));
while (u128(rt + 1) * (rt + 1) <= val) {
rt++;
}
while (u128(rt) * rt > val) {
rt--;
}
return rt;
}
// Offsets count odd values. Return the next multiple's offset in the next block.
inline std::uint64_t strike(std::vector<unsigned char> &com, std::size_t cnt,
std::uint64_t p, std::uint64_t off) {
for (; off < cnt; off += p) com[off] = true;
return off - cnt;
}
} // namespace segmented_sieve_internal
/// @brief Call callback(index, prime) for every prime not exceeding n and
/// return pi(n). The callback is invoked in increasing prime order.
/// bs counts odd values per block; zero is treated as one.
template <class Callback>
std::uint64_t for_each_prime(std::uint64_t n, Callback &&cal,
std::size_t bs = 1U << 20) {
if (n < 2) {
return 0;
}
bs = std::max<std::size_t>(bs, 1);
std::uint64_t lim = segmented_sieve_internal::integer_sqrt(n);
std::vector<std::uint64_t> ps;
{
std::vector<unsigned char> cmp(lim + 1);
for (std::uint64_t val = 2; val <= lim; val++) {
if (cmp[val]) continue;
if (val > 2) ps.push_back(val);
if (val <= lim / val) {
for (std::uint64_t mul = val * val; mul <= lim; mul += val)
cmp[mul] = true;
}
}
}
std::uint64_t np = 0;
cal(np++, std::uint64_t(2));
std::vector<std::uint64_t> nxt(ps.size());
std::size_t active = 0;
std::vector<unsigned char> com(std::min<std::uint64_t>(bs, (n - 1) / 2));
for (std::uint64_t l = 3; l <= n;) {
std::uint64_t ava = (n - l) / 2 + 1;
std::size_t cnt =
std::size_t(std::min<std::uint64_t>(ava, std::uint64_t(bs)));
std::uint64_t r = l + 2 * (cnt - 1);
std::fill_n(com.begin(), cnt, 0);
while (active < ps.size() && ps[active] <= r / ps[active]) {
// A prime becomes active in the block containing its square.
nxt[active] = (ps[active] * ps[active] - l) / 2;
++active;
}
for (std::size_t j = 0; j < active; ++j)
nxt[j] = segmented_sieve_internal::strike(com, cnt, ps[j], nxt[j]);
for (std::size_t off = 0; off < cnt; off++) {
if (!com[off]) {
cal(np++, l + 2 * off);
}
}
if (r == n || r + 2 < r) {
break;
}
l = r + 2;
}
return np;
}
} // namespace noya
#endif // NOYA_SEGMENTED_SIEVE_HPP
/// @complexity Time: O(n log log n + ceil(n/bs) pi(sqrt(n))).
/// Space: O(sqrt(n) + bs).
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <vector>
namespace noya {
namespace segmented_sieve_internal {
inline std::uint64_t integer_sqrt(std::uint64_t val) {
using u128 = unsigned __int128;
std::uint64_t rt = std::uint64_t(std::sqrt(static_cast<long double>(val)));
while (u128(rt + 1) * (rt + 1) <= val) {
rt++;
}
while (u128(rt) * rt > val) {
rt--;
}
return rt;
}
// Offsets count odd values. Return the next multiple's offset in the next block.
inline std::uint64_t strike(std::vector<unsigned char> &com, std::size_t cnt,
std::uint64_t p, std::uint64_t off) {
for (; off < cnt; off += p) com[off] = true;
return off - cnt;
}
} // namespace segmented_sieve_internal
/// @brief Call callback(index, prime) for every prime not exceeding n and
/// return pi(n). The callback is invoked in increasing prime order.
/// bs counts odd values per block; zero is treated as one.
template <class Callback>
std::uint64_t for_each_prime(std::uint64_t n, Callback &&cal,
std::size_t bs = 1U << 20) {
if (n < 2) {
return 0;
}
bs = std::max<std::size_t>(bs, 1);
std::uint64_t lim = segmented_sieve_internal::integer_sqrt(n);
std::vector<std::uint64_t> ps;
{
std::vector<unsigned char> cmp(lim + 1);
for (std::uint64_t val = 2; val <= lim; val++) {
if (cmp[val]) continue;
if (val > 2) ps.push_back(val);
if (val <= lim / val) {
for (std::uint64_t mul = val * val; mul <= lim; mul += val)
cmp[mul] = true;
}
}
}
std::uint64_t np = 0;
cal(np++, std::uint64_t(2));
std::vector<std::uint64_t> nxt(ps.size());
std::size_t active = 0;
std::vector<unsigned char> com(std::min<std::uint64_t>(bs, (n - 1) / 2));
for (std::uint64_t l = 3; l <= n;) {
std::uint64_t ava = (n - l) / 2 + 1;
std::size_t cnt =
std::size_t(std::min<std::uint64_t>(ava, std::uint64_t(bs)));
std::uint64_t r = l + 2 * (cnt - 1);
std::fill_n(com.begin(), cnt, 0);
while (active < ps.size() && ps[active] <= r / ps[active]) {
// A prime becomes active in the block containing its square.
nxt[active] = (ps[active] * ps[active] - l) / 2;
++active;
}
for (std::size_t j = 0; j < active; ++j)
nxt[j] = segmented_sieve_internal::strike(com, cnt, ps[j], nxt[j]);
for (std::size_t off = 0; off < cnt; off++) {
if (!com[off]) {
cal(np++, l + 2 * off);
}
}
if (r == n || r + 2 < r) {
break;
}
l = r + 2;
}
return np;
}
} // namespace noya
// NOYA_SEGMENTED_SIEVE_HPP