萌新刚学OI 114ms，求助

#include<bits/stdc++.h>
using namespace std;
const int MAXN=1505;
const long long mod=1000000007;
long long f[MAXN],A[MAXN][MAXN],m1,m2;
long long T,n,p,m,k;
long long q_power(long long a,long long b){
	long long ans=1;
	while(b){
		if(b&1)ans*=a,ans%=mod;
		a=a*a,a%=mod;
		b>>=1;
	}
	return ans;
}
long long get_m(long long a,long long b){
	return (a*q_power(b,mod-2))%mod;
}//return a/b % mod
void init_f(){
	f[0]=q_power(get_m(m,m+1),k);
	//cout<<f[0]; 
	long long mm=get_m(1,m);
	for(int i=1;i<=min(n+1,k);i++){
		f[i]=get_m(k-i+1,i)*mm%mod*f[i-1]%mod;
		f[i]%=mod;//cout<<f[i]<<" ";
	}
	//for(int i=1;i<=k;i++)
	//cout<<endl;
}
void make_square(){
	for(int i=1;i<n;i++){ 
		for(int j=1;j<=i;j++)A[i][j]=(f[i-j+1]%mod+m*f[i-j]%mod)%mod*m1%mod;
		A[i][i+1]=f[0]*m1%mod;
		for(int j=i+2;j<=n;j++)A[i][j]=0;
	}
	for(int i=1;i<=n;i++)A[n][i]=f[n-i]%mod;
	for(int i=1;i<=n;i++)A[i][i]=(A[i][i]-1+mod)%mod;
	for(int i=1;i<=n;i++)A[i][n+1]=mod-1;
}
void gauss(){
/*	for(int i=1;i<=n;i++){
		for(int j=1;j<=n+1;j++)cout<<A[i][j]<<" ";
		cout<<endl;
	}*/
	for(int i=1;i<=n;i++){
		long long k=get_m(1,A[i][i]);
		A[i][i]=1;
		if(i!=n)A[i][i+1]=(A[i][i+1]*k)%mod;
		A[i][n+1]=(A[i][n+1]*k)%mod;
		for(int j=i+1;j<=n;j++){		
			A[j][i+1]=(A[j][i+1]-(A[i][i+1]*A[j][i])%mod+mod)%mod;
			A[j][n+1]=(A[j][n+1]-(A[i][n+1]*A[j][i])%mod+mod)%mod;
			A[j][i]=0;
		//	cout<<A[j][i]<<" "<<A[j][i+1]<<" "<<A[j][n+1]<<endl;
		}
	}
/*	for(int i=1;i<=n;i++){
		for(int j=1;j<=n+1;j++)cout<<A[i][j]<<" ";
		cout<<endl;
	}*/
	for(int i=n;i>1;i--){
		A[i-1][n+1]=(A[i-1][n+1]-(A[i][n+1]*A[i-1][i])%mod+mod)%mod;
		A[i-1][i]=0;
	}
	/*(for(int i=1;i<=n;i++){
		for(int j=1;j<=n+1;j++)cout<<A[i][j]<<" ";
		cout<<endl;
	}*/
}
int main(){
	cin>>T;
	while(T--){
		cin>>n>>p>>m>>k;
		if(k==0){
			cout<<-1<<endl;
			continue;
		}else if(m==0){
			if(k==1)cout<<-1<<endl;
			else{
				long long ans=0;
				while(p>0){
					if(p<n)p++;
					p-=k;
					ans++;
				}
				cout<<ans<<endl;
			}
			continue;
		}
		m1=get_m(1,m+1),m2=get_m(m,m+1);
		init_f();
		make_square();
		gauss();
		cout<<A[p][n+1]<<endl;
	}  
	return 0;
}