从k=1的特殊情况，在这种情况下，矩阵乘法时要将元素再加一次，因为这个WA了挺久。
以下是本人代码：
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define ll long long
using namespace std;
ll n,r,k,p;
struct matrix {ll a[50][50];} M[50],A;
matrix mult(matrix A,matrix B)
{
matrix S;
memset(S.a,0,sizeof(S.a));
for(int i=0;i<k;i++)
for(int j=0;j<k;j++)
for(int l=0;l<k;l++)
S.a[i][j]=(S.a[i][j]+A.a[i][l]*B.a[l][j])%p;
if (k==1) S.a[0][0]=(S.a[0][0]+A.a[0][0]*B.a[0][0])%p;
return S;
}
matrix power(ll x)
{
matrix S;
memset(S.a,0,sizeof(S.a));
for(int i=0;i<k;i++) S.a[i][i]=1;
int j=0;
while(x)
{
if (x&1) S=mult(S,M[j]);
x>>=1;j++;
}
return S;
}
int main()
{
scanf("%lld%lld%lld%lld",&n,&p,&k,&r);
memset(M[0].a,0,sizeof(M[0].a));
for(int i=0;i<k;i++)
{
M[0].a[i][(i-1+k)%k]=1;
M[0].a[i][i]=1;
}
for(int i=1;i<=45;i++) M[i]=mult(M[i-1],M[i-1]);
A=power(n*k);
printf("%lld",A.a[r][0]);
return 0;
}
