結果
| 問題 | No.5025 Forest Fire |
| コンテスト | |
| ユーザー |
 Shki
|
| 提出日時 | 2026-10-20 19:13:57 |
| 言語 | C++23(gcc16) (gcc 16.2.0 + boost 1.92.0 + ACL) |
| 結果 |
AC
不安定
|
| 実行時間 | 1,579 ms / 2,000 ms |
| + 55µs | |
| コード長 | 46,041 bytes |
| 記録 | |
| コンパイル時間 | 7,702 ms |
| コンパイル使用メモリ | 423,836 KB |
| 実行使用メモリ | 9,932 KB |
| スコア | 168,955,440 |
| 最終ジャッジ日時 | 2026-10-24 14:40:04 |
| 合計ジャッジ時間 | 144,441 ms |
|
ジャッジサーバーID （参考情報） |
judge5_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,
// exact re-ranked swaps, annealed 1-15-point exchange,
// route-preserving shortest-path rerouting, reverse DSU restoration greedy,
// and three cooling cycles with two reheats on common samples.
// GUIDED_MODE=0 original; 1 guided 1..7 only; 2 60/35/5 mix; 3 adaptive mix.
// Judges should compile with: g++ -std=c++17 -O2
#ifndef GUIDED_MODE
#define GUIDED_MODE 2
#endif
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;
}
}
}
// Reverse construction: start with every tree removed, then restore
// V-K vertices with minimum marginal increase in sampled fire damage.
// A lazy heap avoids rescanning every candidate after each insertion.
// In cyclic graphs marginal costs are NOT necessarily monotone (an
// alternative bridge can lower a vertex's cost). Periodic full heap
// refreshes are therefore essential; lazy steps are approximations.
bool reverseRestoration(vector<uint8_t>& answer,
int samplesCount, double deadline) {
const int nt=(int)trees.size();
const int need=nt-K;
if(need<=0) return false;
const int nslots=samplesCount*V;
vector<int> parent(nslots,-1);
vector<int16_t> componentSize(nslots,0);
vector<uint8_t> retained(V,0);
auto findRoot=[&](int p)->int {
while(parent[p]!=p){
int q=parent[p];
parent[p]=parent[q];
p=parent[p];
}
return p;
};
auto roots=[&](int sample,int vertex,int *out)->int {
int n=0,base=sample*V;
for(int j=0;j<degree[vertex];j++){
int v=adj[vertex][j];
if(parent[base+v]<0)continue;
int root=findRoot(base+v);
bool dup=false;
for(int k=0;k<n;k++)if(out[k]==root){dup=true;break;}
if(!dup)out[n++]=root;
}
return n;
};
auto calculate=[&](int vertex)->ll {
ll total=0;
int nearby[4];
for(int s=0;s<samplesCount;s++){
if(!dry[s][vertex])continue;
int count=roots(s,vertex,nearby);
ll sum=1,prior=0;
for(int j=0;j<count;j++){
ll z=componentSize[nearby[j]];
sum+=z;prior+=z*z;
}
total+=sum*sum-prior;
}
return total;
};
auto add=[&](int vertex) {
retained[vertex]=1;
int nearby[4];
for(int s=0;s<samplesCount;s++){
if(!dry[s][vertex])continue;
int base=s*V;
int count=roots(s,vertex,nearby);
int r=base+vertex;
parent[r]=r;
int total=1;
for(int j=0;j<count;j++){
int q=nearby[j];
total+=componentSize[q];
parent[q]=r;
}
componentSize[r]=total;
}
};
struct Entry { ll cost; uint64_t tie; int vertex; };
struct Greater {bool operator()(const Entry&a,const Entry&b)const {
if(a.cost!=b.cost)return a.cost>b.cost;
return a.tie>b.tie;
}};
priority_queue<Entry,vector<Entry>,Greater> pq;
// Keep the existing SA's RNG trajectory unchanged whenever this
// additional initial candidate is rejected.
Random reverseRnd(0x9f4a7c159e3779b9ULL);
// Precompute singleton cost exactly, avoiding 2,500 DSU scans.
for(int v:trees){
int singleton=0;
for(int s=0;s<samplesCount;s++)singleton+=dry[s][v];
pq.push({singleton,reverseRnd.next(),v});
}
int restored=0,checks=0,lastRefresh=-1;
while(restored<need && !pq.empty()){
if((++checks&511)==0 && !timeOkay(deadline))return false;
if(restored>0 && restored%96==0 && restored!=lastRefresh) {
// Fresh priorities also account for marginal costs that
// DECREASE when other bridges connect adjacent components.
priority_queue<Entry,vector<Entry>,Greater> refreshed;
for(int v:trees)if(!retained[v])
refreshed.push({calculate(v),reverseRnd.next(),v});
pq.swap(refreshed);
lastRefresh=restored;
}
Entry e=pq.top();pq.pop();
if(retained[e.vertex])continue;
ll actual=calculate(e.vertex);
// Compare the refreshed cost to the best remaining *cached*
// cost. Periodic global refresh keeps stale priorities in check.
if(!pq.empty() &&
(actual>pq.top().cost ||
(actual==pq.top().cost && e.tie>pq.top().tie))){
pq.push({actual,e.tie,e.vertex});
continue;
}
add(e.vertex);
++restored;
}
answer.assign(V,0);
for(int v:trees)answer[v]=!retained[v];
return restored==need;
}
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;
}
// 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;
}
// At the beginning of EACH cooling cycle, sample all swap sizes 1..15
// uniformly. Gradually concentrate the distribution towards 1-point swaps.
// P(k|p) = exp(-2*p*(k-1)) / sum_j exp(-2*p*(j-1)), where p in [0,1].
// Reheating resets p to 0, restoring the uniform distribution.
int swapSize(double progress) {
const int maxK = min({15, K, (int)trees.size() - K});
if (maxK == 0) return 0;
const double lambda = 2.0 * max(0.0, min(1.0, progress));
double weightSum = 0.0;
double weights[16];
for (int k=1; k<=maxK; ++k) {
weights[k] = exp(-lambda * (k - 1));
weightSum += weights[k];
}
double choice = rnd.real() * weightSum;
for (int k=1; k<=maxK; ++k) {
choice -= weights[k];
if (choice < 0) return k;
}
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;
}
// Guided 1..7-point swaps. At the start of every SA cooling cycle,
// each k is equiprobable; close to the end, prefer k=1.
// Guided destinations are all sampled WITHOUT replacement from the
// top 20% of the latest Lowlink estimated per-tree removal gains.
// Original cut locations are uniform as requested.
int guidedSwapSize(double progress) {
const int m = min({7, K, (int)trees.size()-K});
if (m <= 0) return 0;
const double lambda = 2.0 * clamp(progress, 0.0, 1.0);
// One exp() instead of m calls, using a geometric series.
double decay = exp(-lambda);
double sum = 0.0, w = 1.0;
for (int j=0; j<m; ++j) {sum+=w; w*=decay;}
double r = rnd.real()*sum;
w=1.0;
for (int k=1; k<=m; ++k) {
r-=w;
if (r < 0.0) return k;
w*=decay;
}
return m;
}
bool guidedMultiSwap(vector<uint8_t>& out,
const vector<uint8_t>& cut,
const vector<int>& cutList,
const vector<int>& topGain,
double progress) {
if (cutList.empty() || topGain.empty()) return false;
int k=guidedSwapSize(progress);
if (k<=0)return false;
// The ranking can become stale after accepted moves. Avoid selecting
// an already cut vertex, even if it is still in the old top-20%.
int live=0;
for(int v:topGain) live+=(!cut[v] && d[v]>0);
k=min(k,live);
if (k<=0)return false;
out=cut;
// First remove k distinct old cuts. out[a] makes duplicates impossible.
for(int j=0;j<k;j++) {
int a=-1;
for(int trial=0;trial<24;trial++) {
int v=cutList[rnd.integer((int)cutList.size())];
if(out[v]) {a=v;break;}
}
if(a<0) for(int v:cutList) if(out[v]) {a=v;break;}
if(a<0) return false;
out[a]=0;
}
// Add k distinct new cuts. `cut[v]` (not out[v]) ensures a restored
// old cut is never accidentally selected again as a destination.
for(int j=0;j<k;j++) {
int b=-1;
for(int trial=0;trial<24;trial++) {
int v=topGain[rnd.integer((int)topGain.size())];
if(!cut[v] && !out[v] && d[v]>0) {b=v;break;}
}
if(b<0) for(int v:topGain)
if(!cut[v] && !out[v] && d[v]>0) {b=v;break;}
if(b<0)return false;
out[b]=1;
}
return true;
}
bool unbiasedSwap(vector<uint8_t>& out,
const vector<uint8_t>& cut,
const vector<int>& cutList) {
if(cutList.empty() || (int)trees.size() == (int)cutList.size()) return false;
out=cut;
int a=cutList[rnd.integer((int)cutList.size())];
int b=-1;
for(int k=0;k<35;k++){
int v=trees[rnd.integer((int)trees.size())];
if(!cut[v]) { b=v; break; }
}
if(b<0)for(int v:trees)if(!cut[v]){b=v;break;}
if(b<0)return false;
out[a]=0;
out[b]=1;
return true;
}
// 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;
#if GUIDED_MODE == 1
// Pure guided 1..7-point neighbor, no LNS.
return guidedMultiSwap(out,cut,cutList,highGain,progress);
#elif GUIDED_MODE == 2
// 60% top-20% guided 1..7 swap; 35% original 1..15 swap,
// 5% completely random single swap.
int modeRoll=rnd.integer(100);
if(modeRoll<60) return guidedMultiSwap(out,cut,cutList,highGain,progress);
if(modeRoll<65) return unbiasedSwap(out,cut,cutList);
#elif GUIDED_MODE == 3
// Temperature-dependent mixture: broad early, guided late.
// p=0: 20% guided, 75% broad, 5% uniform.
// p=1: 85% guided, 10% broad, 5% uniform.
int modeRoll=rnd.integer(10000);
int guidedThreshold=(int)(2000+6500*progress);
if(modeRoll<guidedThreshold) return guidedMultiSwap(out,cut,cutList,highGain,progress);
if(modeRoll<guidedThreshold+500) return unbiasedSwap(out,cut,cutList);
#endif
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.55);
fillOnce(greed);
registerCandidate(greed);
#ifdef PROFILE
cerr<<"greed "<<elapsed()<<" score="<<bestValidation<<"\n";
#endif
// Independent reverse-DSU initial solution: leave the 100 most
// expensive-to-restore vertices as an adaptive firebreak.
if(timeOkay(0.61) && (int)trees.size()>K+20) {
vector<uint8_t> reverse;
if(reverseRestoration(reverse,24,0.75)) {
ll reverseValidation = evaluate(reverse,96,96);
if(reverseValidation < bestValidation) {
bestValidation = reverseValidation;
best = reverse;
}
#ifdef PROFILE
cerr<<"reverse "<<elapsed()<<" score="<<reverseValidation
<<" best="<<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.70) && 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.81)) 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.81)) 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.93)) 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 GUIDED_MODE == 0
if(highGain.size()>160)highGain.resize(160);
#else
// All uncut trees have been ranked. Retain the top 20%,
// with a minimum of 1 tree for sparse instances.
size_t shortlist=max<size_t>(1,(highGain.size()+4)/5);
if(highGain.size()>shortlist)highGain.resize(shortlist);
#endif
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 run exact-guided one-exchange alongside guided 1..7 and 1..15-point swaps.
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 (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;
}