結果
| 問題 | No.5025 Forest Fire |
| コンテスト | |
| ユーザー |
 mech_39
|
| 提出日時 | 2026-10-19 20:13:20 |
| 言語 | C++23(gcc16) (gcc 16.2.0 + boost 1.92.0 + ACL) |
| 結果 |
AC
不安定
|
| 実行時間 | 1,754 ms / 2,000 ms |
| + 421µs | |
| コード長 | 47,285 bytes |
| 記録 | |
| コンパイル時間 | 8,694 ms |
| コンパイル使用メモリ | 426,188 KB |
| 実行使用メモリ | 9,928 KB |
| スコア | 184,563,539 |
| 最終ジャッジ日時 | 2026-10-24 14:32:42 |
| 合計ジャッジ時間 | 154,478 ms |
|
ジャッジサーバーID （参考情報） |
judge2_0 / <nil> |
| 純コード判定しない問題か言語 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 50 |
ソースコード
#define MODE 3
#define GUIDE_DEF 20
#define H_DEAD1 0.35
#define H_DEAD2 0.75
#ifndef MODE
#define MODE 0
#endif
#ifndef GUIDE_DEF
#define GUIDE_DEF 20
#endif
#ifndef F_DEAD
#define F_DEAD 0.9
#endif
#ifndef H_DEAD1
#define H_DEAD1 0.55
#endif
#ifndef H_DEAD2
#define H_DEAD2 1.05
#endif
#include <bits/stdc++.h>
using namespace std;
int N, K, V;
vector<int> dry10;
vector<char> tree0;
chrono::steady_clock::time_point T0;
double elapsed() { return chrono::duration<double>(chrono::steady_clock::now() - T0).count(); }
const int DX[4] = {0, 0, 1, -1}, DY[4] = {1, -1, 0, 0};
double TIME_LIMIT = 1.80;
struct RNG
{
uint64_t s;
explicit RNG(uint64_t seed) : s(seed) {}
inline uint64_t next()
{
uint64_t z = (s += 0x9E3779B97F4A7C15ULL);
z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
return z ^ (z >> 31);
}
inline int nextInt(int n) { return (int)(next() % (uint64_t)n); }
inline double nextDouble() { return (next() >> 11) * (1.0 / 9007199254740992.0); }
};
vector<array<int, 4>> NB; // grid neighbours (-1 if outside)
struct Samples
{
int S;
vector<vector<uint8_t>> dry;
void gen(int s_, uint64_t seed)
{
S = s_;
RNG rng(seed);
dry.assign(S, vector<uint8_t>(V));
for (int s = 0; s < S; s++)
for (int v = 0; v < V; v++)
dry[s][v] = (rng.next() % 10) < (uint64_t)dry10[v];
}
};
vector<int> ufp, ufs;
int findr(int x)
{
while (ufp[x] != x)
{
ufp[x] = ufp[ufp[x]];
x = ufp[x];
}
return x;
}
double evaluate(const vector<char> &tr, const Samples &sm)
{
ufp.resize(V);
ufs.resize(V);
double tot = 0;
vector<char> act(V);
for (int s = 0; s < sm.S; s++)
{
const auto &dr = sm.dry[s];
for (int v = 0; v < V; v++)
{
ufp[v] = v;
ufs[v] = 1;
act[v] = tr[v] && dr[v];
}
for (int v = 0; v < V; v++)
{
if (!act[v])
continue;
if (v % N + 1 < N && act[v + 1])
{
int a = findr(v), b = findr(v + 1);
if (a != b)
{
if (ufs[a] < ufs[b])
swap(a, b);
ufp[b] = a;
ufs[a] += ufs[b];
}
}
if (v + N < V && act[v + N])
{
int a = findr(v), b = findr(v + N);
if (a != b)
{
if (ufs[a] < ufs[b])
swap(a, b);
ufp[b] = a;
ufs[a] += ufs[b];
}
}
}
long long sq = 0;
for (int v = 0; v < V; v++)
if (act[v] && ufp[v] == v)
sq += (long long)ufs[v] * ufs[v];
tot += sq;
}
return tot / sm.S / ((double)N * N);
}
// ---------------- incremental articulation greedy ----------------
struct APGreedy
{
const Samples *sm;
vector<char> tr;
vector<vector<long long>> contrib;
vector<double> gain;
vector<int> disc, low, sub, parent, stv, sti, order, mark;
vector<long long> sepSq, sepCnt;
int stamp = 0;
inline int nb(int s, int v, int d) const
{
int u = NB[v][d];
return (u >= 0 && tr[u] && sm->dry[s][u]) ? u : -1;
}
void tarjan(int s, int r)
{
order.clear();
int top = 0, timer = 0;
stv[top] = r;
sti[top] = 0;
top++;
mark[r] = stamp;
disc[r] = low[r] = timer++;
parent[r] = -1;
sub[r] = 1;
sepSq[r] = 0;
sepCnt[r] = 0;
order.push_back(r);
while (top > 0)
{
int v = stv[top - 1];
int &it = sti[top - 1];
if (it < 4)
{
int u = nb(s, v, it++);
if (u < 0)
continue;
if (mark[u] != stamp)
{
mark[u] = stamp;
disc[u] = low[u] = timer++;
parent[u] = v;
sub[u] = 1;
sepSq[u] = 0;
sepCnt[u] = 0;
order.push_back(u);
stv[top] = u;
sti[top] = 0;
top++;
}
else if (u != parent[v])
low[v] = min(low[v], disc[u]);
}
else
{
top--;
int pv = parent[v];
if (pv >= 0)
{
low[pv] = min(low[pv], low[v]);
sub[pv] += sub[v];
if (low[v] >= disc[pv])
{
sepSq[pv] += (long long)sub[v] * sub[v];
sepCnt[pv] += sub[v];
}
}
}
}
long long c = sub[r], c2 = c * c;
for (int v : order)
{
long long rest = c - 1 - sepCnt[v];
long long val = c2 - (sepSq[v] + rest * rest);
contrib[s][v] = val;
gain[v] += val;
}
}
void init(const Samples *sm_, const vector<char> &tr_)
{
sm = sm_;
tr = tr_;
contrib.assign(sm->S, vector<long long>(V, 0));
gain.assign(V, 0);
disc.assign(V, 0);
low.assign(V, 0);
sub.assign(V, 0);
parent.assign(V, 0);
stv.assign(V, 0);
sti.assign(V, 0);
mark.assign(V, -1);
sepSq.assign(V, 0);
sepCnt.assign(V, 0);
for (int s = 0; s < sm->S; s++)
{
stamp++;
for (int v = 0; v < V; v++)
if (tr[v] && sm->dry[s][v] && mark[v] != stamp)
tarjan(s, v);
}
}
vector<int> bfsq;
void cut(int x)
{
for (int s = 0; s < sm->S; s++)
{
if (!sm->dry[s][x])
continue;
stamp++;
bfsq.clear();
bfsq.push_back(x);
mark[x] = stamp;
for (size_t h = 0; h < bfsq.size(); h++)
{
int v = bfsq[h];
for (int d = 0; d < 4; d++)
{
int u = nb(s, v, d);
if (u >= 0 && mark[u] != stamp)
{
mark[u] = stamp;
bfsq.push_back(u);
}
}
}
for (int v : bfsq)
{
gain[v] -= contrib[s][v];
contrib[s][v] = 0;
}
tr[x] = 0;
stamp++;
for (int v : bfsq)
if (v != x && mark[v] != stamp)
tarjan(s, v);
tr[x] = 1;
}
tr[x] = 0;
gain[x] = -1e18;
}
int best() const
{
int b = -1;
double bg = -1e17;
for (int v = 0; v < V; v++)
if (tr[v] && gain[v] > bg)
{
bg = gain[v];
b = v;
}
return b;
}
};
vector<int> greedyFill(const Samples &sm, vector<char> tr, vector<int> cuts, double tEnd)
{
APGreedy g;
g.init(&sm, tr);
while ((int)cuts.size() < K && elapsed() < tEnd)
{
int b = g.best();
if (b < 0)
break;
g.cut(b);
tr[b] = 0;
cuts.push_back(b);
}
if ((int)cuts.size() < K)
{
vector<int> cand;
for (int v = 0; v < V; v++)
if (tr[v])
cand.push_back(v);
sort(cand.begin(), cand.end(), [&](int a, int b)
{ return g.gain[a] > g.gain[b]; });
for (int v : cand)
{
if ((int)cuts.size() >= K)
break;
tr[v] = 0;
cuts.push_back(v);
}
}
return cuts;
}
// ---------------- walls via min vertex cut ----------------
struct Dinic
{
struct E
{
int to, cap;
};
vector<E> es;
vector<vector<int>> g;
vector<int> level, it;
void init(int n)
{
es.clear();
g.assign(n, {});
}
void add(int a, int b, int c)
{
g[a].push_back(es.size());
es.push_back({b, c});
g[b].push_back(es.size());
es.push_back({a, 0});
}
bool bfs(int s, int t)
{
level.assign(g.size(), -1);
queue<int> q;
q.push(s);
level[s] = 0;
while (!q.empty())
{
int v = q.front();
q.pop();
for (int id : g[v])
if (es[id].cap > 0 && level[es[id].to] < 0)
{
level[es[id].to] = level[v] + 1;
q.push(es[id].to);
}
}
return level[t] >= 0;
}
int dfs(int v, int t, int f)
{
if (v == t)
return f;
for (int &i = it[v]; i < (int)g[v].size(); i++)
{
E &e = es[g[v][i]];
if (e.cap > 0 && level[e.to] == level[v] + 1)
{
int d = dfs(e.to, t, min(f, e.cap));
if (d > 0)
{
e.cap -= d;
es[g[v][i] ^ 1].cap += d;
return d;
}
}
}
return 0;
}
int maxflow(int s, int t, int limit)
{
int flow = 0;
while (flow < limit && bfs(s, t))
{
it.assign(g.size(), 0);
int f;
while ((f = dfs(s, t, INT_MAX)) > 0)
{
flow += f;
if (flow >= limit)
break;
}
}
return flow;
}
};
// ---------------- SA state with per-sample component labels ----------------
struct SAState
{
const Samples *sm;
int S;
vector<char> tr; // tree and not cut
vector<vector<int>> comp; // [s][v] component id (valid if active)
vector<vector<int>> csz; // [s][id] size
vector<vector<int>> freeIds; // [s] free id pool
long long total = 0; // sum over samples of sum |C|^2
// scratch
vector<int> mark;
int stamp = 0;
vector<int> owner; // group owner for interleaved BFS
vector<int> q[4];
inline bool act(int s, int v) const { return tr[v] && sm->dry[s][v]; }
int newId(int s)
{
int id = freeIds[s].back();
freeIds[s].pop_back();
return id;
}
void init(const Samples *sm_, const vector<char> &tr_)
{
sm = sm_;
S = sm->S;
tr = tr_;
comp.assign(S, vector<int>(V, -1));
csz.assign(S, vector<int>(V + 1, 0));
freeIds.assign(S, {});
mark.assign(V, 0);
owner.assign(V, -1);
total = 0;
for (int s = 0; s < S; s++)
{
for (int id = V; id >= 0; id--)
freeIds[s].push_back(id);
for (int v = 0; v < V; v++)
{
if (!act(s, v) || comp[s][v] >= 0)
continue;
int id = newId(s);
vector<int> st = {v};
comp[s][v] = id;
int cnt = 0;
while (!st.empty())
{
int x = st.back();
st.pop_back();
cnt++;
for (int d = 0; d < 4; d++)
{
int u = NB[x][d];
if (u >= 0 && act(s, u) && comp[s][u] < 0)
{
comp[s][u] = id;
st.push_back(u);
}
}
}
csz[s][id] = cnt;
total += (long long)cnt * cnt;
}
}
}
// ---- delta of restoring a (a is currently cut) in sample s: returns delta and fills merged info
inline long long restoreDelta(int s, int a, int *ids, int &nid, long long &msz) const
{
nid = 0;
msz = 1;
long long old = 0;
for (int d = 0; d < 4; d++)
{
int u = NB[a][d];
if (u < 0 || !act(s, u))
continue;
int id = comp[s][u];
bool dup = false;
for (int k = 0; k < nid; k++)
if (ids[k] == id)
dup = true;
if (dup)
continue;
ids[nid++] = id;
msz += csz[s][id];
old += (long long)csz[s][id] * csz[s][id];
}
return msz * msz - old;
}
// ---- split computation: remove b from its component (total size csize) with optional extra active cell
// `extra` (restored a, treated active) and the component considered is the one containing b.
// Interleaved BFS from b's neighbours; returns sum of squares of pieces (excluding b).
// If `pieces` != nullptr, stores the finished (small) pieces' cell lists for relabeling.
long long splitSq(int s, int b, int extra, long long csize, vector<vector<int>> *pieces)
{
int starts[4], ns = 0;
for (int d = 0; d < 4; d++)
{
int u = NB[b][d];
if (u >= 0 && (u == extra || act(s, u)))
starts[ns++] = u;
}
if (ns == 0)
return 0;
if (ns == 1)
{
long long r = csize - 1;
return r * r;
}
stamp++;
int gpar[4];
long long gsize[4];
bool done[4];
size_t head[4];
for (int k = 0; k < ns; k++)
{
gpar[k] = k;
gsize[k] = 1;
done[k] = false;
head[k] = 0;
q[k].clear();
q[k].push_back(starts[k]);
mark[starts[k]] = stamp;
owner[starts[k]] = k;
}
auto gfind = [&](int x)
{ while (gpar[x] != x) x = gpar[x]; return x; };
long long finSq = 0, finCnt = 0;
int alive = ns;
while (alive > 1)
{
for (int k = 0; k < ns; k++)
{
int r = gfind(k);
if (done[r] || head[k] >= q[k].size())
continue;
int x = q[k][head[k]++];
for (int d = 0; d < 4; d++)
{
int u = NB[x][d];
if (u < 0 || u == b)
continue;
if (!(u == extra || act(s, u)))
continue;
if (mark[u] == stamp)
{
int r1 = gfind(owner[u]), r2 = gfind(k);
if (r1 != r2)
{
gpar[r2] = r1;
gsize[r1] += gsize[r2];
alive--;
}
}
else
{
mark[u] = stamp;
owner[u] = k;
q[k].push_back(u);
gsize[gfind(k)]++;
}
}
}
if (alive <= 1)
break;
for (int r = 0; r < ns; r++)
{
if (gpar[r] != r || done[r])
continue;
bool empty = true;
for (int t = 0; t < ns; t++)
if (gfind(t) == r && head[t] < q[t].size())
{
empty = false;
break;
}
if (!empty)
continue;
if (alive <= 1)
break;
done[r] = true;
alive--;
finSq += gsize[r] * gsize[r];
finCnt += gsize[r];
if (pieces)
{
vector<int> cells;
for (int t = 0; t < ns; t++)
if (gfind(t) == r)
cells.insert(cells.end(), q[t].begin(), q[t].end());
pieces->push_back(move(cells));
}
}
}
long long rest = csize - 1 - finCnt;
return finSq + rest * rest;
}
// delta of moving cut a -> b (a cut now, b tree now) in sample s
long long moveDelta(int s, int a, int b)
{
bool da = sm->dry[s][a], db = sm->dry[s][b];
if (!da && !db)
return 0;
int ids[4], nid;
long long msz;
if (da && !db)
return restoreDelta(s, a, ids, nid, msz);
if (!da && db)
{
long long c = csz[s][comp[s][b]];
return splitSq(s, b, -1, c, nullptr) - c * c;
}
long long dR = restoreDelta(s, a, ids, nid, msz);
int cb = comp[s][b];
bool inMerged = false;
for (int k = 0; k < nid; k++)
if (ids[k] == cb)
inMerged = true;
if (!inMerged)
{
long long c = csz[s][cb];
return dR + splitSq(s, b, -1, c, nullptr) - c * c;
}
return dR + splitSq(s, b, a, msz, nullptr) - msz * msz;
}
void commitRestore(int s, int a)
{
if (!sm->dry[s][a])
return;
int ids[4], nid;
long long msz;
long long d = restoreDelta(s, a, ids, nid, msz);
total += d;
if (nid == 0)
{
int id = newId(s);
comp[s][a] = id;
csz[s][id] = 1;
return;
}
// keep the largest id, relabel others
int keep = ids[0];
for (int k = 1; k < nid; k++)
if (csz[s][ids[k]] > csz[s][keep])
keep = ids[k];
comp[s][a] = keep;
for (int k = 0; k < nid; k++)
{
int id = ids[k];
if (id == keep)
continue;
// relabel component id -> keep via DFS from a neighbour having id
int start = -1;
for (int d2 = 0; d2 < 4; d2++)
{
int u = NB[a][d2];
if (u >= 0 && act(s, u) && comp[s][u] == id)
{
start = u;
break;
}
}
vector<int> st = {start};
comp[s][start] = keep;
while (!st.empty())
{
int x = st.back();
st.pop_back();
for (int d2 = 0; d2 < 4; d2++)
{
int u = NB[x][d2];
if (u >= 0 && act(s, u) && comp[s][u] == id)
{
comp[s][u] = keep;
st.push_back(u);
}
}
}
csz[s][id] = 0;
freeIds[s].push_back(id);
}
csz[s][keep] = (int)msz;
}
void commitCut(int s, int b)
{
if (!sm->dry[s][b])
return;
int id = comp[s][b];
long long c = csz[s][id];
vector<vector<int>> pieces;
long long sq = splitSq(s, b, -1, c, &pieces);
total += sq - c * c;
comp[s][b] = -1;
long long moved = 0;
for (auto &cells : pieces)
{
int nid2 = newId(s);
for (int x : cells)
comp[s][x] = nid2;
csz[s][nid2] = cells.size();
moved += cells.size();
}
csz[s][id] = (int)(c - 1 - moved);
if (csz[s][id] == 0)
freeIds[s].push_back(id);
}
void applyMove(int a, int b)
{
tr[a] = 1;
for (int s = 0; s < S; s++)
commitRestore(s, a);
for (int s = 0; s < S; s++)
{
if (!sm->dry[s][b])
continue;
commitCut(s, b);
}
tr[b] = 0;
}
};
// ---------------- wall candidates: min vertex cuts on dryness-thresholded regions ----------------
struct Cand
{
vector<int> cells;
double proxy = 0;
double gain = 0;
};
struct WallGen
{
vector<int> loc;
WallGen() { loc.assign(V, -1); }
vector<int> bfsDist(int src, const vector<int> &label, int lab)
{
vector<int> d(V, -1);
static vector<int> q;
q.clear();
q.push_back(src);
d[src] = 0;
for (size_t h = 0; h < q.size(); h++)
{
int v = q[h];
for (int k = 0; k < 4; k++)
{
int u = NB[v][k];
if (u < 0 || label[u] != lab || d[u] >= 0)
continue;
d[u] = d[v] + 1;
q.push_back(u);
}
}
return d;
}
// try min cut between A = {ds <= ra} and B = {dt <= rb}; returns cut cells and proxy gain/cell
bool tryCut(const vector<int> &cells, const vector<int> &label, int lab, const vector<int> &ds, const vector<int> &dt,
int ra, int rb, int budget, Cand &out)
{
int n = cells.size();
for (int k = 0; k < n; k++)
loc[cells[k]] = k;
Dinic din;
din.init(2 * n + 2);
int SRC = 2 * n, SNK = 2 * n + 1;
const int INF = 1 << 29;
bool ok = true;
for (int k = 0; k < n && ok; k++)
{
int v = cells[k];
bool inA = ds[v] <= ra, inB = dt[v] <= rb;
if (inA && inB)
{
ok = false;
break;
}
din.add(2 * k, 2 * k + 1, (inA || inB) ? INF : 1);
if (inA)
din.add(SRC, 2 * k, INF);
if (inB)
din.add(2 * k + 1, SNK, INF);
for (int d = 0; d < 4; d++)
{
int u = NB[v][d];
if (u < 0 || label[u] != lab)
continue;
din.add(2 * k + 1, 2 * loc[u], INF);
}
}
if (!ok)
{
for (int v : cells)
loc[v] = -1;
return false;
}
int f = din.maxflow(SRC, SNK, budget + 1);
if (f > budget || f == 0)
{
for (int v : cells)
loc[v] = -1;
return false;
}
vector<char> vis(2 * n + 2, 0);
vector<int> q = {SRC};
vis[SRC] = 1;
for (size_t h = 0; h < q.size(); h++)
{
int v = q[h];
for (int id : din.g[v])
if (din.es[id].cap > 0 && !vis[din.es[id].to])
{
vis[din.es[id].to] = 1;
q.push_back(din.es[id].to);
}
}
vector<int> cut;
for (int k = 0; k < n; k++)
if (vis[2 * k] && !vis[2 * k + 1])
cut.push_back(cells[k]);
// proxy: sizes of parts within the region
vector<char> isCut(n, 0);
for (int v : cut)
isCut[loc[v]] = 1;
vector<char> seen(n, 0);
long long sq = 0;
for (int k0 = 0; k0 < n; k0++)
{
if (isCut[k0] || seen[k0])
continue;
long long cnt = 0;
vector<int> st = {k0};
seen[k0] = 1;
while (!st.empty())
{
int k = st.back();
st.pop_back();
cnt++;
int v = cells[k];
for (int d = 0; d < 4; d++)
{
int u = NB[v][d];
if (u < 0 || label[u] != lab)
continue;
int ku = loc[u];
if (isCut[ku] || seen[ku])
continue;
seen[ku] = 1;
st.push_back(ku);
}
}
sq += cnt * cnt;
}
for (int v : cells)
loc[v] = -1;
if ((int)cut.size() != f)
return false;
out.cells = cut;
out.proxy = (double)((long long)n * n - sq) / cut.size();
return true;
}
void regionCands(const vector<int> &cells, const vector<int> &label, int lab, int budget, mt19937 &rng, vector<Cand> &out)
{
vector<pair<int, int>> pairs;
{
vector<int> d0 = bfsDist(cells[rng() % cells.size()], label, lab);
int s = cells[0];
for (int v : cells)
if (d0[v] > d0[s])
s = v;
vector<int> d1 = bfsDist(s, label, lab);
int t = cells[0];
for (int v : cells)
if (d1[v] > d1[t])
t = v;
pairs.push_back({s, t});
int xl = cells[0], xr = cells[0], yl = cells[0], yr = cells[0];
for (int v : cells)
{
if (v % N < xl % N)
xl = v;
if (v % N > xr % N)
xr = v;
if (v / N < yl / N)
yl = v;
if (v / N > yr / N)
yr = v;
}
pairs.push_back({xl, xr});
pairs.push_back({yl, yr});
}
vector<Cand> local;
for (auto [s, t] : pairs)
{
if (s == t)
continue;
vector<int> ds = bfsDist(s, label, lab), dt = bfsDist(t, label, lab);
vector<int> sd, td;
for (int v : cells)
{
sd.push_back(ds[v]);
td.push_back(dt[v]);
}
sort(sd.begin(), sd.end());
sort(td.begin(), td.end());
const double fr[][2] = {{0.4, 0.4}, {0.25, 0.25}, {0.15, 0.6}, {0.6, 0.15}, {0.1, 0.75}, {0.75, 0.1}};
for (auto &f : fr)
{
int ra = sd[min((int)sd.size() - 1, (int)(f[0] * sd.size()))];
int rb = td[min((int)td.size() - 1, (int)(f[1] * td.size()))];
Cand c;
if (tryCut(cells, label, lab, ds, dt, ra, rb, budget, c))
{
sort(c.cells.begin(), c.cells.end());
bool dup = false;
for (auto &o : local)
if (o.cells == c.cells)
dup = true;
if (!dup)
local.push_back(c);
}
}
}
sort(local.begin(), local.end(), [](const Cand &a, const Cand &b)
{ return a.proxy > b.proxy; });
for (int i = 0; i < (int)local.size() && i < 3; i++)
out.push_back(local[i]);
}
// candidates for all thresholds
map<uint64_t, vector<Cand>> cache;
vector<Cand> generate(const vector<char> &tr, int budget, const vector<int> &thetas, mt19937 &rng)
{
vector<Cand> out;
for (int th : thetas)
{
vector<int> label(V, -1);
int c = 0;
for (int s = 0; s < V; s++)
{
if (!tr[s] || dry10[s] < th || label[s] >= 0)
continue;
vector<int> cells = {s};
label[s] = c;
for (size_t h = 0; h < cells.size(); h++)
{
int v = cells[h];
for (int d = 0; d < 4; d++)
{
int u = NB[v][d];
if (u < 0 || !tr[u] || dry10[u] < th || label[u] >= 0)
continue;
label[u] = c;
cells.push_back(u);
}
}
if (cells.size() >= 40)
{
uint64_t h = 1469598103934665603ULL ^ (uint64_t)th;
vector<int> sorted = cells;
sort(sorted.begin(), sorted.end());
for (int v : sorted)
{
h ^= (uint64_t)v + 0x9e37;
h *= 1099511628211ULL;
}
auto it = cache.find(h);
if (it == cache.end())
{
vector<Cand> rc;
regionCands(cells, label, c, budget, rng, rc);
it = cache.emplace(h, rc).first;
}
for (auto &cd : it->second)
if ((int)cd.cells.size() <= budget)
out.push_back(cd);
}
c++;
}
}
return out;
}
};
long long sumSq(const vector<char> &tr, const Samples &sm)
{
ufp.resize(V);
ufs.resize(V);
long long tot = 0;
static vector<char> act;
act.resize(V);
for (int s = 0; s < sm.S; s++)
{
const auto &dr = sm.dry[s];
for (int v = 0; v < V; v++)
{
ufp[v] = v;
ufs[v] = 1;
act[v] = tr[v] && dr[v];
}
for (int v = 0; v < V; v++)
{
if (!act[v])
continue;
if (v % N + 1 < N && act[v + 1])
{
int a = findr(v), b = findr(v + 1);
if (a != b)
{
if (ufs[a] < ufs[b])
swap(a, b);
ufp[b] = a;
ufs[a] += ufs[b];
}
}
if (v + N < V && act[v + N])
{
int a = findr(v), b = findr(v + N);
if (a != b)
{
if (ufs[a] < ufs[b])
swap(a, b);
ufp[b] = a;
ufs[a] += ufs[b];
}
}
}
for (int v = 0; v < V; v++)
if (act[v] && ufp[v] == v)
tot += (long long)ufs[v] * ufs[v];
}
return tot;
}
struct Walls
{
vector<char> tr;
vector<int> bfsDist(int src, const vector<int> &label, int lab)
{
vector<int> d(V, -1);
queue<int> q;
q.push(src);
d[src] = 0;
while (!q.empty())
{
int v = q.front();
q.pop();
for (int k = 0; k < 4; k++)
{
int u = NB[v][k];
if (u < 0 || label[u] != lab || d[u] >= 0)
continue;
d[u] = d[v] + 1;
q.push(u);
}
}
return d;
}
int labelComponents(vector<int> &label, vector<int> &sizes)
{
label.assign(V, -1);
sizes.clear();
int c = 0;
for (int s = 0; s < V; s++)
{
if (!tr[s] || label[s] >= 0)
continue;
int cnt = 0;
queue<int> q;
q.push(s);
label[s] = c;
while (!q.empty())
{
int v = q.front();
q.pop();
cnt++;
for (int k = 0; k < 4; k++)
{
int u = NB[v][k];
if (u < 0 || !tr[u] || label[u] >= 0)
continue;
label[u] = c;
q.push(u);
}
}
sizes.push_back(cnt);
c++;
}
return c;
}
struct Split
{
double gain = -1;
vector<int> cut;
};
map<uint64_t, Split> splitCache;
Split tryCut(const vector<int> &cells, const vector<int> &label, int lab, const vector<int> &ds, const vector<int> &dt,
int ra, int rb, int budget)
{
Split res;
int n = cells.size();
vector<int> loc(V, -1);
for (int k = 0; k < n; k++)
loc[cells[k]] = k;
Dinic din;
din.init(2 * n + 2);
int SRC = 2 * n, SNK = 2 * n + 1;
const int INF = 1 << 29;
for (int k = 0; k < n; k++)
{
int v = cells[k];
bool inA = ds[v] <= ra, inB = dt[v] <= rb;
if (inA && inB)
return res;
din.add(2 * k, 2 * k + 1, (inA || inB) ? INF : 1);
if (inA)
din.add(SRC, 2 * k, INF);
if (inB)
din.add(2 * k + 1, SNK, INF);
for (int d = 0; d < 4; d++)
{
int u = NB[v][d];
if (u < 0 || label[u] != lab)
continue;
din.add(2 * k + 1, 2 * loc[u], INF);
}
}
int f = din.maxflow(SRC, SNK, budget + 1);
if (f > budget || f == 0)
return res;
vector<char> vis(2 * n + 2, 0);
queue<int> q;
q.push(SRC);
vis[SRC] = 1;
while (!q.empty())
{
int v = q.front();
q.pop();
for (int id : din.g[v])
if (din.es[id].cap > 0 && !vis[din.es[id].to])
{
vis[din.es[id].to] = 1;
q.push(din.es[id].to);
}
}
vector<int> cut;
for (int k = 0; k < n; k++)
if (vis[2 * k] && !vis[2 * k + 1])
cut.push_back(cells[k]);
if ((int)cut.size() != f)
return res;
for (int v : cut)
tr[v] = 0;
vector<char> seen(V, 0);
long long sq = 0;
for (int v0 : cells)
{
if (!tr[v0] || seen[v0])
continue;
long long cnt = 0;
queue<int> qq;
qq.push(v0);
seen[v0] = 1;
while (!qq.empty())
{
int v = qq.front();
qq.pop();
cnt++;
for (int d = 0; d < 4; d++)
{
int u = NB[v][d];
if (u < 0 || !tr[u] || seen[u] || label[u] != lab)
continue;
seen[u] = 1;
qq.push(u);
}
}
sq += cnt * cnt;
}
for (int v : cut)
tr[v] = 1;
res.gain = (double)((long long)n * n - sq) / cut.size();
res.cut = cut;
return res;
}
Split bestSplitOfRegion(int lab, const vector<int> &label, int budget, mt19937 &rng)
{
vector<int> cells;
for (int v = 0; v < V; v++)
if (label[v] == lab)
cells.push_back(v);
Split best;
if (cells.size() < 8)
return best;
vector<pair<int, int>> pairs;
{
vector<int> d0 = bfsDist(cells[rng() % cells.size()], label, lab);
int s = cells[0];
for (int v : cells)
if (d0[v] > d0[s])
s = v;
vector<int> d1 = bfsDist(s, label, lab);
int t = cells[0];
for (int v : cells)
if (d1[v] > d1[t])
t = v;
pairs.push_back({s, t});
int xl = cells[0], xr = cells[0], yl = cells[0], yr = cells[0];
for (int v : cells)
{
if (v % N < xl % N)
xl = v;
if (v % N > xr % N)
xr = v;
if (v / N < yl / N)
yl = v;
if (v / N > yr / N)
yr = v;
}
pairs.push_back({xl, xr});
pairs.push_back({yl, yr});
}
for (auto [s, t] : pairs)
{
if (s == t)
continue;
vector<int> ds = bfsDist(s, label, lab), dt = bfsDist(t, label, lab);
vector<int> sd, td;
for (int v : cells)
{
sd.push_back(ds[v]);
td.push_back(dt[v]);
}
sort(sd.begin(), sd.end());
sort(td.begin(), td.end());
const double fr[][2] = {{0.4, 0.4}, {0.25, 0.25}, {0.15, 0.6}, {0.6, 0.15}, {0.1, 0.75}, {0.75, 0.1}};
for (auto &f : fr)
{
int ra = sd[min((int)sd.size() - 1, (int)(f[0] * sd.size()))];
int rb = td[min((int)td.size() - 1, (int)(f[1] * td.size()))];
Split sp = tryCut(cells, label, lab, ds, dt, ra, rb, budget);
if (sp.gain > best.gain)
best = sp;
}
}
return best;
}
vector<vector<int>> build(const vector<char> &tr0, int budget, double tEnd)
{
tr = tr0;
mt19937 rng(1);
vector<vector<int>> seq;
while (budget > 0 && elapsed() < tEnd)
{
vector<int> label, sizes;
int c = labelComponents(label, sizes);
Split best;
vector<uint64_t> hs(c, 1469598103934665603ULL);
for (int v = 0; v < V; v++)
if (label[v] >= 0)
{
uint64_t &h = hs[label[v]];
h ^= (uint64_t)v + 0x9e37;
h *= 1099511628211ULL;
}
for (int lab = 0; lab < c; lab++)
{
if (sizes[lab] < 8)
continue;
auto it = splitCache.find(hs[lab]);
if (it == splitCache.end())
it = splitCache.emplace(hs[lab], bestSplitOfRegion(lab, label, K, rng)).first;
const Split &sp = it->second;
if ((int)sp.cut.size() > budget)
continue;
if (sp.gain > best.gain)
best = sp;
}
if (best.gain <= 0)
break;
for (int v : best.cut)
tr[v] = 0;
budget -= best.cut.size();
seq.push_back(best.cut);
}
return seq;
}
};
// ---------------- SA with periodic resampling ----------------
vector<int> runSA(vector<int> cuts, double tEnd, uint64_t seedBase, int S, int segments, const Samples *valid = nullptr)
{
if (cuts.empty())
return cuts;
vector<int> bestValid = cuts;
double bestValidE = 1e18;
if (valid)
{
vector<char> t0 = tree0;
for (int v : cuts)
t0[v] = 0;
bestValidE = evaluate(t0, *valid);
}
vector<int> treeList;
for (int v = 0; v < V; v++)
if (tree0[v])
treeList.push_back(v);
RNG rr(seedBase * 7 + 1);
double tStart = elapsed();
if (tEnd - tStart < 0.02)
return cuts;
vector<char> tr = tree0;
for (int v : cuts)
tr[v] = 0;
Samples train;
SAState st;
double T0sa = -1, T1sa = -1, temp = 0;
long long iters = 0, acc = 0;
for (int seg = 0; seg < segments; seg++)
{
double segEnd = tStart + (tEnd - tStart) * (seg + 1) / segments;
train.gen(S, seedBase + 1000 * seg);
st.init(&train, tr);
vector<int> topCells;
int sinceRefresh = 1 << 30;
const int GUIDE = GUIDE_DEF;
if (T0sa < 0)
{
T0sa = max(1.0, (double)st.total / S * 0.002);
T1sa = T0sa * 0.01;
}
while (true)
{
if ((iters & 63) == 0)
{
double t = elapsed();
if (t > segEnd)
break;
double prog = (t - tStart) / (tEnd - tStart);
temp = T0sa * pow(T1sa / T0sa, prog);
}
iters++;
if (GUIDE && sinceRefresh >= GUIDE)
{
// exact single-cut gains on the current training samples via articulation points
APGreedy g;
g.init(&train, st.tr);
vector<int> cand;
for (int v = 0; v < V; v++)
if (st.tr[v] && g.gain[v] > 0)
cand.push_back(v);
int top = min<int>(cand.size(), 64);
partial_sort(cand.begin(), cand.begin() + top, cand.end(), [&](int x, int y)
{ return g.gain[x] > g.gain[y]; });
topCells.assign(cand.begin(), cand.begin() + top);
sinceRefresh = 0;
}
int ka = rr.nextInt(cuts.size());
int a = cuts[ka];
int b;
if (GUIDE && !topCells.empty() && rr.nextInt(10) < 3)
{
b = topCells[rr.nextInt(topCells.size())];
}
else if (rr.nextInt(10) < 7)
{
int r = 1 + rr.nextInt(3);
int ai = a / N, aj = a % N;
int bi = ai + rr.nextInt(2 * r + 1) - r, bj = aj + rr.nextInt(2 * r + 1) - r;
if (bi < 0 || bi >= N || bj < 0 || bj >= N)
continue;
b = bi * N + bj;
}
else
{
b = treeList[rr.nextInt(treeList.size())];
}
if (!st.tr[b] || b == a)
continue;
long long delta = 0;
for (int s = 0; s < S; s++)
delta += st.moveDelta(s, a, b);
double dd = (double)delta / S;
if (dd <= 0 || rr.nextDouble() < exp(-dd / temp))
{
st.applyMove(a, b);
cuts[ka] = b;
acc++;
sinceRefresh++;
}
}
tr = st.tr;
if (valid)
{
double e = evaluate(tr, *valid);
if (e < bestValidE)
{
bestValidE = e;
bestValid = cuts;
}
}
}
cerr << "SA iters=" << iters << " acc=" << acc << "\n";
return valid ? bestValid : cuts;
}
vector<int> hybridInit(const Samples &train, const Samples &test, double tEnd)
{
Walls W;
vector<vector<int>> seq = W.build(tree0, K, elapsed() + (tEnd - elapsed()) * 0.5);
const double fracs[] = {1.0, 0.75, 0.5, 0.25, 0.0};
double bestE = 1e18;
vector<int> best;
for (int fi = 0; fi < 5; fi++)
{
double f = fracs[fi];
double tE = elapsed() + (tEnd - elapsed()) / (5 - fi);
vector<char> tr = tree0;
vector<int> cs;
int used = 0;
for (auto &w : seq)
{
if (used + (int)w.size() > f * K + 1e-9)
break;
for (int v : w)
{
tr[v] = 0;
cs.push_back(v);
}
used += w.size();
}
cs = greedyFill(train, tr, cs, tE);
vector<char> tr2 = tree0;
for (int v : cs)
tr2[v] = 0;
double E = evaluate(tr2, test);
if (E < bestE)
{
bestE = E;
best = cs;
}
}
return best;
}
vector<int> unifiedInit(const Samples &train, double tEnd)
{
vector<char> tr = tree0;
vector<int> cuts;
int budget = K;
APGreedy ap;
ap.init(&train, tr);
WallGen wg;
mt19937 rng(7);
const vector<int> thetas = {1, 5, 6, 7};
vector<Cand> cands;
bool dirty = true;
int singles = 0;
while (budget > 0 && elapsed() < tEnd)
{
if (dirty || singles >= 6)
{
cands = wg.generate(tr, budget, thetas, rng);
long long base = sumSq(tr, train);
for (auto &c : cands)
{
vector<char> t2 = tr;
int cost = 0;
for (int v : c.cells)
if (t2[v])
{
t2[v] = 0;
cost++;
}
c.gain = cost > 0 ? (double)(base - sumSq(t2, train)) / cost : -1;
if (elapsed() > tEnd)
break;
}
dirty = false;
singles = 0;
}
int bs = ap.best();
double gs = bs >= 0 ? ap.gain[bs] : -1;
int bw = -1;
double gw = -1;
for (int i = 0; i < (int)cands.size(); i++)
{
int cost = 0;
for (int v : cands[i].cells)
if (tr[v])
cost++;
if (cost == 0 || cost > budget)
continue;
if (cands[i].gain > gw)
{
gw = cands[i].gain;
bw = i;
}
}
if (bw >= 0 && gw > gs)
{
for (int v : cands[bw].cells)
if (tr[v])
{
ap.cut(v);
tr[v] = 0;
cuts.push_back(v);
budget--;
}
dirty = true;
}
else
{
if (bs < 0)
break;
ap.cut(bs);
tr[bs] = 0;
cuts.push_back(bs);
budget--;
singles++;
}
}
if (budget > 0)
cuts = greedyFill(train, tr, cuts, tEnd + 0.1);
return cuts;
}
int main()
{
T0 = chrono::steady_clock::now();
cin >> N >> K;
V = N * N;
vector<string> G(N);
for (auto &s : G)
cin >> s;
tree0.assign(V, 0);
dry10.assign(V, 0);
NB.assign(V, {-1, -1, -1, -1});
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
{
int v = i * N + j;
dry10[v] = G[i][j] - '0';
tree0[v] = G[i][j] != '0';
for (int d = 0; d < 4; d++)
{
int a = i + DX[d], b = j + DY[d];
if (a >= 0 && a < N && b >= 0 && b < N)
NB[v][d] = a * N + b;
}
}
int nTrees = count(tree0.begin(), tree0.end(), 1);
K = min(K, nTrees);
Samples train, test;
train.gen(48, 777);
test.gen(160, 999);
auto evalCuts = [&](const vector<int> &cs)
{ vector<char> t = tree0; for (int v : cs) t[v] = 0; return evaluate(t, test); };
vector<int> cuts;
if (MODE == 0)
cuts = greedyFill(train, tree0, {}, 0.6);
else if (MODE == 1)
cuts = hybridInit(train, test, F_DEAD);
else if (MODE == 2)
cuts = unifiedInit(train, 1.0);
else if (MODE == 4)
{
vector<int> c0 = greedyFill(train, tree0, {}, 0.25);
vector<int> c1 = hybridInit(train, test, 0.55);
vector<int> c2 = unifiedInit(train, 0.85);
double e0 = evalCuts(c0), e1 = evalCuts(c1), e2 = evalCuts(c2);
cerr << "greedy=" << e0 << " hybrid=" << e1 << " unified=" << e2 << "\n";
cuts = c0;
if (e1 < min(e0, e2))
cuts = c1;
else if (e2 < e0)
cuts = c2;
}
else
{
vector<int> c1 = hybridInit(train, test, H_DEAD1);
vector<int> c2 = unifiedInit(train, H_DEAD2);
double e1 = evalCuts(c1), e2 = evalCuts(c2);
cerr << "hybrid=" << e1 << " unified=" << e2 << "\n";
cuts = (e1 <= e2) ? c1 : c2;
}
double E0 = evalCuts(cuts);
cerr << "init E(test)=" << E0 << " t=" << elapsed() << "\n";
int SS = 48;
int SEG = 4;
vector<int> out = runSA(cuts, TIME_LIMIT - 0.05, 4242, SS, SEG, &test);
cerr << "final E(test)=" << evalCuts(out) << " t=" << elapsed() << "\n";
cout << out.size() << "\n";
for (int v : out)
cout << v / N << " " << v % N << "\n";
}