CF1151F Sonya and Informatics

矩阵乘法，dp

题意

给定一个长度为 \(n\) 的 01 数组 \(a\)，进行 \(k\) 次操作，每次操作等概率随机选择两个不同位置并交换它们的值。求 \(k\) 次操作后数组变为非递减（即所有 0 在前，所有 1 在后）的概率，结果对 \(10^9+7\) 取模。

• 数组长度 \(n\)：\(2 \leq n \leq 100\)
• 操作次数 \(k\)：\(1 \leq k \leq 10^9\)
• 数组元素 \(a_i \in \{0, 1\}\)

思路

设 \(0\) 的数量为 \(m\)，设计 dp 状态为操作了 \(k\) 次，前 \(m\) 个数中有 \(i\) 个 \(0\) 的方案数。

有转移：\(dp_{k+1}[i]=dp_{k}[i]\times (\dfrac{n(n-1)}{2}-i\times (n-2m+i)-(m-i)^2)+dp_k[i-1]\times (m-i)^2+dp_k[i+1]\times i\times (n-2m+i)\).

具体来说，前 \(m\) 个数有 \(i\) 个零，如果这次交换在内部进行，那么不会改变零的数量，对应 \(\dfrac{n(n-1)}{2}-i\times (n-2m+i)-(m-i)^2\);

内部的 \(0\) 与外部的 \(1\) 交换，会让零的数量减一，对应 \(i\times (n-2m+i)\);

另外一种就是内部的 \(1\) 换成了外部的 \(0\) 的零，会让零的数量加一，对应 \((m-i)^2\)。

最后用矩阵加速一下就做完了。

代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=1e9+7;
int n,m,k,cnt,sum;
int a[110];
struct MAT
{
    int n,m;
    int a[110][110];
    MAT(int n=0,int m=0):n(n),m(m) {memset(a,0,sizeof(a));}
};
void output(MAT x)
{
    for(int i=0;i<=x.n;i++)
    {
        for(int j=0;j<=x.m;j++) cout<<x.a[i][j]<<" ";
        cout<<endl;
    }
    cout<<endl;
}
MAT operator*(MAT x,MAT y)
{
    MAT re=MAT(x.n,y.m);
    for(int i=0;i<=x.n;i++)
        for(int j=0;j<=y.m;j++)
            for(int k=0;k<=x.m;k++)
                (re.a[i][j]+=x.a[i][k]*y.a[k][j]%mod)%=mod;
    return re;
}
MAT mat_ksm(MAT a,int b)
{
    if(b==1) return a;
    MAT re=mat_ksm(a,b>>1);
    re=re*re;
    if(b&1) re=re*a;
    return re;
}
int ksm(int a,int b)
{
    if(b==1) return a;
    int re=ksm(a,b>>1);
    re=re*re%mod;
    if(b&1) re=re*a%mod;
    return re;
}
int add(int i) {return max(0ll,(m-i)*(m-i));}
int del(int i) {return max(0ll,i*(n-2*m+i));}
int keep(int i) {return max(0ll,n*(n-1)/2-add(i)-del(i));}
signed main()
{
    cin>>n>>k;
    for(int i=1;i<=n;i++)
    {
        cin>>a[i];
        m+=!a[i];
    }
    MAT m1=MAT(m,m),m2=MAT(m,1);
    for(int i=0;i<=m;i++)
    {
        m1.a[i][i]=keep(i);
        if(i-1>=0) m1.a[i][i-1]=add(i-1);
        if(i+1<=m) m1.a[i][i+1]=del(i+1);
    }
    m1=mat_ksm(m1,k);
    for(int i=1;i<=m;i++) cnt+=(a[i]==0);
    m2.a[cnt][1]=1;
    m2=m1*m2;
    for(int i=0;i<=m;i++) (sum+=m2.a[i][1])%=mod;
    cout<<m2.a[m][1]*ksm(sum,mod-2)%mod;
    return 0;
}