CF1151F Sonya and Informatics（概率期望，DP，矩阵快速幂）

明明是水题结果没切掉……降智了……

首先令 $c$ 为序列中 $0$ 的个数，那么排序后序列肯定是前面 $c$ 个 $0$，后面 $n-c$ 个 $1$。

那么就能上 DP 了。（居然卡在这里……）

$f[i][j]$ 表示经过 $i$ 次操作后，前 $c$ 个数中有 $j$ 个 $0$ 的方案数。答案就是 $dfrac{f[k][c]}{sum f[k][i]}$。

这个状态的好处就是可以直接求出以下这些值：

• 前 $c$ 个数中 $1$ 的个数为 $c-j$
• 后 $c$ 个数中 $0$ 的个数为 $c-j$
• 后 $c$ 个数中 $1$ 的个数为 $n-2c+j$（所以 $jge 2c-n$）

初始状态：令初始序列前 $c$ 个数中 $0$ 的个数为 $cnt$，那么 $f[0][cnt]=1$，其它的 $f[0][i]=0$。

转移：

$$f[i][j]+=f[i-1][j](dfrac{c(c-1)}{2}+dfrac{(n-c)(n-c-1)}{2}+j(c-j)+(c-j)(n-2c+j))$$

括号中第一个是前 $c$ 个中交换，第二个是后 $n-c$ 个中交换，第三个是前面的 $0$ 和后面的 $0$ 交换，第四个是前面的 $1$ 和后面的 $1$ 交换。

$$f[i][j]+=f[i-1][j-1](c-j+1)^2$$

前 $c$ 个中的 $0$ 和后 $n-c$ 个中的 $1$ 交换。

$$f[i][j]+=f[i-1][j+1](j+1)(n-2c+j+1)$$

前 $c$ 个中的 $1$ 和后 $n-c$ 个中的 $0$ 交换。

然后发现第 $i$ 层之和第 $i-1$ 层有关，那么可以矩阵快速幂。

时间复杂度 $O(n^3log k)$。

#include<bits/stdc++.h>
using namespace std;
const int maxn=101,mod=1000000007;
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
#define MEM(x,v) memset(x,v,sizeof(x))
inline int read(){
    char ch=getchar();int x=0,f=0;
    while(ch<'0' || ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar();
    return f?-x:x;
}
int n,k,a[maxn],c,cnt;
inline int add(int a,int b){return a+b<mod?a+b:a+b-mod;}
inline int sub(int a,int b){return a<b?a-b+mod:a-b;}
inline int mul(int a,int b){return 1ll*a*b%mod;}
inline int qpow(int a,int b){
    int ans=1;
    for(;b;b>>=1,a=mul(a,a)) if(b&1) ans=mul(ans,a);
    return ans;
}
struct matrix{
    int a[maxn][maxn];
    matrix(){MEM(a,0);}
    matrix operator*(const matrix &t)const{
        matrix ans;
        FOR(i,0,c) FOR(k,0,c) FOR(j,0,c) ans.a[i][j]=add(ans.a[i][j],mul(a[i][k],t.a[k][j]));
        return ans;
    }
}beg,fac,ans;
matrix qpow(matrix a,int b){
    matrix ans;
    FOR(i,0,c) ans.a[i][i]=1;
    for(;b;b>>=1,a=a*a) if(b&1) ans=ans*a;
    return ans;
}
int main(){
    n=read();k=read();
    FOR(i,1,n) a[i]=read(),c+=!a[i];
    FOR(i,1,c) cnt+=!a[i];
    beg.a[cnt][0]=1;
    FOR(i,0,c){
        fac.a[i][i]=(1ll*c*(c-1)/2+1ll*(n-c)*(n-c-1)/2+1ll*i*(c-i))%mod;
        if(i>=2*c-n) fac.a[i][i]=add(fac.a[i][i],mul(c-i,n-2*c+i));
        if(i) fac.a[i][i-1]=mul(c-i+1,c-i+1);
        if(i!=c && i+1>=2*c-n) fac.a[i][i+1]=mul(i+1,n-2*c+i+1);
    }
    ans=qpow(fac,k)*beg;
    int s=0;
    FOR(i,0,c) s=add(s,ans.a[i][0]);
    s=mul(qpow(s,mod-2),ans.a[c][0]);
    printf("%d
",s);
}