@desription@
求：
[sum_{i=0}^{infty}C_{nk}^{ik + r}mod p
]
@solution@
这不是单位根反演裸题吗
注意到模数并没有特殊性质，可以想到从组合数的递推公式入手。
当 k = 1 实际上就是快速幂，所以联想到说不定可以矩阵幂优化。
记 (f_n(r) = sum_{i=0}^{infty}C_{nk}^{ik + r})，则：
[f_n(r) = sum_{j=0}^{k}C_{k}^{j} imes f_{n-1}(r - j mod k)
]
从组合意义上不难理解上面的式子。
那么就可以构造矩阵 A 使得 (a_{ij} = C_{k}^{i - j mod k} + [i - j = 0 mod k])，跑个 (O(k^3log n)) 的矩阵加速。
虽然这道题数据范围看上去好像可以过，不过还有更快的优化方法。
注意到 A 是一个循环矩阵，即 (a_{ij} = g(i - j mod k))。同时可以证明 A 的幂也一定是循环矩阵。
因此，只需要维护 A 对应的 g 就可以还原出 A，而 g 的维护是 (O(k^2)) 的卷积，因此总时间复杂度降为 (O(k^2log n))。
如果你觉得上述算法不好理解，事实上对于本题还有另一种理解方式：
[f_{a+b}(r) = sum_{j=0}^{k}f_a(j) imes f_b(r - j mod k)
]
一样地，从组合意义不难理解上式。
也就是 dp 状态之间可以作 “复合” 运算。那么直接对 dp 做倍增，然后把需要的 2 的幂 “复合” 起来即可。
嘛，不过这种理解方法和上面那个是等价的就是了。
@accepted code@
#include <cstdio>
#include <algorithm>
int n, p, k, r;
inline int add(int x, int y) {return (x + y >= p ? x + y - p : x + y);}
inline int mul(int x, int y) {return 1LL * x * y % p;}
int t[55];
void merge(int *A, int *B, int *C) {
for(int i=0;i<k;i++) t[i] = 0;
for(int i=0;i<k;i++)
for(int j=0;j<k;j++) {
int p = (i + j >= k ? i + j - k : i + j);
t[p] = add(t[p], mul(A[i], B[j]));
}
for(int i=0;i<k;i++) C[i] = t[i];
}
int f[55], g[55];
void get(int n) {
for(int i=0;i<k;i++) g[i] = (i == 0);
for(int i=n;i;i>>=1,merge(f, f, f))
if( i & 1 ) merge(f, g, g);
}
int c[55][55];
int main() {
scanf("%d%d%d%d", &n, &p, &k, &r);
for(int i=0;i<=k;i++) {
c[i][0] = 1;
for(int j=1;j<=i;j++)
c[i][j] = add(c[i-1][j], c[i-1][j-1]);
}
for(int i=0;i<k;i++)
f[i] = c[k][i];
f[0] = add(f[0], c[k][k]);
get(n), printf("%d
", g[r]);
}
@details@
事实上，这种思路还可以延伸到 (a_{ij} = g(i - j*A)) 之类的情况，不过和上面的略有些不同就是了。