[BZOJ4770]图样(概率期望、二进制数位dp)
题面
https://darkbzoj.tk/problem/4770
题解
设\(f(n,m)\)表示n个点,每一个点的权值从\([0,2^{m+1})\)随机生成,这n个点的最小生成树期望权值和。
\(E(i,j,m)\)表示我们有两堆数,分别有i,j个,都从\([0,2^{m+1})\)随机生成,从两堆数中分别取一个,异或和的最小值的期望。
\(S(i,j,m,x)\)表示我们有两堆数,分别有i,j个,都从\([0,2^{m+1})\)随机生成,从两堆数中分别取一个,异或和的最小值\(\geq x\)的概率。
-
转移\(f(n,m)\):
如果n个点的权值的最高位(即\(2^m\)位)上的值都是0或都是1,即为\(f(n,m-1)\);否则,一定是最高位0、最高位1的分别先形成两个连通块,然后有且仅有一条边连接这两个连通块。这条边的大小可以用E表示。
\[f(n,m)={\frac{f(n,m-1)}{2^{n-1}}}+\sum\limits_{i=1}^{n}[E(i,n-i,m-1)+f(i,m-1)+f(n-i,m-1)+2^m] \] -
转移\(E(i,j,m)\):
由期望的定义,显然有
\[E(i,j,m)=\sum\limits_{x=0}^{2^{m+1}-1}x(S(i,j,m,x)-S(i,j,m,x+1)) \] -
转移\(S(i,j,m,x)\):
这是最重点的部分,要分两种情况考虑。下设\(lg[x]\)表示x的最高位,如\(lg[255]=7\)。显然当\(lg[x]>m\)时S为0。
-
\(lg[x]=m\):首先两堆数的\(2^m\)位上的值一定是一堆全是0,另一堆全是1。随后便可转化为$${\frac{1}{2{i+j-1}}}S(i,j,m-1,x-2)$$。
-
\(lg[x]<m\):枚举第一堆中\(2^m\)位上为0的数的个数\(i_0\),第二堆的个数\(j_0\)。如果两堆中的两个数的\(2^m\)位上的数字不同,那么它们的异或和一定大于m。所以只用考虑相同的情况。这部分的总贡献是
\[\sum\limits_{i_0=0}^{i} \sum\limits_{j_0=0}^{j} \frac{\tbinom{i}{i_0}\tbinom{j}{j_0}}{2^{i+j}}S(i_0,j_0,m-1,x)S(i-i_0,j-j_0,m-1,x) \]
-
最后求的就是\(f(n,m-1)\)。
现在考虑时间复杂度。\(f\)是\(O(nm)\)状态,\(O(n)\)转移。总计\(O(n^2m)\)。
\(E\)是\(O(n^2m)\)状态,\(O(2^m)\)转移。总计\(O(n^2m2^m)\)。
\(S\)则比较大,\(O(n^2m2^m)\)状态,\(O(n^2)\)转移,看上去要T飞了~
其实不然。这些状态根本跑不满;\(i_0\)与\(j_0\)都只需要枚举到i,j,立即砍掉一半常数;再加上\(S(i,j,m,x)\)限制\(i<j\)(否则交换),以及各种边界的特判,我的程序在最大数据上,开/不开O2是691/216ms。
况且还有一法,n,m乘起来不过400,打表不香吗www
代码
#include<bits/stdc++.h>
using namespace std;
#define rg register
#define In inline
#define ll long long
const int N = 50;
const int M = 8;
const int W = 256;
const ll mod = 258280327;
const ll iv2 = 129140164;
namespace ModCalc{
In void Inc(ll &x,ll y){
x += y;if(x >= mod)x -= mod;
}
In void Dec(ll &x,ll y){
x -= y;if(x < 0)x += mod;
}
In ll Add(ll x,ll y){
Inc(x,y);return x;
}
In ll Sub(ll x,ll y){
Dec(x,y);return x;
}
}
using namespace ModCalc;
ll power(ll a,ll n){
ll s = 1,x = a;
while(n){
if(n & 1)s = s * x % mod;
x = x * x % mod;
n >>= 1;
}
return s;
}
ll jc[N+5],iv[N+5],C[N+5][N+5];
ll pv[800+5],lg[W+5];
void prepro(){
jc[0] = 1;
for(rg int i = 1;i <= N;i++)jc[i] = jc[i-1] * i % mod;
iv[N] = power(jc[N],mod - 2);
for(rg int i = N - 1;i >= 0;i--)iv[i] = iv[i+1] * (i + 1) % mod;
for(rg int i = 0;i <= N;i++)
for(rg int j = 0;j <= i;j++)C[i][j] = jc[i] * iv[j] % mod * iv[i-j] % mod;
pv[0] = 1;
for(rg int i = 1;i <= 800;i++)pv[i] = pv[i-1] * iv2 % mod;
lg[0] = -1;
for(rg int i = 1;i <= W;i++)lg[i] = lg[i>>1] + 1;
}
ll S[N+5][N+5][M+5][W+5];
ll E[N+5][N+5][M+5];
ll f[N+5][M+5];
ll calcS(int i,int j,int m,int x){
if(!x || !i || !j)return 1;
if(lg[x] > m)return 0;
if(i > j)return calcS(j,i,m,x);
if(S[i][j][m][x] != -1)return S[i][j][m][x];
ll ans = 0;
if(lg[x] == m){
ll Pr = pv[i+j-1];
ans = Pr * calcS(i,j,m-1,x^(1<<lg[x])) % mod;
}
else{
for(rg int i0 = 0;i0 <= i;i0++)
for(rg int j0 = 0;j0 <= j;j0++){
ll Pr = pv[i+j] * C[i][i0] % mod * C[j][j0] % mod;
Inc(ans,Pr * calcS(i0,j0,m-1,x) % mod * calcS(i-i0,j-j0,m-1,x) % mod);
}
}
return S[i][j][m][x] = ans;
}
ll calcE(int i,int j,int m){
if(m == -1)return 0;
if(i > j)return calcE(j,i,m);
if(E[i][j][m] != -1)return E[i][j][m];
ll ans = 0;
for(rg int x = 0;x < (1<<(m+1));x++){
Inc(ans,1ll * x * Sub(calcS(i,j,m,x),calcS(i,j,m,x+1)) % mod);
}
return E[i][j][m] = ans;
}
ll calcf(int n,int m){
if(f[n][m] != -1)return f[n][m];
if(n == 1 || m == -1)return 0;
ll ans = calcf(n,m - 1) * pv[n-1] % mod;
for(rg int i = 1;i < n;i++){
ll Pr = C[n][i] * pv[n] % mod;
ll curf = Add(Add(1ll<<m,calcE(i,n-i,m-1)),Add(calcf(i,m-1),calcf(n-i,m-1)));
Inc(ans,Pr * curf % mod);
}
return f[n][m] = ans;
}
ll n,m;
int main(){
// freopen("B4770.in","r",stdin);
// freopen("B4770.out","w",stdout);
cin >> n >> m;
prepro();
memset(S,-1,sizeof(S));
memset(E,-1,sizeof(E));
memset(f,-1,sizeof(f));
cout << calcf(n,m - 1) << endl;
return 0;
}