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