看了题解才悟了，我还是太菜了。
solution
要求
\[\left( \sum_{i = 0}^\infty C_{nk}^{ik + r} \right) \bmod p
\]
这个形式很像生成函数吧。我们套用生成函数：
\[G(x) = \sum_{i=0}^{\infty}\begin{pmatrix}nk \\ i\end{pmatrix}x^i
\]
所求即为
\[\sum_{i \bmod k = r}[x^i](1 + x)^{nk}
\]
这里有个小技巧：循环卷积。
我们卷积时，显然有：
\[x^{ik+r}x^{jk+l} \Rightarrow x^{(i+j+[l+r>k])k + (l+r) \bmod k}
\]
\(ik+r\) 即模意义下为 \(r\)，\(jk+l\) 即模意义下为 \(l\)。
我们发现他们两个的乘积放到 \(l + r \bmod k\) 了。
所以卷积时可以直接套模意义。
就是循环卷积。
然后初始化：
显然：初始时生成函数多项式为 \((1+x)\)，即 \([x^0](1+x) = [x^1](1+x) = 1\)。
当 \(k=1\) 时较特殊，所有项放到 \([x^0]\) 了，所以为 \(2\)。
记得取模。
#include<bits/stdc++.h>
#define int long long
using namespace std;
int p,k;
vector<int> operator*(vector<int> a,vector<int> b)
{
vector<int> res(k);
for(int i=0;i<k;i++)
{
for(int j=0;j<k;j++)
{
res[(i + j) % k] = (res[(i + j) % k] + a[i] * b[j] % p) % p;
}
}
return res;
}
vector<int> fpow(vector<int> a,int b)
{
vector<int> r(k);
r[0] = 1;
while(b)
{
if(b & 1)r = r * a;
a=a*a;
b>>=1;
}
return r;
}
int n,r;
signed main()
{
cin>>n>>p>>k>>r;
vector<int> a(k);
if(k == 1)a[0] = 2 % p;
else a[0] = a[1] = 1;
vector<int> ans(k);
ans = fpow(a,n * k);
cout<<ans[r];
return 0;
}