题面

式子看起来很吓人哈

不过我们来考虑它的实际意义

就是在 nk 个物品中选择 模k余数为r 的方案数和

设 \(f[i,j]\) 表示选 在 i 个物品中选择 模k余数为j 的方案数和

容易得到 \(f[i,j]=f[i-1,j-1]+f[i-1,j]\)

当然，为了避免负数，我们写成这种模样 \(f[i,j]=f[i-1,(j+k-1)\%k]+f[i-1,j]\)

这种形式可以写成矩阵

然后就矩阵快速幂 啦啦啦啦

---

当然，我们可以继续优化

发现这个式子也可以不由 i-1转移过来,而是由两个状态合并而来

\(f[x+y,k]=f[x,0]f[y,k]+f[x,1]f[y,k-1]...f[x,k]f[y,0]\)

然后可以dp快速幂了

---

然后说如果p特殊，可以用单位根反演？(%%%boshi)

---

就只上第一种做法的代码吧（准确来说我只会这一种）

注意当k=1时，矩阵会退化为单个变量，用赋值的方法会挂；应该自加

#include<bits/stdc++.h>
using namespace std;
#define int long long
int const MAXN=60;
int n,p,k,r;
struct Mat{
	int m[MAXN][MAXN];
	Mat(){
		memset(m,0,sizeof(m));
	}
	Mat operator*(const Mat& b) const{
		Mat x;
		for(int i=0;i<k;i++){
			for(int l=0;l<k;l++){
				int a=m[i][l];if(!a)continue;
				for(int j=0;j<k;j++){
					x.m[i][j]+=1ll*a*b.m[l][j]%p;
					if(x.m[i][j]>=p)x.m[i][j]-=p;
				}
			}
		}
		return x;
	}
	void init(){
		for(int i=0;i<k;i++)m[i][i]=1;
	}
	void print(){
		for(int i=0;i<k;i++){
			for(int j=0;j<k;j++){
				printf("%d ",m[i][j]);
			}
			printf("\n");
		}
	}
}a,b;
void mpow(Mat &x,int b){
	Mat ans;ans.init();
	while(b){
		if(b&1)ans=x*ans;
		x=x*x;
		b>>=1;
	}
	x=ans;
}
signed main(){
	scanf("%lld%lld%lld%lld",&n,&p,&k,&r);
	a.m[0][0]=1;
	for(int i=0;i<k;i++){
		b.m[i][i]++;
		b.m[i][(i+k-1)%k]++;
	}
	mpow(b,n*k);
	a=b*a;
	printf("%lld\n",a.m[r][0]); 
	return 0;
}