CP-Algorithms Library
This documentation is automatically generated by competitive-verifier/competitive-verifier
cp-algo/math/sieve.hpp
- View this file on GitHub
- Last update: 2025-12-13 22:15:46+01:00
- Include:
#include "cp-algo/math/sieve.hpp"
Depends on
- cp-algo/structures/bit_array.hpp
- cp-algo/structures/bit_array_util.hpp
- cp-algo/util/big_alloc.hpp
- cp-algo/util/bit.hpp
- cp-algo/util/checkpoint.hpp
- cp-algo/util/simd.hpp
Verified with
Code
#ifndef CP_ALGO_MATH_SIEVE_HPP
#define CP_ALGO_MATH_SIEVE_HPP
#include "../structures/bit_array.hpp"
#include "../structures/bit_array_util.hpp"
#include "../util/bit.hpp"
#include "../util/checkpoint.hpp"
#include <cstdint>
#include <cstddef>
#include <vector>
#include <span>
#include <algorithm>
#include <cassert>
CP_ALGO_BIT_PRAGMA_PUSH
namespace cp_algo::math {
using cp_algo::structures::dynamic_bit_array;
using cp_algo::structures::bit_array;
constexpr auto wheel_primes = std::array{2u, 3u, 5u, 7u};
constexpr uint8_t period = std::ranges::fold_left(wheel_primes, 1u, std::multiplies{});
constexpr uint8_t coprime = std::ranges::fold_left(wheel_primes, 1u, [](auto a, auto b){ return a * (b - 1); });
constexpr auto coprime_wheel = [](auto x) {
return std::ranges::all_of(wheel_primes, [x](auto p){ return x % p; });
};
// Residues coprime to period
constexpr auto res_wheel = []() {
std::array<uint8_t, coprime> res;
int idx = 0;
for(uint8_t i = 1; i < period; i += 2) {
if (coprime_wheel(i)) {
res[idx++] = i;
}
}
return res;
}();
// Maps residue mod period to pre-upper_bound index in res_wheel
constexpr auto state_wheel = []() {
std::array<uint8_t, period> state;
uint8_t idx = 0;
for(uint8_t i = 0; i < period; i++) {
state[i] = idx;
idx += coprime_wheel(i);
}
return state;
}();
// Add to reach next coprime residue
constexpr auto add_wheel = []() {
std::array<uint8_t, period> add;
for(uint8_t i = 0; i < period; i++) {
add[i] = 1;
while (!coprime_wheel(i + add[i])) {
add[i]++;
}
}
return add;
}();
constexpr auto gap_wheel = []() {
std::array<uint8_t, coprime> gap;
for(uint8_t i = 0; i < coprime; i++) {
gap[i] = add_wheel[res_wheel[i]];
}
return gap;
}();
// Convert value to ordinal (index in compressed bit array)
constexpr uint32_t to_ord(uint32_t x) {
return (x / period) * coprime + state_wheel[x % period];
}
// Convert ordinal to value
constexpr uint32_t to_val(uint32_t x) {
return (x / coprime) * period + res_wheel[x % coprime];
}
constexpr size_t sqrt_threshold = 1 << 16;
constexpr auto sqrt_prime_bits = []() {
const int size = sqrt_threshold / 2;
bit_array<size> prime;
prime.set_all();
prime.reset(to_ord(1));
for(uint32_t i = res_wheel[1]; to_ord(i * i) < size; i += add_wheel[i % period]) {
if (prime[to_ord(i)]) {
for(uint32_t k = i; to_ord(i * k) < size; k += add_wheel[k % period]) {
prime.reset(to_ord(i * k));
}
}
}
return prime;
}();
constexpr size_t num_primes = []() {
size_t cnt = 0;
for(uint32_t i = res_wheel[1]; i < sqrt_threshold; i += add_wheel[i % period]) {
cnt += sqrt_prime_bits[to_ord(i)];
}
return cnt;
}();
constexpr auto sqrt_primes = []() {
std::array<uint32_t, num_primes> primes;
size_t j = 0;
for(uint32_t i = res_wheel[1]; i < sqrt_threshold; i += add_wheel[i % period]) {
if (sqrt_prime_bits[to_ord(i)]) {
primes[j++] = i;
}
}
return primes;
}();
struct wheel_t {
dynamic_bit_array mask;
uint32_t product;
};
constexpr auto make_wheel(big_vector<uint32_t> primes, uint32_t product) {
assert(product % period == 0 && product / period * coprime % dynamic_bit_array::width == 0);
wheel_t wheel;
wheel.product = product;
wheel.mask.resize(product / period * coprime);
wheel.mask.set_all();
for(auto p: primes) {
for (uint32_t k = 1; p * k < product; k += add_wheel[k % period]) {
wheel.mask.reset(to_ord(p * k));
}
}
return wheel;
}
constexpr void sieve_dense(auto &prime, uint32_t l, uint32_t r, wheel_t const& wheel) {
if (l >= r) return;
const uint32_t width = (uint32_t)dynamic_bit_array::width; // 64
uint32_t wl = l / width;
uint32_t wr = (r + width - 1) / width;
uint32_t N = (uint32_t)wheel.mask.words;
auto loop = [&](uint32_t i, uint32_t block) {
auto p_ptr = std::assume_aligned<32>(&prime.word(i));
auto m_ptr = std::assume_aligned<32>(&wheel.mask.word(0));
#pragma GCC unroll coprime
for (uint32_t j = 0; j < block; j++) {
p_ptr[j] &= m_ptr[j];
}
};
while (wl + N <= wr) {
loop(wl, N);
wl += N;
}
loop(wl, wr - wl);
}
template <class BitArray>
constexpr auto sieve_wheel(BitArray& prime, uint32_t l, uint32_t r, size_t i, int state) {
static const auto ord_step = []() {
big_vector<std::array<uint32_t, 2 * coprime>> ord_steps(num_primes);
for (uint32_t i = 0; i < size(sqrt_primes); i++) {
auto p = sqrt_primes[i];
auto &ords = ord_steps[i];
auto last = to_ord(p);
for(uint32_t j = 0; j < coprime; j++) {
auto next = to_ord(p * (res_wheel[j] + gap_wheel[j]));
ords[j] = ords[j + coprime] = next - last;
last = next;
}
}
return ord_steps;
}();
auto &ords = ord_step[i];
auto advance = [&]() {
prime.reset(std::exchange(l, l + ords[state++]));
};
uint32_t p = sqrt_primes[i];
while (l + p * coprime <= r) {
#pragma GCC unroll coprime
for (size_t j = 0; j < coprime; j++) {
advance();
}
state -= coprime;
}
while (l < r) {
advance();
}
state = state >= coprime ? state - coprime : state;
return std::pair{l, state};
}
// Primes smaller or equal than N
constexpr dynamic_bit_array sieve_wheel(uint32_t N) {
N++;
dynamic_bit_array prime(to_ord(N));
prime.set_all();
static const auto [wheels, medium_primes_begin] = []() {
constexpr size_t max_wheel_size = 1 << 20;
uint32_t product = period * dynamic_bit_array::width >> (size(wheel_primes) - 2);
big_vector<uint32_t> current;
big_vector<wheel_t> wheels;
for(size_t i = 0; i < size(sqrt_primes); i++) {
uint32_t p = sqrt_primes[i];
if (product * p > max_wheel_size) {
wheels.push_back(make_wheel(current, product));
current = {p};
product = (period * dynamic_bit_array::width >> (size(wheel_primes) - 2)) * p;
if (product > max_wheel_size) {
checkpoint("make wheels");
return std::pair{wheels, i};
}
} else {
current.push_back(p);
product *= p;
}
}
assert(false);
}();
static constexpr uint32_t dense_block = 1 << 25;
for(uint32_t start = 0; start < N; start += dense_block) {
uint32_t r = std::min(start + dense_block, N);
for(auto const& wheel: wheels) {
auto l = start / wheel.product * wheel.product;
sieve_dense(prime, to_ord(l), to_ord(r), wheel);
}
}
checkpoint("dense sieve");
static constexpr uint32_t sparse_block = 1 << 22;
auto [pos, state] = []() {
big_vector<uint32_t> pos(num_primes);
big_vector<uint8_t> state(num_primes);
for (auto [i, p]: sqrt_primes | std::views::enumerate) {
pos[i] = to_ord(p * p);
state[i] = state_wheel[p % period];
}
return std::pair{pos, state};
}();
for(uint32_t start = 0; start < N; start += sparse_block) {
uint32_t r = to_ord(std::min(start + sparse_block, N));
for(size_t i = medium_primes_begin; i < size(sqrt_primes); i++) {
if(state[i] >= r) break;
std::tie(pos[i], state[i]) = sieve_wheel(prime, pos[i], r, i, state[i]);
}
}
checkpoint("sparse sieve");
for(size_t i = 0; i < std::min(prime.words, sqrt_prime_bits.words); i++) {
prime.word(i) = sqrt_prime_bits.word(i);
}
return prime;
}
}
#pragma GCC pop_options
#endif // CP_ALGO_MATH_SIEVE_HPP
#line 1 "cp-algo/math/sieve.hpp"
#line 1 "cp-algo/structures/bit_array.hpp"
#line 1 "cp-algo/util/bit.hpp"
#line 1 "cp-algo/util/simd.hpp"
#include <experimental/simd>
#include <cstdint>
#include <cstddef>
#include <memory>
#if defined(__x86_64__) && !defined(CP_ALGO_DISABLE_AVX2)
#define CP_ALGO_SIMD_AVX2_TARGET _Pragma("GCC target(\"avx2\")")
#else
#define CP_ALGO_SIMD_AVX2_TARGET
#endif
#define CP_ALGO_SIMD_PRAGMA_PUSH \
_Pragma("GCC push_options") \
CP_ALGO_SIMD_AVX2_TARGET
CP_ALGO_SIMD_PRAGMA_PUSH
namespace cp_algo {
template<typename T, size_t len>
using simd [[gnu::vector_size(len * sizeof(T))]] = T;
using i64x4 = simd<int64_t, 4>;
using u64x4 = simd<uint64_t, 4>;
using u32x8 = simd<uint32_t, 8>;
using i32x4 = simd<int32_t, 4>;
using u32x4 = simd<uint32_t, 4>;
using i16x4 = simd<int16_t, 4>;
using u8x32 = simd<uint8_t, 32>;
using dx4 = simd<double, 4>;
dx4 abs(dx4 a) {
return dx4{
std::abs(a[0]),
std::abs(a[1]),
std::abs(a[2]),
std::abs(a[3])
};
}
// https://stackoverflow.com/a/77376595
// works for ints in (-2^51, 2^51)
static constexpr dx4 magic = dx4() + (3ULL << 51);
inline i64x4 lround(dx4 x) {
return i64x4(x + magic) - i64x4(magic);
}
inline dx4 to_double(i64x4 x) {
return dx4(x + i64x4(magic)) - magic;
}
inline dx4 round(dx4 a) {
return dx4{
std::nearbyint(a[0]),
std::nearbyint(a[1]),
std::nearbyint(a[2]),
std::nearbyint(a[3])
};
}
inline u64x4 low32(u64x4 x) {
return x & uint32_t(-1);
}
inline auto swap_bytes(auto x) {
return decltype(x)(__builtin_shufflevector(u32x8(x), u32x8(x), 1, 0, 3, 2, 5, 4, 7, 6));
}
inline u64x4 montgomery_reduce(u64x4 x, uint32_t mod, uint32_t imod) {
#ifdef __AVX2__
auto x_ninv = u64x4(_mm256_mul_epu32(__m256i(x), __m256i() + imod));
x += u64x4(_mm256_mul_epu32(__m256i(x_ninv), __m256i() + mod));
#else
auto x_ninv = u64x4(u32x8(low32(x)) * imod);
x += x_ninv * uint64_t(mod);
#endif
return swap_bytes(x);
}
inline u64x4 montgomery_mul(u64x4 x, u64x4 y, uint32_t mod, uint32_t imod) {
#ifdef __AVX2__
return montgomery_reduce(u64x4(_mm256_mul_epu32(__m256i(x), __m256i(y))), mod, imod);
#else
return montgomery_reduce(x * y, mod, imod);
#endif
}
inline u32x8 montgomery_mul(u32x8 x, u32x8 y, uint32_t mod, uint32_t imod) {
return u32x8(montgomery_mul(u64x4(x), u64x4(y), mod, imod)) |
u32x8(swap_bytes(montgomery_mul(u64x4(swap_bytes(x)), u64x4(swap_bytes(y)), mod, imod)));
}
inline dx4 rotate_right(dx4 x) {
static constexpr u64x4 shuffler = {3, 0, 1, 2};
return __builtin_shuffle(x, shuffler);
}
template<std::size_t Align = 32>
inline bool is_aligned(const auto* p) noexcept {
return (reinterpret_cast<std::uintptr_t>(p) % Align) == 0;
}
template<class Target>
inline Target& vector_cast(auto &&p) {
return *reinterpret_cast<Target*>(std::assume_aligned<alignof(Target)>(&p));
}
}
#pragma GCC pop_options
#line 6 "cp-algo/util/bit.hpp"
#include <array>
#include <bit>
#if defined(__x86_64__) && !defined(CP_ALGO_DISABLE_AVX2)
#define CP_ALGO_BIT_OPS_TARGET _Pragma("GCC target(\"avx2,bmi,bmi2,lzcnt,popcnt\")")
#else
#define CP_ALGO_BIT_OPS_TARGET _Pragma("GCC target(\"bmi,bmi2,lzcnt,popcnt\")")
#endif
#define CP_ALGO_BIT_PRAGMA_PUSH \
_Pragma("GCC push_options") \
CP_ALGO_BIT_OPS_TARGET
CP_ALGO_BIT_PRAGMA_PUSH
namespace cp_algo {
template<typename Uint>
constexpr size_t bit_width = sizeof(Uint) * 8;
// n < 64
uint64_t mask(size_t n) {
return (1ULL << n) - 1;
}
size_t order_of_bit(auto x, size_t k) {
return k ? std::popcount(x << (bit_width<decltype(x)> - k)) : 0;
}
inline size_t kth_set_bit(uint64_t x, size_t k) {
return std::countr_zero(_pdep_u64(1ULL << k, x));
}
template<int fl = 0>
void with_bit_floor(size_t n, auto &&callback) {
if constexpr (fl >= 63) {
return;
} else if (n >> (fl + 1)) {
with_bit_floor<fl + 1>(n, callback);
} else {
callback.template operator()<1ULL << fl>();
}
}
void with_bit_ceil(size_t n, auto &&callback) {
with_bit_floor(n, [&]<size_t N>() {
if(N == n) {
callback.template operator()<N>();
} else {
callback.template operator()<N << 1>();
}
});
}
inline uint32_t read_bits(char const* p) {
return _mm256_movemask_epi8(__m256i(vector_cast<u8x32 const>(p[0]) + (127 - '0')));
}
inline uint64_t read_bits64(char const* p) {
return read_bits(p) | (uint64_t(read_bits(p + 32)) << 32);
}
inline void write_bits(char *p, uint32_t bits) {
static constexpr u8x32 shuffler = {
0, 0, 0, 0, 0, 0, 0, 0,
1, 1, 1, 1, 1, 1, 1, 1,
2, 2, 2, 2, 2, 2, 2, 2,
3, 3, 3, 3, 3, 3, 3, 3
};
auto shuffled = u8x32(_mm256_shuffle_epi8(__m256i() + bits, __m256i(shuffler)));
static constexpr u8x32 mask = {
1, 2, 4, 8, 16, 32, 64, 128,
1, 2, 4, 8, 16, 32, 64, 128,
1, 2, 4, 8, 16, 32, 64, 128,
1, 2, 4, 8, 16, 32, 64, 128
};
for(int z = 0; z < 32; z++) {
p[z] = shuffled[z] & mask[z] ? '1' : '0';
}
}
inline void write_bits64(char *p, uint64_t bits) {
write_bits(p, uint32_t(bits));
write_bits(p + 32, uint32_t(bits >> 32));
}
}
#pragma GCC pop_options
#line 1 "cp-algo/util/big_alloc.hpp"
#include <set>
#include <map>
#include <deque>
#include <stack>
#include <queue>
#include <vector>
#include <string>
#line 12 "cp-algo/util/big_alloc.hpp"
#include <iostream>
#include <generator>
#include <forward_list>
// Single macro to detect POSIX platforms (Linux, Unix, macOS)
#if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
# define CP_ALGO_USE_MMAP 1
# include <sys/mman.h>
#else
# define CP_ALGO_USE_MMAP 0
#endif
namespace cp_algo {
template <typename T, size_t Align = 32>
class big_alloc {
static_assert( Align >= alignof(void*), "Align must be at least pointer-size");
static_assert(std::popcount(Align) == 1, "Align must be a power of two");
public:
using value_type = T;
template <class U> struct rebind { using other = big_alloc<U, Align>; };
constexpr bool operator==(const big_alloc&) const = default;
constexpr bool operator!=(const big_alloc&) const = default;
big_alloc() noexcept = default;
template <typename U, std::size_t A>
big_alloc(const big_alloc<U, A>&) noexcept {}
[[nodiscard]] T* allocate(std::size_t n) {
std::size_t padded = round_up(n * sizeof(T));
std::size_t align = std::max<std::size_t>(alignof(T), Align);
#if CP_ALGO_USE_MMAP
if (padded >= MEGABYTE) {
void* raw = mmap(nullptr, padded,
PROT_READ | PROT_WRITE,
MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
madvise(raw, padded, MADV_HUGEPAGE);
madvise(raw, padded, MADV_POPULATE_WRITE);
return static_cast<T*>(raw);
}
#endif
return static_cast<T*>(::operator new(padded, std::align_val_t(align)));
}
void deallocate(T* p, std::size_t n) noexcept {
if (!p) return;
std::size_t padded = round_up(n * sizeof(T));
std::size_t align = std::max<std::size_t>(alignof(T), Align);
#if CP_ALGO_USE_MMAP
if (padded >= MEGABYTE) { munmap(p, padded); return; }
#endif
::operator delete(p, padded, std::align_val_t(align));
}
private:
static constexpr std::size_t MEGABYTE = 1 << 20;
static constexpr std::size_t round_up(std::size_t x) noexcept {
return (x + Align - 1) / Align * Align;
}
};
template<typename T> using big_vector = std::vector<T, big_alloc<T>>;
template<typename T> using big_basic_string = std::basic_string<T, std::char_traits<T>, big_alloc<T>>;
template<typename T> using big_deque = std::deque<T, big_alloc<T>>;
template<typename T> using big_stack = std::stack<T, big_deque<T>>;
template<typename T> using big_queue = std::queue<T, big_deque<T>>;
template<typename T> using big_priority_queue = std::priority_queue<T, big_vector<T>>;
template<typename T> using big_forward_list = std::forward_list<T, big_alloc<T>>;
using big_string = big_basic_string<char>;
template<typename Key, typename Value, typename Compare = std::less<Key>>
using big_map = std::map<Key, Value, Compare, big_alloc<std::pair<const Key, Value>>>;
template<typename T, typename Compare = std::less<T>>
using big_multiset = std::multiset<T, Compare, big_alloc<T>>;
template<typename T, typename Compare = std::less<T>>
using big_set = std::set<T, Compare, big_alloc<T>>;
template<typename Ref, typename V = void>
using big_generator = std::generator<Ref, V, big_alloc<std::byte>>;
}
// Deduction guide to make elements_of with big_generator default to big_alloc
namespace std::ranges {
template<typename Ref, typename V>
elements_of(cp_algo::big_generator<Ref, V>&&) -> elements_of<cp_algo::big_generator<Ref, V>&&, cp_algo::big_alloc<std::byte>>;
}
#line 5 "cp-algo/structures/bit_array.hpp"
#include <cassert>
CP_ALGO_SIMD_PRAGMA_PUSH
namespace cp_algo::structures {
template<typename C>
concept Resizable = requires(C& c, std::size_t n) { c.resize(n); };
template<class Cont>
struct _bit_array {
using word_t = typename Cont::value_type;
static constexpr size_t width = bit_width<word_t>;
size_t words, n;
alignas(32) Cont data;
constexpr void resize(size_t N) {
n = N;
words = (n + width - 1) / width;
if constexpr (Resizable<Cont>) {
data.resize(words);
} else {
assert(std::size(data) >= words);
}
}
constexpr _bit_array(): data() {
if constexpr (!Resizable<Cont>) {
resize(std::size(data) * width);
} else {
resize(0);
}
}
constexpr _bit_array(size_t N): data() {
resize(N);
}
constexpr auto&& word(this auto&& self, size_t x) {
return self.data[x];
}
constexpr void set_all(word_t val = -1) {
for(size_t i = 0; i < words; i++) {
data[i] = val;
}
}
constexpr void reset() {
set_all(0);
}
constexpr void set(size_t x) {
word(x / width) |= 1ULL << (x % width);
}
constexpr void reset(size_t x) {
word(x / width) &= ~(1ULL << (x % width));
}
constexpr void flip(size_t x) {
word(x / width) ^= 1ULL << (x % width);
}
constexpr bool test(size_t x) const {
return (word(x / width) >> (x % width)) & 1;
}
constexpr bool operator[](size_t x) const {
return test(x);
}
constexpr size_t size() const {
return n;
}
auto operator <=> (_bit_array const& t) const = default;
constexpr _bit_array& xor_hint(_bit_array const& t, size_t hint) {
for(size_t i = hint / width; i < words; i++) {
data[i] ^= t.data[i];
}
return *this;
}
constexpr _bit_array& operator ^= (_bit_array const& t) {
return xor_hint(t, 0);
}
constexpr _bit_array operator ^ (_bit_array const& t) const {
return _bit_array(*this) ^= t;
}
};
template<size_t N>
using bit_array = _bit_array<std::array<uint64_t, (N + 63) / 64>>;
using dynamic_bit_array = _bit_array<big_vector<uint64_t>>;
}
#pragma GCC pop_options
#line 1 "cp-algo/structures/bit_array_util.hpp"
#line 9 "cp-algo/structures/bit_array_util.hpp"
#include <span>
CP_ALGO_BIT_PRAGMA_PUSH
namespace cp_algo::structures {
template<typename BitArray>
constexpr void from_string(BitArray& arr, std::span<char> bits) {
arr.resize(std::size(bits));
int64_t i = 0;
constexpr int width = BitArray::width;
for(; i + width <= (int64_t)std::size(bits); i += width) {
arr.word(i / width) = read_bits64(bits.data() + i);
}
for(; i < (int64_t)std::size(bits); i++) {
if(bits[i] & 1) {
arr.set(i);
}
}
}
template<typename BitArray>
constexpr big_string to_string(BitArray const& arr) {
big_string res(arr.words * BitArray::width, '0');
for(size_t i = 0; i < arr.words; i++) {
write_bits64(res.data() + i * BitArray::width, arr.word(i));
}
res.resize(arr.n);
return res;
}
template<typename BitArray>
constexpr size_t count(BitArray const& arr, size_t n) {
size_t res = 0;
for(size_t i = 0; i < n / BitArray::width; i++) {
res += std::popcount(arr.word(i));
}
if (n % BitArray::width) {
res += std::popcount(arr.word(n / BitArray::width) & mask(n % BitArray::width));
}
return res;
}
template<typename BitArray>
constexpr size_t count(BitArray const& arr) {
return count(arr, arr.n);
}
template<typename BitArray>
constexpr size_t ctz(BitArray const& arr) {
size_t res = 0;
size_t i = 0;
while(i < arr.words && arr.word(i) == 0) {
res += BitArray::width;
i++;
}
if(i < arr.words) {
res += std::countr_zero(arr.word(i));
}
return std::min(res, arr.n);
}
template<typename BitArray>
constexpr size_t skip(BitArray const& arr, size_t pos = 0, int k = 0) {
if (pos >= arr.n) return arr.n;
size_t word_idx = pos / BitArray::width;
auto w = arr.word(word_idx) >> (pos % BitArray::width);
auto popcnt = std::popcount(w);
if (popcnt > k) {
return pos + cp_algo::kth_set_bit(w, k);
}
k -= popcnt;
while (++word_idx < arr.words) {
w = arr.word(word_idx);
auto popcnt = std::popcount(w);
if (popcnt > k) [[unlikely]] {
return word_idx * BitArray::width + cp_algo::kth_set_bit(w, k);
}
k -= popcnt;
}
return arr.n;
}
}
#pragma GCC pop_options
#line 1 "cp-algo/util/checkpoint.hpp"
#line 5 "cp-algo/util/checkpoint.hpp"
#include <chrono>
#line 8 "cp-algo/util/checkpoint.hpp"
namespace cp_algo {
#ifdef CP_ALGO_CHECKPOINT
big_map<big_string, double> checkpoints;
double last;
#endif
template<bool final = false>
void checkpoint([[maybe_unused]] auto const& _msg) {
#ifdef CP_ALGO_CHECKPOINT
big_string msg = _msg;
double now = (double)clock() / CLOCKS_PER_SEC;
double delta = now - last;
last = now;
if(msg.size() && !final) {
checkpoints[msg] += delta;
}
if(final) {
for(auto const& [key, value] : checkpoints) {
std::cerr << key << ": " << value * 1000 << " ms\n";
}
std::cerr << "Total: " << now * 1000 << " ms\n";
}
#endif
}
template<bool final = false>
void checkpoint() {
checkpoint<final>("");
}
}
#line 11 "cp-algo/math/sieve.hpp"
#include <algorithm>
#line 13 "cp-algo/math/sieve.hpp"
CP_ALGO_BIT_PRAGMA_PUSH
namespace cp_algo::math {
using cp_algo::structures::dynamic_bit_array;
using cp_algo::structures::bit_array;
constexpr auto wheel_primes = std::array{2u, 3u, 5u, 7u};
constexpr uint8_t period = std::ranges::fold_left(wheel_primes, 1u, std::multiplies{});
constexpr uint8_t coprime = std::ranges::fold_left(wheel_primes, 1u, [](auto a, auto b){ return a * (b - 1); });
constexpr auto coprime_wheel = [](auto x) {
return std::ranges::all_of(wheel_primes, [x](auto p){ return x % p; });
};
// Residues coprime to period
constexpr auto res_wheel = []() {
std::array<uint8_t, coprime> res;
int idx = 0;
for(uint8_t i = 1; i < period; i += 2) {
if (coprime_wheel(i)) {
res[idx++] = i;
}
}
return res;
}();
// Maps residue mod period to pre-upper_bound index in res_wheel
constexpr auto state_wheel = []() {
std::array<uint8_t, period> state;
uint8_t idx = 0;
for(uint8_t i = 0; i < period; i++) {
state[i] = idx;
idx += coprime_wheel(i);
}
return state;
}();
// Add to reach next coprime residue
constexpr auto add_wheel = []() {
std::array<uint8_t, period> add;
for(uint8_t i = 0; i < period; i++) {
add[i] = 1;
while (!coprime_wheel(i + add[i])) {
add[i]++;
}
}
return add;
}();
constexpr auto gap_wheel = []() {
std::array<uint8_t, coprime> gap;
for(uint8_t i = 0; i < coprime; i++) {
gap[i] = add_wheel[res_wheel[i]];
}
return gap;
}();
// Convert value to ordinal (index in compressed bit array)
constexpr uint32_t to_ord(uint32_t x) {
return (x / period) * coprime + state_wheel[x % period];
}
// Convert ordinal to value
constexpr uint32_t to_val(uint32_t x) {
return (x / coprime) * period + res_wheel[x % coprime];
}
constexpr size_t sqrt_threshold = 1 << 16;
constexpr auto sqrt_prime_bits = []() {
const int size = sqrt_threshold / 2;
bit_array<size> prime;
prime.set_all();
prime.reset(to_ord(1));
for(uint32_t i = res_wheel[1]; to_ord(i * i) < size; i += add_wheel[i % period]) {
if (prime[to_ord(i)]) {
for(uint32_t k = i; to_ord(i * k) < size; k += add_wheel[k % period]) {
prime.reset(to_ord(i * k));
}
}
}
return prime;
}();
constexpr size_t num_primes = []() {
size_t cnt = 0;
for(uint32_t i = res_wheel[1]; i < sqrt_threshold; i += add_wheel[i % period]) {
cnt += sqrt_prime_bits[to_ord(i)];
}
return cnt;
}();
constexpr auto sqrt_primes = []() {
std::array<uint32_t, num_primes> primes;
size_t j = 0;
for(uint32_t i = res_wheel[1]; i < sqrt_threshold; i += add_wheel[i % period]) {
if (sqrt_prime_bits[to_ord(i)]) {
primes[j++] = i;
}
}
return primes;
}();
struct wheel_t {
dynamic_bit_array mask;
uint32_t product;
};
constexpr auto make_wheel(big_vector<uint32_t> primes, uint32_t product) {
assert(product % period == 0 && product / period * coprime % dynamic_bit_array::width == 0);
wheel_t wheel;
wheel.product = product;
wheel.mask.resize(product / period * coprime);
wheel.mask.set_all();
for(auto p: primes) {
for (uint32_t k = 1; p * k < product; k += add_wheel[k % period]) {
wheel.mask.reset(to_ord(p * k));
}
}
return wheel;
}
constexpr void sieve_dense(auto &prime, uint32_t l, uint32_t r, wheel_t const& wheel) {
if (l >= r) return;
const uint32_t width = (uint32_t)dynamic_bit_array::width; // 64
uint32_t wl = l / width;
uint32_t wr = (r + width - 1) / width;
uint32_t N = (uint32_t)wheel.mask.words;
auto loop = [&](uint32_t i, uint32_t block) {
auto p_ptr = std::assume_aligned<32>(&prime.word(i));
auto m_ptr = std::assume_aligned<32>(&wheel.mask.word(0));
#pragma GCC unroll coprime
for (uint32_t j = 0; j < block; j++) {
p_ptr[j] &= m_ptr[j];
}
};
while (wl + N <= wr) {
loop(wl, N);
wl += N;
}
loop(wl, wr - wl);
}
template <class BitArray>
constexpr auto sieve_wheel(BitArray& prime, uint32_t l, uint32_t r, size_t i, int state) {
static const auto ord_step = []() {
big_vector<std::array<uint32_t, 2 * coprime>> ord_steps(num_primes);
for (uint32_t i = 0; i < size(sqrt_primes); i++) {
auto p = sqrt_primes[i];
auto &ords = ord_steps[i];
auto last = to_ord(p);
for(uint32_t j = 0; j < coprime; j++) {
auto next = to_ord(p * (res_wheel[j] + gap_wheel[j]));
ords[j] = ords[j + coprime] = next - last;
last = next;
}
}
return ord_steps;
}();
auto &ords = ord_step[i];
auto advance = [&]() {
prime.reset(std::exchange(l, l + ords[state++]));
};
uint32_t p = sqrt_primes[i];
while (l + p * coprime <= r) {
#pragma GCC unroll coprime
for (size_t j = 0; j < coprime; j++) {
advance();
}
state -= coprime;
}
while (l < r) {
advance();
}
state = state >= coprime ? state - coprime : state;
return std::pair{l, state};
}
// Primes smaller or equal than N
constexpr dynamic_bit_array sieve_wheel(uint32_t N) {
N++;
dynamic_bit_array prime(to_ord(N));
prime.set_all();
static const auto [wheels, medium_primes_begin] = []() {
constexpr size_t max_wheel_size = 1 << 20;
uint32_t product = period * dynamic_bit_array::width >> (size(wheel_primes) - 2);
big_vector<uint32_t> current;
big_vector<wheel_t> wheels;
for(size_t i = 0; i < size(sqrt_primes); i++) {
uint32_t p = sqrt_primes[i];
if (product * p > max_wheel_size) {
wheels.push_back(make_wheel(current, product));
current = {p};
product = (period * dynamic_bit_array::width >> (size(wheel_primes) - 2)) * p;
if (product > max_wheel_size) {
checkpoint("make wheels");
return std::pair{wheels, i};
}
} else {
current.push_back(p);
product *= p;
}
}
assert(false);
}();
static constexpr uint32_t dense_block = 1 << 25;
for(uint32_t start = 0; start < N; start += dense_block) {
uint32_t r = std::min(start + dense_block, N);
for(auto const& wheel: wheels) {
auto l = start / wheel.product * wheel.product;
sieve_dense(prime, to_ord(l), to_ord(r), wheel);
}
}
checkpoint("dense sieve");
static constexpr uint32_t sparse_block = 1 << 22;
auto [pos, state] = []() {
big_vector<uint32_t> pos(num_primes);
big_vector<uint8_t> state(num_primes);
for (auto [i, p]: sqrt_primes | std::views::enumerate) {
pos[i] = to_ord(p * p);
state[i] = state_wheel[p % period];
}
return std::pair{pos, state};
}();
for(uint32_t start = 0; start < N; start += sparse_block) {
uint32_t r = to_ord(std::min(start + sparse_block, N));
for(size_t i = medium_primes_begin; i < size(sqrt_primes); i++) {
if(state[i] >= r) break;
std::tie(pos[i], state[i]) = sieve_wheel(prime, pos[i], r, i, state[i]);
}
}
checkpoint("sparse sieve");
for(size_t i = 0; i < std::min(prime.words, sqrt_prime_bits.words); i++) {
prime.word(i) = sqrt_prime_bits.word(i);
}
return prime;
}
}
#pragma GCC pop_options
#ifndef CP_ALGO_MATH_SIEVE_HPP
#define CP_ALGO_MATH_SIEVE_HPP
#include "../structures/bit_array.hpp"
#include "../structures/bit_array_util.hpp"
#include "../util/bit.hpp"
#include "../util/checkpoint.hpp"
#include <cstdint>
#include <cstddef>
#include <vector>
#include <span>
#include <algorithm>
#include <cassert>
CP_ALGO_BIT_PRAGMA_PUSHnamespace cp_algo::math{using cp_algo::structures::dynamic_bit_array;using cp_algo::structures::bit_array;constexpr auto wheel_primes=std::array{2u,3u,5u,7u};constexpr uint8_t period=std::ranges::fold_left(wheel_primes,1u,std::multiplies{});constexpr uint8_t coprime=std::ranges::fold_left(wheel_primes,1u,[](auto a,auto b){return a*(b-1);});constexpr auto coprime_wheel=[](auto x){return std::ranges::all_of(wheel_primes,[x](auto p){return x%p;});};constexpr auto res_wheel=[](){std::array<uint8_t,coprime>res;int idx=0;for(uint8_t i=1;i<period;i+=2){if(coprime_wheel(i)){res[idx++]=i;}}return res;}();constexpr auto state_wheel=[](){std::array<uint8_t,period>state;uint8_t idx=0;for(uint8_t i=0;i<period;i++){state[i]=idx;idx+=coprime_wheel(i);}return state;}();constexpr auto add_wheel=[](){std::array<uint8_t,period>add;for(uint8_t i=0;i<period;i++){add[i]=1;while(!coprime_wheel(i+add[i])){add[i]++;}}return add;}();constexpr auto gap_wheel=[](){std::array<uint8_t,coprime>gap;for(uint8_t i=0;i<coprime;i++){gap[i]=add_wheel[res_wheel[i]];}return gap;}();constexpr uint32_t to_ord(uint32_t x){return(x/period)*coprime+state_wheel[x%period];}constexpr uint32_t to_val(uint32_t x){return(x/coprime)*period+res_wheel[x%coprime];}constexpr size_t sqrt_threshold=1<<16;constexpr auto sqrt_prime_bits=[](){const int size=sqrt_threshold/2;bit_array<size>prime;prime.set_all();prime.reset(to_ord(1));for(uint32_t i=res_wheel[1];to_ord(i*i)<size;i+=add_wheel[i%period]){if(prime[to_ord(i)]){for(uint32_t k=i;to_ord(i*k)<size;k+=add_wheel[k%period]){prime.reset(to_ord(i*k));}}}return prime;}();constexpr size_t num_primes=[](){size_t cnt=0;for(uint32_t i=res_wheel[1];i<sqrt_threshold;i+=add_wheel[i%period]){cnt+=sqrt_prime_bits[to_ord(i)];}return cnt;}();constexpr auto sqrt_primes=[](){std::array<uint32_t,num_primes>primes;size_t j=0;for(uint32_t i=res_wheel[1];i<sqrt_threshold;i+=add_wheel[i%period]){if(sqrt_prime_bits[to_ord(i)]){primes[j++]=i;}}return primes;}();struct wheel_t{dynamic_bit_array mask;uint32_t product;};constexpr auto make_wheel(big_vector<uint32_t>primes,uint32_t product){assert(product%period==0&&product/period*coprime%dynamic_bit_array::width==0);wheel_t wheel;wheel.product=product;wheel.mask.resize(product/period*coprime);wheel.mask.set_all();for(auto p:primes){for(uint32_t k=1;p*k<product;k+=add_wheel[k%period]){wheel.mask.reset(to_ord(p*k));}}return wheel;}constexpr void sieve_dense(auto&prime,uint32_t l,uint32_t r,wheel_t const&wheel){if(l>=r)return;const uint32_t width=(uint32_t)dynamic_bit_array::width;uint32_t wl=l/width;uint32_t wr=(r+width-1)/width;uint32_t N=(uint32_t)wheel.mask.words;auto loop=[&](uint32_t i,uint32_t block){auto p_ptr=std::assume_aligned<32>(&prime.word(i));auto m_ptr=std::assume_aligned<32>(&wheel.mask.word(0));
#pragma GCC unroll coprime
for(uint32_t j=0;j<block;j++){p_ptr[j]&=m_ptr[j];}};while(wl+N<=wr){loop(wl,N);wl+=N;}loop(wl,wr-wl);}template<class BitArray>constexpr auto sieve_wheel(BitArray&prime,uint32_t l,uint32_t r,size_t i,int state){static const auto ord_step=[](){big_vector<std::array<uint32_t,2*coprime>>ord_steps(num_primes);for(uint32_t i=0;i<size(sqrt_primes);i++){auto p=sqrt_primes[i];auto&ords=ord_steps[i];auto last=to_ord(p);for(uint32_t j=0;j<coprime;j++){auto next=to_ord(p*(res_wheel[j]+gap_wheel[j]));ords[j]=ords[j+coprime]=next-last;last=next;}}return ord_steps;}();auto&ords=ord_step[i];auto advance=[&](){prime.reset(std::exchange(l,l+ords[state++]));};uint32_t p=sqrt_primes[i];while(l+p*coprime<=r){
#pragma GCC unroll coprime
for(size_t j=0;j<coprime;j++){advance();}state-=coprime;}while(l<r){advance();}state=state>=coprime?state-coprime:state;return std::pair{l,state};}constexpr dynamic_bit_array sieve_wheel(uint32_t N){N++;dynamic_bit_array prime(to_ord(N));prime.set_all();static const auto[wheels,medium_primes_begin]=[](){constexpr size_t max_wheel_size=1<<20;uint32_t product=period*dynamic_bit_array::width>>(size(wheel_primes)-2);big_vector<uint32_t>current;big_vector<wheel_t>wheels;for(size_t i=0;i<size(sqrt_primes);i++){uint32_t p=sqrt_primes[i];if(product*p>max_wheel_size){wheels.push_back(make_wheel(current,product));current={p};product=(period*dynamic_bit_array::width>>(size(wheel_primes)-2))*p;if(product>max_wheel_size){checkpoint("make wheels");return std::pair{wheels,i};}}else{current.push_back(p);product*=p;}}assert(false);}();static constexpr uint32_t dense_block=1<<25;for(uint32_t start=0;start<N;start+=dense_block){uint32_t r=std::min(start+dense_block,N);for(auto const&wheel:wheels){auto l=start/wheel.product*wheel.product;sieve_dense(prime,to_ord(l),to_ord(r),wheel);}}checkpoint("dense sieve");static constexpr uint32_t sparse_block=1<<22;auto[pos,state]=[](){big_vector<uint32_t>pos(num_primes);big_vector<uint8_t>state(num_primes);for(auto[i,p]:sqrt_primes|std::views::enumerate){pos[i]=to_ord(p*p);state[i]=state_wheel[p%period];}return std::pair{pos,state};}();for(uint32_t start=0;start<N;start+=sparse_block){uint32_t r=to_ord(std::min(start+sparse_block,N));for(size_t i=medium_primes_begin;i<size(sqrt_primes);i++){if(state[i]>=r)break;std::tie(pos[i],state[i])=sieve_wheel(prime,pos[i],r,i,state[i]);}}checkpoint("sparse sieve");for(size_t i=0;i<std::min(prime.words,sqrt_prime_bits.words);i++){prime.word(i)=sqrt_prime_bits.word(i);}return prime;}}
#pragma GCC pop_options
#endif
#line 1 "cp-algo/math/sieve.hpp"
#line 1 "cp-algo/structures/bit_array.hpp"
#line 1 "cp-algo/util/bit.hpp"
#line 1 "cp-algo/util/simd.hpp"
#include <experimental/simd>
#include <cstdint>
#include <cstddef>
#include <memory>
#if defined(__x86_64__) && !defined(CP_ALGO_DISABLE_AVX2)
#define CP_ALGO_SIMD_AVX2_TARGET _Pragma("GCC target(\"avx2\")")
#else
#define CP_ALGO_SIMD_AVX2_TARGET
#endif
#define CP_ALGO_SIMD_PRAGMA_PUSH \
_Pragma("GCC push_options")\CP_ALGO_SIMD_AVX2_TARGETCP_ALGO_SIMD_PRAGMA_PUSHnamespace cp_algo{template<typename T,size_t len>using simd[[gnu::vector_size(len*sizeof(T))]]=T;using i64x4=simd<int64_t,4>;using u64x4=simd<uint64_t,4>;using u32x8=simd<uint32_t,8>;using i32x4=simd<int32_t,4>;using u32x4=simd<uint32_t,4>;using i16x4=simd<int16_t,4>;using u8x32=simd<uint8_t,32>;using dx4=simd<double,4>;dx4 abs(dx4 a){return dx4{std::abs(a[0]),std::abs(a[1]),std::abs(a[2]),std::abs(a[3])};}static constexpr dx4 magic=dx4()+(3ULL<<51);inline i64x4 lround(dx4 x){return i64x4(x+magic)-i64x4(magic);}inline dx4 to_double(i64x4 x){return dx4(x+i64x4(magic))-magic;}inline dx4 round(dx4 a){return dx4{std::nearbyint(a[0]),std::nearbyint(a[1]),std::nearbyint(a[2]),std::nearbyint(a[3])};}inline u64x4 low32(u64x4 x){return x&uint32_t(-1);}inline auto swap_bytes(auto x){return decltype(x)(__builtin_shufflevector(u32x8(x),u32x8(x),1,0,3,2,5,4,7,6));}inline u64x4 montgomery_reduce(u64x4 x,uint32_t mod,uint32_t imod){
#ifdef __AVX2__
auto x_ninv=u64x4(_mm256_mul_epu32(__m256i(x),__m256i()+imod));x+=u64x4(_mm256_mul_epu32(__m256i(x_ninv),__m256i()+mod));
#else
auto x_ninv=u64x4(u32x8(low32(x))*imod);x+=x_ninv*uint64_t(mod);
#endif
return swap_bytes(x);}inline u64x4 montgomery_mul(u64x4 x,u64x4 y,uint32_t mod,uint32_t imod){
#ifdef __AVX2__
return montgomery_reduce(u64x4(_mm256_mul_epu32(__m256i(x),__m256i(y))),mod,imod);
#else
return montgomery_reduce(x*y,mod,imod);
#endif
}inline u32x8 montgomery_mul(u32x8 x,u32x8 y,uint32_t mod,uint32_t imod){return u32x8(montgomery_mul(u64x4(x),u64x4(y),mod,imod))|u32x8(swap_bytes(montgomery_mul(u64x4(swap_bytes(x)),u64x4(swap_bytes(y)),mod,imod)));}inline dx4 rotate_right(dx4 x){static constexpr u64x4 shuffler={3,0,1,2};return __builtin_shuffle(x,shuffler);}template<std::size_t Align=32>inline bool is_aligned(const auto*p)noexcept{return(reinterpret_cast<std::uintptr_t>(p)%Align)==0;}template<class Target>inline Target&vector_cast(auto&&p){return*reinterpret_cast<Target*>(std::assume_aligned<alignof(Target)>(&p));}}
#pragma GCC pop_options
#line 6 "cp-algo/util/bit.hpp"
#include <array>
#include <bit>
#if defined(__x86_64__) && !defined(CP_ALGO_DISABLE_AVX2)
#define CP_ALGO_BIT_OPS_TARGET _Pragma("GCC target(\"avx2,bmi,bmi2,lzcnt,popcnt\")")
#else
#define CP_ALGO_BIT_OPS_TARGET _Pragma("GCC target(\"bmi,bmi2,lzcnt,popcnt\")")
#endif
#define CP_ALGO_BIT_PRAGMA_PUSH \
_Pragma("GCC push_options")\CP_ALGO_BIT_OPS_TARGETCP_ALGO_BIT_PRAGMA_PUSHnamespace cp_algo{template<typename Uint>constexpr size_t bit_width=sizeof(Uint)*8;uint64_t mask(size_t n){return(1ULL<<n)-1;}size_t order_of_bit(auto x,size_t k){return k?std::popcount(x<<(bit_width<decltype(x)>-k)):0;}inline size_t kth_set_bit(uint64_t x,size_t k){return std::countr_zero(_pdep_u64(1ULL<<k,x));}template<int fl=0>void with_bit_floor(size_t n,auto&&callback){if constexpr(fl>=63){return;}else if(n>>(fl+1)){with_bit_floor<fl+1>(n,callback);}else{callback.template operator()<1ULL<<fl>();}}void with_bit_ceil(size_t n,auto&&callback){with_bit_floor(n,[&]<size_t N>(){if(N==n){callback.template operator()<N>();}else{callback.template operator()<N<<1>();}});}inline uint32_t read_bits(char const*p){return _mm256_movemask_epi8(__m256i(vector_cast<u8x32 const>(p[0])+(127-'0')));}inline uint64_t read_bits64(char const*p){return read_bits(p)|(uint64_t(read_bits(p+32))<<32);}inline void write_bits(char*p,uint32_t bits){static constexpr u8x32 shuffler={0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3};auto shuffled=u8x32(_mm256_shuffle_epi8(__m256i()+bits,__m256i(shuffler)));static constexpr u8x32 mask={1,2,4,8,16,32,64,128,1,2,4,8,16,32,64,128,1,2,4,8,16,32,64,128,1,2,4,8,16,32,64,128};for(int z=0;z<32;z++){p[z]=shuffled[z]&mask[z]?'1':'0';}}inline void write_bits64(char*p,uint64_t bits){write_bits(p,uint32_t(bits));write_bits(p+32,uint32_t(bits>>32));}}
#pragma GCC pop_options
#line 1 "cp-algo/util/big_alloc.hpp"
#include <set>
#include <map>
#include <deque>
#include <stack>
#include <queue>
#include <vector>
#include <string>
#line 12 "cp-algo/util/big_alloc.hpp"
#include <iostream>
#include <generator>
#include <forward_list>
#if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
# define CP_ALGO_USE_MMAP 1
# include <sys/mman.h>
#else
# define CP_ALGO_USE_MMAP 0
#endif
namespace cp_algo{template<typename T,size_t Align=32>class big_alloc{static_assert(Align>=alignof(void*),"Align must be at least pointer-size");static_assert(std::popcount(Align)==1,"Align must be a power of two");public:using value_type=T;template<class U>struct rebind{using other=big_alloc<U,Align>;};constexpr bool operator==(const big_alloc&)const=default;constexpr bool operator!=(const big_alloc&)const=default;big_alloc()noexcept=default;template<typename U,std::size_t A>big_alloc(const big_alloc<U,A>&)noexcept{}[[nodiscard]]T*allocate(std::size_t n){std::size_t padded=round_up(n*sizeof(T));std::size_t align=std::max<std::size_t>(alignof(T),Align);
#if CP_ALGO_USE_MMAP
if(padded>=MEGABYTE){void*raw=mmap(nullptr,padded,PROT_READ|PROT_WRITE,MAP_PRIVATE|MAP_ANONYMOUS,-1,0);madvise(raw,padded,MADV_HUGEPAGE);madvise(raw,padded,MADV_POPULATE_WRITE);return static_cast<T*>(raw);}
#endif
return static_cast<T*>(::operator new(padded,std::align_val_t(align)));}void deallocate(T*p,std::size_t n)noexcept{if(!p)return;std::size_t padded=round_up(n*sizeof(T));std::size_t align=std::max<std::size_t>(alignof(T),Align);
#if CP_ALGO_USE_MMAP
if(padded>=MEGABYTE){munmap(p,padded);return;}
#endif
::operator delete(p,padded,std::align_val_t(align));}private:static constexpr std::size_t MEGABYTE=1<<20;static constexpr std::size_t round_up(std::size_t x)noexcept{return(x+Align-1)/Align*Align;}};template<typename T>using big_vector=std::vector<T,big_alloc<T>>;template<typename T>using big_basic_string=std::basic_string<T,std::char_traits<T>,big_alloc<T>>;template<typename T>using big_deque=std::deque<T,big_alloc<T>>;template<typename T>using big_stack=std::stack<T,big_deque<T>>;template<typename T>using big_queue=std::queue<T,big_deque<T>>;template<typename T>using big_priority_queue=std::priority_queue<T,big_vector<T>>;template<typename T>using big_forward_list=std::forward_list<T,big_alloc<T>>;using big_string=big_basic_string<char>;template<typename Key,typename Value,typename Compare=std::less<Key>>using big_map=std::map<Key,Value,Compare,big_alloc<std::pair<const Key,Value>>>;template<typename T,typename Compare=std::less<T>>using big_multiset=std::multiset<T,Compare,big_alloc<T>>;template<typename T,typename Compare=std::less<T>>using big_set=std::set<T,Compare,big_alloc<T>>;template<typename Ref,typename V=void>using big_generator=std::generator<Ref,V,big_alloc<std::byte>>;}namespace std::ranges{template<typename Ref,typename V>elements_of(cp_algo::big_generator<Ref,V>&&)->elements_of<cp_algo::big_generator<Ref,V>&&,cp_algo::big_alloc<std::byte>>;}
#line 5 "cp-algo/structures/bit_array.hpp"
#include <cassert>
CP_ALGO_SIMD_PRAGMA_PUSHnamespace cp_algo::structures{template<typename C>concept Resizable=requires(C&c,std::size_t n){c.resize(n);};template<class Cont>struct _bit_array{using word_t=typename Cont::value_type;static constexpr size_t width=bit_width<word_t>;size_t words,n;alignas(32)Cont data;constexpr void resize(size_t N){n=N;words=(n+width-1)/width;if constexpr(Resizable<Cont>){data.resize(words);}else{assert(std::size(data)>=words);}}constexpr _bit_array():data(){if constexpr(!Resizable<Cont>){resize(std::size(data)*width);}else{resize(0);}}constexpr _bit_array(size_t N):data(){resize(N);}constexpr auto&&word(this auto&&self,size_t x){return self.data[x];}constexpr void set_all(word_t val=-1){for(size_t i=0;i<words;i++){data[i]=val;}}constexpr void reset(){set_all(0);}constexpr void set(size_t x){word(x/width)|=1ULL<<(x%width);}constexpr void reset(size_t x){word(x/width)&=~(1ULL<<(x%width));}constexpr void flip(size_t x){word(x/width)^=1ULL<<(x%width);}constexpr bool test(size_t x)const{return(word(x/width)>>(x%width))&1;}constexpr bool operator[](size_t x)const{return test(x);}constexpr size_t size()const{return n;}auto operator<=>(_bit_array const&t)const=default;constexpr _bit_array&xor_hint(_bit_array const&t,size_t hint){for(size_t i=hint/width;i<words;i++){data[i]^=t.data[i];}return*this;}constexpr _bit_array&operator^=(_bit_array const&t){return xor_hint(t,0);}constexpr _bit_array operator^(_bit_array const&t)const{return _bit_array(*this)^=t;}};template<size_t N>using bit_array=_bit_array<std::array<uint64_t,(N+63)/64>>;using dynamic_bit_array=_bit_array<big_vector<uint64_t>>;}
#pragma GCC pop_options
#line 1 "cp-algo/structures/bit_array_util.hpp"
#line 9 "cp-algo/structures/bit_array_util.hpp"
#include <span>
CP_ALGO_BIT_PRAGMA_PUSHnamespace cp_algo::structures{template<typename BitArray>constexpr void from_string(BitArray&arr,std::span<char>bits){arr.resize(std::size(bits));int64_t i=0;constexpr int width=BitArray::width;for(;i+width<=(int64_t)std::size(bits);i+=width){arr.word(i/width)=read_bits64(bits.data()+i);}for(;i<(int64_t)std::size(bits);i++){if(bits[i]&1){arr.set(i);}}}template<typename BitArray>constexpr big_string to_string(BitArray const&arr){big_string res(arr.words*BitArray::width,'0');for(size_t i=0;i<arr.words;i++){write_bits64(res.data()+i*BitArray::width,arr.word(i));}res.resize(arr.n);return res;}template<typename BitArray>constexpr size_t count(BitArray const&arr,size_t n){size_t res=0;for(size_t i=0;i<n/BitArray::width;i++){res+=std::popcount(arr.word(i));}if(n%BitArray::width){res+=std::popcount(arr.word(n/BitArray::width)&mask(n%BitArray::width));}return res;}template<typename BitArray>constexpr size_t count(BitArray const&arr){return count(arr,arr.n);}template<typename BitArray>constexpr size_t ctz(BitArray const&arr){size_t res=0;size_t i=0;while(i<arr.words&&arr.word(i)==0){res+=BitArray::width;i++;}if(i<arr.words){res+=std::countr_zero(arr.word(i));}return std::min(res,arr.n);}template<typename BitArray>constexpr size_t skip(BitArray const&arr,size_t pos=0,int k=0){if(pos>=arr.n)return arr.n;size_t word_idx=pos/BitArray::width;auto w=arr.word(word_idx)>>(pos%BitArray::width);auto popcnt=std::popcount(w);if(popcnt>k){return pos+cp_algo::kth_set_bit(w,k);}k-=popcnt;while(++word_idx<arr.words){w=arr.word(word_idx);auto popcnt=std::popcount(w);if(popcnt>k)[[unlikely]]{return word_idx*BitArray::width+cp_algo::kth_set_bit(w,k);}k-=popcnt;}return arr.n;}}
#pragma GCC pop_options
#line 1 "cp-algo/util/checkpoint.hpp"
#line 5 "cp-algo/util/checkpoint.hpp"
#include <chrono>
#line 8 "cp-algo/util/checkpoint.hpp"
namespace cp_algo{
#ifdef CP_ALGO_CHECKPOINT
big_map<big_string,double>checkpoints;double last;
#endif
template<bool final=false>void checkpoint([[maybe_unused]]auto const&_msg){
#ifdef CP_ALGO_CHECKPOINT
big_string msg=_msg;double now=(double)clock()/CLOCKS_PER_SEC;double delta=now-last;last=now;if(msg.size()&&!final){checkpoints[msg]+=delta;}if(final){for(auto const&[key,value]:checkpoints){std::cerr<<key<<": "<<value*1000<<" ms\n";}std::cerr<<"Total: "<<now*1000<<" ms\n";}
#endif
}template<bool final=false>void checkpoint(){checkpoint<final>("");}}
#line 11 "cp-algo/math/sieve.hpp"
#include <algorithm>
#line 13 "cp-algo/math/sieve.hpp"
CP_ALGO_BIT_PRAGMA_PUSHnamespace cp_algo::math{using cp_algo::structures::dynamic_bit_array;using cp_algo::structures::bit_array;constexpr auto wheel_primes=std::array{2u,3u,5u,7u};constexpr uint8_t period=std::ranges::fold_left(wheel_primes,1u,std::multiplies{});constexpr uint8_t coprime=std::ranges::fold_left(wheel_primes,1u,[](auto a,auto b){return a*(b-1);});constexpr auto coprime_wheel=[](auto x){return std::ranges::all_of(wheel_primes,[x](auto p){return x%p;});};constexpr auto res_wheel=[](){std::array<uint8_t,coprime>res;int idx=0;for(uint8_t i=1;i<period;i+=2){if(coprime_wheel(i)){res[idx++]=i;}}return res;}();constexpr auto state_wheel=[](){std::array<uint8_t,period>state;uint8_t idx=0;for(uint8_t i=0;i<period;i++){state[i]=idx;idx+=coprime_wheel(i);}return state;}();constexpr auto add_wheel=[](){std::array<uint8_t,period>add;for(uint8_t i=0;i<period;i++){add[i]=1;while(!coprime_wheel(i+add[i])){add[i]++;}}return add;}();constexpr auto gap_wheel=[](){std::array<uint8_t,coprime>gap;for(uint8_t i=0;i<coprime;i++){gap[i]=add_wheel[res_wheel[i]];}return gap;}();constexpr uint32_t to_ord(uint32_t x){return(x/period)*coprime+state_wheel[x%period];}constexpr uint32_t to_val(uint32_t x){return(x/coprime)*period+res_wheel[x%coprime];}constexpr size_t sqrt_threshold=1<<16;constexpr auto sqrt_prime_bits=[](){const int size=sqrt_threshold/2;bit_array<size>prime;prime.set_all();prime.reset(to_ord(1));for(uint32_t i=res_wheel[1];to_ord(i*i)<size;i+=add_wheel[i%period]){if(prime[to_ord(i)]){for(uint32_t k=i;to_ord(i*k)<size;k+=add_wheel[k%period]){prime.reset(to_ord(i*k));}}}return prime;}();constexpr size_t num_primes=[](){size_t cnt=0;for(uint32_t i=res_wheel[1];i<sqrt_threshold;i+=add_wheel[i%period]){cnt+=sqrt_prime_bits[to_ord(i)];}return cnt;}();constexpr auto sqrt_primes=[](){std::array<uint32_t,num_primes>primes;size_t j=0;for(uint32_t i=res_wheel[1];i<sqrt_threshold;i+=add_wheel[i%period]){if(sqrt_prime_bits[to_ord(i)]){primes[j++]=i;}}return primes;}();struct wheel_t{dynamic_bit_array mask;uint32_t product;};constexpr auto make_wheel(big_vector<uint32_t>primes,uint32_t product){assert(product%period==0&&product/period*coprime%dynamic_bit_array::width==0);wheel_t wheel;wheel.product=product;wheel.mask.resize(product/period*coprime);wheel.mask.set_all();for(auto p:primes){for(uint32_t k=1;p*k<product;k+=add_wheel[k%period]){wheel.mask.reset(to_ord(p*k));}}return wheel;}constexpr void sieve_dense(auto&prime,uint32_t l,uint32_t r,wheel_t const&wheel){if(l>=r)return;const uint32_t width=(uint32_t)dynamic_bit_array::width;uint32_t wl=l/width;uint32_t wr=(r+width-1)/width;uint32_t N=(uint32_t)wheel.mask.words;auto loop=[&](uint32_t i,uint32_t block){auto p_ptr=std::assume_aligned<32>(&prime.word(i));auto m_ptr=std::assume_aligned<32>(&wheel.mask.word(0));
#pragma GCC unroll coprime
for(uint32_t j=0;j<block;j++){p_ptr[j]&=m_ptr[j];}};while(wl+N<=wr){loop(wl,N);wl+=N;}loop(wl,wr-wl);}template<class BitArray>constexpr auto sieve_wheel(BitArray&prime,uint32_t l,uint32_t r,size_t i,int state){static const auto ord_step=[](){big_vector<std::array<uint32_t,2*coprime>>ord_steps(num_primes);for(uint32_t i=0;i<size(sqrt_primes);i++){auto p=sqrt_primes[i];auto&ords=ord_steps[i];auto last=to_ord(p);for(uint32_t j=0;j<coprime;j++){auto next=to_ord(p*(res_wheel[j]+gap_wheel[j]));ords[j]=ords[j+coprime]=next-last;last=next;}}return ord_steps;}();auto&ords=ord_step[i];auto advance=[&](){prime.reset(std::exchange(l,l+ords[state++]));};uint32_t p=sqrt_primes[i];while(l+p*coprime<=r){
#pragma GCC unroll coprime
for(size_t j=0;j<coprime;j++){advance();}state-=coprime;}while(l<r){advance();}state=state>=coprime?state-coprime:state;return std::pair{l,state};}constexpr dynamic_bit_array sieve_wheel(uint32_t N){N++;dynamic_bit_array prime(to_ord(N));prime.set_all();static const auto[wheels,medium_primes_begin]=[](){constexpr size_t max_wheel_size=1<<20;uint32_t product=period*dynamic_bit_array::width>>(size(wheel_primes)-2);big_vector<uint32_t>current;big_vector<wheel_t>wheels;for(size_t i=0;i<size(sqrt_primes);i++){uint32_t p=sqrt_primes[i];if(product*p>max_wheel_size){wheels.push_back(make_wheel(current,product));current={p};product=(period*dynamic_bit_array::width>>(size(wheel_primes)-2))*p;if(product>max_wheel_size){checkpoint("make wheels");return std::pair{wheels,i};}}else{current.push_back(p);product*=p;}}assert(false);}();static constexpr uint32_t dense_block=1<<25;for(uint32_t start=0;start<N;start+=dense_block){uint32_t r=std::min(start+dense_block,N);for(auto const&wheel:wheels){auto l=start/wheel.product*wheel.product;sieve_dense(prime,to_ord(l),to_ord(r),wheel);}}checkpoint("dense sieve");static constexpr uint32_t sparse_block=1<<22;auto[pos,state]=[](){big_vector<uint32_t>pos(num_primes);big_vector<uint8_t>state(num_primes);for(auto[i,p]:sqrt_primes|std::views::enumerate){pos[i]=to_ord(p*p);state[i]=state_wheel[p%period];}return std::pair{pos,state};}();for(uint32_t start=0;start<N;start+=sparse_block){uint32_t r=to_ord(std::min(start+sparse_block,N));for(size_t i=medium_primes_begin;i<size(sqrt_primes);i++){if(state[i]>=r)break;std::tie(pos[i],state[i])=sieve_wheel(prime,pos[i],r,i,state[i]);}}checkpoint("sparse sieve");for(size_t i=0;i<std::min(prime.words,sqrt_prime_bits.words);i++){prime.word(i)=sqrt_prime_bits.word(i);}return prime;}}
#pragma GCC pop_options