結果
| 問題 | No.5025 Forest Fire |
| コンテスト | |
| ユーザー |
 Shki
|
| 提出日時 | 2026-10-20 18:06:04 |
| 言語 | C++23(gcc16) (gcc 16.2.0 + boost 1.92.0 + ACL) |
| 結果 |
AC
不安定
|
| 実行時間 | 1,579 ms / 2,000 ms |
| + 142µs | |
| コード長 | 38,846 bytes |
| 記録 | |
| コンパイル時間 | 5,867 ms |
| コンパイル使用メモリ | 422,992 KB |
| 実行使用メモリ | 9,932 KB |
| スコア | 157,491,501 |
| 最終ジャッジ日時 | 2026-10-24 14:38:34 |
| 合計ジャッジ時間 | 142,578 ms |
|
ジャッジサーバーID （参考情報） |
judge4_0 / <nil> |
| 純コード判定しない問題か言語 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 50 |
ソースコード
#include <bits/stdc++.h>
using namespace std;
// C++17, no external libraries. N = 50, K = 100.
// This solver minimizes the expected squared connected component sizes.
// Additional search: weighted 8-neighbor minimum-cost separators,
// articulation-aware greedy reconstruction (LNS), exact re-ranked swaps,
// route-preserving shortest-path rerouting, and annealing on common samples.
// Judges should compile with: g++ -std=c++17 -O2
using ll = long long;
using Clock = chrono::steady_clock;
struct Random {
uint64_t state;
explicit Random(uint64_t s) : state(s) {}
uint64_t next() {
uint64_t z = (state += 0x9e3779b97f4a7c15ULL);
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
return z ^ (z >> 31);
}
int integer(int n) { return (int)(next() % (uint64_t)n); }
double real() { return (next() >> 11) * 0x1.0p-53; }
};
class Solver {
static constexpr int NS = 192;
int N, K, V;
bool sparseEnough = false;
vector<int> d, trees;
vector<array<int,4>> adj;
vector<uint8_t> degree;
vector<vector<uint8_t>> dry;
Random rnd;
Clock::time_point start;
const double TIME_LIMIT = 1.62; // Leave margin for slower judges.
static constexpr int COOLING_CYCLES = 3; // Two reheats between 3 SA cooling cycles.
vector<int> tin, low, sub, separateSum, separateSq, dfsOrder;
int timer = 0;
const vector<uint8_t>* currentCut = nullptr;
const vector<uint8_t>* currentDry = nullptr;
vector<int> mark, queueBuffer, comp, sizes;
int stamp = 0;
double elapsed() const {
return chrono::duration<double>(Clock::now() - start).count();
}
bool timeOkay(double deadline) const { return elapsed() < deadline; }
ll evalOne(const vector<uint8_t>& cut, const vector<uint8_t>& sample) {
if (++stamp == INT_MAX) { fill(mark.begin(), mark.end(), 0); stamp = 1; }
ll result = 0;
for (int root = 0; root < V; ++root) {
if (!sample[root] || cut[root] || mark[root] == stamp) continue;
int h = 0, t = 1;
queueBuffer[0] = root;
mark[root] = stamp;
while (h < t) {
int u = queueBuffer[h++];
for (int z = 0; z < degree[u]; ++z) {
int w = adj[u][z];
if (sample[w] && !cut[w] && mark[w] != stamp) {
mark[w] = stamp;
queueBuffer[t++] = w;
}
}
}
result += (ll)t*t;
}
return result;
}
ll evaluate(const vector<uint8_t>& cut, int from, int count) {
ll val = 0;
for (int s = from; s < from + count; ++s) val += evalOne(cut, dry[s]);
return val;
}
// Exact local objective updates for the fixed SA sample set.
// Label the currently retained dry components once after every accepted move.
// For a trial only components touching changed cells can change their sizes.
static constexpr int SA_SAMPLES = 32;
vector<vector<int16_t>> cachedId;
vector<vector<int16_t>> cachedSize;
ll cachedScore[SA_SAMPLES] = {};
vector<int> touchedComp;
int touchedGeneration=0;
bool saCacheReady=false;
void invalidateSACache() { saCacheReady=false; }
void buildSACache(const vector<uint8_t>& cut) {
if((int)cachedId.size()!=SA_SAMPLES) {
cachedId.assign(SA_SAMPLES,vector<int16_t>(V));
cachedSize.resize(SA_SAMPLES);
touchedComp.assign(V,0);
}
for(int s=0;s<SA_SAMPLES;s++) {
auto &ids=cachedId[s];
auto &sizesHere=cachedSize[s];
fill(ids.begin(),ids.end(),int16_t(-1));
sizesHere.clear();
ll total=0;
const auto &sample=dry[s];
for(int r=0;r<V;r++) {
if(!sample[r] || cut[r] || ids[r]>=0)continue;
int id=(int)sizesHere.size();
int h=0,t=1;queueBuffer[0]=r;ids[r]=id;
while(h<t) {
int u=queueBuffer[h++];
for(int j=0;j<degree[u];j++) {
int v=adj[u][j];
if(ids[v]<0 && sample[v] && !cut[v]) {
ids[v]=id;
queueBuffer[t++]=v;
}
}
}
sizesHere.push_back((int16_t)t);
total+=(ll)t*t;
}
cachedScore[s]=total;
}
saCacheReady=true;
}
ll evaluateSATrial(const vector<uint8_t>& current,
const vector<uint8_t>& trial) {
if(!saCacheReady)buildSACache(current);
int changes[2500], nc=0;
for(int v:trees)if(current[v]!=trial[v])changes[nc++]=v;
if(nc==0){ll sum=0;for(ll val:cachedScore)sum+=val;return sum;}
// For very large modifications, a full scan is faster than
// building an affected-component frontier.
if(nc>45)return evaluate(trial,0,SA_SAMPLES);
ll total=0;
for(int s=0;s<SA_SAMPLES;s++) {
const auto &sample=dry[s];
const auto &ids=cachedId[s];
const auto &sizesHere=cachedSize[s];
++touchedGeneration;
ll oldSq=0;
auto include=[&](int u) {
int id=ids[u];
if(id<0 || touchedComp[id]==touchedGeneration)return;
touchedComp[id]=touchedGeneration;
ll sz=sizesHere[id];
oldSq+=sz*sz;
};
for(int i=0;i<nc;i++) {
int u=changes[i];
if(!sample[u])continue;
if(trial[u])include(u); // a previously retained dry vertex is cut
else for(int j=0;j<degree[u];j++)include(adj[u][j]);
}
if(++stamp==INT_MAX){fill(mark.begin(),mark.end(),0);stamp=1;}
ll newSq=0;
auto bfs=[&](int root) {
if(!sample[root] || trial[root] || mark[root]==stamp)return;
int h=0,t=1;
queueBuffer[0]=root;mark[root]=stamp;
while(h<t) {
int u=queueBuffer[h++];
for(int j=0;j<degree[u];j++) {
int v=adj[u][j];
if(sample[v] && !trial[v] && mark[v]!=stamp) {
mark[v]=stamp;queueBuffer[t++]=v;
}
}
}
newSq+=(ll)t*t;
};
for(int i=0;i<nc;i++) {
int u=changes[i];
if(!sample[u])continue;
if(!trial[u])bfs(u);
else for(int j=0;j<degree[u];j++)bfs(adj[u][j]);
}
total+=cachedScore[s]-oldSq+newSq;
}
#ifdef CHECK_INCREMENTAL
ll full=evaluate(trial,0,SA_SAMPLES);
if(full!=total){cerr<<"INVALID incremental="<<total<<" full="<<full<<"\n";abort();}
#endif
return total;
}
void dfs(int u, int parent) {
tin[u] = low[u] = timer++;
sub[u] = 1;
separateSum[u] = separateSq[u] = 0;
dfsOrder.push_back(u);
for (int z = 0; z < degree[u]; ++z) {
int v = adj[u][z];
if ((*currentCut)[v] || !(*currentDry)[v]) continue;
if (tin[v] == -1) {
dfs(v, u);
sub[u] += sub[v];
low[u] = min(low[u], low[v]);
if (low[v] >= tin[u]) {
separateSum[u] += sub[v];
separateSq[u] += sub[v]*sub[v];
}
} else if (v != parent) {
low[u] = min(low[u], tin[v]);
}
}
}
void removalGains(const vector<uint8_t>& cut, int S, vector<ll>& gain, int startSample=0) {
fill(gain.begin(), gain.end(), 0);
currentCut = &cut;
for (int s = 0; s < S; ++s) {
currentDry = &dry[(s+startSample)%NS];
fill(tin.begin(), tin.end(), -1);
dfsOrder.clear();
timer = 0;
for (int r = 0; r < V; ++r) {
if (cut[r] || !(*currentDry)[r] || tin[r] != -1) continue;
int first = (int)dfsOrder.size();
dfs(r, -1);
int componentSize = (int)dfsOrder.size() - first;
int totalSq = componentSize*componentSize;
for (int i = first; i < (int)dfsOrder.size(); ++i) {
int u = dfsOrder[i];
int leftover = componentSize - 1 - separateSum[u];
gain[u] += (ll)totalSq - separateSq[u] - leftover*leftover;
}
}
}
}
void restoreCosts(const vector<uint8_t>& cut, int S, vector<ll>& cost) {
fill(cost.begin(), cost.end(), 0LL);
for (int s = 0; s < S; ++s) {
fill(comp.begin(), comp.end(), -1);
sizes.clear();
for (int r = 0; r < V; ++r) {
if (cut[r] || !dry[s][r] || comp[r] != -1) continue;
int id = (int)sizes.size();
int h = 0, t = 1;
queueBuffer[0] = r;
comp[r] = id;
while (h < t) {
int u = queueBuffer[h++];
for (int z = 0; z < degree[u]; ++z) {
int w = adj[u][z];
if (!cut[w] && dry[s][w] && comp[w] == -1) {
comp[w] = id;
queueBuffer[t++] = w;
}
}
}
sizes.push_back(t);
}
for (int v : trees) {
if (!cut[v] || !dry[s][v]) continue;
int uniqueIds[4], count = 0;
ll joined = 1, previous = 0;
for (int z = 0; z < degree[v]; ++z) {
int id = comp[adj[v][z]];
if (id == -1) continue;
bool duplicate = false;
for (int j = 0; j < count; ++j) duplicate |= (uniqueIds[j] == id);
if (duplicate) continue;
uniqueIds[count++] = id;
ll q = sizes[id];
joined += q;
previous += q*q;
}
cost[v] += joined*joined - previous;
}
}
}
int numberCut(const vector<uint8_t>& cut) const {
int c = 0;
for (int v : trees) c += cut[v];
return c;
}
void greedyComplete(vector<uint8_t>& cut, int samplesCount, double deadline) {
vector<ll> gain(V);
int count = numberCut(cut);
while (count < K && timeOkay(deadline)) {
removalGains(cut, samplesCount, gain);
ll bestVal = -1;
int bestVertex = -1;
for (int v : trees) {
if (!cut[v] && gain[v] > bestVal) {
bestVal = gain[v];
bestVertex = v;
}
}
if (bestVertex < 0) break;
cut[bestVertex] = 1;
++count;
}
}
// Cheap single pass to fill a geometric candidate to exactly K cuts.
void fillOnce(vector<uint8_t>& cut) {
if (numberCut(cut) >= K) return;
vector<ll> gains(V);
removalGains(cut, 12, gains);
vector<int> remaining;
for (int v : trees) if (!cut[v]) remaining.push_back(v);
sort(remaining.begin(), remaining.end(), [&](int a, int b) {
return gains[a] > gains[b];
});
int current = numberCut(cut);
for (int v : remaining) {
if (current >= K) break;
cut[v] = 1;
++current;
}
}
// Water-aware thin blocking barriers. The resulting path can bend.
vector<int> getBarrier(bool horizontal, int target, int seed) {
const int INF = 1000000000;
vector<int> dp(V, INF), parent(V, -1);
auto at = [&](int progress, int pos) {
return horizontal ? pos*N + progress : progress*N + pos;
};
for (int pos = 0; pos < N; ++pos) {
int v = at(0, pos);
dp[pos] = 500*!!d[v] + 6*d[v] + 9*abs(pos-target);
}
for (int x = 1; x < N; ++x) {
for (int p = 0; p < N; ++p) {
int best = INF, idx = -1;
for (int delta = -1; delta <= 1; ++delta) {
int prev = p+delta;
if (prev < 0 || prev >= N) continue;
int trial = dp[(x-1)*N+prev];
if (trial < best) { best = trial; idx = prev; }
}
int v = at(x,p);
int waviness = (seed&1) ? 5 : 12;
dp[x*N+p] = best + 500*!!d[v] + 6*d[v] + waviness*abs(p-target);
parent[x*N+p] = idx;
}
}
int p = 0, best = INF;
for (int j = 0; j < N; ++j) {
int score = dp[(N-1)*N+j] + 9*abs(j-target);
if (score < best) best = score, p = j;
}
vector<int> result;
result.reserve(N);
for (int x = N-1; x >= 0; --x) {
result.push_back(at(x,p));
if (x) p = parent[x*N+p];
}
return result;
}
// Node-weighted shortest 8-connected separator path. Unlike the old
// monotone DP path, this can detour backwards through large lakes.
vector<int> minCostBarrier(bool horizontal, int target, int kind) {
const int INF = 1000000000;
vector<int> dist(V,INF),previous(V,-1);
// We split the grid into sides by removing an 8-connected path.
// On a 4-neighbor primal grid, 8-connectivity is the correct dual.
priority_queue<pair<int,int>, vector<pair<int,int>>, greater<pair<int,int>>> q;
auto cost=[&](int v)->int {
int y=v/N,x=v%N;
int perp=horizontal?y:x;
int off=abs(perp-target);
int width=kind==0?2:(kind==1?4:7);
int bias=off*off/width;
if(d[v]==0)return 2+bias;
int w=110+(kind==2?0:d[v]*2);
return w+bias;
};
for(int v=0;v<V;v++) {
int y=v/N,x=v%N;
if((horizontal?x:y)!=0)continue;
int edgePerp=horizontal?y:x;
int startCost=cost(v)+40*abs(edgePerp-target);
dist[v]=startCost;
q.push({startCost,v});
}
int dest=-1, destCost=INF;
while(!q.empty()){
auto [ds,u]=q.top();q.pop();
if(ds!=dist[u])continue;
int y=u/N,x=u%N;
int progress=horizontal?x:y;
if(progress==N-1){
int finalCost=ds+40*abs((horizontal?y:x)-target);
if(finalCost<destCost) destCost=finalCost,dest=u;
// still explore; water may provide a better target end
if(ds>destCost)break;
}
for(int dy=-1;dy<=1;dy++)for(int dx=-1;dx<=1;dx++){
if(!dy&&!dx)continue;
int ny=y+dy,nx=x+dx;
if(ny<0||ny>=N||nx<0||nx>=N)continue;
int v=ny*N+nx;
int nd=ds+cost(v)+((dy&&dx)?5:0);
if(nd<dist[v]) {dist[v]=nd;previous[v]=u;q.push({nd,v});}
}
}
vector<int> path;
for(int u=dest;u>=0;u=previous[u]) path.push_back(u);
return path;
}
// Ruin-and-recreate large neighborhood search. Restore a whole
// collection of cuts, then re-cut the best vertices using recomputed
// articulation-based marginal benefits after every insertion.
bool destroyRepair(const vector<uint8_t>& initial, vector<uint8_t>& trial,
int amount, int sampleOffset, bool local,
const vector<int>& lowRest) {
trial=initial;
vector<int> cutList;
for(int v:trees)if(trial[v])cutList.push_back(v);
if((int)cutList.size()<amount)return false;
int anchor=cutList[rnd.integer(cutList.size())];
int removed=0;
if(local){
sort(cutList.begin(),cutList.end(),[&](int a,int b){
int da=abs(a/N-anchor/N)+abs(a%N-anchor%N);
int db=abs(b/N-anchor/N)+abs(b%N-anchor%N);
return da!=db?da<db:a<b;
});
} else if(!lowRest.empty() && rnd.integer(100)<65) {
// Recover low-value currently selected cuts first.
cutList=lowRest;
} else {
for(int j=(int)cutList.size()-1;j>0;j--)swap(cutList[j],cutList[rnd.integer(j+1)]);
}
for(int i=0;i<(int)cutList.size()&&removed<amount;i++){
int v=cutList[i];
if(!trial[v])continue;
trial[v]=0;
removed++;
}
if(!removed)return false;
vector<ll> gain(V);
for(int i=0;i<removed;i++) {
removalGains(trial,8,gain,sampleOffset);
int bestV=-1;
// Occasionally sample among the five best candidates, rather
// than deterministically rebuilding the destroyed solution.
pair<ll,int> top[5];
for(auto& t:top)t={LLONG_MIN,-1};
for(int v:trees){
if(trial[v])continue;
ll value=gain[v];
for(int j=0;j<5;j++)if(value>top[j].first){
for(int k=4;k>j;k--)top[k]=top[k-1];
top[j]={value,v};
break;
}
}
int take = (rnd.integer(100)<75)?0:rnd.integer(5);
if(top[take].second<0)take=0;
bestV=top[take].second;
if(bestV<0)return false;
trial[bestV]=1;
}
return true;
}
// Steepest-exchange search: approximate every single swap with exact
// marginal restoration/removal costs, check top-ranked swaps with the
// full stochastic objective. This is not a random SA neighbor.
bool steepestExchange(vector<uint8_t>& current, ll& currentScore,
vector<uint8_t>& best, ll& bestValidation,
double deadline, int sampleOffset) {
vector<ll> gain(V), cost(V);
removalGains(current,12,gain,sampleOffset);
restoreCosts(current,12,cost);
vector<int> a,b;
for(int v:trees){if(current[v])a.push_back(v);else b.push_back(v);}
if(a.empty()||b.empty())return false;
int ac=min(14,(int)a.size()),bc=min(24,(int)b.size());
partial_sort(a.begin(),a.begin()+ac,a.end(),[&](int x,int y){return cost[x]<cost[y];});
partial_sort(b.begin(),b.begin()+bc,b.end(),[&](int x,int y){return gain[x]>gain[y];});
struct Move{ll approx;int a,b;};
vector<Move> cand;
for(int i=0;i<ac;i++)for(int j=0;j<bc;j++)cand.push_back({cost[a[i]]-gain[b[j]],a[i],b[j]});
int sz=min(42,(int)cand.size());
partial_sort(cand.begin(),cand.begin()+sz,cand.end(),[&](auto&x,auto&y){return x.approx<y.approx;});
vector<uint8_t> trial=current,chosen;
ll minScore=currentScore;
for(int i=0;i<sz&&timeOkay(deadline);i++){
int x=cand[i].a,y=cand[i].b;
trial[x]=0;trial[y]=1;
ll score=evaluate(trial,0,32);
if(score<minScore){minScore=score;chosen=trial;}
trial[x]=1;trial[y]=0;
}
if(chosen.empty())return false;
current.swap(chosen);currentScore=minScore;
ll val=evaluate(current,96,96);
if(val<bestValidation){bestValidation=val;best=current;}
return true;
}
vector<uint8_t> makeCut(const vector<int>& a, const vector<int>& b) {
vector<uint8_t> cut(V, 0);
for (int v : a) if (d[v]) cut[v] = 1;
for (int v : b) if (d[v]) cut[v] = 1;
return cut;
}
int swapSize(double progress) {
int maxK = min({7,K,(int)trees.size()-K});
if (maxK == 0) return 0;
double decay = 1.0 + 1.4*progress, sum = 0;
double weight[8];
for (int i = 1; i <= maxK; ++i) {
weight[i] = 0.96*exp(-decay*(i-1)) + 0.04;
sum += weight[i];
}
double r = rnd.real()*sum;
for (int i = 1; i <= maxK; ++i) {
r -= weight[i];
if (r < 0) return i;
}
return maxK;
}
// Path-preserving nonlocal move. Reroute an 8-connected segment of an
// existing firebreak through an alternative narrow winding corridor.
// Restrict to a local box around both endpoints for speed and locality.
bool reroute(const vector<uint8_t>& cut, vector<uint8_t>& result,
const vector<int>& selected, const vector<int>& lowRest) {
if(selected.size()<12)return false;
int src=selected[rnd.integer((int)selected.size())];
// Find nearby 8-connected cut vertices.
vector<int> seen(V,-1),q(V),pre(V,-1);
int h=0,t=1;seen[src]=0;q[0]=src;
while(h<t&&t<200){
int u=q[h++],y=u/N,x=u%N;
for(int dy=-1;dy<=1;dy++)for(int dx=-1;dx<=1;dx++){
if(!dx&&!dy)continue;
int ny=y+dy,nx=x+dx;
if(ny<0||ny>=N||nx<0||nx>=N)continue;
int v=ny*N+nx;
if(cut[v]&&seen[v]<0){seen[v]=seen[u]+1;pre[v]=u;q[t++]=v;}
}
}
vector<int> ends;
for(int i=0;i<t;i++)if(seen[q[i]]>=7&&seen[q[i]]<=20)ends.push_back(q[i]);
if(ends.empty())return false;
int dst=ends[rnd.integer((int)ends.size())];
vector<int> oldPath;
for(int u=dst;u!=src&&u>=0;u=pre[u])oldPath.push_back(u);
oldPath.push_back(src);
if(oldPath.size()<7)return false;
int minY=N,maxY=0,minX=N,maxX=0;
for(int v:oldPath){minY=min(minY,v/N);maxY=max(maxY,v/N);minX=min(minX,v%N);maxX=max(maxX,v%N);}
minY=max(0,minY-3);maxY=min(N-1,maxY+3);
minX=max(0,minX-3);maxX=min(N-1,maxX+3);
// Early anchors are harder to remove; unused trees in the narrow
// corridor have a small penalty to produce a true alternative path.
result=cut;
for(int v:oldPath)if(v!=src&&v!=dst) result[v]=0;
const int INF=1000000000;
vector<int> dist(V,INF),par(V,-1);
priority_queue<pair<int,int>,vector<pair<int,int>>,greater<pair<int,int>>> pq;
dist[src]=0;pq.push({0,src});
int wander=rnd.integer(13);
while(!pq.empty()){
auto [cd,u]=pq.top();pq.pop();
if(cd!=dist[u])continue;
if(u==dst)break;
int y=u/N,x=u%N;
for(int dy=-1;dy<=1;dy++)for(int dx=-1;dx<=1;dx++){
if(!dx&&!dy)continue;
int ny=y+dy,nx=x+dx;
if(ny<minY||ny>maxY||nx<minX||nx>maxX)continue;
int v=ny*N+nx;
int cost=result[v]?0:(d[v]?90+3*d[v]:1);
// perturb path tie breaking and discourage the original path
if(!result[v]&&d[v])cost+=int((uint64_t(v)*1836311903ULL+uint64_t(wander)*5039ULL)%17);
int nd=cd+cost;
if(nd<dist[v]){dist[v]=nd;par[v]=u;pq.push({nd,v});}
}
}
if(dist[dst]>=INF)return false;
for(int u=dst;u>=0;u=par[u])if(d[u])result[u]=1;
int delta=numberCut(result)-K;
if(delta>0){
// Remove noncritical cuts, preferring those with the lowest
// restore marginal gain. Avoid destroying the new path.
vector<int> pathFlag(V,0);
for(int u=dst;u>=0;u=par[u])pathFlag[u]=1;
for(int v:lowRest)if(delta>0&&result[v]&&!pathFlag[v])result[v]=0,--delta;
if(delta>0)return false;
} else if(delta<0){
vector<ll> gain(V);
removalGains(result,8,gain,rnd.integer(NS));
vector<int> remaining;
for(int v:trees)if(!result[v])remaining.push_back(v);
int k=min((int)remaining.size(),-delta);
if(k>0){
partial_sort(remaining.begin(),remaining.begin()+k,remaining.end(),
[&](int a,int b){return gain[a]>gain[b];});
for(int i=0;i<k;i++)result[remaining[i]]=1;
}
}
return numberCut(result)==K && result!=cut;
}
// Every successful neighbor retains exactly the same number of cuts.
bool propose(vector<uint8_t>& out, const vector<uint8_t>& cut,
const vector<int>& cutList, const vector<int>& lowRest,
const vector<int>& highGain, double progress) {
if (cutList.empty()) return false;
out = cut;
int roll = rnd.integer(100);
if(roll<8 && reroute(cut,out,cutList,lowRest))return true;
if (roll>=8 && roll<14) {
// Move all cut sites on one (popular) row or column by one cell.
bool row = (rnd.integer(2) == 0);
int counts[50] = {};
for (int u : cutList) ++counts[row ? u/N : u%N];
int maxCount = *max_element(counts, counts+N);
if (maxCount >= 5) {
vector<int> lines;
for (int i=0;i<N;i++) if (counts[i] >= max(5,maxCount-3)) lines.push_back(i);
int line = lines[rnd.integer((int)lines.size())];
int off = rnd.integer(2) ? 1 : -1;
if (line+off >= 0 && line+off < N) {
vector<pair<int,int>> changes;
for (int j = 0; j < N; ++j) {
int old = row ? line*N+j : j*N+line;
int nxt = row ? (line+off)*N+j : j*N+line+off;
if (cut[old] && d[nxt] && !cut[nxt]) changes.push_back({old,nxt});
}
if (changes.size() >= 3) {
for (auto [a,b]:changes) { out[a]=0; out[b]=1; }
return true;
}
}
}
}
int k = swapSize(progress);
if (k == 0) return false;
vector<int> removed;
removed.reserve(k);
bool cluster = k>=2 && rnd.integer(100)<45;
int anchor = cutList[rnd.integer((int)cutList.size())];
for (int j=0;j<k;j++) {
int a = -1;
if (cluster) {
for (int trial=0;trial<30;trial++) {
int q=cutList[rnd.integer((int)cutList.size())];
int distance=abs(q/N-anchor/N)+abs(q%N-anchor%N);
if (out[q] && distance <= k+2) {a=q;break;}
}
}
if (a < 0 && !lowRest.empty() && rnd.integer(100)<60) {
int top = min((int)lowRest.size(), 30+7*k);
for (int trial=0;trial<20;trial++) {
int q=lowRest[rnd.integer(top)];
if (out[q]) { a=q;break; }
}
}
if (a < 0) {
for (int trial=0;trial<20;trial++) {
int q=cutList[rnd.integer((int)cutList.size())];
if (out[q]) {a=q;break;}
}
}
if (a < 0) for (int q:cutList) if (out[q]) {a=q;break;}
if (a < 0) return false;
removed.push_back(a);
out[a]=0;
}
for (int j=0;j<k;j++) {
int a=removed[j], b=-1;
if (rnd.integer(100)<75) {
int radius=rnd.integer(100)<80?2:5;
for (int tries=0;tries<35;tries++) {
int y=a/N + rnd.integer(2*radius+1)-radius;
int x=a%N + rnd.integer(2*radius+1)-radius;
if (x<0||x>=N||y<0||y>=N) continue;
int q=y*N+x;
if (d[q] && !cut[q] && !out[q]) {b=q;break;}
}
}
if (b < 0 && !highGain.empty()) {
int top=min((int)highGain.size(), 100+15*k);
for (int tries=0;tries<35;tries++) {
int q=highGain[rnd.integer(top)];
if (!cut[q] && !out[q]) {b=q;break;}
}
}
if (b<0) {
for (int tries=0;tries<35;tries++) {
int q=trees[rnd.integer((int)trees.size())];
if (!cut[q] && !out[q]) {b=q;break;}
}
}
if (b<0) for (int q:trees) if (!cut[q] && !out[q]) {b=q;break;}
if (b<0) return false;
out[b]=1;
}
return true;
}
public:
Solver(int n, int k, vector<int> dd)
: N(n), K(k), V(n*n), d(move(dd)), adj(V), degree(V,0),
rnd(0x31aa707f682e1cb9ULL), start(Clock::now()),
tin(V), low(V), sub(V), separateSum(V), separateSq(V),
mark(V), queueBuffer(V), comp(V) {
for (int v=0;v<V;v++) {
if (d[v]>0) trees.push_back(v);
int y=v/N,x=v%N;
if (y) adj[v][degree[v]++]=v-N;
if (y+1<N) adj[v][degree[v]++]=v+N;
if (x) adj[v][degree[v]++]=v-1;
if (x+1<N) adj[v][degree[v]++]=v+1;
}
// Guard against log(0) in scorers that don't handle E==0.
// If <=K trees exist, leaving one is always legal and ensures E>0.
if (!trees.empty() && (int)trees.size()<=K) {
sparseEnough = true;
K=(int)trees.size()-1;
}
else K=min(K,(int)trees.size());
dry.assign(NS,vector<uint8_t>(V));
for (int s=0;s<NS;s++) {
for (int v=0;v<V;v++) dry[s][v]=(d[v] && rnd.integer(10)<d[v]);
}
}
vector<uint8_t> solve() {
vector<uint8_t> empty(V,0);
if (sparseEnough) {
// Keep exactly one lowest-d tree: E remains positive, so log(E)
// is defined even if the checker has no special E==0 handling.
int survivor=trees[0];
for (int v:trees) if (d[v]<d[survivor]) survivor=v;
for (int v:trees) if (v!=survivor) empty[v]=1;
return empty;
}
if (K==0) return empty;
if (K==(int)trees.size()) {
for(int v:trees) empty[v]=1;
return empty;
}
vector<uint8_t> best=empty;
ll bestValidation = LLONG_MAX;
auto registerCandidate = [&](const vector<uint8_t>& cand) {
ll value = evaluate(cand,96,96);
if (value < bestValidation) {bestValidation=value;best=cand;}
};
// Initial greedy: always valid even if interrupted.
vector<uint8_t> greed(V,0);
greedyComplete(greed,32,0.57);
fillOnce(greed);
registerCandidate(greed);
#ifdef PROFILE
cerr<<"greed "<<elapsed()<<" score="<<bestValidation<<"\n";
#endif
// Only invoke expensive geometric search when the greedy solution
// leaves a substantial fraction of the baseline burn risk.
double fraction = double(evaluate(greed,96,64)) / max(1LL,evaluate(empty,96,64));
if (timeOkay(0.65) && K>=50 && fraction>0.40) {
vector<vector<int>> horizontal, vertical;
for (int t : {10,15,20,25,30,35,40}) {
horizontal.push_back(getBarrier(true,t,0));
vertical.push_back(getBarrier(false,t,0));
}
// Bidirectional shortest paths can exploit winding waterways.
for(int t:{12,20,25,30,38})for(int kind=0;kind<2;kind++){
horizontal.push_back(minCostBarrier(true,t,kind));
vertical.push_back(minCostBarrier(false,t,kind));
}
vector<pair<ll,vector<uint8_t>>> candidates;
auto addCandidate=[&](const vector<int>&a,const vector<int>&b) {
vector<uint8_t> c=makeCut(a,b);
if (numberCut(c)>K) return;
candidates.emplace_back(evaluate(c,0,6),move(c));
};
for (const auto &h:horizontal) for (const auto &v:vertical) {
if (!timeOkay(0.73)) break;
addCandidate(h,v);
}
for(int i=0;i<(int)horizontal.size();i++) for(int j=i+1;j<(int)horizontal.size();j++) {
if (!timeOkay(0.73)) break;
addCandidate(horizontal[i],horizontal[j]);
addCandidate(vertical[i],vertical[j]);
}
sort(candidates.begin(),candidates.end(),
[](const auto&a,const auto&b){return a.first<b.first;});
if(candidates.size()>5) candidates.resize(5);
for(auto& item:candidates) {
if(!timeOkay(0.85)) break;
fillOnce(item.second);
registerCandidate(item.second);
}
}
#ifdef PROFILE
cerr<<"init "<<elapsed()<<" score="<<bestValidation<<"\n";
#endif
// Simulated annealing on common random numbers, with independent validation.
if(!timeOkay(1.02)) return best;
vector<uint8_t> current=best, trial;
const int SA_S=SA_SAMPLES;
ll currentScore=evaluate(current,0,SA_S);
ll bestSA=currentScore;
double saStart=elapsed(), saFinish=TIME_LIMIT-0.04;
double T0=max(10.0,(double)currentScore/SA_S*0.010);
double T1=max(0.3,(double)currentScore/SA_S*0.000010);
int accepted=0, lastRankUpdate=-10000, iter=0;
int coolingCycle=0, reheats=0;
int cycleIterations[COOLING_CYCLES] = {};
vector<int> cutList, lowRest, highGain;
vector<ll> gains(V), costs(V);
auto refreshCutList=[&]() {
cutList.clear();
for(int v:trees) if(current[v]) cutList.push_back(v);
};
auto refreshRanks=[&]() {
removalGains(current,12,gains);
restoreCosts(current,12,costs);
lowRest=cutList;
sort(lowRest.begin(),lowRest.end(),[&](int a,int b){return costs[a]<costs[b];});
highGain.clear();
for(int v:trees) if(!current[v]) highGain.push_back(v);
sort(highGain.begin(),highGain.end(),[&](int a,int b){return gains[a]>gains[b];});
if(highGain.size()>160)highGain.resize(160);
lastRankUpdate=accepted;
};
refreshCutList();
while(timeOkay(saFinish)) {
++iter;
// Divide the *available* SA time into three complete cooling cycles.
// At t=1/3 and t=2/3, temperature immediately jumps back to T0.
// Retain the current solution and the best independently validated one.
double globalProgress = (elapsed()-saStart)/max(0.001,saFinish-saStart);
globalProgress = max(0.0,min(1.0,globalProgress));
const int targetCycle = min(COOLING_CYCLES-1, int(globalProgress*COOLING_CYCLES));
if (targetCycle != coolingCycle) {
reheats += targetCycle-coolingCycle;
coolingCycle = targetCycle;
// Restart the guided-neighbor ranking at each new heating cycle.
lastRankUpdate = -10000;
}
++cycleIterations[coolingCycle];
double localProgress = globalProgress*COOLING_CYCLES-coolingCycle;
localProgress = max(0.0,min(1.0,localProgress));
double temp = max(1e-9, T0 * pow(T1/T0,localProgress) * SA_S);
// Periodically visit an entirely different neighborhood:
// exact-guided one-exchange and 8-24-point greedy destroy/repair.
if(iter%97==0 && timeOkay(saFinish-0.12)){
if(steepestExchange(current,currentScore,best,bestValidation,
min(saFinish-0.10,elapsed()+0.045),rnd.integer(NS))){
accepted++;
refreshCutList();
invalidateSACache();
}
continue;
}
if(iter%27==0 && timeOkay(saFinish-0.11)) {
int amount=7+rnd.integer(20);
bool local=(rnd.integer(100)<65);
if(destroyRepair(current,trial,amount,rnd.integer(NS),local,lowRest)){
ll score=evaluateSATrial(current,trial);
if(score<=currentScore || rnd.real()<exp(-double(score-currentScore)/temp)){
current.swap(trial);
currentScore=score;
invalidateSACache();
accepted++;
refreshCutList();
}
if(currentScore<bestSA){bestSA=currentScore;registerCandidate(current);}
}
continue;
}
if (lastRankUpdate < 0 || (accepted-lastRankUpdate>=24 && timeOkay(saFinish-0.08))) {
refreshRanks();
}
if(!propose(trial,current,cutList,lowRest,highGain,localProgress))continue;
ll newVal=evaluateSATrial(current,trial);
ll delta=newVal-currentScore;
if(delta<=0 || (rnd.real()<exp(-double(delta)/max(1e-9,temp)))) {
current.swap(trial);
currentScore=newVal;
invalidateSACache();
++accepted;
refreshCutList();
}
if(currentScore<bestSA) {
bestSA=currentScore;
if(timeOkay(saFinish-0.005)) registerCandidate(current);
}
}
// Validate the final state if possible.
if(timeOkay(TIME_LIMIT-0.006)) registerCandidate(current);
#ifdef PROFILE
cerr<<"done "<<elapsed()<<" iter="<<iter<<" accepted="<<accepted
<<" reheats="<<reheats<<" cycles="<<cycleIterations[0]<<","<<cycleIterations[1]
<<","<<cycleIterations[2]<<" score="<<bestValidation<<"\n";
#endif
return best;
}
};
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
int N, K;
if(!(cin >> N >> K)) return 0;
if(N <= 0 || N > 50 || K < 0) { cout << "0\n"; return 0; }
// Each dry level is a SINGLE DIGIT (0..9).
// Accept both contiguous digit rows and whitespace-separated cells.
// Read characters rather than extracting multi-digit ints.
vector<int> d(N*N);
int filled = 0;
char ch;
while (filled < N*N && cin.get(ch)) {
if (ch >= '0' && ch <= '9') d[filled++] = ch - '0';
}
if (filled != N*N) {
cout << "0\n";
return 0;
}
Solver solver(N,K,d);
auto cut=solver.solve();
// Output sanitizer: never emit an invalid cell, repeated coordinate or >K cuts.
vector<int> output;
for (int v=0;v<N*N;v++) if (cut[v] && d[v]>0 && (int)output.size()<K)
output.push_back(v);
cout << output.size() << '\n';
for (int v : output) cout << v/N << ' ' << v%N << '\n';
return 0;
}