Codeforces Round #553 (Div. 2) F. Sonya and Informatics （dp 矩阵快速幂）

题目链接：

https://codeforces.com/contest/1151/problem/F

题目大意：

给定一个数组包含$01$两种字符，每次操作可以选其中两个调换位置，问经过$k$次操作之后变为00000011111这样非严格升序的概率，用逆元的形式输出。

题解思路：

分别记$0$的数量为$zero$，$1$的数量为$one$在$[1-zero]$中的$0$的数量为$now$

先不考虑时空复杂度的问题：
设$dp[i][j]$表示经过i轮(次调换)后，前$zero$个数中有$j$个数为$0$
$res=dp[k][zero]\ast inv((\frac{n\ast (n-1)}{2}{)}^k)$

具体的转移过程见下图：

转移过程中 $dp[i][x]$只和$dp[i-1][y]$有关，也就是说转移方程中不会出现$i$作为系数转移，就可以构建转移矩阵加速
构建过程如下图：

由于$dp$数组记录的是对应的情况数量，最后的结果就是$dp[k][zero]$*$inv({\sum }_{i=1}^{zero}dp[k][i])$ ，$inv$里面也可以直接表示为$(\frac{n\ast (n-1)}{2}{)}^k$

#include <bits/stdc++.h>
using namespace std;
#define IOS ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define ll long long
#define int ll
#define debug cout<<"fuck"<<endl;
#define pb  push_back
#define endl '\n'

const int maxn=(int)100+5;
const int mod=(int)1e9+7;
struct mat
{
    int a[maxn][maxn];
    mat()
    {
        memset(a,0,sizeof(a));
    }
};
int n,k;
int zero=0,one=0;
int now=0;
ll quick_pow_mod(ll a,ll b,ll c)
{
    ll res=1;
    while(b)
    {
        if(b & 1) {res=(res*a)%c;}
        a=(a*a)%c;b=b>>1;
    }
    return res;
}

mat mul(mat a,mat b)
{
    mat res;
    for (int (i) = 0; (i) <=zero; ++(i)) //枚举res的行
    {
        for (int j = 0; j <=zero; ++j)   //枚举res的列
        {
            for (int k = 0; k <=zero; ++k)
            {
                res.a[i][j]=(res.a[i][j]+(a.a[i][k]*b.a[k][j])%mod)%mod;
            }
        }
    }
    return res;
}
mat qpow(mat a,int k)
{
    mat res;
    for(int i=0;i<=zero;i++)
    {
        res.a[i][i]=1;
    }
    while(k)
    {
        if(k&1)
        {
            res=mul(res,a);
        }
        a=mul(a,a);
        k>>=1;
    }
    return res;
}

int a[maxn];
mat dp;//dp[i][j]表示共有i个0，有j个0在应该在的位置上的状态概率
signed  main()
{
    IOS
    cin>>n>>k;
    for(int i=1;i<=n;i++)
    {
        cin>>a[i];
        if(a[i]==0)zero++;
        else one++;
    }
    for(int i=1;i<=zero;i++)
    {
        if(a[i]==0)now++;
    }

    //构造转移矩阵
    for(int i=0;i<=zero;i++)//i对应now
    {
        dp.a[i][i]=(zero*(zero-1)/2+one*(one-1)/2+ (zero-i)*(one-zero+i) + i*(zero-i))%mod;
        if(i-1>=0)
        {
            dp.a[i][i-1]=((zero-i+1)*(zero-i+1))%mod;
        }
        if(i+1<=zero)
        {
            dp.a[i][i+1]=((i+1)*(one-zero+i+1))%mod;
        }
    }
    mat res=qpow(dp,k);
    ll ans=0;
    ans=(res.a[zero][now]*quick_pow_mod(quick_pow_mod(n*(n-1)/2,k,mod),mod-2,mod))%mod;
    cout<<ans<<endl;

    return 0;
}