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\n",ans);
    return 0;
}