[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}}}+sumlimits_{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)=sumlimits_{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]})$$。

  • (lg[x]<m)：枚举第一堆中(2^m)位上为0的数的个数(i_0)，第二堆的个数(j_0)。如果两堆中的两个数的(2^m)位上的数字不同，那么它们的异或和一定大于m。所以只用考虑相同的情况。这部分的总贡献是

    [sumlimits_{i_0=0}^{i} sumlimits_{j_0=0}^{j} frac{ binom{i}{i_0} binom{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;
}