P7450 [THUSCH2017]巧克力

· · 题解

主要讲讲这类算法的正确率怎么分析，这样以后就不用瞎猜了。

顺便提一下，这题用的随机化 trick 在学术界叫作 color-coding，是很正经的算法，大家可以去搜一下论文。

color-coding 的流程是这样的：将每种图案随机映射到 [0,k-1] 中的一个整数，从而将问题归约到一个新问题：找出一个联通块，使得 0\sim k-1 在其中至少出现一次，最小化联通块的大小，如果联通块大小相同，则最小化美味度的中位数大小。

在介绍新问题的解法之前，我们先简单分析一下这个算法的正确率。

假设存在一个包含恰好 k 种图案的联通块，且它是最优的。考虑这个联通块中的图案，所有可能的染色方案总共 k^k 种。其中，恰好是一个 0\sim k-1 排列的染色方案总共 k! 种。

所以，运行一次上述算法，找出这个联通块的概率就是 P=\frac{k!}{k^k}。k=5 时，P≈\frac{1}{26}。

btw，如果 k 任意，根据斯特林公式，有如下估计：P≈ \frac{\sqrt{2\pi k}}{e^k}，虽然并不 practical。

如果最优联通块包含的图案数大于 k 种，找出它的概率一定大于 \frac{k!}{k^k}，所以我们只需考虑 P=\frac{k!}{k^k} 的情况。

假设我们运行了 T 次，那么错误率 \epsilon=(1-P)^T。

这个东西不容易分析，我们先对它进行一些处理。

令 \lambda=\frac{1}{P}，则

\epsilon=& \left(1-\frac{1}{\lambda}\right)^T \\=&\left(\frac{\lambda-1}{\lambda}\right)^T \\=&\left(\frac{1}{1+\frac{1}{\lambda-1}}\right)^T \\=&\left(\left(\frac{1}{1+\frac{1}{\lambda-1}}\right)^{\lambda-1}\right)^{\frac{T}{\lambda-1}} \\=&\left(\frac{1}{\left(1+\frac{1}{\lambda-1}\right)^{\lambda-1}}\right)^{\frac{T}{\lambda-1}} \end{aligned}

众所周知 \lim_{n\to \infty}{\left(1+\frac{1}{n}\right)}^n=e，所以 \epsilon ≈\left(\frac{1}{e}\right)^{\frac{T}{\lambda-1}}。

也就是说，我们每运行 \lambda-1 次算法，错误率就会变为原来的 \frac{1}{e}，这个结论可以记一下。

代入数据算一下，取 T=150\sim 200 ，得到的答案就非常准确了。

接下来考虑新问题怎么做。

第一个限制是好处理的，仿照 Steiner Tree 的思路，设 f_{S,i,j} 代表包含了颜色集合 S，当前树根位于 (i,j)，最小的联通块大小。

唯一不同的是贡献在点上，所以每次枚举子集的时候要减去根节点的贡献。

对于第二个限制，我们考虑常见的套路：二分 mid，将美味度 \le mid 的巧克力的权值设为 -1，否则设为 1。

如果存在权值和小于 0 的最小联通块，说明最小中位数 \le mid，二分左边界；否则二分右边界。

如何实现 check 函数呢？

我们将 f_{S,i,j} 设为一个二元组 (cnt,sum)，分别表示联通块大小和权值和。当 cnt 相等时，最小化 sum 即可。

时间复杂度 O(T \log nm (3^k nm+2^knm\log nm))。（写完才发现这里的运行次数 T 好像和数据组数重复了，反正你们理解就行）

这里还有一个剪枝：因为每次二分求出的联通块大小都是相等的，如果这个联通块大小大于当前最优解 ans，可以不用二分。

亲测非常有效，对于我这种 STL 选手，运行时间直接从 5s 降到了 1.3s。

代码很丑。

#include<bits/stdc++.h>
#define vi vector<int>
#define vvi vector<vi>
#define pII pair<int,int>
#define fi first
#define se second
#define rep(i,x,y) for(int i=(x);i<=(y);++i)
using namespace std;
const int inf=0x3f3f3f3f;
mt19937 gen(time(0));
int dx[4]={0,1,0,-1},dy[4]={1,0,-1,0};
int T,n,m,k,V,C,clr[250],num[250];
pII ans=make_pair(inf,inf);
vvi c,a;
vector<vector<pII>>f[1<<5];
vector<vector<bool>>vs;
pII operator+(const pII&a,const pII&b){return make_pair(a.fi+b.fi,a.se+b.se);}
pII operator-(const pII&a,const pII&b){return make_pair(a.fi-b.fi,a.se-b.se);}
pII operator-(const pII&b){return make_pair(-b.fi,-b.se);}
int id(int x){return lower_bound(num+1,num+V+1,x)-num;}
int sgn(int x){return x&1?-1:1;}
bool valid(int x,int y){return x>=1&&x<=n&&y>=1&&y<=m&&~c[x][y];}
int now=inf,mid;
pII add(int x,int y){return make_pair(1,sgn(a[x][y]<=mid));}
void dij(int p){
    priority_queue<pair<pII,pII>>q;
    vs=vector<vector<bool>>(n+1,vector<bool>(m+1));
    rep(i,1,n)rep(j,1,m)if(valid(i,j))q.emplace(-f[p][i][j],make_pair(i,j));
    while(q.size()){
        pII now=q.top().se;q.pop();
        int x=now.fi,y=now.se;
        if(vs[x][y])continue;
        vs[x][y]=1;
        rep(k,0,3){
            int nx=x+dx[k],ny=y+dy[k];
            if(!valid(nx,ny))continue;
            if(f[p][nx][ny]>f[p][x][y]+add(nx,ny)){
                f[p][nx][ny]=f[p][x][y]+add(nx,ny);
                q.emplace(-f[p][nx][ny],make_pair(nx,ny));
            }
        }
    }
}
bool check(){
    rep(p,0,(1<<k)-1)f[p]=vector<vector<pII>>(n+1,vector<pII>(m+1,make_pair(inf,inf)));
    rep(i,1,n)rep(j,1,m)if(valid(i,j))
        f[1<<clr[c[i][j]]][i][j]=add(i,j);
    rep(p,1,(1<<k)-1){
        for(int q=p&p-1;q;q=q-1&p)rep(i,1,n)rep(j,1,m)if(valid(i,j)){ 
            f[p][i][j]=min(f[p][i][j],f[q][i][j]+f[p-q][i][j]-add(i,j));
        }
        dij(p);
    }
    pII ans=make_pair(inf,inf);
    rep(i,1,n)rep(j,1,m)ans=min(ans,f[(1<<k)-1][i][j]);
    now=ans.fi;
    return ans.se<=0;
}
int main(){
    cin>>T;
    while(T--){
        cin>>n>>m>>k;
        ans=make_pair(inf,inf);V=C=0;
        c=a=vvi(n+1,vi(m+1));
        rep(i,1,n)rep(j,1,m)scanf("%d",&c[i][j]),C=max(C,c[i][j]);
        rep(i,1,n)rep(j,1,m)scanf("%d",&a[i][j]),num[++V]=a[i][j];
        sort(num+1,num+V+1);V=unique(num+1,num+V+1)-(num+1);
        rep(i,1,n)rep(j,1,m)a[i][j]=id(a[i][j]);
        int tt=150;
        while(tt--){
            rep(i,1,C)clr[i]=gen()%k;
            now=inf;
            int l=1,r=V;
            while(l<r){
                mid=l+r>>1;
                if(check())r=mid;
                else l=mid+1;
                if(now>ans.fi)break;
            }
            ans=min(ans,make_pair(now,num[l]));
        }
        if(ans.fi==inf)puts("-1 -1");
        else cout<<ans.fi<<" "<<ans.se<<endl;
    }
    return 0;
}