#pragma GCC optimize("O3,unroll-loops") #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace std; using lint = long long; using pint = pair; using plint = pair; struct fast_ios { fast_ios(){ cin.tie(nullptr), ios::sync_with_stdio(false), cout << fixed << setprecision(20); }; } fast_ios_; #define ALL(x) (x).begin(), (x).end() #define FOR(i, begin, end) for(int i=(begin),i##_end_=(end);i=i##_begin_;i--) #define REP(i, n) FOR(i,0,n) #define IREP(i, n) IFOR(i,0,n) template bool chmax(T &m, const T q) { return m < q ? (m = q, true) : false; } template bool chmin(T &m, const T q) { return m > q ? (m = q, true) : false; } const std::vector> grid_dxs{{1, 0}, {-1, 0}, {0, 1}, {0, -1}}; int floor_lg(long long x) { return x <= 0 ? -1 : 63 - __builtin_clzll(x); } template T1 floor_div(T1 num, T2 den) { return (num > 0 ? num / den : -((-num + den - 1) / den)); } template std::pair operator+(const std::pair &l, const std::pair &r) { return std::make_pair(l.first + r.first, l.second + r.second); } template std::pair operator-(const std::pair &l, const std::pair &r) { return std::make_pair(l.first - r.first, l.second - r.second); } template std::vector sort_unique(std::vector vec) { sort(vec.begin(), vec.end()), vec.erase(unique(vec.begin(), vec.end()), vec.end()); return vec; } template int arglb(const std::vector &v, const T &x) { return std::distance(v.begin(), std::lower_bound(v.begin(), v.end(), x)); } template int argub(const std::vector &v, const T &x) { return std::distance(v.begin(), std::upper_bound(v.begin(), v.end(), x)); } template IStream &operator>>(IStream &is, std::vector &vec) { for (auto &v : vec) is >> v; return is; } template OStream &operator<<(OStream &os, const std::vector &vec); template OStream &operator<<(OStream &os, const std::array &arr); template OStream &operator<<(OStream &os, const std::unordered_set &vec); template OStream &operator<<(OStream &os, const pair &pa); template OStream &operator<<(OStream &os, const std::deque &vec); template OStream &operator<<(OStream &os, const std::set &vec); template OStream &operator<<(OStream &os, const std::multiset &vec); template OStream &operator<<(OStream &os, const std::unordered_multiset &vec); template OStream &operator<<(OStream &os, const std::pair &pa); template OStream &operator<<(OStream &os, const std::map &mp); template OStream &operator<<(OStream &os, const std::unordered_map &mp); template OStream &operator<<(OStream &os, const std::tuple &tpl); template OStream &operator<<(OStream &os, const std::vector &vec) { os << '['; for (auto v : vec) os << v << ','; os << ']'; return os; } template OStream &operator<<(OStream &os, const std::array &arr) { os << '['; for (auto v : arr) os << v << ','; os << ']'; return os; } template std::istream &operator>>(std::istream &is, std::tuple &tpl) { std::apply([&is](auto &&... args) { ((is >> args), ...);}, tpl); return is; } template OStream &operator<<(OStream &os, const std::tuple &tpl) { os << '('; std::apply([&os](auto &&... args) { ((os << args << ','), ...);}, tpl); return os << ')'; } template OStream &operator<<(OStream &os, const std::unordered_set &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; } template OStream &operator<<(OStream &os, const std::deque &vec) { os << "deq["; for (auto v : vec) os << v << ','; os << ']'; return os; } template OStream &operator<<(OStream &os, const std::set &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; } template OStream &operator<<(OStream &os, const std::multiset &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; } template OStream &operator<<(OStream &os, const std::unordered_multiset &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; } template OStream &operator<<(OStream &os, const std::pair &pa) { return os << '(' << pa.first << ',' << pa.second << ')'; } template OStream &operator<<(OStream &os, const std::map &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; } template OStream &operator<<(OStream &os, const std::unordered_map &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; } #ifdef HITONANODE_LOCAL const string COLOR_RESET = "\033[0m", BRIGHT_CYAN = "\033[1;36m"; #define dbgif(cond, x) ((cond) ? std::cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << " (L" << __LINE__ << ") " << __FILE__ << std::endl : std::cerr) #define dbg(x) dbgif(true, x) #else #define dbg(x) ((void)0) #define dbgif(cond, x) ((void)0) #endif // https://hitonanode.github.io/cplib-cpp/unionfind/unionfind.hpp #include #include #include #include // UnionFind Tree (0-indexed), based on size of each disjoint set struct UnionFind { std::vector par, cou; UnionFind(int N = 0) : par(N), cou(N, 1) { iota(par.begin(), par.end(), 0); } int find(int x) { return (par[x] == x) ? x : (par[x] = find(par[x])); } bool unite(int x, int y) { x = find(x), y = find(y); if (x == y) return false; if (cou[x] < cou[y]) std::swap(x, y); par[y] = x, cou[x] += cou[y]; return true; } int count(int x) { return cou[find(x)]; } bool same(int x, int y) { return find(x) == find(y); } std::vector> groups() { std::vector> ret(par.size()); for (int i = 0; i < int(par.size()); ++i) ret[find(i)].push_back(i); ret.erase(std::remove_if(ret.begin(), ret.end(), [&](const std::vector &v) { return v.empty(); }), ret.end()); return ret; } }; // https://hitonanode.github.io/cplib-cpp/random/rand_nondeterministic.hpp #include #include #include #include struct RNGMt19937 { using lint = long long; std::mt19937 mt; RNGMt19937(long long seed) : mt(seed) {} lint operator()(lint x) { // [0, x) assert(x > 0); return this->operator()(0, x); } lint operator()(lint l, lint r) { assert(l < r); std::uniform_int_distribution d(l, r - 1); return d(mt); } }; // RNGMt19937 rnd(123456789); // RNGMt19937 rnd(std::chrono::steady_clock::now().time_since_epoch().count()); // Fast random number generator based on xoshiro256++. struct FastRNG { using ull = unsigned long long; ull s[4]; static ull rotl(ull x, int k) { return (x << k) | (x >> (64 - k)); } static ull splitmix64(ull &x) { ull z = (x += 0x9e3779b97f4a7c15ULL); z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL; z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL; return z ^ (z >> 31); } FastRNG(ull seed) { for (auto &x : s) x = splitmix64(seed); } ull next_u64() { ull res = rotl(s[0] + s[3], 23) + s[0]; ull t = s[1] << 17; s[2] ^= s[0]; s[3] ^= s[1]; s[1] ^= s[2]; s[0] ^= s[3]; s[2] ^= t; s[3] = rotl(s[3], 45); return res; } uint32_t next_u32() { return next_u64() >> 32; } // [0, n) uint32_t operator()(uint32_t n) { assert(n > 0); return (uint64_t(next_u32()) * n) >> 32; } // [l, r) int operator()(int l, int r) { assert(l < r); uint32_t width = uint64_t(int64_t(r) - int64_t(l)); return int(int64_t(l) + (*this)(width)); } uint16_t next_u16() { return next_u64() >> 48; } }; inline FastRNG rnd(123456789); // FastRNG rnd(std::chrono::steady_clock::now().time_since_epoch().count()); // https://hitonanode.github.io/cplib-cpp/heuristic/exponential_dist_sampler.hpp #include #include #include #include template struct ExponentialDistSampler { std::array minuslogps; constexpr ExponentialDistSampler() { for (int i = 0; i < (1 << D); ++i) minuslogps.at(i) = -log((0.5 + i) / (1 << D)); } double sample(uint32_t random_mask) const { return minuslogps.at(random_mask & ((1 << D) - 1)); } // p ~ U(0, 1) => -log(p) ~ Ex(1) // P[exp(-|dx| / T) >= p] = P[|dx| <= -log(p) * T] bool check_sa(double abs_dx, double T, uint32_t random_mask) const { return abs_dx <= minuslogps.at(random_mask & ((1 << D) - 1)) * T; } }; const ExponentialDistSampler<16> log_ps; template struct Annealer { static_assert(UPDATE_INTERVAL > 0); double temp; double ratio; int next_update = UPDATE_INTERVAL; Annealer(double start_temp, double end_temp, int iterations) : temp(start_temp) { assert(iterations > 0); const int updates = (iterations - 1) / UPDATE_INTERVAL; ratio = updates == 0 ? 1.0 : pow(end_temp / start_temp, 1.0 / updates); // ratio = updates == 0 ? 0.0 : (end_temp - start_temp) / updates; } void update(int iter) { if (iter >= next_update) { temp *= ratio; // temp += ratio; next_update += UPDATE_INTERVAL; } } template bool accept_score(int delta_score, RNG &rng, const Sampler &sampler) const { if (delta_score >= 0) return true; return sampler.check_sa(-delta_score, temp, rng.next_u16()); } template bool accept_cost(int delta_cost, RNG &rng, const Sampler &sampler) const { if (delta_cost <= 0) return true; return sampler.check_sa(delta_cost, temp, rng.next_u16()); } }; constexpr int N = 50, K = 100; using InputState = array; using RandState = array; using BanState = array; RandState Gen() { array ret{}; for (auto &c : ret) c = '0' + rnd(10); return ret; } int f(int x, int y) { return x * N + y; } bool isin(int x, int y) { return 0 <= x and x < N and 0 <= y and y < N; } int QuadSum(const InputState &input, const RandState &rand_state, const BanState &ban_state) { int ret = 0; vector visited(N * N); auto IsAlive = [&](int x, int y) { if (!isin(x, y)) return false; const int idx = f(x, y); if (ban_state[idx]) return false; return input[idx] > rand_state[idx]; }; REP(i, N) REP(j, N) { if (visited[f(i, j)]) continue; int sum = 0; auto rec = [&](auto &&self, int x, int y) -> void { if (!IsAlive(x, y)) return; const int idx = f(x, y); if (visited[idx]) return; visited[idx] = 1; sum++; self(self, x + 1, y); self(self, x - 1, y); self(self, x, y + 1); self(self, x, y - 1); }; rec(rec, i, j); ret += sum * sum; } return ret; } void BurnProb(const InputState &input, const RandState &rand_state, const BanState &ban_state, vector &burn_weight) { assert((int)burn_weight.size() == N * N); UnionFind uf(N * N); auto IsAlive = [&](int x, int y) { const int idx = f(x, y); if (ban_state[idx]) return false; return !ban_state[idx] and (input[idx] > rand_state[idx]); }; REP(i, N) REP(j, N) { if (!IsAlive(i, j)) continue; if (i + 1 < N and IsAlive(i + 1, j)) uf.unite(f(i, j), f(i + 1, j)); if (j + 1 < N and IsAlive(i, j + 1)) uf.unite(f(i, j), f(i, j + 1)); } REP(i, N) REP(j, N) { burn_weight[f(i, j)] += uf.count(f(i, j)); } } double CalcE(const InputState &input, const vector &rand_states, const BanState &ban_state) { lint quad_sum = 0; for (const auto &rand_state : rand_states) { quad_sum += QuadSum(input, rand_state, ban_state); } return double(quad_sum) / ((lint)rand_states.size() * N * N); } vector BurnProbMap(const InputState &input, const vector &rand_states, const BanState &ban_state) { vector ret(N * N); for(const auto &rand_state : rand_states) { BurnProb(input, rand_state, ban_state, ret); } for (auto &x : ret) x /= (lint)rand_states.size() * N * N; return ret; } const vector dxs8{1, -1, 0, 0, 1, 1, -1, -1}; const vector dys8{0, 0, 1, -1, 1, -1, 1, -1}; struct Line { int x0, y0; int dir; template friend OStream &operator<<(OStream &os, const Line &line) { os << "Line(" << line.x0 << ", " << line.y0 << ", " << line.dir << ")"; return os; } }; void Apply(BanState &state, const BanState &ban_candidate, const Line &line) { const int dx = dxs8[line.dir]; const int dy = dys8[line.dir]; int x = line.x0, y = line.y0; assert(isin(x, y)); while (true) { if (!ban_candidate[f(x, y)]) break; if (state[f(x, y)]) break; state[f(x, y)] = true; x += dx; y += dy; if (!isin(x, y)) break; } } BanState Apply(const vector &lines, const BanState &ban_candidate) { BanState ret{}; for (const auto &line : lines) Apply(ret, ban_candidate, line); return ret; } int main() { InputState input_state{}; vector S(N); { int N_, K_; cin >> N_ >> K_; assert(N == N_); assert(K == K_); for (auto &s : S) cin >> s; REP(i, N) REP(j, N) input_state[f(i, j)] = S[i][j]; } vector rand_states; REP(_, 50) rand_states.push_back(Gen()); dbg(rand_states.size()); // REP(_, 1000) { // const auto e1 = CalcE(input_state, rand_states, ban0); // // dbg(e1); // } // const int score = round(1e6 * (1 + log2(e0 / e1))); // dbg(score); // dbg(BurnProbMap(input_state, rand_states, ban0)); auto IsTarget = [&](int x, int y) { const int idx = f(x, y); if (input_state[idx] >= '5') return true; if (input_state[idx] < '4') return false; for (auto [dx, dy] : grid_dxs) { const int nx = x + dx; const int ny = y + dy; if (isin(nx, ny) && input_state[f(nx, ny)] >= '8') return true; } return false; }; BanState ban_candidate{}; REP(x, N) REP(y, N) { if (IsTarget(x, y)) ban_candidate[f(x, y)] = true; } auto FastEval = [&](const BanState &ban_state) { UnionFind uf(N * N); REP(x, N) REP(y, N) { if (ban_state[f(x, y)]) continue; if (!IsTarget(x, y)) continue; if (x + 1 < N && !ban_state[f(x + 1, y)] && IsTarget(x + 1, y)) uf.unite(f(x, y), f(x + 1, y)); if (y + 1 < N && !ban_state[f(x, y + 1)] && IsTarget(x, y + 1)) uf.unite(f(x, y), f(x, y + 1)); } lint ret = 0; REP(idx, N * N) { if (ban_state[idx]) continue; if (!IsTarget(idx / N, idx % N)) continue; if (uf.find(idx) == idx) { const lint c = uf.count(idx); ret += c * c * c; } } return ret; }; BanState state{}; vector best_lines; while (true) { const int k = count(ALL(state), true); const int rem = K - k; if (rem <= 0) break; double max_eval = -1e200; Line best_line{-1, -1, -1}; const auto orig_e = FastEval(state); REP(x0, N) REP(y0, N) { if (!IsTarget(x0, y0)) continue; if (state[f(x0, y0)]) continue; if ((x0 + y0) % 2) continue; FOR(dir, 4, 8) { const int dx = dxs8[dir]; const int dy = dys8[dir]; int nban = 0; BanState next_state = state; FOR(len, 1, N + 1) { const int nx = x0 + dx * (len - 1); const int ny = y0 + dy * (len - 1); if (!isin(nx, ny)) break; if (!IsTarget(nx, ny)) break; if (chmax(next_state[f(nx, ny)], true)) { ++nban; if (nban > 30 or nban > rem) break; const int nnx = nx + dx; const int nny = ny + dy; if (!isin(nnx, nny) or !IsTarget(nnx, nny) or state[f(nnx, nny)]) { const double imp = orig_e - FastEval(next_state); if (chmax(max_eval, imp * imp * imp / nban / nban)) { best_line = Line{x0, y0, dir}; } break; } } else { break; } } } } dbg(make_tuple(best_line, max_eval)); if (best_line.x0 != -1) { Apply(state, ban_candidate, best_line); best_lines.push_back(best_line); } else { break; } } dbg(best_lines); auto best_eval = CalcE(input_state, rand_states, Apply(best_lines, ban_candidate)); dbg(best_eval); { // auto state = best_state; auto current_lines = best_lines; auto state_eval = best_eval; Annealer sa(0.1, 0.01, 800); REP(_, 800) { sa.update(_); auto lines = current_lines; const int idx = rnd(lines.size()); const int d = rnd(8); lines[idx].x0 += dxs8[d]; lines[idx].y0 += dys8[d]; if (!isin(lines[idx].x0, lines[idx].y0) or !IsTarget(lines[idx].x0, lines[idx].y0)) continue; auto new_state = Apply(lines, ban_candidate); if (count(ALL(new_state), true) > K) continue; const auto new_eval = CalcE(input_state, rand_states, new_state); if (chmin(best_eval, new_eval)) { best_eval = new_eval; best_lines = lines; dbg(make_tuple(_, new_eval)); } if (sa.accept_cost(new_eval - state_eval, rnd, log_ps)) { dbg(make_tuple(_, new_eval, state_eval)); state_eval = new_eval; current_lines = lines; } } } { auto best_state = Apply(best_lines, ban_candidate); cout << count(ALL(best_state), true) << "\n"; REP(i, N) REP(j, N) { if (best_state[f(i, j)]) cout << i << " " << j << "\n"; } } }