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

不難找到DP方程與輔助DP方程，發現DP方程具有後效性，於是高斯消元即可。

但樸素消元顯然無法通過，註意到f[i]的方程至多與f[i+1]有關，於是從下往上依次消去最後一個數，剩下的就是一個下三角，直接求解即可。

註意中間與指數有關的計算能預處理的就不用快速冪，以及階乘等值可以在程序開頭預處理。

復雜度$O(n^2)$，不知道為什麽和別人的代碼相比常數巨大。

 1 #include<cstdio>
 2 #include<algorithm>
 3 #define
 
 rep(i,l,r) for (int i=(l); i<=(r); i++)
 4 using namespace std;
 5 
 6 const int N=3010,mod=1e9+7;
 7 int n,m,p,k,T,d[N],pw[N],fac[N],inv[N],P[N][N],a[N][N];
 8 
 9 int ksm(int a,int b){
10     int res=1;
11     for (; b; a=1ll*a*a%mod,b>>=1)
12         if (b & 1) res=1ll*res*a%mod;
 
13     return res;
14 }
15 
16 bool Gauss(){
17     for (int i=n; i; i--){
18         if (!a[i][i]) return 0;
19         int t=1ll*a[i-1][i]*ksm(a[i][i],mod-2)%mod;
20         rep(j,0,i) a[i-1][j]=(a[i-1][j]-1ll*t*a[i][j]%mod+mod)%mod;
21         a[i-1][n+1]=(a[i-1][n+1]-1ll*t*a[i][n+1]%mod+mod)%mod;
 
22     }
23     rep(i,1,n){
24         rep(j,0,i-1) a[i][n+1]=(a[i][n+1]-1ll*a[i][j]*a[j][n+1]%mod+mod)%mod;
25         a[i][n+1]=1ll*a[i][n+1]*ksm(a[i][i],mod-2)%mod;
26     }
27     return 1;
28 }
29 
30 int main(){
31     freopen("heal.in","r",stdin);
32     freopen("heal.out","w",stdout);
33     n=1500;
34     fac[0]=1; rep(i,1,n) fac[i]=1ll*fac[i-1]*i%mod;
35     inv[n]=ksm(fac[n],mod-2);
36     for (int i=n-1; ~i; i--) inv[i]=1ll*inv[i+1]*(i+1)%mod;
37     for (scanf("%d",&T); T--; ){
38         scanf("%d%d%d%d",&n,&p,&m,&k);
39         rep(i,0,n+1) rep(j,0,n+1) a[i][j]=P[i][j]=0;
40         d[0]=1; rep(i,1,min(n,k)) d[i]=1ll*d[i-1]*(k-i+1)%mod;
41         pw[0]=ksm(m,k); int t=ksm(m,mod-2);
42         rep(i,1,min(n,k)) pw[i]=1ll*pw[i-1]*t%mod;
43         if (k<=n) pw[k]=1;
44         if (!k || (!m && k==1)){ puts("-1"); continue; }
45         t=ksm(ksm(m+1,k),mod-2);
46         rep(i,1,n){
47             int sm=0;
48             rep(j,0,min(i,k))
49                 P[i][j]=(i==j)?(1-sm+mod)%mod:1ll*d[j]*inv[j]%mod*pw[j]%mod*t%mod,sm+=P[i][j];
50         }
51         a[0][0]=1; int inv=ksm(m+1,mod-2);
52         rep(i,1,n-1){
53             a[i][n+1]=a[i][i]=mod-1;
54             rep(j,0,i+1){
55                 a[i][i-j+1]=(a[i][i-j+1]+1ll*P[i+1][j]*inv)%mod;
56                 a[i][i-j]=(a[i][i-j]+1ll*P[i][j]*inv%mod*m)%mod;
57             }
58         }
59         a[n][n+1]=a[n][n]=mod-1;
60         rep(j,0,n) a[n][n-j]=(a[n][n-j]+P[n][j])%mod;
61         if (Gauss()) printf("%d\n",a[p][n+1]); else puts("-1");
62     }
63     return 0;
64 }

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