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$（所以 $j\ge 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^3\log k)$。

#include<bits/stdc++.h>

using namespace std;

const int maxn=,mod=;

#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=,f=;

    while(ch<'' || ch>'') f|=ch=='-',ch=getchar();

    while(ch>='' && ch<='') x=x*+ch-'',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=;

    for(;b;b>>=,a=mul(a,a)) if(b&) ans=mul(ans,a);

    return ans;

}

struct matrix{

    int a[maxn][maxn];

    matrix(){MEM(a,);}

    matrix operator*(const matrix &t)const{

        matrix ans;

        FOR(i,,c) FOR(k,,c) FOR(j,,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,,c) ans.a[i][i]=;

    for(;b;b>>=,a=a*a) if(b&) ans=ans*a;

    return ans;

}

int main(){

    n=read();k=read();

    FOR(i,,n) a[i]=read(),c+=!a[i];

    FOR(i,,c) cnt+=!a[i];

    beg.a[cnt][]=;

    FOR(i,,c){

        fac.a[i][i]=(1ll*c*(c-)/+1ll*(n-c)*(n-c-)/+1ll*i*(c-i))%mod;

        if(i>=*c-n) fac.a[i][i]=add(fac.a[i][i],mul(c-i,n-*c+i));

        if(i) fac.a[i][i-]=mul(c-i+,c-i+);

        if(i!=c && i+>=*c-n) fac.a[i][i+]=mul(i+,n-*c+i+);

    }

    ans=qpow(fac,k)*beg;

    int s=;

    FOR(i,,c) s=add(s,ans.a[i][]);

    s=mul(qpow(s,mod-),ans.a[c][]);

    printf("%d\n",s);

}