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

                                                F. Sonya and Informatics

题意：

          输入n(1e2),k(1e9),接下来n个数ai（0或1），而k表示k个操作，每个操作随机取i,j(1<=i<j<=n)并交换，问最后ai是不下降序列的概率为多少（对1e9+7取模）

题解：

          一眼概率dp，接下来找状态，考虑最终的状态，0在前面，1在后面，那么设0的个数cnt0个，1的个数cnt1个，dp[i][j]表示j次操作后前cnt0个位置中有i个0的操作方案，最后所求为dp[cnt0][k]/(总方案=(n*(n-1)/2)^k)

考虑转移：dp[i][j-1]->dp[i][j]    不改变前cnt0个数中0的个数 0和0交换,1和1交换，前cnt0个数中01交换，后n-cnt0个数中01交换。

                 dp[i-1][j-1]->dp[i][j]   前cnt0个数中0的个数加1    前cnt0个数中1和后面的0交换。

                 dp[i+1][j-1]->dp[i][j]  前cnt0个数中0的个数减1    前cnt0个数中0和后面的1交换。

发现转移和j没有关系，可以用矩阵快速幂加速，时间复杂度O(n^3log(k))

代码：

#include<bits/stdc++.h>
using namespace std;
#define ll lonng long
const int N=109;
const int mod=1e9+7;
int n,a[N],k;
int cnt0,cnt1;
int b[N][N];
void jzc(int a[][N],int b[][N],int c){
    int ans[N][N];
    for(int i=0;i<=cnt0;i++)for(int j=0;j<=cnt0;j++)ans[i][j]=0;
    for(int i=0;i<=cnt0;i++)
      for(int j=0;j<=cnt0;j++)
        for(int k=0;k<=cnt0;k++)ans[i][j]=(ans[i][j]+1ll*a[i][k]*b[k][j])%mod;
    for(int i=0;i<=cnt0;i++)for(int j=0;j<=cnt0;j++)a[i][j]=ans[i][j];
}
void jzksm(int a[][N],int b,int c){
    int ans[N][N],base[N][N];
    for(int i=0;i<=cnt0;i++){
        for(int j=0;j<=cnt0;j++){
            if(i==j)ans[i][j]=1;else ans[i][j]=0;
            base[i][j]=a[i][j];
        }
    }
    while(b){
        if(b&1)jzc(ans,base,c);
        jzc(base,base,c);
        b>>=1;
    }
    for(int i=0;i<=cnt0;i++)
      for(int j=0;j<=cnt0;j++)a[i][j]=ans[i][j];
}
void print(){
    for(int i=0;i<=cnt0;i++){
        for(int j=0;j<=cnt0;j++)cout<<b[i][j]<<" ";
        cout<<endl;
    }
    cout<<endl;
}
int poww(int a,int b,int c){
    int ans=1,base=a;
    while(b){
        if(b&1)ans=1ll*ans*base%c;
        base=1ll*base*base%c;
        b>>=1;
    }
    return ans;
}
int main(){
   // freopen("tt.in","r",stdin),freopen("tt.out","w",stdout);
    ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
    cin>>n>>k;
    for(int i=1;i<=n;i++)cin>>a[i],cnt0+=(a[i]==0),cnt1+=(a[i]==1);
    int st=0;for(int i=1;i<=cnt0;i++)st+=(a[i]==0);
    for(int i=0;i<=cnt0;i++){
        if(i)b[i-1][i]=(cnt0-i+1)*(cnt0-i+1);
        b[i][i]=cnt0*(cnt0-1)/2+(cnt1-1)*(cnt1)/2+i*(cnt0-i)+(cnt0-i)*(cnt1-cnt0+i);
        if(i<cnt0)b[i+1][i]=(i+1)*(cnt1-cnt0+i+1);
    }
   // print();
    jzksm(b,k,mod);
   // print();
    int sum=n*(n-1)/2;sum=poww(sum,k,mod);
    cout<<b[st][cnt0]*1ll*poww(sum,mod-2,mod)%mod<<endl;
    return 0;
}