@bzoj - 4870@ 组合数问题

目录

• @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}\times f_{n-1}(r - j \mod k) \]

从组合意义上不难理解上面的式子。
那么就可以构造矩阵 A 使得 \(a_{ij} = C_{k}^{i - j \mod k} + [i - j = 0 \mod k]\)，跑个 \(O(k^3\log n)\) 的矩阵加速。

虽然这道题数据范围看上去好像可以过，不过还有更快的优化方法。
注意到 A 是一个循环矩阵，即 \(a_{ij} = g(i - j \mod k)\)。同时可以证明 A 的幂也一定是循环矩阵。
因此，只需要维护 A 对应的 g 就可以还原出 A，而 g 的维护是 \(O(k^2)\) 的卷积，因此总时间复杂度降为 \(O(k^2\log n)\)。

如果你觉得上述算法不好理解，事实上对于本题还有另一种理解方式：

\[f_{a+b}(r) = \sum_{j=0}^{k}f_a(j)\times 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\n", g[r]);
}

@details@

事实上，这种思路还可以延伸到 \(a_{ij} = g(i - j*A)\) 之类的情况，不过和上面的略有些不同就是了。