CF1151F Sonya and Informatics

矩阵优化dp超好题，dp状态十分巧妙

题解：

DP部分：

题目要求排成一个非降序列，由于只有0和1，这意味着要将所有0放在1的前面，又因为0和1的数量是确定的，所以哪些位置给0，哪些位置给1也是确定的

记$m$表示序列中0的个数，我们知道1~m中的0和m+1~n中的1是合法的，反之不合法

状态很巧妙，设$f[i][j]$表示交换$i$次后，前$m$个位置有$j$个0

光凭状态中的一个$j$，我们就可以知道1~m中0和1的个数，m+1~n中0和1的个数，换句话说，知道了合法的数的个数

记$cnt$表示数列中前$m$个数中0的个数，初始化$f[0][cnt]=1$，表示现有状况

考虑转移，容易发现每次交换两个数，$f[i][j]$只会由$f[i-1][j-1]$，$f[i-1][j+1]$，$f[i-1][j]$得到，可以进行如下分类讨论：

  • $f[i-1][j-1](m-(j-1))(m-(j-1))$

表示将前$m$个数中的0与后$n-m$个数中的1做交换，其中第一个$(m-(j-1))$的意思是前m个数中1的个数，第二个的意思是后$n-m$个数中0的个数

  • $f[i-1][j+1](j+1)(n-m-(m-(j+1)))$

与第一中转移相反，即将前$m$个数中的1与后$n-m$个数中的0做交换

  • $f[i-1][j]j(m-j)$

在前$m$个数中选两个数交换
（注意是一个1一个0！）

  • $f[i-1][j](m-j)(n-m-(m-j))$

在后$n-m$个数中选两个数交换
（注意是一个1一个0！）

  • $f[i-1][j]m(m-1)/2$

表示将两个0做交换

  • $f[i-1][j](n-m)(n-m-1)/2$

表示将两个1做交换

最后$f[i][j]$的值就是以上6部分转移之和

而根据题目的要求，$ans=\frac{f[k][m]}{\sum\limits_{i=0}^{m}f[k][i]}$，表示$k$次交换后的前$m$个数都是0的方案数（合法方案）除以总的方案数

这样你就得到了一个$O(nk)$的dp，下面按照分类讨论的顺序dp代码：

int cnt=0;
for(int i=1;i<=m;i++)
    cnt+=(a[i]==0);
f[0][cnt]=1;
for(int i=1;i<=k;i++)
    for(int j=0;j<=m;j++){
        if(j!=0) f[i][j]=(f[i][j]+f[i-1][j-1]*(m-(j-1))*(m-(j-1))%mod)%mod;
        if(j!=m) f[i][j]=(f[i][j]+f[i-1][j+1]*(j+1)*(n-m-(m-(j+1)))%mod)%mod;
        f[i][j]=(f[i][j]+f[i-1][j]*j*(m-j)%mod)%mod;
        f[i][j]=(f[i][j]+f[i-1][j]*(m-j)*(n-m-(m-j))%mod)%mod;
        f[i][j]=(f[i][j]+f[i-1][j]*m*(m-1)/2%mod)%mod;
        f[i][j]=(f[i][j]+f[i-1][j]*(n-m)*(n-m-1)/2%mod)%mod;
    }

矩乘部分：

注意到$k$非常大，而$n$却很小，这意味着还有优化的空间

重新观察转移方程，发现$f[i][j]$可以表示为$f[i-1][j-1]k_1+f[i-1][j]k_2+f[i-1][j+1]*k_3$的形式，其中$k_1$，$k_2$，$k_3$是三个只与$j$相关的常数

这个系数恒定不变的性质，就是矩乘优化的标志

矩乘由两部分组成，$1m$状态矩阵$A$，和$mm$转移矩阵$B$

$A$描述的是每一阶段的$f$数组，下标与$f$数组的$j$意义相同，因此初始化$A_{cnt}=1$

$B$是用来将$A$从旧状态转移到新状态的，$B_{i,j}$表示从$A_i$转移到$A_j$的系数

具体地，只要仔细地根据转移方程将系数填进$B$矩阵就好了

这样答案所需要的状态$A’$就等于$A*B^k$，意思是转移了$k$次

而又因为矩乘的结合律，可以用快速幂

复杂度被优化到了$O(n^3\log k)$，其中$n^3$是矩乘的，$\log k$是快速幂的

代码：

#include <bits/stdc++.h>
using namespace std;
template<class t> inline t read(t &x){
    x=0;char c=getchar();bool f=0;
    while(!isdigit(c)) f|=c=='-',c=getchar();
    while(isdigit(c)) x=(x<<1)+(x<<3)+(c^48),c=getchar();
    if(f) x=-x;return  x;
}
template<class t> inline void write(t x){
    if(x<0){putchar('-'),write(-x);}
    else{if(x>9)write(x/10);putchar('0'+x%10);}
}

const int mod=1e9+7;
const int N=105;
int n,k,m,cnt,a[N];

struct mat{
    int a[N][N];
    mat(){memset(a,0,sizeof a);}
    inline mat operator * (const mat &nt) const { //重载矩阵乘法
        mat res;
        for(int i=0;i<=m;i++) for(int j=0;j<=m;j++) for(int k=0;k<=m;k++)
            res.a[i][j]=(1ll*res.a[i][j]+1ll*a[i][k]*nt.a[k][j]%mod)%mod;
        return res;
    }
}A,B;

mat mpow(mat x,int y){ //矩阵快速幂
    mat res;
    for(int i=0;i<=m;i++)
        res.a[i][i]=1;
    for(;y;y>>=1,x=x*x) if(y&1) res=res*x;
    return res;
}

int fpow(int x,int y){ //普通快速幂（用来求逆元）
    int res=1;
    for(;y;y>>=1,x=1ll*x*x%mod) if(y&1) res=1ll*res*x%mod;
    return res;
}

signed main(){
    read(n);read(k);
    for(int i=1;i<=n;i++){
        read(a[i]);
        m+=(a[i]==0);
    }
    for(int i=1;i<=m;i++)
        cnt+=(a[i]==0);
    A.a[0][cnt]=1;
    for(int i=0;i<=m;i++){ //枚举，把转移方程塞到转移矩阵里
        if(i!=0) B.a[i-1][i]=1ll*(m-(i-1))*(m-(i-1))%mod;
        B.a[i][i]=1ll*i*(m-i)%mod;
        B.a[i][i]=(1ll*B.a[i][i]+1ll*(m-i)*(n-m-(m-i))%mod)%mod;
        B.a[i][i]=(1ll*B.a[i][i]+1ll*m*(m-1)/2%mod)%mod;
        B.a[i][i]=(1ll*B.a[i][i]+1ll*(n-m)*(n-m-1)/2%mod)%mod;
        if(i!=m) B.a[i+1][i]=1ll*(i+1)*(n-m-(m-(i+1)))%mod;
    }
    B=mpow(B,k);
    A=A*B;
    int sum=0;
    for(int i=0;i<=m;i++)
        sum=(1ll*sum+A.a[0][i])%mod; //统计总方案数（答案的分母）
    write(1ll*A.a[0][m]*fpow(sum,mod-2)%mod);
}
赏

fighter