發現不能直接記憶化搜尋，因為當前可能狀態可能轉移到血量+1的狀態，但可以發現如果血量為上限值時就沒有這種情況了，因此先寫出每一個狀態的轉移式，為一串其他狀態*係數+常數，然後從第n位下來把每一個狀態轉移式中>=當前狀態血量的狀態消去，再dfs即可。

程式碼：

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const ll mod=1e9+7;
const int N=1510;

ll n,p,m,k,f[N],xs[N][N],cs[N],C[N],jc[N],ijc[N],X,Y,g[N];

void Ad(ll &x,ll y)
{x+=y;if(x>=mod)x-=mod;}
void Dw(ll &x,ll y)
{x=(x>=y)?x-y:x-y+mod;}
ll qpow(ll x,ll y)
{
    ll res=1;
    while(y)
    {
        if(y&1)res=res*x%mod;
        x=x*x%mod,y>>=1;
    }
    return res;
}

void f__k(int x)
{
    ll tp;
    if(xs[x][x+1])
    {
        tp=xs[x][x+1];
        xs[x][x+1]=0;
        for(int i=x;i>=0;i--)
            Ad(xs[x][i],tp*xs[x+1][i]%mod);
        Ad(cs[x],tp*cs[x+1]%mod);
    }
    if(xs[x][x])
    {
        tp=1;
        Dw(tp,xs[x][x]);
        tp=qpow(tp,mod-2);
        xs[x][x]=0;
        for(int i=0;i<x;i++)
            xs[x][i]=xs[x][i]*tp%mod;
        cs[x]=cs[x]*tp%mod;
    }
    
}
void get_xs(ll x)
{
    for(int i=0;i<=min(k,x);i++)
        xs[x][x+1-i]=X*Y%mod*C[i]%mod*g[i]%mod;
    for(int i=0;i<=min(k,x);i++)
        Ad(xs[x][x-i],m*Y%mod*X%mod*C[i]%mod*g[i]%mod);
    cs[x]=1;
}
ll dfs(int nw)
{
    if(f[nw]!=-1)return f[nw];
    f[nw]=cs[nw];
    for(int i=nw-1;i>=0;i--)
        Ad(f[nw],xs[nw][i]*dfs(i)%mod);
    return f[nw];
}
bool check()
{
    if(k==0)return 1;
    if(m==0)
    {
        if(k==1)
        {
            if(n>1)return 1;
        }
    }
    return 0;
}
int main()
{
    int T;ll tp;
    jc[0]=1;
    for(int i=1;i<N;i++)
        jc[i]=jc[i-1]*i%mod;
    ijc[N-1]=qpow(jc[N-1],mod-2);
    for(int i=N-2;i>=0;i--)
        ijc[i]=ijc[i+1]*(i+1)%mod;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%lld%lld%lld",&n,&p,&m,&k);
        if(check()){puts("-1");continue;}
        memset(xs,0,sizeof xs),memset(cs,0,sizeof cs);
        for(int i=0;i<=p;i++)f[i]=-1;f[0]=0,tp=1;
        for(int i=0;i<=n&&i<=k;i++)
        {
            C[i]=tp*ijc[i]%mod;
            tp=tp*(k-i)%mod;
        }
        X=qpow(m+1,k),X=qpow(X,mod-2),Y=qpow(m+1,mod-2);
        for(int i=0;i<=n&&i<=k;i++)
            g[i]=qpow(m,k-i);
        for(int j=0;j<=min(k,n);j++)
            xs[n][n-j]=X*C[j]%mod*g[j]%mod;
        cs[n]=1,f__k(n);
        for(int i=n-1;i>=1;i--)
            get_xs(i),f__k(i);
        printf("%lld\n",dfs(p));
    }
}