问题分析

我们需要计算：

$$\left( \sum_{i=0}^\infty \mathrm{C}_{nk}^{ik + r} \right) \bmod p $$

其中 $n \leq 10^9$，$k \leq 50$，不能直接计算组合数。

核心思想：组合意义 + 矩阵快速幂

1. 组合意义转化

设 $m = nk$，原式等价于：

• 从 $m$ 个物品中选取若干物品
• 要求选取的物品数模 $k$ 等于 $r$

即：

$$\sum_{j \equiv r \ (\text{mod } k)} \mathrm{C}_{m}^{j} $$

2. 动态规划模型

定义状态：

$$f[i][x] = \text{从前 } i \text{ 个物品中选取，物品数模 } k = x \text{ 的方案数} $$

状态转移：

• 不选第 $i$ 个物品：$f[i][x] \leftarrow f[i-1][x]$
• 选第 $i$ 个物品：$f[i][x] \leftarrow f[i-1][(x-1) \bmod k]$

即：

$$f[i][x] = f[i-1][x] + f[i-1][(x-1) \bmod k]$$

3. 矩阵形式

将状态向量记为：

$$\mathbf{A}_i = \begin{pmatrix} f[i][0] \\ f[i][1] \\ \vdots \\ f[i][k-1] \end{pmatrix} $$

则：

$$\mathbf{A}_i = M \cdot \mathbf{A}_{i-1}$$

其中 $M$ 是 $k \times k$ 的转移矩阵：

• $M_{x,x} = 1$（不选）
• $M_{x,(x-1) \bmod k} = 1$（选）

初始状态：

$$\mathbf{A}_0 = \begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix} $$

最终答案：

$$\text{ans} = \mathbf{A}_m[r] = (M^m \cdot \mathbf{A}_0)[r] $$

算法详解

数据结构

typedef long long ll1;
ll1 mod;
const int N = 105;
ll1 n, k, r;
ll1 m;  // m = n*k
ll1 tmp[N];
ll1 A[N];  // 当前结果向量
ll1 B[N];  // 转移矩阵的幂

算法步骤

步骤1：初始化

void init() {
    scanf("%lld%lld%lld%lld", &n, &mod, &k, &r);
    m = n * k;      // 总物品数
    n = k;          // 复用n作为模数k
    
    // B: 转移矩阵对应的多项式（循环卷积）
    B[0] = B[1] = 1;  // 对应 x^0 + x^1
    
    // A: 初始状态向量
    A[0] = 1;          // f[0][0] = 1
    
    if (n == 1)        // k=1的特殊情况
        B[0]++;        // x^0系数变为2
}

为什么 B[0]=B[1]=1？

• 转移方程：$f[i][x] = f[i-1][x] + f[i-1][(x-1) \bmod k]$
• 用多项式表示：$A_i(x) = (1 + x) \cdot A_{i-1}(x) \pmod{x^k - 1}$
• 这里 $1+x$ 对应系数 [1, 1, 0, ..., 0]

步骤2：循环卷积乘法

void Mul(ll1 *a, ll1 *b, ll1 *c) {
    for (int i = 0; i < n; i++) {
        tmp[i] = 0;
        for (int j = 0; j < n; j++) {
            (tmp[i] += a[j] * b[(n + (i - j)) % n] % mod) %= mod;
        }
    }
    for (int i = 0; i < n; i++)
        c[i] = tmp[i];
}

循环卷积原理：

• 在模 $x^k - 1$ 意义下，$x^k \equiv 1$
• 乘法：$c_i = \sum_{j=0}^{k-1} a_j \cdot b_{(i-j) \bmod k}$
• 这等价于 $k \times k$ 矩阵乘法，但用循环卷积实现更高效

步骤3：快速幂计算

void work() {
    for (; m; m >>= 1) {
        if (m & 1) {
            Mul(A, B, A);  // A = A * B
        }
        Mul(B, B, B);      // B = B * B
    }
    
    ll1 ans = A[r];
    printf("%lld\n", ans);
}

快速幂原理：

• 计算 $(1+x)^m \pmod{x^k - 1}$
• $A$ 初始为 1（对应 $m=0$）
• 每次将 $B$ 平方，需要时乘到 $A$ 上
• 时间复杂度：$O(k^2 \log m)$

正确性证明

定理1：多项式表示

设 $A_m(x) = \sum_{x=0}^{k-1} f[m][x] \cdot x^x$，则：

$$A_m(x) = (1+x)^m \pmod{x^k - 1}$$

证明：

• 归纳法：$A_0(x) = 1$
• $A_i(x) = (1+x) \cdot A_{i-1}(x) \pmod{x^k - 1}$
• 因为 $x^k \equiv 1$，所以是循环卷积

定理2：答案提取

$f[m][r]$ 是 $A_m(x)$ 中 $x^r$ 的系数，即：

$$\text{ans} = [x^r] (1+x)^m \pmod{x^k - 1}$$

复杂度分析

时间复杂度

• 循环卷积乘法：$O(k^2)$
• 快速幂：$O(\log m)$ 次乘法
• 总复杂度：$O(k^2 \log m)$
• 由于 $k \leq 50$，$\log m \leq 60$，完全可行

空间复杂度

• $O(k)$ 存储多项式和临时数组

示例分析

样例1

n=2, p=10007, k=2, r=0

• $m = nk = 4$
• 计算 $(1+x)^4 \pmod{x^2 - 1}$
• $(1+x)^4 = 1 + 4x + 6x^2 + 4x^3 + x^4$
• 模 $x^2-1$：$x^2 \equiv 1, x^3 \equiv x, x^4 \equiv 1$
• 化简：$(1+6+1) + (4+4)x = 8 + 8x$
• $r=0$ 的系数为 8

样例2

n=20, p=10007, k=20, r=0

• $m = 400$
• 计算 $(1+x)^{400} \pmod{x^{20} - 1}$
• 用快速幂计算，得到系数和模 10007 为 176

算法优势

1. 高效：$O(k^2 \log nk)$ 处理 $n=10^9$
2. 巧妙：将组合数求和转化为多项式问题
3. 简洁：代码简短，核心仅两个函数
4. 通用：适用于任意模数 $p$

特殊情况处理

k=1 的情况

当 $k=1$ 时：

• 转移矩阵变为 $2 \times 2$？实际上 $k=1$ 时只有模 1 的余数 0
• 转移：$f[i][0] = 2 \cdot f[i-1][0]$
• 对应多项式：$(1+x) \equiv 2 \pmod{x-1}$
• 代码中特殊处理：if (n == 1) B[0]++

总结

本题的核心技巧：

1. 组合转化：将组合数求和转化为模 $k$ 计数问题
2. 多项式表示：用循环卷积表示模 $k$ 的递推
3. 快速幂优化：$O(k^2 \log n)$ 计算 $(1+x)^m$

通过将组合问题转化为代数问题，利用循环卷积和快速幂，高效解决了大数据范围下的组合数求和问题。