@bzoj

目录

• @desription@
• @solution@
• @accepted code@
• @details@

---

@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)) 之类的情况，不过和上面的略有些不同就是了。