考虑组合数递推公式(inom{n}{m}=inom{n-1}{m}+inom{n-1}{m-1}),然后代入原式,有
(sum_{i=0}^{infty} inom{nk}{ik+r}=sum_{i=0}^{infty} inom{nk-1}{ik+r}+sum_{i=0}^{infty} inom{nk-1}{ik+r-1})
设(f_{i,j}=sum_{i=0}^{infty} inom{i}{ik+j}),那么就会有(f_{i,j}=f_{i-1,j}+f_{i-1,(j-1)mod k}),就可以矩乘转移了.或者可以发现(f_{i,j}=[x^j]((x+1)^ipmod{x^k-1})),使用倍增+FFT可以做到(O(klog nlog k))
#include<bits/stdc++.h>
#define LL long long
#define db double
using namespace std;
const int N=1e6+10;
int rd()
{
int x=0,w=1;char ch=0;
while(ch<'0'||ch>'9'){if(ch=='-') w=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){x=x*10+(ch^48);ch=getchar();}
return x*w;
}
int n,mod,k,r;
int fpow(int a,int b,int mod){int an=1;while(b){if(b&1) an=1ll*an*a%mod;a=1ll*a*a%mod,b>>=1;}return an;}
struct matrix
{
int a[50][50];
matrix(int az=0){memset(a,0,sizeof(a));for(int i=0;i<50;++i) a[i][i]=az;}
matrix operator * (const matrix &bb) const
{
matrix an;
for(int i=0;i<k;++i)
for(int l=0;l<k;++l)
for(int j=0;j<k;++j)
an.a[i][j]=(1ll*a[i][l]*bb.a[l][j]+an.a[i][j])%mod;
return an;
}
}aa,bb;
int main()
{
n=rd(),mod=rd(),k=rd(),r=rd();
for(int i=0;i<k;++i) ++bb.a[i][i],++bb.a[i][(i+1)%k];
aa.a[0][0]=1;
LL rs=1ll*n*k;
while(rs)
{
if(rs&1) aa=aa*bb;
bb=bb*bb,rs>>=1;
}
printf("%d
",aa.a[0][r]);
return 0;
}