結果
| 問題 | No.5025 Forest Fire |
| コンテスト | |
| ユーザー |
 mech_39
|
| 提出日時 | 2026-10-21 09:26:19 |
| 言語 | C++23(gcc16) (gcc 16.2.0 + boost 1.92.0 + ACL) |
| 結果 |
AC
不安定
|
| 実行時間 | 1,897 ms / 2,000 ms |
| + 96µs | |
| コード長 | 105,651 bytes |
| 記録 | |
| コンパイル時間 | 8,462 ms |
| コンパイル使用メモリ | 489,908 KB |
| 実行使用メモリ | 16,172 KB |
| スコア | 199,604,765 |
| 最終ジャッジ日時 | 2026-10-24 14:41:17 |
| 合計ジャッジ時間 | 162,453 ms |
|
ジャッジサーバーID （参考情報） |
judge2_0 / <nil> |
| 純コード判定しない問題か言語 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 50 |
ソースコード
#include <bits/stdc++.h>
using namespace std;
#define WSA_S 12
#define S_END 32
#define WSA_SEG 10
#define GREEDY_M 64
#define GREEDY_MC 8
#define P_SNAP 10
#define REVERT 1.03
#define SLACK 4
#define LAMBDA 0.01
#define P_MINCUT 20
#define PCACHE 32
#define MEMO 200000
#define P_LNS 20
#define LNS_M 4
#define TWOSTAGE 0.03
#define WSA_T 0.80
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};
#define TL_DEF 1.95
double TIME_LIMIT = TL_DEF;
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;
struct Cand
{
vector<int> cells;
double proxy = 0;
double gain = 0;
};
namespace FE
{
int W, PV;
int OFF[4];
const int OPP[4] = {1, 0, 3, 2};
vector<int> p2o, o2p;
vector<int> pd10;
vector<double> pp;
vector<int> prow, pcol;
void setup()
{
W = N + 2;
PV = W * W;
OFF[0] = 1;
OFF[1] = -1;
OFF[2] = W;
OFF[3] = -W;
p2o.assign(PV, -1);
o2p.assign(V, -1);
pd10.assign(PV, 0);
pp.assign(PV, 0.0);
prow.assign(PV, 0);
pcol.assign(PV, 0);
for (int P = 0; P < PV; P++)
{
prow[P] = P / W;
pcol[P] = P % W;
}
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
{
int P = (i + 1) * W + (j + 1), o = i * N + j;
p2o[P] = o;
o2p[o] = P;
pd10[P] = tree0[o] ? dry10[o] : 0;
pp[P] = pd10[P] / 10.0;
}
}
}
constexpr uint16_t F_ACT = 0x8000, F_SEP = 0x4000, F_LAB = 0x0FFF;
#define VL 4
typedef int32_t v8i __attribute__((vector_size(4 * VL)));
typedef uint16_t v8u16 __attribute__((vector_size(2 * VL), aligned(2)));
static inline v8i ld8(const uint16_t *p)
{
v8u16 x;
memcpy(&x, p, 2 * VL);
return __builtin_convertvector(x, v8i);
}
struct Engine
{
int S = 0, SW = 0;
vector<uint64_t> dryb;
vector<uint16_t> F, TIN, SUB, SZ;
vector<uint64_t> rows;
vector<char> tr;
vector<long long> sq;
long long total = 0;
long long visitedCells = 0;
struct Tmp
{
int vis;
uint16_t tin, low, sub;
uint8_t rem, pdir;
};
vector<Tmp> tmp;
uint64_t visRows[64];
bool useT1 = false;
int stamp = 0;
vector<int> stk, ord;
inline bool dry(int v, int s) const { return (dryb[(size_t)v * SW + (s >> 6)] >> (s & 63)) & 1; }
void genSamples(int S_, uint64_t seed)
{
using namespace FE;
S = S_;
SW = (S + 63) / 64;
dryb.assign((size_t)PV * SW, 0);
RNG rng(seed);
vector<int> perm(S);
for (int v = 0; v < PV; v++)
{
if (pd10[v] == 0)
continue;
double ex = pd10[v] * S / 10.0;
int k = (int)floor(ex);
if (rng.nextDouble() < ex - k)
k++;
for (int i = 0; i < S; i++)
perm[i] = i;
for (int i = 0; i < k; i++)
{
int j = i + rng.nextInt(S - i);
swap(perm[i], perm[j]);
dryb[(size_t)v * SW + (perm[i] >> 6)] |= 1ULL << (perm[i] & 63);
}
}
}
void build(const vector<char> &tr_)
{
using namespace FE;
tr = tr_;
const size_t n = (size_t)PV * S;
F.assign(n, 0);
TIN.assign(n, 0);
SUB.assign(n, 0);
SZ.assign(n, 0);
zrow.assign(S + 16, 0);
rows.assign((size_t)S * W, 0);
for (int v = 0; v < PV; v++)
if (tr[v])
for (int s = 0; s < S; s++)
if (dry(v, s))
{
F[(size_t)v * S + s] = F_ACT;
rows[(size_t)s * W + prow[v]] |= 1ULL << pcol[v];
}
sq.assign(S, 0);
total = 0;
tmp.assign(PV, Tmp{0, 0, 0, 0, 0, 0});
stamp = 0;
stk.assign(PV + 1, 0);
ord.reserve(PV);
for (int s = 0; s < S; s++)
{
stamp++;
for (int v = 0; v < PV; v++)
if ((F[(size_t)v * S + s] & F_ACT) && tmp[v].vis != stamp)
{
long long c = tarjan(s, v);
sq[s] += c * c;
}
total += sq[s];
}
}
inline unsigned nbrMask(const uint64_t *R, int v) const
{
using namespace FE;
const int r = prow[v], c = pcol[v];
const uint64_t mid = R[r] >> (c - 1);
return (unsigned)(((mid >> 2) & 1) | ((mid & 1) << 1) | (((R[r + 1] >> c) & 1) << 2) | (((R[r - 1] >> c) & 1) << 3));
}
int tarjan(int s, int r)
{
using namespace FE;
const uint64_t *R = rows.data() + (size_t)s * W;
Tmp *T = tmp.data();
int t = 0, top = 0;
ord.clear();
{
Tmp &X = T[r];
X.vis = stamp;
X.tin = 0;
X.low = 0;
X.sub = 1;
X.pdir = 0;
X.rem = (uint8_t)nbrMask(R, r);
}
stk[top++] = r;
ord.push_back(r);
while (top)
{
const int v = stk[top - 1];
Tmp &X = T[v];
if (X.rem)
{
const int d = __builtin_ctz(X.rem);
X.rem &= X.rem - 1;
const int u = v + OFF[d];
Tmp &U = T[u];
if (U.vis != stamp)
{
++t;
U.vis = stamp;
U.tin = (uint16_t)t;
U.low = (uint16_t)t;
U.sub = 1;
U.pdir = (uint8_t)OPP[d];
U.rem = (uint8_t)(nbrMask(R, u) & ~(1u << OPP[d]));
stk[top++] = u;
ord.push_back(u);
}
else if (U.tin < X.low)
X.low = U.tin;
}
else
{
top--;
if (top)
{
Tmp &P = T[stk[top - 1]];
P.sub += X.sub;
if (X.low < P.low)
P.low = X.low;
}
}
}
const int c = t + 1;
const size_t SS = S;
const uint16_t lab = (uint16_t)r;
for (int v : ord)
{
const Tmp &X = T[v];
uint16_t f = F_ACT | lab;
if (v != r)
{
f |= (uint16_t)(X.pdir << 12);
const Tmp &P = T[v + OFF[X.pdir]];
if (X.low >= P.tin)
f |= F_SEP;
}
const size_t idx = (size_t)v * SS + s;
F[idx] = f;
TIN[idx] = X.tin;
SUB[idx] = X.sub;
SZ[idx] = (uint16_t)c;
}
visitedCells += c;
return c;
}
int tarjan2(int s, int r)
{
using namespace FE;
const uint64_t *R = rows.data() + (size_t)s * W;
Tmp *T = tmp.data();
const size_t SS = S;
uint16_t *Fp = F.data() + s;
uint16_t *Tp = TIN.data() + s;
uint16_t *Up = SUB.data() + s;
int t = 0, top = 0;
const int ordStart = (int)ord.size();
const uint16_t lab = (uint16_t)r;
visRows[prow[r]] |= 1ULL << pcol[r];
T[r].tin = 0;
T[r].low = 0;
T[r].sub = 1;
T[r].pdir = 4;
stk[top++] = r;
ord.push_back(r);
while (top)
{
const int v = stk[top - 1];
const unsigned m = nbrMask(R, v) & ~nbrMask(visRows, v);
if (m)
{
const int d = __builtin_ctz(m);
const int u = v + OFF[d];
visRows[prow[u]] |= 1ULL << pcol[u];
++t;
Tmp &U = T[u];
U.tin = (uint16_t)t;
U.sub = 1;
U.pdir = (uint8_t)OPP[d];
unsigned bm = nbrMask(R, u) & nbrMask(visRows, u);
uint16_t lo = (uint16_t)t;
while (bm)
{
int e = __builtin_ctz(bm);
bm &= bm - 1;
uint16_t x = T[u + OFF[e]].tin;
if (x < lo)
lo = x;
}
U.low = lo;
stk[top++] = u;
ord.push_back(u);
}
else
{
top--;
Tmp &X = T[v];
uint16_t f = F_ACT | lab;
if (top)
{
Tmp &P = T[stk[top - 1]];
P.sub += X.sub;
if (X.low < P.low)
P.low = X.low;
f |= (uint16_t)(X.pdir << 12);
if (X.low >= P.tin)
f |= F_SEP;
}
const size_t idx = (size_t)v * SS;
Fp[idx] = f;
Tp[idx] = X.tin;
Up[idx] = X.sub;
}
}
const int c = t + 1;
uint16_t *Zp = SZ.data() + s;
for (int k = ordStart; k < (int)ord.size(); k++)
Zp[(size_t)ord[k] * SS] = (uint16_t)c;
visitedCells += c;
return c;
}
vector<uint16_t> zrow;
void swapParts(int a, int b, long long &SMa, long long &SMb) const
{
using namespace FE;
const size_t SS = S;
const uint16_t *Z = zrow.data();
const uint16_t *bF = &F[b * SS], *bT = &TIN[b * SS], *bZ = &SZ[b * SS];
const uint16_t *nF[4], *nT[4], *nU[4], *nZ[4];
int bo[4];
for (int d = 0; d < 4; d++)
{
int u = b + OFF[d];
bo[d] = 0x8000 | 0x4000 | (OPP[d] << 12);
if (tr[u])
{
nF[d] = &F[u * SS];
nT[d] = &TIN[u * SS];
nU[d] = &SUB[u * SS];
nZ[d] = &SZ[u * SS];
}
else
{
nF[d] = nT[d] = nU[d] = nZ[d] = Z;
}
}
const uint16_t *aF[4], *aT[4], *aZ[4];
for (int d = 0; d < 4; d++)
{
int u = a + OFF[d];
if (a >= 0 && tr[u] && u != b)
{
aF[d] = &F[u * SS];
aT[d] = &TIN[u * SS];
aZ[d] = &SZ[u * SS];
}
else
{
aF[d] = aT[d] = aZ[d] = Z;
}
}
const v8i ACTv = (v8i){} + 0x8000, LABv = (v8i){} + 0x0FFF, ONE = (v8i){} + 1, ZERO = (v8i){};
const v8i BO0 = (v8i){} + bo[0], BO1 = (v8i){} + bo[1], BO2 = (v8i){} + bo[2], BO3 = (v8i){} + bo[3];
v8i accA = ZERO, accB = ZERO;
for (int s = 0; s < S; s += VL)
{
const v8i Bf = ld8(bF + s), Bt = ld8(bT + s), Bz = ld8(bZ + s);
const v8i bd = (Bf & ACTv) != ZERO;
const v8i lb = Bf & LABv;
v8i lo[4], hi[4];
v8i sepSum = ZERO, sep2 = ZERO;
v8i wk[4], wz[4];
#define CHILD(j, BOj) \
{ \
const v8i f = ld8(nF[j] + s), tn = ld8(nT[j] + s), su = ld8(nU[j] + s); \
const v8i isC = (f == (BOj | lb)) & (tn > Bt) & bd; \
const v8i sm = su & isC; \
lo[j] = tn & isC; \
hi[j] = (tn + su) & isC; \
sepSum += sm; \
sep2 += sm * sm; \
const v8i act = (f & ACTv) != ZERO; \
wk[j] = (f & LABv & act) | (~act & ((v8i){} + (0x2000 + j))); \
wz[j] = ld8(nZ[j] + s) & act; \
}
CHILD(0, BO0)
CHILD(1, BO1)
CHILD(2, BO2)
CHILD(3, BO3)
#undef CHILD
const v8i rest = Bz - ONE - sepSum;
const v8i mbDry = Bz * Bz - sep2 - rest * rest;
v8i sum = wz[0], sum2 = wz[0] * wz[0];
{
v8i z1 = wz[1] & (wk[1] != wk[0]);
sum += z1;
sum2 += z1 * z1;
v8i z2 = wz[2] & (wk[2] != wk[0]) & (wk[2] != wk[1]);
sum += z2;
sum2 += z2 * z2;
v8i z3 = wz[3] & (wk[3] != wk[0]) & (wk[3] != wk[1]) & (wk[3] != wk[2]);
sum += z3;
sum2 += z3 * z3;
}
const v8i mbWet = (ONE + sum) * (ONE + sum) - sum2;
accB += (mbDry & bd) | (mbWet & ~bd);
v8i key[4], zs[4];
#define ANB(k) \
{ \
const v8i f = ld8(aF[k] + s), tn = ld8(aT[k] + s), zz = ld8(aZ[k] + s); \
const v8i act = (f & ACTv) != ZERO; \
const v8i lab = f & LABv; \
const v8i inB = (lab == lb) & bd; \
v8i pk = (v8i){} + 4100, pz = rest; \
for (int j = 0; j < 4; j++) \
{ \
const v8i in = (tn >= lo[j]) & (tn < hi[j]); \
pk = (in & ((v8i){} + (4096 + j))) | (~in & pk); \
pz = (in & (hi[j] - lo[j])) | (~in & pz); \
} \
const v8i kk = (inB & pk) | (~inB & lab); \
const v8i z = (inB & pz) | (~inB & zz); \
key[k] = (kk & act) | (~act & ((v8i){} + (0x3000 + k))); \
zs[k] = z & act; \
}
ANB(0)
ANB(1)
ANB(2)
ANB(3)
#undef ANB
sum = zs[0];
sum2 = zs[0] * zs[0];
{
v8i z1 = zs[1] & (key[1] != key[0]);
sum += z1;
sum2 += z1 * z1;
v8i z2 = zs[2] & (key[2] != key[0]) & (key[2] != key[1]);
sum += z2;
sum2 += z2 * z2;
v8i z3 = zs[3] & (key[3] != key[0]) & (key[3] != key[1]) & (key[3] != key[2]);
sum += z3;
sum2 += z3 * z3;
}
accA += (ONE + sum) * (ONE + sum) - sum2;
}
long long ra = 0, rb = 0;
for (int i = 0; i < VL; i++)
{
ra += accA[i];
rb += accB[i];
}
SMa = ra;
SMb = rb;
}
inline double swapDelta(int a, int b) const
{
long long ma, mb;
swapParts(a, b, ma, mb);
return (FE::pp[a] * ma - FE::pp[b] * mb) / S;
}
void commit(int a, int b)
{
using namespace FE;
const size_t SS = S;
if (a >= 0)
tr[a] = 1;
if (b >= 0)
tr[b] = 0;
for (int s = 0; s < S; s++)
{
const bool da = a >= 0 && dry(a, s);
const bool bb = b >= 0 && (F[b * SS + s] & F_ACT);
if (!da && !bb)
continue;
int labs[5], nl = 0;
long long old = 0;
auto addOld = [&](int u)
{
const uint16_t f = F[u * SS + s];
if (!(f & F_ACT))
return;
int l = f & F_LAB;
for (int k = 0; k < nl; k++)
if (labs[k] == l)
return;
labs[nl++] = l;
long long z = SZ[u * SS + s];
old += z * z;
};
if (da)
for (int d = 0; d < 4; d++)
addOld(a + OFF[d]);
if (bb)
addOld(b);
uint64_t *R = rows.data() + (size_t)s * W;
if (bb)
{
const size_t ib = b * SS + s;
F[ib] = 0;
TIN[ib] = 0;
SUB[ib] = 0;
SZ[ib] = 0;
R[prow[b]] &= ~(1ULL << pcol[b]);
}
if (da)
R[prow[a]] |= 1ULL << pcol[a];
long long nw = 0;
stamp++;
if (da)
{
long long c = tarjan(s, a);
nw += c * c;
}
if (bb)
for (int d = 0; d < 4; d++)
{
int u = b + OFF[d];
if ((F[u * SS + s] & F_ACT) && tmp[u].vis != stamp)
{
long long c = tarjan(s, u);
nw += c * c;
}
}
sq[s] += nw - old;
total += nw - old;
}
}
double E() const { return (double)total / S / ((double)N * N); }
void cutGains(vector<double> &g) const
{
using namespace FE;
g.assign(PV, -1e18);
for (int v = 0; v < PV; v++)
{
if (!tr[v])
continue;
long long ma, mb;
swapParts(-1, v, ma, mb);
g[v] = pp[v] * mb / S;
}
}
double cutGain(int v) const
{
long long ma, mb;
swapParts(-1, v, ma, mb);
return FE::pp[v] * mb / S;
}
double restoreCost(int a) const
{
using namespace FE;
const size_t SS = S;
int an[4], na = 0;
for (int d = 0; d < 4; d++)
{
int u = a + OFF[d];
if (tr[u])
an[na++] = u;
}
long long sm = 0;
for (int s = 0; s < S; s++)
{
int keys[4], nk = 0;
long long sum = 0, sum2 = 0;
for (int j = 0; j < na; j++)
{
const uint16_t f = F[an[j] * SS + s];
if (!(f & F_ACT))
continue;
int key = f & F_LAB;
bool dup = false;
for (int q = 0; q < nk; q++)
if (keys[q] == key)
dup = true;
if (dup)
continue;
keys[nk++] = key;
long long z = SZ[an[j] * SS + s];
sum += z;
sum2 += z * z;
}
sm += (1 + sum) * (1 + sum) - sum2;
}
return pp[a] * sm / S;
}
};
struct Validator
{
int S = 0;
vector<uint64_t> dryRows;
void gen(int S_, uint64_t seed)
{
using namespace FE;
S = S_;
dryRows.assign((size_t)S * W, 0);
RNG rng(seed);
vector<int> perm(S);
for (int v = 0; v < PV; v++)
{
if (pd10[v] == 0)
continue;
double ex = pd10[v] * S / 10.0;
int k = (int)floor(ex);
if (rng.nextDouble() < ex - k)
k++;
for (int i = 0; i < S; i++)
perm[i] = i;
for (int i = 0; i < k; i++)
{
int j = i + rng.nextInt(S - i);
swap(perm[i], perm[j]);
dryRows[(size_t)perm[i] * W + prow[v]] |= 1ULL << pcol[v];
}
}
}
double meanSq(const vector<char> &trP) const
{
using namespace FE;
vector<uint64_t> trRow(W, 0);
for (int v = 0; v < PV; v++)
if (trP[v])
trRow[prow[v]] |= 1ULL << pcol[v];
static vector<int> par, sz;
static vector<int> rs, re, rid;
par.resize(64 * W);
sz.resize(64 * W);
long long tot = 0;
for (int s = 0; s < S; s++)
{
const uint64_t *D = dryRows.data() + (size_t)s * W;
int nNodes = 0;
int pS[64], pE[64], pI[64], np = 0;
int cS[64], cE[64], cI[64];
for (int r = 1; r <= N; r++)
{
uint64_t m = D[r] & trRow[r];
int nc = 0;
while (m)
{
int st = __builtin_ctzll(m);
uint64_t x = ~(m >> st);
int len = x ? __builtin_ctzll(x) : 64 - st;
int id = nNodes++;
par[id] = id;
sz[id] = len;
cS[nc] = st;
cE[nc] = st + len;
cI[nc] = id;
nc++;
m &= ~(((len >= 64) ? ~0ULL : ((1ULL << len) - 1)) << st);
}
int p = 0;
for (int k = 0; k < nc; k++)
{
while (p < np && pE[p] <= cS[k])
p++;
int q = p;
while (q < np && pS[q] < cE[k])
{
int x = pI[q], y = cI[k];
while (par[x] != x)
{
par[x] = par[par[x]];
x = par[x];
}
while (par[y] != y)
{
par[y] = par[par[y]];
y = par[y];
}
if (x != y)
{
if (sz[x] < sz[y])
swap(x, y);
par[y] = x;
sz[x] += sz[y];
}
q++;
}
if (q > p)
p = q - 1;
}
np = nc;
for (int k = 0; k < nc; k++)
{
pS[k] = cS[k];
pE[k] = cE[k];
pI[k] = cI[k];
}
}
long long q = 0;
for (int i = 0; i < nNodes; i++)
if (par[i] == i)
q += (long long)sz[i] * sz[i];
tot += q;
}
return (double)tot / S;
}
double E(const vector<char> &trP) const { return meanSq(trP) / ((double)N * N); }
};
struct RegionGain
{
vector<int> isW, vis, stk;
int wst = 0, vst = 0;
double gain(const Engine &E, const vector<int> &cells, int s0 = 0, int s1 = -1)
{
using namespace FE;
if (s1 < 0)
s1 = E.S;
if (isW.empty())
{
isW.assign(PV, 0);
vis.assign(PV, 0);
stk.assign(PV + 1, 0);
}
++wst;
int wc[64], nw = 0;
for (int c : cells)
if (E.tr[c] && isW[c] != wst && nw < 64)
{
isW[c] = wst;
wc[nw++] = c;
}
if (nw == 0)
return 0;
const size_t SS = E.S;
long long tot = 0;
for (int s = s0; s < s1; s++)
{
const uint64_t *R = E.rows.data() + (size_t)s * W;
int labs[64], nl = 0;
long long old = 0;
for (int k = 0; k < nw; k++)
{
const uint16_t f = E.F[wc[k] * SS + s];
if (!(f & F_ACT))
continue;
int l = f & F_LAB;
bool dup = false;
for (int q = 0; q < nl; q++)
if (labs[q] == l)
{
dup = true;
break;
}
if (dup)
continue;
labs[nl++] = l;
long long z = E.SZ[wc[k] * SS + s];
old += z * z;
}
if (nl == 0)
continue;
++vst;
long long nwq = 0;
for (int k = 0; k < nw; k++)
{
const int c = wc[k];
if (!(E.F[c * SS + s] & F_ACT))
continue;
for (int d = 0; d < 4; d++)
{
const int u0 = c + OFF[d];
if (!((R[prow[u0]] >> pcol[u0]) & 1) || isW[u0] == wst || vis[u0] == vst)
continue;
long long cnt = 0;
int top = 0;
stk[top++] = u0;
vis[u0] = vst;
while (top)
{
const int x = stk[--top];
cnt++;
for (int e = 0; e < 4; e++)
{
const int u = x + OFF[e];
if (((R[prow[u]] >> pcol[u]) & 1) && isW[u] != wst && vis[u] != vst)
{
vis[u] = vst;
stk[top++] = u;
}
}
}
nwq += cnt * cnt;
}
}
tot += old - nwq;
}
return (double)tot / (s1 - s0);
}
};
struct PolishStats
{
int rounds = 0, applied = 0;
};
int polishSwaps(Engine &E, vector<int> &cuts, double tEnd, int topA, int topB, double eps, PolishStats *ps = nullptr)
{
using namespace FE;
vector<double> g, rc(cuts.size());
vector<int> ia(cuts.size()), cand;
int applied = 0;
while (elapsed() < tEnd)
{
E.cutGains(g);
for (size_t k = 0; k < cuts.size(); k++)
rc[k] = E.restoreCost(cuts[k]);
for (size_t k = 0; k < cuts.size(); k++)
ia[k] = k;
int na = min<int>(topA, cuts.size());
partial_sort(ia.begin(), ia.begin() + na, ia.end(), [&](int x, int y)
{ return rc[x] < rc[y]; });
cand.clear();
for (int v = 0; v < PV; v++)
if (E.tr[v] && g[v] > 0)
cand.push_back(v);
int nb = min<int>(topB, cand.size());
partial_sort(cand.begin(), cand.begin() + nb, cand.end(), [&](int x, int y)
{ return g[x] > g[y]; });
double best = -eps;
int bk = -1, bb = -1;
for (int i = 0; i < na; i++)
for (int j = 0; j < nb; j++)
{
const int a = cuts[ia[i]], b = cand[j];
double d = E.swapDelta(a, b);
if (d < best)
{
best = d;
bk = ia[i];
bb = b;
}
}
if (ps)
ps->rounds++;
if (bk < 0)
break;
E.commit(cuts[bk], bb);
cuts[bk] = bb;
applied++;
if (ps)
ps->applied++;
}
return applied;
}
struct FastFlow
{
int nN = 0, nE = 0;
vector<int> head, nxt, to, cap, level, it, q;
void init(int n)
{
nN = n;
nE = 0;
head.assign(n, -1);
if ((int)level.size() < n)
{
level.resize(n);
it.resize(n);
q.resize(n);
}
}
inline void add(int a, int b, int c)
{
if (nE + 2 > (int)to.size())
{
int m = max(1024, (int)to.size() * 2);
to.resize(m);
nxt.resize(m);
cap.resize(m);
}
to[nE] = b;
cap[nE] = c;
nxt[nE] = head[a];
head[a] = nE++;
to[nE] = a;
cap[nE] = 0;
nxt[nE] = head[b];
head[b] = nE++;
}
bool bfs(int s, int t)
{
for (int i = 0; i < nN; i++)
level[i] = -1;
int qh = 0, qt = 0;
q[qt++] = s;
level[s] = 0;
while (qh < qt)
{
int v = q[qh++];
for (int e = head[v]; e >= 0; e = nxt[e])
if (cap[e] > 0 && level[to[e]] < 0)
{
level[to[e]] = level[v] + 1;
q[qt++] = to[e];
}
}
return level[t] >= 0;
}
int dfs(int v, int t, int f)
{
if (v == t)
return f;
for (int &e = it[v]; e >= 0; e = nxt[e])
{
int u = to[e];
if (cap[e] > 0 && level[u] == level[v] + 1)
{
int d = dfs(u, t, min(f, cap[e]));
if (d > 0)
{
cap[e] -= d;
cap[e ^ 1] += d;
return d;
}
}
}
return 0;
}
int maxflow(int s, int t, int limit)
{
int flow = 0;
while (flow < limit && bfs(s, t))
{
for (int i = 0; i < nN; i++)
it[i] = head[i];
int f;
while ((f = dfs(s, t, INT_MAX)) > 0)
{
flow += f;
if (flow >= limit)
break;
}
}
return flow;
}
void reach(int s, vector<char> &vis)
{
vis.assign(nN, 0);
int qh = 0, qt = 0;
q[qt++] = s;
vis[s] = 1;
while (qh < qt)
{
int v = q[qh++];
for (int e = head[v]; e >= 0; e = nxt[e])
if (cap[e] > 0 && !vis[to[e]])
{
vis[to[e]] = 1;
q[qt++] = to[e];
}
}
}
};
struct WallGen2
{
vector<int> loc;
FastFlow din;
vector<char> vis, isCut, seen;
vector<int> st, bq;
WallGen2() { loc.assign(V, -1); }
vector<int> bfsDist(int src, const vector<int> &label, int lab)
{
vector<int> d(V, -1);
bq.clear();
bq.push_back(src);
d[src] = 0;
for (size_t h = 0; h < bq.size(); h++)
{
int v = bq[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;
bq.push_back(u);
}
}
return d;
}
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;
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;
}
din.reach(SRC, vis);
vector<int> cut;
for (int k = 0; k < n; k++)
if (vis[2 * k] && !vis[2 * k + 1])
cut.push_back(cells[k]);
isCut.assign(n, 0);
for (int v : cut)
isCut[loc[v]] = 1;
seen.assign(n, 0);
long long sq = 0;
for (int k0 = 0; k0 < n; k0++)
{
if (isCut[k0] || seen[k0])
continue;
long long cnt = 0;
st.clear();
st.push_back(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; });
#define WG_KEEP 3
for (int i = 0; i < (int)local.size() && i < WG_KEEP; i++)
out.push_back(local[i]);
}
map<uint64_t, vector<Cand>> cache;
vector<Cand> generate(const vector<char> &tr, int budget, const vector<int> &thetas, mt19937 &rng)
{
vector<Cand> out;
vector<int> label(V, -1);
for (int th : thetas)
{
fill(label.begin(), label.end(), -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);
}
}
#define WG_MINREG 40
if (cells.size() >= WG_MINREG)
{
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;
}
};
#define INIT_REGEN 6
#define INIT_SCREEN_DIV 4
#define INIT_SCREEN 6
#define INIT_FULLG 8
long long nWallEval = 0;
vector<int> unifiedInit2(Engine &eng, double tEnd, int &nWalls)
{
using namespace FE;
vector<char> trO = tree0;
vector<int> cuts;
int budget = K;
static WallGen2 wg;
mt19937 rng(7);
#define INIT_THETAS 1, 5, 6, 7
const vector<int> thetas = {INIT_THETAS};
vector<Cand> cands;
vector<double> wgain;
bool dirty = true, regen = true;
int singles = 0;
RegionGain rg;
vector<double> g;
eng.cutGains(g);
nWalls = 0;
auto cutCell = [&](int o)
{
int P = o2p[o];
eng.commit(-1, P);
trO[o] = 0;
cuts.push_back(o);
budget--;
};
vector<char> fresh;
int sinceFull = 0, lastCutFrom = 0;
while (budget > 0 && elapsed() < tEnd)
{
if (dirty || singles >= INIT_REGEN)
{
if (regen)
{
cands = wg.generate(trO, budget, thetas, rng);
wgain.assign(cands.size(), 1e18);
fresh.assign(cands.size(), 0);
}
regen = false;
{
const int s1 = max(8, eng.S / INIT_SCREEN_DIV);
vector<pair<double, int>> sc;
for (size_t i = 0; i < cands.size(); i++)
{
vector<int> cp;
for (int v : cands[i].cells)
if (trO[v])
cp.push_back(o2p[v]);
int cost = cp.size();
if (cost == 0 || cost > budget)
{
wgain[i] = -1;
fresh[i] = 1;
continue;
}
double gq = rg.gain(eng, cp, 0, s1) / cost;
sc.push_back({gq, (int)i});
wgain[i] = -1;
fresh[i] = 1;
}
int nt = min<int>(INIT_SCREEN, sc.size());
partial_sort(sc.begin(), sc.begin() + nt, sc.end(), greater<pair<double, int>>());
for (int k = 0; k < nt; k++)
{
int i = sc[k].second;
vector<int> cp;
for (int v : cands[i].cells)
if (trO[v])
cp.push_back(o2p[v]);
wgain[i] = rg.gain(eng, cp) / cp.size();
}
}
dirty = false;
singles = 0;
}
int bs = -1;
double gs = -1e18;
for (int P = 0; P < PV; P++)
if (eng.tr[P] && g[P] > gs)
{
gs = g[P];
bs = P;
}
int bw = -1;
double gw = -1;
while (true)
{
bw = -1;
gw = -1;
for (int i = 0; i < (int)cands.size(); i++)
{
if (wgain[i] <= gw)
continue;
int cost = 0;
for (int v : cands[i].cells)
if (trO[v])
cost++;
if (cost == 0 || cost > budget)
{
wgain[i] = -1;
fresh[i] = 1;
continue;
}
gw = wgain[i];
bw = i;
}
if (bw < 0 || fresh[bw])
break;
vector<int> cp;
for (int v : cands[bw].cells)
if (trO[v])
cp.push_back(o2p[v]);
wgain[bw] = rg.gain(eng, cp) / cp.size();
fresh[bw] = 1;
nWallEval++;
}
#define INIT_WALL_BONUS 2.0
const double wbonus = INIT_WALL_BONUS;
if (bw >= 0 && gw * wbonus > gs)
{
for (int v : cands[bw].cells)
if (trO[v])
cutCell(v);
dirty = true;
regen = true;
nWalls++;
}
else
{
if (bs < 0)
break;
cutCell(p2o[bs]);
singles++;
}
if (++sinceFull >= INIT_FULLG || dirty)
{
eng.cutGains(g);
sinceFull = 0;
}
else
{
static vector<int> idx;
idx.clear();
for (int k = lastCutFrom; k < (int)cuts.size(); k++)
{
const int c = o2p[cuts[k]];
for (int di = -3; di <= 3; di++)
for (int dj = -3; dj <= 3; dj++)
{
int r = prow[c] + di, q = pcol[c] + dj;
if (r < 1 || r > N || q < 1 || q > N)
continue;
int P = r * W + q;
if (eng.tr[P])
idx.push_back(P);
}
}
vector<int> top;
for (int P = 0; P < PV; P++)
if (eng.tr[P] && g[P] > 0)
top.push_back(P);
int nt = min<int>(48, top.size());
partial_sort(top.begin(), top.begin() + nt, top.end(), [&](int x, int y)
{ return g[x] > g[y]; });
for (int k = 0; k < nt; k++)
idx.push_back(top[k]);
for (int P : idx)
g[P] = eng.cutGain(P);
for (int k = lastCutFrom; k < (int)cuts.size(); k++)
g[o2p[cuts[k]]] = -1e18;
}
lastCutFrom = cuts.size();
}
while (budget > 0)
{
int bs = -1;
double gs = -1e18;
for (int P = 0; P < PV; P++)
if (eng.tr[P] && g[P] > gs)
{
gs = g[P];
bs = P;
}
if (bs < 0)
break;
cutCell(p2o[bs]);
if (elapsed() < tEnd + 0.1)
eng.cutGains(g);
else
g[bs] = -1e18;
}
return cuts;
}
#define SA_SEG 8
#define SA_S_LIST 64, 64, 96, 96, 128, 192, 256, 384
#define SA_T0 0.002
#define SA_T1R 0.01
#define SA_GUIDE 20
#define VALID_S 512
double evalOn(Engine &E, const vector<char> &trP)
{
E.build(trP);
return E.E();
}
struct SAStats
{
long long iters = 0, acc = 0, polished = 0;
double validE = 0;
};
vector<int> runSA2(vector<int> cutsO, double tEnd, uint64_t seedBase, SAStats &st)
{
using namespace FE;
if (cutsO.empty())
return cutsO;
vector<int> cuts;
for (int v : cutsO)
cuts.push_back(o2p[v]);
vector<char> tr(PV, 0);
for (int v = 0; v < PV; v++)
tr[v] = pd10[v] > 0;
for (int v : cuts)
tr[v] = 0;
vector<int> treeList;
for (int v = 0; v < PV; v++)
if (pd10[v] > 0)
treeList.push_back(v);
Validator valid;
valid.gen(VALID_S, seedBase * 31 + 7);
double tv0 = elapsed();
double bestValidE = valid.E(tr);
double tValid = elapsed() - tv0;
tEnd -= tValid * 1.5;
vector<int> bestCuts = cuts;
RNG rr(seedBase * 7 + 1);
const double tStart = elapsed();
Engine eng;
double T0sa = -1, T1sa = -1, temp = 0;
const int segments = SA_SEG;
for (int seg = 0; seg < segments; seg++)
{
const double segEnd = tStart + (tEnd - tStart) * (seg + 1) / segments;
static const int SLIST[] = {SA_S_LIST};
const int Sseg = SLIST[min<int>(seg, sizeof(SLIST) / sizeof(int) - 1)];
eng.genSamples(Sseg, seedBase + 1000 * seg + 17);
eng.build(tr);
if (T0sa < 0)
{
T0sa = max(1.0, (double)eng.total / eng.S * SA_T0);
T1sa = T0sa * SA_T1R;
}
vector<int> topCells, lowCuts;
int sinceRefresh = 1 << 30;
vector<double> g;
#define SA_POLISH_FRAC 0.1
#define SA_PAIR 0
const double saEnd = segEnd - (segEnd - (seg == 0 ? tStart : tStart + (tEnd - tStart) * seg / segments)) * SA_POLISH_FRAC;
while (true)
{
if ((st.iters & 63) == 0)
{
double t = elapsed();
if (t > saEnd)
break;
double prog = (t - tStart) / (tEnd - tStart);
temp = T0sa * pow(T1sa / T0sa, prog);
}
st.iters++;
if (SA_GUIDE && sinceRefresh >= SA_GUIDE)
{
eng.cutGains(g);
vector<int> cand;
for (int v : treeList)
if (eng.tr[v] && g[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[x] > g[y]; });
topCells.assign(cand.begin(), cand.begin() + top);
if (SA_PAIR)
{
vector<double> rc(cuts.size());
lowCuts.resize(cuts.size());
for (size_t k = 0; k < cuts.size(); k++)
{
rc[k] = eng.restoreCost(cuts[k]);
lowCuts[k] = k;
}
int nl = min<int>(16, cuts.size());
partial_sort(lowCuts.begin(), lowCuts.begin() + nl, lowCuts.end(), [&](int x, int y)
{ return rc[x] < rc[y]; });
lowCuts.resize(nl);
}
sinceRefresh = 0;
}
int ka = rr.nextInt(cuts.size());
int b;
const int rsel = rr.nextInt(10);
if (SA_PAIR && !lowCuts.empty() && !topCells.empty() && rsel < SA_PAIR)
{
ka = lowCuts[rr.nextInt(lowCuts.size())];
b = topCells[rr.nextInt(topCells.size())];
}
else if (SA_GUIDE && !topCells.empty() && rsel < 3 + SA_PAIR)
{
b = topCells[rr.nextInt(topCells.size())];
}
else if (rr.nextInt(10) < 7)
{
#define SA_NEAR_R 3
int r = 1 + rr.nextInt(SA_NEAR_R);
int bi = prow[cuts[ka]] + rr.nextInt(2 * r + 1) - r, bj = pcol[cuts[ka]] + rr.nextInt(2 * r + 1) - r;
if (bi < 1 || bi > N || bj < 1 || bj > N)
continue;
b = bi * W + bj;
}
else
{
b = treeList[rr.nextInt(treeList.size())];
}
const int a = cuts[ka];
if (!eng.tr[b] || b == a)
continue;
const double dd = eng.swapDelta(a, b);
if (dd <= 0 || rr.nextDouble() < exp(-dd / temp))
{
eng.commit(a, b);
cuts[ka] = b;
st.acc++;
sinceRefresh++;
}
}
{
PolishStats ps;
polishSwaps(eng, cuts, segEnd, 16, 48, 0.0, &ps);
st.polished += ps.applied;
}
tr = eng.tr;
double e = valid.E(tr);
if (e < bestValidE)
{
bestValidE = e;
bestCuts = cuts;
}
}
st.validE = bestValidE;
vector<int> out;
for (int v : bestCuts)
out.push_back(p2o[v]);
return out;
}
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];
}
};
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;
}
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;
}
};
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;
}
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]);
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]);
}
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;
}
};
void genBits(int S, uint64_t seed, vector<uint64_t> &bits)
{
bits.assign((size_t)S * N, 0);
RNG rng(seed);
for (int s = 0; s < S; s++)
for (int v = 0; v < V; v++)
if ((rng.next() % 10) < (uint64_t)dry10[v])
bits[(size_t)s * N + v / N] |= 1ULL << (v % N);
}
struct FastMC
{
int S = 0;
vector<uint64_t> bits;
void gen(int S_, uint64_t seed)
{
S = S_;
genBits(S, seed, bits);
}
int par[1600], sz[1600];
inline int fnd(int x)
{
while (par[x] != x)
{
par[x] = par[par[x]];
x = par[x];
}
return x;
}
long long sampleSq(int s, const uint64_t *T)
{
const uint64_t *B = &bits[(size_t)s * N];
int nr = 0, ps = 0;
uint64_t pm = 0, pst = 0;
long long q = 0;
for (int r = 0; r < N; r++)
{
uint64_t m = B[r] & T[r];
uint64_t st = m & ~(m << 1), en = m & ~(m >> 1);
int s0 = nr;
uint64_t st2 = st;
while (st2)
{
int l = __builtin_ctzll(st2), h = __builtin_ctzll(en);
st2 &= st2 - 1;
en &= en - 1;
int len = h - l + 1;
par[nr] = nr;
sz[nr] = len;
q += (long long)len * len;
nr++;
}
uint64_t ov = m & pm;
uint64_t os = ov & ~(ov << 1);
while (os)
{
int c = __builtin_ctzll(os);
os &= os - 1;
uint64_t lowm = (c == 63) ? ~0ULL : ((2ULL << c) - 1);
int i = ps + __builtin_popcountll(pst & lowm) - 1;
int j = s0 + __builtin_popcountll(st & lowm) - 1;
int a = fnd(i), b = fnd(j);
if (a != b)
{
if (sz[a] < sz[b])
swap(a, b);
q += 2LL * sz[a] * sz[b];
par[b] = a;
sz[a] += sz[b];
}
}
ps = s0;
pm = m;
pst = st;
}
return q;
}
long long sumSq(const uint64_t *T)
{
long long tot = 0;
for (int s = 0; s < S; s++)
tot += sampleSq(s, T);
return tot;
}
double E(const uint64_t *T) { return (double)sumSq(T) / S / ((double)N * N); }
};
void masksFromTree(const vector<char> &tr, uint64_t *T)
{
for (int r = 0; r < N; r++)
{
uint64_t m = 0;
for (int c = 0; c < N; c++)
if (tr[r * N + c])
m |= 1ULL << c;
T[r] = m;
}
}
struct Wall
{
int a, b, f;
vector<int> cut;
};
#define PATH_MARGIN 4
struct PathFinder
{
vector<int> dist, prv, seen;
int stamp = 0;
PathFinder()
{
dist.assign(V, 0);
prv.assign(V, -1);
seen.assign(V, 0);
}
inline int cellCost(int v, int f) const
{
if (!tree0[v])
return 0;
if (dry10[v] <= f)
return 12;
return 100 + (9 - dry10[v]);
}
vector<int> bk[128];
vector<int> done, allowed;
void compute(Wall &w, int margin = PATH_MARGIN)
{
if (done.empty())
done.assign(V, 0);
w.cut.clear();
int ra = w.a / N, ca = w.a % N, rb = w.b / N, cb = w.b % N;
int r0 = max(0, min(ra, rb) - margin), r1 = min(N - 1, max(ra, rb) + margin);
int c0 = max(0, min(ca, cb) - margin), c1 = min(N - 1, max(ca, cb) + margin);
stamp++;
dist[w.a] = cellCost(w.a, w.f);
prv[w.a] = -1;
seen[w.a] = stamp;
int curD = dist[w.a];
bk[curD & 127].push_back(w.a);
int pending = 1;
while (pending > 0)
{
auto &B = bk[curD & 127];
if (B.empty())
{
curD++;
continue;
}
int v = B.back();
B.pop_back();
pending--;
if (dist[v] != curD || done[v] == stamp)
continue;
done[v] = stamp;
if (v == w.b)
break;
int r = v / N, c = v % N;
for (int dr = -1; dr <= 1; dr++)
{
int nr = r + dr;
if (nr < r0 || nr > r1)
continue;
for (int dc = -1; dc <= 1; dc++)
{
if (!dr && !dc)
continue;
int nc = c + dc;
if (nc < c0 || nc > c1)
continue;
int u = nr * N + nc;
int nd = curD + cellCost(u, w.f);
if (seen[u] != stamp || nd < dist[u])
{
seen[u] = stamp;
dist[u] = nd;
prv[u] = v;
bk[nd & 127].push_back(u);
pending++;
}
}
}
}
for (int i = 0; i < 128; i++)
bk[i].clear();
for (int v = w.b; v >= 0; v = prv[v])
{
if (tree0[v] && dry10[v] > w.f)
w.cut.push_back(v);
if (v == w.a)
break;
}
}
};
struct WallSA
{
FastMC mc, valid;
uint64_t Tc[64];
long long curTot = 0, lastTot = 0;
PathFinder pf;
vector<Wall> walls;
vector<int> cnt;
int cost = 0;
uint64_t T[64];
vector<int> dryCells;
#define WSEED 12345
RNG rng{WSEED};
long long iters = 0, accepted = 0;
long long mvTry[12] = {0}, mvAcc[12] = {0};
void addCells(const Wall &w, int sgn)
{
for (int v : w.cut)
{
if (sgn > 0)
{
if (cnt[v]++ == 0)
{
cost++;
T[v / N] &= ~(1ULL << (v % N));
}
}
else
{
if (--cnt[v] == 0)
{
cost--;
T[v / N] |= 1ULL << (v % N);
}
}
}
}
int rayEnd(int X, int dy, int dx, int stopD)
{
int r = X / N, c = X % N;
while (true)
{
int nr = r + dy, nc = c + dx;
if (nr < 0 || nr >= N || nc < 0 || nc >= N)
return r * N + c;
int u = nr * N + nc;
if (!tree0[u] || cnt[u] > 0 || dry10[u] <= stopD)
return u;
r = nr;
c = nc;
}
}
bool randomWall(Wall &w)
{
int X = dryCells[rng.nextInt(dryCells.size())];
if (cnt[X] > 0)
return false;
static const int DIRS[4][2] = {{1, 0}, {0, 1}, {1, 1}, {1, -1}};
int k = rng.nextInt(4);
int stopD = rng.nextInt(2) ? 3 : 4;
w.a = rayEnd(X, DIRS[k][0], DIRS[k][1], stopD);
w.b = rayEnd(X, -DIRS[k][0], -DIRS[k][1], stopD);
w.f = 0;
computePath(w);
return !w.cut.empty();
}
WallGen wg;
vector<int> mlabel, mq;
vector<int> bfsIn(int src, int lab)
{
vector<int> d(V, -1);
mq.clear();
mq.push_back(src);
d[src] = 0;
for (size_t h = 0; h < mq.size(); h++)
{
int v = mq[h];
for (int k = 0; k < 4; k++)
{
int u = NB[v][k];
if (u < 0 || mlabel[u] != lab || d[u] >= 0)
continue;
d[u] = d[v] + 1;
mq.push_back(u);
}
}
return d;
}
bool proposeMinCut(vector<Wall> &out)
{
static const int TH[3] = {5, 6, 7};
int th = TH[rng.nextInt(3)];
mlabel.assign(V, -1);
vector<vector<int>> comps;
for (int v0 = 0; v0 < V; v0++)
{
if (!tree0[v0] || cnt[v0] > 0 || dry10[v0] < th || mlabel[v0] >= 0)
continue;
int lab = comps.size();
comps.push_back({v0});
mlabel[v0] = lab;
auto &cl = comps.back();
for (size_t h = 0; h < cl.size(); h++)
{
int v = cl[h];
for (int k = 0; k < 4; k++)
{
int u = NB[v][k];
if (u < 0 || !tree0[u] || cnt[u] > 0 || dry10[u] < th || mlabel[u] >= 0)
continue;
mlabel[u] = lab;
cl.push_back(u);
}
}
}
double totw = 0;
for (auto &c : comps)
if (c.size() >= 40)
totw += (double)c.size() * c.size();
if (totw <= 0)
return false;
double x = rng.nextDouble() * totw;
int lab = -1;
for (int i = 0; i < (int)comps.size(); i++)
if (comps[i].size() >= 40)
{
x -= (double)comps[i].size() * comps[i].size();
if (x <= 0)
{
lab = i;
break;
}
}
if (lab < 0)
return false;
const auto &cells = comps[lab];
vector<int> d0 = bfsIn(cells[rng.nextInt(cells.size())], lab);
int s = cells[0];
for (int v : cells)
if (d0[v] > d0[s])
s = v;
vector<int> ds = bfsIn(s, lab);
int mx = 0;
for (int v : cells)
mx = max(mx, ds[v]);
vector<int> far;
for (int v : cells)
if (ds[v] >= mx * 6 / 10)
far.push_back(v);
int t = far[rng.nextInt(far.size())];
vector<int> dt = bfsIn(t, 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());
static const double FR[3] = {0.15, 0.25, 0.4};
double fa = FR[rng.nextInt(3)], fb = FR[rng.nextInt(3)];
int ra = sd[min((int)sd.size() - 1, (int)(fa * sd.size()))];
int rb = td[min((int)td.size() - 1, (int)(fb * td.size()))];
Cand c;
if (!wg.tryCut(cells, mlabel, lab, ds, dt, ra, rb, 45, c))
return false;
vector<int> inC(V, 0);
for (int v : c.cells)
inC[v] = 1;
vector<int> seen(V, 0);
for (int v0 : c.cells)
{
if (seen[v0])
continue;
vector<int> grp = {v0};
seen[v0] = 1;
for (size_t h = 0; h < grp.size(); h++)
{
int v = grp[h], r = v / N, cc = v % N;
for (int dr = -1; dr <= 1; dr++)
for (int dc = -1; dc <= 1; dc++)
{
int nr = r + dr, nc = cc + dc;
if (nr < 0 || nr >= N || nc < 0 || nc >= N)
continue;
int u = nr * N + nc;
if (inC[u] && !seen[u])
{
seen[u] = 1;
grp.push_back(u);
}
}
}
auto farthest = [&](int src)
{
map<int, int> dd;
vector<int> q = {src};
dd[src] = 0;
int best = src;
for (size_t h = 0; h < q.size(); h++)
{
int v = q[h], r = v / N, cc = v % N;
if (dd[v] > dd[best])
best = v;
for (int dr = -1; dr <= 1; dr++)
for (int dc = -1; dc <= 1; dc++)
{
int nr = r + dr, nc = cc + dc;
if (nr < 0 || nr >= N || nc < 0 || nc >= N)
continue;
int u = nr * N + nc;
if (inC[u] && !dd.count(u))
{
dd[u] = dd[v] + 1;
q.push_back(u);
}
}
}
return best;
};
int a = farthest(grp[0]), b = farthest(a);
Wall w;
w.a = a;
w.b = b;
w.f = th - 1;
computePath(w);
{
int m = -1, md = 1;
for (int v : grp)
{
int dmin = 99;
for (int u : w.cut)
dmin = min(dmin, max(abs(u / N - v / N), abs(u % N - v % N)));
if (dmin > md)
{
md = dmin;
m = v;
}
}
if (m >= 0)
{
Wall w1 = w, w2 = w;
w1.b = m;
w2.a = m;
computePath(w1);
computePath(w2);
if (!w1.cut.empty())
out.push_back(w1);
if (!w2.cut.empty())
out.push_back(w2);
continue;
}
}
if (!w.cut.empty())
out.push_back(w);
}
return !out.empty();
}
int shiftCell(int v, int r)
{
int i = v / N + rng.nextInt(2 * r + 1) - r, j = v % N + rng.nextInt(2 * r + 1) - r;
i = max(0, min(N - 1, i));
j = max(0, min(N - 1, j));
return i * N + j;
}
unordered_map<uint64_t, vector<int>> memo;
long long memoHit = 0, memoMiss = 0;
void computePath(Wall &w)
{
{
uint64_t key = ((uint64_t)w.a * V + w.b) * 16 + w.f;
auto it = memo.find(key);
if (it != memo.end())
{
w.cut = it->second;
memoHit++;
return;
}
memoMiss++;
pf.compute(w);
sort(w.cut.begin(), w.cut.end());
if (memo.size() > MEMO)
memo.clear();
memo.emplace(key, w.cut);
return;
}
pf.compute(w);
sort(w.cut.begin(), w.cut.end());
}
struct Tree
{
int src = -1, f = -1, r0 = 0, r1 = -1, c0 = 0, c1 = -1;
long long lastUse = 0;
vector<int> prv;
};
Tree trees[PCACHE];
long long useClock = 0, treeBuilds = 0, treeHits = 0;
void buildTree(Tree &tr, int s, int f, int r0, int r1, int c0, int c1)
{
tr.src = s;
tr.f = f;
tr.r0 = r0;
tr.r1 = r1;
tr.c0 = c0;
tr.c1 = c1;
if (tr.prv.empty())
tr.prv.assign(V, -1);
auto &dist = pf.dist;
auto &seen = pf.seen;
auto &done = pf.done;
if (done.empty())
done.assign(V, 0);
int stp = ++pf.stamp;
dist[s] = pf.cellCost(s, f);
tr.prv[s] = -1;
seen[s] = stp;
int curD = dist[s];
pf.bk[curD & 127].push_back(s);
int pending = 1;
while (pending > 0)
{
auto &B = pf.bk[curD & 127];
if (B.empty())
{
curD++;
continue;
}
int v = B.back();
B.pop_back();
pending--;
if (dist[v] != curD || done[v] == stp)
continue;
done[v] = stp;
int r = v / N, c = v % N;
for (int dr = -1; dr <= 1; dr++)
{
int nr = r + dr;
if (nr < r0 || nr > r1)
continue;
for (int dc = -1; dc <= 1; dc++)
{
if (!dr && !dc)
continue;
int nc = c + dc;
if (nc < c0 || nc > c1)
continue;
int u = nr * N + nc;
int nd = curD + pf.cellCost(u, f);
if (seen[u] != stp || nd < dist[u])
{
seen[u] = stp;
dist[u] = nd;
tr.prv[u] = v;
pf.bk[nd & 127].push_back(u);
pending++;
}
}
}
}
}
void cachedPath(Wall &w, int s, int t)
{
uint64_t mkey = ((uint64_t)w.a * V + w.b) * 16 + w.f;
{
auto it = memo.find(mkey);
if (it != memo.end())
{
w.cut = it->second;
memoHit++;
return;
}
memoMiss++;
}
useClock++;
int tr_ = t / N, tc_ = t % N;
Tree *hit = nullptr;
for (auto &tr : trees)
if (tr.src == s && tr.f == w.f && tr_ >= tr.r0 && tr_ <= tr.r1 && tc_ >= tr.c0 && tc_ <= tr.c1)
{
hit = &tr;
break;
}
if (!hit)
{
Tree *lru = &trees[0];
for (auto &tr : trees)
if (tr.lastUse < lru->lastUse)
lru = &tr;
int M = PATH_MARGIN + 3;
int sr = s / N, sc = s % N;
buildTree(*lru, s, w.f, max(0, min(sr, tr_) - M), min(N - 1, max(sr, tr_) + M), max(0, min(sc, tc_) - M), min(N - 1, max(sc, tc_) + M));
hit = lru;
treeBuilds++;
}
else
treeHits++;
hit->lastUse = useClock;
w.cut.clear();
for (int v = t; v >= 0; v = hit->prv[v])
{
if (tree0[v] && dry10[v] > w.f)
w.cut.push_back(v);
if (v == s)
break;
}
sort(w.cut.begin(), w.cut.end());
if (memo.size() > MEMO)
memo.clear();
memo.emplace(mkey, w.cut);
}
vector<long long> qcur, qnew;
long long earlyRej = 0, fullEval = 0;
double energy()
{
bool same = true;
for (int r = 0; r < N; r++)
if (T[r] != Tc[r])
{
same = false;
break;
}
if (same)
{
lastTot = curTot;
qnew = qcur;
return log((double)curTot / mc.S + 1e-9);
}
long long tot = 0;
for (int s = 0; s < mc.S; s++)
tot += (qnew[s] = mc.sampleSq(s, T));
lastTot = tot;
return log((double)tot / mc.S + 1e-9);
}
double energyT(double allowed)
{
if (TWOSTAGE <= 0)
return energy();
bool same = true;
for (int r = 0; r < N; r++)
if (T[r] != Tc[r])
{
same = false;
break;
}
if (same)
{
lastTot = curTot;
qnew = qcur;
return log((double)curTot / mc.S + 1e-9);
}
int S = mc.S, k1 = S / 2;
long long pn = 0, pc = 0;
for (int s = 0; s < k1; s++)
{
pn += (qnew[s] = mc.sampleSq(s, T));
pc += qcur[s];
}
double est = (double)curTot * pn / max(1LL, pc);
if (log(est / S + 1e-9) > allowed + TWOSTAGE)
{
earlyRej++;
return 1e18;
}
long long tot = pn;
for (int s = k1; s < S; s++)
tot += (qnew[s] = mc.sampleSq(s, T));
fullEval++;
lastTot = tot;
return log((double)tot / mc.S + 1e-9);
}
void commit()
{
for (int r = 0; r < N; r++)
Tc[r] = T[r];
curTot = lastTot;
qcur.swap(qnew);
qnew.assign(mc.S, 0);
}
void resample(int S, uint64_t seed)
{
mc.gen(S, seed);
for (int r = 0; r < N; r++)
Tc[r] = T[r];
qcur.assign(S, 0);
qnew.assign(S, 0);
curTot = 0;
for (int s = 0; s < S; s++)
curTot += (qcur[s] = mc.sampleSq(s, T));
}
void greedyInit(int M, double tEnd)
{
double cur = energy();
while (elapsed() < tEnd)
{
double bestScore = 0;
vector<Wall> best;
bool found = false;
for (int m = 0; m < M + GREEDY_MC; m++)
{
vector<Wall> ws;
if (m < M)
{
Wall w;
if (!randomWall(w))
continue;
ws.push_back(w);
}
else if (!proposeMinCut(ws))
continue;
for (auto &w : ws)
addCells(w, +1);
if (cost <= KW)
{
int before = cost;
for (auto &w : ws)
addCells(w, -1);
int added = before - cost;
for (auto &w : ws)
addCells(w, +1);
double e = energy();
double sc = (cur - e) / max(1, added);
if (sc > bestScore)
{
bestScore = sc;
best = ws;
found = true;
}
}
for (auto &w : ws)
addCells(w, -1);
}
if (!found)
break;
for (auto &w : best)
{
addCells(w, +1);
walls.push_back(w);
}
energy();
commit();
cur = log((double)curTot / mc.S + 1e-9);
}
}
int KW = 0;
double lam = 0;
int ccost = 0;
inline double pen(int c) const { return lam * max(0, c - KW); }
vector<pair<double, vector<Wall>>> snaps;
vector<int> run(double tEnd, int S, int segs, double Tm0, double Tm1)
{
KW = K;
cnt.assign(V, 0);
cost = 0;
masksFromTree(tree0, T);
for (int r = 0; r < N; r++)
Tc[r] = T[r];
for (int v = 0; v < V; v++)
if (tree0[v] && dry10[v] >= 6)
dryCells.push_back(v);
if (dryCells.empty())
for (int v = 0; v < V; v++)
if (tree0[v])
dryCells.push_back(v);
valid.gen(160, 31337);
double tStart = elapsed();
resample(S, 5555);
#define NINIT 1
if (GREEDY_M > 0)
{
vector<Wall> bestInit;
double bestInitE = 1e18;
for (int it = 0; it < NINIT; it++)
{
for (auto &w : walls)
addCells(w, -1);
walls.clear();
for (int r = 0; r < N; r++)
Tc[r] = T[r];
curTot = mc.sumSq(T);
greedyInit(GREEDY_M, tStart + (tEnd - tStart) * 0.25 * (it + 1) / NINIT);
double ve = valid.E(T);
if (ve < bestInitE)
{
bestInitE = ve;
bestInit = walls;
}
}
for (auto &w : walls)
addCells(w, -1);
walls = bestInit;
for (auto &w : walls)
addCells(w, +1);
for (int r = 0; r < N; r++)
Tc[r] = T[r];
curTot = mc.sumSq(T);
}
cerr << " greedy walls=" << walls.size() << " cost=" << cost << " validE=" << valid.E(T) << " t=" << elapsed() << "\n";
double tS2 = elapsed();
double bestValidE = valid.E(T);
vector<Wall> bestWalls = walls;
for (int seg = 0; seg < segs; seg++)
{
double segEnd = tS2 + (tEnd - tS2) * (seg + 1) / segs;
int Sseg = (int)lround(S * pow((double)S_END / S, segs > 1 ? (double)seg / (segs - 1) : 1.0));
if (seg > 0 || Sseg != S)
resample(Sseg, 5555 + seg * 7919);
double cur = log((double)curTot / mc.S + 1e-9);
ccost = cost;
double temp = Tm0;
while (true)
{
if ((iters & 15) == 0)
{
double t = elapsed();
if (t > segEnd)
break;
double prog = (t - tS2) / (tEnd - tS2);
temp = Tm0 * pow(Tm1 / Tm0, prog);
double segStart = tS2 + (tEnd - tS2) * seg / segs;
double sp = (t - segStart) / max(1e-9, segEnd - segStart);
lam = LAMBDA * (1.0 + 30.0 * max(0.0, (sp - 0.8) / 0.2));
}
iters++;
int mv = rng.nextInt(100);
int nw = walls.size();
if (nw >= 2 && rng.nextInt(10000) < P_LNS)
{
mvTry[8]++;
int worst = -1;
double we = 1e18;
for (int k = 0; k < nw; k++)
{
addCells(walls[k], -1);
double e = energy();
addCells(walls[k], +1);
if (e < we)
{
we = e;
worst = k;
}
}
Wall old = walls[worst];
addCells(old, -1);
double be = 1e18;
vector<Wall> bws;
for (int m = 0; m < LNS_M; m++)
{
vector<Wall> nws;
if (!proposeMinCut(nws))
continue;
for (auto &w : nws)
addCells(w, +1);
if (cost <= KW + SLACK)
{
double e = energy() + pen(cost);
if (e < be)
{
be = e;
bws = nws;
}
}
for (auto &w : nws)
addCells(w, -1);
}
if (bws.empty())
{
addCells(old, +1);
continue;
}
for (auto &w : bws)
addCells(w, +1);
double e = energy();
if (e + pen(cost) <= cur + pen(ccost) || rng.nextDouble() < exp((cur + pen(ccost) - e - pen(cost)) / temp))
{
walls.erase(walls.begin() + worst);
for (auto &w : bws)
walls.push_back(w);
cur = e;
accepted++;
mvAcc[8]++;
commit();
ccost = cost;
}
else
{
for (auto &w : bws)
addCells(w, -1);
addCells(old, +1);
}
continue;
}
if (rng.nextInt(1000) < P_MINCUT)
{
mvTry[6]++;
vector<Wall> nws;
if (!proposeMinCut(nws))
continue;
for (auto &w : nws)
addCells(w, +1);
int rem = -1;
Wall old;
if (cost > KW + SLACK && nw > 0)
{
rem = rng.nextInt(nw);
old = walls[rem];
addCells(old, -1);
}
if (cost > KW + SLACK)
{
if (rem >= 0)
addCells(old, +1);
for (auto &w : nws)
addCells(w, -1);
continue;
}
double e = energy();
if (e + pen(cost) <= cur + pen(ccost) || rng.nextDouble() < exp((cur + pen(ccost) - e - pen(cost)) / temp))
{
if (rem >= 0)
walls.erase(walls.begin() + rem);
for (auto &w : nws)
walls.push_back(w);
cur = e;
accepted++;
mvAcc[6]++;
commit();
ccost = cost;
}
else
{
if (rem >= 0)
addCells(old, +1);
for (auto &w : nws)
addCells(w, -1);
}
continue;
}
if (nw == 0 || mv < 6)
{
mvTry[0]++;
Wall w;
if (!randomWall(w))
continue;
addCells(w, +1);
if (cost > KW + SLACK)
{
addCells(w, -1);
continue;
}
double e = energy();
if (e + pen(cost) <= cur + pen(ccost) || rng.nextDouble() < exp((cur + pen(ccost) - e - pen(cost)) / temp))
{
walls.push_back(w);
cur = e;
accepted++;
mvAcc[0]++;
commit();
ccost = cost;
}
else
addCells(w, -1);
}
else if (mv < 14)
{
int k = rng.nextInt(nw);
Wall old = walls[k];
addCells(old, -1);
Wall w;
bool hasNew = false;
int type = rng.nextInt(2);
mvTry[1 + type]++;
#define REPL_M 1
double e = 0;
if (type)
{
double be = 1e18;
Wall bw;
for (int m = 0; m < REPL_M; m++)
{
Wall c;
if (!randomWall(c))
continue;
addCells(c, +1);
if (cost <= KW)
{
double ec = energy();
if (ec < be)
{
be = ec;
bw = c;
hasNew = true;
}
}
addCells(c, -1);
}
if (!hasNew)
{
addCells(old, +1);
continue;
}
w = bw;
addCells(w, +1);
e = be;
}
else
e = energy();
if (e + pen(cost) <= cur + pen(ccost) || rng.nextDouble() < exp((cur + pen(ccost) - e - pen(cost)) / temp))
{
if (hasNew)
walls[k] = w;
else
walls.erase(walls.begin() + k);
cur = e;
accepted++;
mvAcc[1 + type]++;
commit();
ccost = cost;
}
else
{
if (hasNew)
addCells(w, -1);
addCells(old, +1);
}
}
else
{
int k = rng.nextInt(nw);
Wall old = walls[k];
Wall w = old;
int t = rng.nextInt(100);
int r = 1 + rng.nextInt(3);
int type;
int movedEnd = -1;
if (t < P_SNAP)
{
type = 7;
static const int D8[8][2] = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}, {1, 1}, {1, -1}, {-1, 1}, {-1, -1}};
int d = rng.nextInt(8);
bool endA = rng.nextInt(2);
int e = endA ? w.a : w.b;
int y = e / N, x = e % N, res = -1;
for (int step = 0; step < 8; step++)
{
int ny = y + D8[d][0], nx = x + D8[d][1];
if (ny < 0 || ny >= N || nx < 0 || nx >= N)
{
res = y * N + x;
break;
}
int u = ny * N + nx;
bool own = false;
for (int c : old.cut)
if (c == u)
{
own = true;
break;
}
if (!tree0[u] || dry10[u] <= 3 || (cnt[u] > 0 && !own))
{
res = u;
break;
}
y = ny;
x = nx;
}
if (res < 0 || res == e)
continue;
if (endA)
w.a = res;
else
w.b = res;
movedEnd = endA ? 0 : 1;
}
else if (t < P_SNAP + (100 - P_SNAP) * 75 / 100)
{
type = 3;
if (rng.nextInt(2))
{
w.a = shiftCell(w.a, r);
movedEnd = 0;
}
else
{
w.b = shiftCell(w.b, r);
movedEnd = 1;
}
}
else if (t < P_SNAP + (100 - P_SNAP) * 92 / 100)
{
type = 4;
int di = rng.nextInt(2 * r + 1) - r, dj = rng.nextInt(2 * r + 1) - r;
int ai = w.a / N + di, aj = w.a % N + dj, bi = w.b / N + di, bj = w.b % N + dj;
if (ai < 0 || ai >= N || aj < 0 || aj >= N || bi < 0 || bi >= N || bj < 0 || bj >= N)
continue;
w.a = ai * N + aj;
w.b = bi * N + bj;
}
else
{
type = 5;
static const int FS[4] = {0, 2, 3, 4};
w.f = FS[rng.nextInt(4)];
if (w.f == old.f)
continue;
}
mvTry[type]++;
if (movedEnd == 0)
cachedPath(w, w.b, w.a);
else if (movedEnd == 1)
cachedPath(w, w.a, w.b);
else
computePath(w);
if (w.cut == old.cut)
{
walls[k] = w;
mvAcc[11]++;
continue;
}
addCells(old, -1);
addCells(w, +1);
if (cost > KW + SLACK)
{
addCells(w, -1);
addCells(old, +1);
mvAcc[10]++;
continue;
}
double allowed = cur + pen(ccost) - pen(cost) - temp * log(max(1e-300, rng.nextDouble()));
double e = energyT(allowed);
if (e <= allowed)
{
walls[k] = w;
cur = e;
accepted++;
mvAcc[type]++;
commit();
ccost = cost;
}
else
{
addCells(w, -1);
addCells(old, +1);
}
}
}
if (cost > KW)
{
cerr << " wseg " << seg << " infeasible cost=" << cost << "\n";
continue;
}
double ve = valid.E(T);
if (ve < bestValidE)
{
bestValidE = ve;
bestWalls = walls;
}
else if (ve > bestValidE * REVERT)
{
for (auto &w : walls)
addCells(w, -1);
walls = bestWalls;
for (auto &w : walls)
addCells(w, +1);
for (int r = 0; r < N; r++)
Tc[r] = T[r];
}
cerr << " wseg " << seg << " walls=" << walls.size() << " cost=" << cost << " trainE=" << exp(cur) / (N * N)
<< " validE=" << ve << " t=" << elapsed() << "\n";
}
cerr << "WallSA iters=" << iters << " acc=" << accepted << " bestValid=" << bestValidE << " moves(try/acc):";
for (int i = 0; i < 10; i++)
cerr << " " << mvTry[i] << "/" << mvAcc[i];
cerr << " same=" << mvAcc[11] << " overK=" << mvAcc[10] << " early=" << earlyRej << "/" << (earlyRej + fullEval);
cerr << " treeBuilds=" << treeBuilds << " hits=" << treeHits;
cerr << " memoHit=" << memoHit << " miss=" << memoMiss;
cerr << "\n";
vector<int> c2(V, 0);
vector<int> cuts;
for (auto &w : bestWalls)
for (int v : w.cut)
if (!c2[v])
{
c2[v] = 1;
cuts.push_back(v);
}
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);
FE::setup();
#define WSA_T0 0.02
#define WSA_T1 0.001
WallSA wsa;
vector<int> cuts = wsa.run(WSA_T, WSA_S, WSA_SEG, WSA_T0, WSA_T1);
cerr << "walls: cuts=" << cuts.size() << " t=" << elapsed() << "\n";
if ((int)cuts.size() < K)
{
Samples train;
train.gen(48, 777);
vector<char> tr = tree0;
for (int v : cuts)
tr[v] = 0;
cuts = greedyFill(train, tr, cuts, elapsed() + 0.1);
}
SAStats st;
vector<int> out = runSA2(cuts, TIME_LIMIT - 0.05, 4242, st);
cerr << "SA iters=" << st.iters << " acc=" << st.acc << " polished=" << st.polished << " validE=" << st.validE << "\n";
cerr << "t=" << elapsed() << "\n";
cout << out.size() << "\n";
for (int v : out)
cout << v / N << " " << v % N << "\n";
}