CF1151F Sonya and Informatics （计数dp+矩阵优化）

题目地址

Solution

（duyi是我们的红太阳）

（这里说一句：这题看上去是一个概率dp，鉴于这题的概率dp写法看上去不好写，我们其实可以写一个计数dp）

首先拿到这个题目我们要能设出一个普通dp。难点在于状态如何设计。(n<=100)状态压缩不可行。

这里有一个设计状态的套路：因为这是一个01序列，最终不降序的状态不就是0全部在前面，1全部在后面吗？设一共有c个0，把这个序列分成[1,c],[c+1,n]左右两个部分，我们假设当前序列左边有t个0，这样序列左边就有c-t个1，序列右边就有c-t个0，序列右边就有n-c-c+t个1，这样我们只要知道序列左边有多少个0，整个序列的情况就被我们刻画出来了（妙啊！）。

我们设 f[i,j] 表示 i 次变化后左边[1,c]区间有j个零的方案总数。假设开始有左边有t个0

初始化 ： (f[0,t]=1)

最终答案：

[Ans=frac{f[K,0]}{sum_{i=0}^c f[K,i]} ]

f[i,j] 转移可以从 f[i-1,j-1] , f[i-1,j] , f[i-1,j+1] 三个地方转移过来。

• f[i-1,j-1] => f[i,j] 左边的0变多了，说明是左边的1和右边的0交换了，有 （左1 * 右0）种情况。

• f[i-1,j] => f[i,j] 左边的0没有变，可能是左边的数内部交换，右边的数内部交换，左边的0和右边的0交换，左边的1和右边的1交换，情况自己算。

• f[i-1,j+1] => f[i,j] 左边的0变少了，说明是左边的0和右边的1换了，情况自己算。

于是我们写出了转移方程，发现每次转移都是由 f[i-1, ] 转移而来，所以我们可以考虑矩阵优化这个dp

开一个 (1*c+1) 的矩阵 (F) 下标从 ([0,c]) ，(F[1,j]) 表示 (f[i,j])。

写一个状态转移矩阵 (A)，使得 (F * A = F`) ，(F`[i,j]) 表示 (f[i+1,j])

(A)怎么构建 这部分就很简单了，留给读者自主思考。

Code

Talk is cheap.Show me the code.

#include<bits/stdc++.h>
#define int long long
using namespace std;
inline int read() {
	int x=0,f=1; char ch=getchar();
	while(ch<'0' || ch>'9') { if(ch=='-') f=-1; ch=getchar(); }
	while(ch>='0'&&ch<='9') { x=(x<<3)+(x<<1)+(ch^48); ch=getchar(); }
	return x * f;
}
const int N = 107,MOD = 1e9+7;
int n,K,c;
bool s[N];
struct Mat {
	int val[N][N]; int Max_n, Max_m;
	Mat() {
		memset(val, 0, sizeof(val)); Max_n = Max_m = 0;
	}
	Mat operator * (const Mat el) const {
		Mat c;
		for(int i=0;i<=Max_n;++i)	//此题比较特殊 从零开始计数 
			for(int j=0;j<=el.Max_m;++j)
				for(int k=0;k<=Max_m;++k)
					c.val[i][j] = (c.val[i][j]+((val[i][k]*el.val[k][j])%MOD))%MOD;
		c.Max_n = Max_n, c.Max_m = el.Max_m;
		return c;
	}
	inline void I(int n) {
		Max_n = Max_m = n;
		for(int i=0;i<=n;++i) val[i][i] = 1;
	}
};
inline int power(int x,int y) {
	int res = 1, base = x;
	while(y) {
		if(y&1) res = (res*base)%MOD; base = (base*base)%MOD; y >>= 1;
	}
	return res;
}
signed main()
{
	n = read(), K = read();
	for(int i=1;i<=n;++i) s[i] = read(), c += (s[i]==0);
	Mat A,F; int t = 0;	//记录c区间以左有多少个0 
	for(int i=1;i<=c;++i) if(s[i]==0) ++t;
	F.val[0][t] = 1; F.Max_n = 1, F.Max_m = c;	// f[j] 表示 [0,j] 的总次数 
	for(int i=0;i<=c;++i) {
		if(i != 0) A.val[i-1][i] = (c-i+1)*(c-i+1);
		A.val[i][i] = (c*(c-1)/2) + ((n-c)*(n-c-1)/2) + (i*(c-i)) + (c-i)*(n-c-c+i);
		if(i != c) A.val[i+1][i] = (i+1)*(n-c-c+i+1);
	}
	A.Max_n = A.Max_m = c;
	int y = K; Mat res; res.I(c);
	while(y) {
		if(y&1) res = res*A; A = A*A; y >>= 1;
	}
	F = F * res;
	int ans = F.val[0][c], sum = 0;
	for(int i=0;i<=c;++i) sum = (sum+F.val[0][i])%MOD;
	ans = (ans * power(sum,MOD-2))%MOD;
	printf("%lld
",ans);
	return 0;
}

Summary

这个题目精华就在于dp状态设计这里，没有使用状压而是用一个数刻画了整个局面。

这种做法可以推广到01序列的其他操作，只要保证最终状态和目前状态的区别很小，就可以考虑记录0,1的差异，由此可以设计出dp状态。