[BZOJ5292][BJOI2018]治疗之雨(概率DP+高斯消元)

https://blog.csdn.net/xyz32768/article/details/83217209

不难找到DP方程与辅助DP方程，发现DP方程具有后效性，于是高斯消元即可。

但朴素消元显然无法通过，注意到f[i]的方程至多与f[i+1]有关，于是从下往上依次消去最后一个数，剩下的就是一个下三角，直接求解即可。

注意中间与指数有关的计算能预处理的就不用快速幂，以及阶乘等值可以在程序开头预处理。

复杂度$O(n^2)$，不知道为什么和别人的代码相比常数巨大。

#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
using namespace std;

const int N=3010,mod=1e9+7;
int n,m,p,k,T,d[N],pw[N],fac[N],inv[N],P[N][N],a[N][N];

int ksm(int a,int b){
    int res=1;
    for (; b; a=1ll*a*a%mod,b>>=1)
        if (b & 1) res=1ll*res*a%mod;
    return res;
}

bool Gauss(){
    for (int i=n; i; i--){
        if (!a[i][i]) return 0;
        int t=1ll*a[i-1][i]*ksm(a[i][i],mod-2)%mod;
        rep(j,0,i) a[i-1][j]=(a[i-1][j]-1ll*t*a[i][j]%mod+mod)%mod;
        a[i-1][n+1]=(a[i-1][n+1]-1ll*t*a[i][n+1]%mod+mod)%mod;
    }
    rep(i,1,n){
        rep(j,0,i-1) a[i][n+1]=(a[i][n+1]-1ll*a[i][j]*a[j][n+1]%mod+mod)%mod;
        a[i][n+1]=1ll*a[i][n+1]*ksm(a[i][i],mod-2)%mod;
    }
    return 1;
}

int main(){
    freopen("heal.in","r",stdin);
    freopen("heal.out","w",stdout);
    n=1500;
    fac[0]=1; rep(i,1,n) fac[i]=1ll*fac[i-1]*i%mod;
    inv[n]=ksm(fac[n],mod-2);
    for (int i=n-1; ~i; i--) inv[i]=1ll*inv[i+1]*(i+1)%mod;
    for (scanf("%d",&T); T--; ){
        scanf("%d%d%d%d",&n,&p,&m,&k);
        rep(i,0,n+1) rep(j,0,n+1) a[i][j]=P[i][j]=0;
        d[0]=1; rep(i,1,min(n,k)) d[i]=1ll*d[i-1]*(k-i+1)%mod;
        pw[0]=ksm(m,k); int t=ksm(m,mod-2);
        rep(i,1,min(n,k)) pw[i]=1ll*pw[i-1]*t%mod;
        if (k<=n) pw[k]=1;
        if (!k || (!m && k==1)){ puts("-1"); continue; }
        t=ksm(ksm(m+1,k),mod-2);
        rep(i,1,n){
            int sm=0;
            rep(j,0,min(i,k))
                P[i][j]=(i==j)?(1-sm+mod)%mod:1ll*d[j]*inv[j]%mod*pw[j]%mod*t%mod,sm+=P[i][j];
        }
        a[0][0]=1; int inv=ksm(m+1,mod-2);
        rep(i,1,n-1){
            a[i][n+1]=a[i][i]=mod-1;
            rep(j,0,i+1){
                a[i][i-j+1]=(a[i][i-j+1]+1ll*P[i+1][j]*inv)%mod;
                a[i][i-j]=(a[i][i-j]+1ll*P[i][j]*inv%mod*m)%mod;
            }
        }
        a[n][n+1]=a[n][n]=mod-1;
        rep(j,0,n) a[n][n-j]=(a[n][n-j]+P[n][j])%mod;
        if (Gauss()) printf("%d\n",a[p][n+1]); else puts("-1");
    }
    return 0;
}