NTT与多项式快速幂：解决大规模序列统计问题的核心技术

1. 从一个计数问题说起：序列统计的挑战

在算法竞赛和日常的算法设计中，我们常常会遇到一类看似简单、实则暗藏玄机的问题：给定一个集合 S，要求统计所有长度为 n 的序列，使得序列中每个元素都属于 S，并且序列中所有元素的乘积（或和）模某个数 m 后，等于一个特定的值 c。

“jzoj4051-序列统计”这个标题，就精准地指向了这类问题。JZOJ是信息学竞赛中一个知名的在线评测系统，而“序列统计”则是这类计数问题的经典描述。标题末尾的“【NTT】”更是点睛之笔，它直接揭示了解决这个问题的核心武器——数论变换。这暗示着，这道题的 n 可能非常大（比如 10^9 级别），集合 S 的大小也可能不小，传统的动态规划（DP）或组合数学方法会因状态空间爆炸而失效。NTT，作为快速傅里叶变换（FFT）在模意义下的高效实现，是处理这类大规模多项式卷积、从而解决组合计数问题的“核武器”。

简单来说，这个问题可以抽象为：已知一个集合 S（元素范围在 0 到 m-1 之间，m 通常是质数），求长度为 n 的序列个数，序列每个位置从 S 中取值，使得所有元素的乘积模 m 等于 c。这里的“乘积”在模 m 下运算，这引入了数论中的循环群结构，是解题的关键。

为什么这个问题有挑战性？假设我们用最朴素的 DP，设dp[i][j]表示长度为 i 的序列，乘积模 m 为 j 的方案数。那么状态数是 O(n * m)，转移需要遍历集合 S，复杂度为 O(n * m * |S|)。当 n 高达 10^9，m 为 10^3 量级时，这个复杂度是完全不可接受的。我们需要一种能够将“长度”这个维度进行指数级压缩的方法。而多项式快速幂，配合 NTT 进行加速，正是这样的方法。其核心思想是，将集合 S 的存在性表示为一个多项式，序列的构造过程转化为这个多项式的 n 次幂，而幂运算可以通过快速幂算法在 O(log n) 次多项式乘法内完成，每次多项式乘法用 NTT 加速到 O(m log m)。这样，总复杂度就从关于 n 的线性级，优化到了对数级。

接下来，我将彻底拆解这个问题，从问题转化、数学原理、到具体的 NTT 实现细节和踩坑经验，手把手带你复现整个解题过程。无论你是第一次接触 NTT，还是对多项式优化计数问题已有了解，相信都能从中获得新的启发。

2. 问题转化：从序列计数到多项式快速幂

我们首先需要把原问题“翻译”成多项式语言。这是整个解决方案的基石，理解这一步，就理解了问题的一半。

2.1 构建初始多项式

设模数为 m（通常题目保证是质数，例如 1004535809，这是一个常用的 NTT 模数）。集合 S 包含一些在 [0, m-1] 范围内的整数。我们构造一个多项式 A(x)：A(x) = a_0 + a_1*x^1 + a_2*x^2 + ... + a_{m-1}*x^{m-1}其中，系数 a_i 的定义为：如果数字 i 在集合 S 中，则 a_i = 1；否则 a_i = 0。

这个多项式 A(x) 的物理意义是什么？x^k的系数为 1，代表数字 k 可以作为序列的一个元素。那么，一个长度为 1 的序列，其乘积模 m 为 k 的方案数，正好就是 a_k。

2.2 卷积与序列拼接

现在考虑长度为 2 的序列。如何由两个长度为 1 的序列拼接而成？假设第一个数字是 i，第二个数字是 j，那么整个序列的乘积模 m 就是(i * j) mod m。 长度为 1、乘积为 i 的方案数是 a_i。 长度为 1、乘积为 j 的方案数是 a_j。 那么，所有满足(i * j) mod m = k的 (i, j) 组合，其对长度为 2、乘积为 k 的方案数的总贡献，就是sum_{i*j ≡ k (mod m)} a_i * a_j。

仔细观察这个求和形式：sum_{i*j ≡ k (mod m)} f(i) * g(j)。这像什么？这非常类似于卷积，但不是普通的加法卷积（下标 i+j），而是乘法卷积（下标 i*j）。在模 m 的剩余系中，如果 m 是质数，那么非零元构成一个乘法群。我们可以利用这个群的性质，将乘法卷积转化为加法卷积。

2.3 利用原根进行下标映射

数论中有一个关键定理：对于质数 m，存在一个原根 g，使得 g^0, g^1, ..., g^{m-2} 恰好遍历 1 到 m-1 的所有整数（模 m 意义下）。0 需要单独处理。

这个性质为我们提供了桥梁。对于任意一个非零数 x (1 <= x < m)，都存在唯一的指数 t (0 <= t < m-1)，使得x ≡ g^t (mod m)。那么，乘法x * y ≡ g^t * g^s ≡ g^{t+s} (mod m)就对应了指数上的加法(t + s) mod (m-1)。

因此，我们可以进行一个下标变换：

1. 对于原多项式 A(x)，我们将其系数 a_i 重新放置。对于 i = 0，我们单独记录（因为 0 没有对数）。对于 i > 0，我们找到它的指数 t，使得i ≡ g^t (mod m)，然后我们构造一个新的多项式 B(x)：B(x) = b_0 + b_1*x^1 + ... + b_{m-2}*x^{m-2}其中，b_t = a_i(i = g^t mod m)。换句话说，B(x) 中x^t的系数，对应原数g^t mod m是否在集合 S 中。
2. 同时，我们还需要处理 0。如果 0 在集合 S 中，它乘以任何数都会得到 0。所以 0 的存在会特殊地影响计数。

经过这个变换后，两个序列的拼接（乘积运算）在新的多项式 B(x) 上，就变成了经典的循环加法卷积：两个长度为 n 的序列，其乘积对应的指数 t 的方案数，就是B(x)与自身卷积后，x^t的系数。注意，这里的卷积是模(x^{m-1} - 1)的循环卷积，因为指数 t 是在模 (m-1) 的意义下相加的。

2.4 推广到长度 n 与多项式快速幂

根据上面的推导，长度为 2 的序列对应的多项式是B(x) * B(x) (mod x^{m-1} - 1)。 那么，长度为 n 的序列对应的多项式就是B(x)^n (mod x^{m-1} - 1)。

这里B(x)^n表示多项式 B(x) 的 n 次幂。计算多项式的高次幂，我们自然想到快速幂算法。在快速幂的每一步中，我们需要做多项式乘法，并且结果要对(x^{m-1} - 1)取模（即进行长度为 m-1 的循环卷积）。

问题的答案（乘积模 m 为 c 的方案数）就藏在这个结果多项式R(x) = B(x)^n mod (x^{m-1}-1)里。

• 如果 c == 0：方案数由两部分组成：(1) 序列中至少有一个元素是 0。这可以通过总方案数减去全为非零元素的方案数得到。总方案数是 |S|^n。全非零的方案数，需要对R(x)的所有系数求和（因为R(x)的系数对应的是乘积为非零数的情况），再乘以（如果 S 中包含 0 时）一些组合因素，具体需要仔细处理 0 的选取。
• 如果 c > 0：找到 c 对应的指数 t_c，使得c ≡ g^{t_c} (mod m)。那么答案就是R(x)中x^{t_c}项的系数。同样需要仔细处理序列中能否包含 0 的情况（如果序列包含 0，乘积必为 0，除非 n=1 且 c=0，否则不可能等于非零的 c）。通常题目会说明 S 是模 m 的子集，且可能包含 0，这就需要我们在计算时，将包含 0 的序列和纯非零序列分开统计，最后根据 c 是否为 0 来合并答案。

至此，我们将一个 O(n * m) 的 DP 问题，转化为了一个计算B(x)^n mod (x^{m-1}-1)的问题，其中 n 以对数复杂度出现在快速幂的迭代次数中。而多项式乘法和取模操作，正是 NTT 大显身手的地方。

3. NTT 核心：为什么是它，以及如何实现循环卷积

快速傅里叶变换（FFT）可以在 O(N log N) 时间内计算两个多项式的卷积，但它涉及浮点数运算，存在精度误差。数论变换（NTT）是 FFT 在有限域（模素数域）上的完美模拟，所有运算在整数模运算下进行，没有精度损失，非常适合这类需要精确计数的算法竞赛题目。

3.1 NTT 模数与原根的选择

NTT 要求模数 p 必须满足p = c * 2^k + 1的形式，并且存在原根 g，使得g^(p-1) ≡ 1 (mod p)，且g^((p-1)/n)可以作为 n 次单位根。常见的 NTT 模数有：

• 998244353 = 119 * 2^23 + 1， 原根 g=3。
• 1004535809 = 479 * 2^21 + 1，原根 g=3。
• 469762049 = 7 * 2^26 + 1， 原根 g=3。

在“序列统计”这类题中，模数 m（题目中乘积模的数）和 NTT 模数 p 通常是两个不同的质数。m 是问题本身的模数（比如 1004535809），而我们需要用 NTT 来计算多项式乘法。这里有一个关键点：我们多项式系数的值（方案数）可能会很大，我们需要对另一个模数 MOD（通常是 1e9+7 之类的）取模作为最终答案。但在 NTT 计算过程中，我们使用的是 NTT 模数 p。我们必须保证，用 p 算出来的卷积结果，在模 MOD 的意义下是正确的。这要求卷积过程中系数的真实值不能超过 p，否则就会因为模 p 丢失信息而导致错误。通常我们通过“三模数NTT”或“拆系数FFT（MTT）”来解决大系数卷积问题。但在本题的典型设定下，集合 S 的元素是 0/1，多项式 B(x) 的系数也是 0/1，卷积过程中系数增长不会太快，如果 n 不是特别大，单模数 NTT（比如就用 1004535809 同时作为 NTT 模数和答案模数）往往是可行的，这大大简化了代码。很多题目也正是这样设计的。

3.2 实现循环卷积

我们需要的运算是：C(x) = A(x) * B(x) mod (x^L - 1)，其中 L = m-1。这就是长度为 L 的循环卷积。 使用 NTT 计算普通线性卷积时，我们需要将长度扩充到至少 2L-1 以避免混叠。对于循环卷积，有一个标准技巧：

1. 将多项式 A(x) 和 B(x) 的长度视为 L（不足补0）。
2. 将它们直接进行长度为 L 的 NTT（而不是 2L）。注意，这里的“长度为 L 的 NTT”要求 L 是 2 的幂，且小于 NTT 模数 p 对应的 2^k。如果 L 不是 2 的幂，需要扩充到大于等于 L 的最小 2 的幂，记为 N。但循环卷积要求模x^L - 1，所以扩充的部分在运算后会被“卷”回来。
3. 更稳妥且通用的做法是：计算长度为 N (N >= 2L-1，且是2的幂) 的普通线性卷积，然后手动实现取模。即：
   • 计算D(x) = A(x) * B(x)(长度为 2L-1 的线性卷积，用长度 N 的 NTT 计算)。
   • 对于结果D(x)，下标为 i 的系数，要加到结果多项式C(x)的下标i mod L上。因为x^i ≡ x^{i mod L} (mod x^L - 1)。

第二种方法虽然多了一步后处理，但思路清晰，且对 L 是否为 2 的幂没有要求，更通用。在快速幂计算B(x)^n的每一步乘法中，我们都这样处理即可。

3.3 NTT 快速幂算法框架

有了以上基础，整个算法的框架就清晰了：

1. 输入：模数 m（质数），集合 S，序列长度 n，目标值 c。
2. 预处理： a. 找到模 m 的一个原根 g。 b. 预处理pow_g[i] = g^i % m和log_g[x] = t (满足 g^t % m == x)的映射表，用于数字和指数之间的快速转换。 c. 构造多项式 B(x)，长度为 L = m-1。遍历 S 中每个非零元素 x，找到t = log_g[x]，令B[t] = 1。如果 S 中包含 0，则标记has_zero = true。
3. 核心计算： a. 使用快速幂算法计算R(x) = B(x)^n mod (x^L - 1)。 - 初始化结果多项式res(x)为单位元，即res[0] = 1，其余为 0。 - 将底数多项式base(x)初始化为B(x)。 - 循环处理 n 的二进制位。如果当前位为 1，则res = multiply(res, base)，然后对结果取模mod (x^L - 1)。 - 无论当前位如何，base = multiply(base, base)，然后取模。 - 其中multiply(P, Q)函数使用 NTT 计算线性卷积，再手动合并为循环卷积。
4. 答案统计： a. 计算总方案数total = pow(|S|, n, MOD)。 b. 计算全非零序列的方案数non_zero_sum = sum(res[i]) % MOD(i from 0 to L-1)。注意，res(x)的系数对应的是乘积（以 g 的指数表示）的方案数，对所有指数求和就得到了所有非零乘积的方案总数。 c. 根据 c 的值输出： - 若c == 0：方案数 =total - non_zero_sum。注意取模防负数。 - 若c > 0：找到t_c = log_g[c]。方案数 =res[t_c]。注意，如果has_zero为真，那么任何包含 0 的序列都不可能得到非零的 c，所以这个结果就是最终答案。

这个框架将时间复杂度从 O(n * m) 优化到了 O(log n * m log m)，足以处理 n 高达 10^9 的规模。

4. 关键实现细节与避坑指南

理论通了，代码实现才是真正的战场。下面我结合自己的踩坑经验，梳理几个关键细节。

4.1 原根的寻找与验证

题目给定的 m 通常是质数，且是常见的 NTT 模数（如 1004535809）。对于这类模数，原根 g=3 几乎总是对的。但为了代码的鲁棒性，实现一个原根查找函数是好的实践。

一个简单的查找方法（适用于质数 p）：

1. 对 p-1 进行质因数分解，得到不同的质因子fac[]。
2. 从 g=2 开始枚举，对于每个质因子f，检查g^((p-1)/f) mod p != 1。如果对所有质因子都满足，则 g 是原根。

// 快速幂函数 long long qpow(long long a, long long b, long long p); // 查找原根 int find_primitive_root(int p) { // 分解 p-1 的质因子 vector<int> factors; int tmp = p - 1; for (int i = 2; i * i <= tmp; ++i) { if (tmp % i == 0) { factors.push_back(i); while (tmp % i == 0) tmp /= i; } } if (tmp > 1) factors.push_back(tmp); // 枚举原根 for (int g = 2; g < p; ++g) { bool ok = true; for (int f : factors) { if (qpow(g, (p - 1) / f, p) == 1) { ok = false; break; } } if (ok) return g; } return -1; // 理论上不会发生 }

在竞赛中，如果已知模数，直接硬编码原根可以节省时间。例如const int G = 3;。

4.2 下标映射表的构建

我们需要在数字和它的离散对数（指数）之间快速转换。由于 m 通常在 10^3 到 10^4 量级，我们可以直接预处理两个数组：

• pow_g[i] = (g^i) % m， i 从 0 到 m-2。
• log_g[x] = i，其中 x = pow_g[i]。对于 x=0，可以特殊处理为 -1。

vector<int> pow_g(m), log_g(m, -1); int cur = 1; // g^0 % m for (int i = 0; i < m-1; ++i) { pow_g[i] = cur; log_g[cur] = i; cur = (1LL * cur * g) % m; }

这样，给定数字 x，log_g[x]就是其指数；给定指数 i，pow_g[i]就是对应的数字。复杂度 O(m)。

4.3 NTT 卷积与循环卷积取模的实现

这是代码的核心部分。我们需要实现一个函数multiply_circular(vector<int> &a, vector<int> &b, int L, int mod)，它返回a和b的循环卷积（模x^L - 1），系数对mod取模。

// 假设已有完整的 NTT 板子，包括 ntt() 函数和逆变换 intt() 函数 // 以及计算 rev 数组的函数 init_rev() const int NTT_MOD = 1004535809; // 常用的NTT模数 const int G = 3; // 对应原根 void ntt(vector<long long> &a, int opt) { /* ... */ } // opt=1 正变换， opt=-1 逆变换 vector<int> multiply_circular(const vector<int> &a, const vector<int> &b, int L, int MOD) { int len_a = a.size(), len_b = b.size(); // 计算线性卷积需要的长度 int len = 1; while (len < len_a + len_b - 1) len <<= 1; vector<long long> A(len, 0), B(len, 0); copy(a.begin(), a.end(), A.begin()); copy(b.begin(), b.end(), B.begin()); // NTT 计算线性卷积 ntt(A, 1); ntt(B, 1); for (int i = 0; i < len; ++i) A[i] = A[i] * B[i] % NTT_MOD; ntt(A, -1); // 此时 A 中存储的是线性卷积结果，长度为 len_a+len_b-1 // 现在将其转换为循环卷积 (模 x^L - 1) vector<int> res(L, 0); for (int i = 0; i < (int)A.size(); ++i) { // 线性卷积的第 i 项贡献到循环卷积的第 i % L 项 res[i % L] = (res[i % L] + A[i]) % MOD; // 注意：A[i] 是模 NTT_MOD 下的值，但 res 是模 MOD 的。 // 这里直接加，是因为我们确信 A[i] 的真实值（在模 MOD 前）不会太大。 // 如果担心，可以在这里用 long long 累加，最后再取模 MOD。 } return res; }

一个极其重要的坑点：NTT_MOD和答案的MOD是不同的。NTT 运算在模NTT_MOD下进行，得到的结果A[i]也是模NTT_MOD后的值。而我们的最终答案需要对MOD取模。我们直接使用A[i]的值参与后续运算，这成立的前提是：多项式系数在真实值（不取任何模）的情况下，不超过NTT_MOD。在本问题中，初始系数是 0/1，经过 n 次幂后，系数的真实值最大可能是C(n+L-1, L-1)量级（可重组合）。当 n 很大时，这个值可能远超NTT_MOD。因此，单模数 NTT 可能溢出。

解决方案：使用“任意模数NTT”（MTT）或者“三模数NTT+CRT”来保证在模MOD意义下的正确性。这是此类题目最易错、也最需要关注的地方。下面简要介绍三模数NTT的思路：

1. 选择三个 NTT 模数p1, p2, p3，例如 998244353, 1004535809, 469762049。
2. 分别用这三个模数做一遍上述的 NTT 卷积和循环卷积取模，得到三组结果res1, res2, res3（每个都是模各自 NTT 模数下的值）。
3. 对于最终结果res的每一位系数r，我们已知：r ≡ res1[i] (mod p1)r ≡ res2[i] (mod p2)r ≡ res3[i] (mod p3)
4. 利用中国剩余定理（CRT）合并这三个同余方程，得到一个模M = p1*p2*p3下的解。由于p1*p2*p3通常远大于我们需要的系数真实值上限和MOD，我们可以得到系数r的真实值（在 long long 范围内）。
5. 最后将真实值r对题目要求的MOD取模，得到最终系数。

实现三模数 NTT 会显著增加代码量，但它是处理大系数多项式乘法的标准且可靠的方法。在竞赛中，如果确认系数不会爆炸，用单模数偷懒是可行的，但理解其风险至关重要。

4.4 零元素的特殊处理

如果集合 S 中包含 0，计数会变得复杂。因为长度为 n 的序列，如果至少有一个位置是 0，整个序列的乘积就是 0。我们需要分情况讨论：

• 目标 c = 0：合法序列包括：(1) 至少有一个 0 的序列。(2) 所有元素非零，但乘积模 m 恰好为 0 的序列？在模质数 m 下，如果所有元素非零，它们的乘积也必然非零（因为模质数的乘法构成群）。所以非零序列的乘积不可能为 0。因此，c=0的方案数就是总序列数减去全非零序列数。ans = (total - non_zero_total + MOD) % MOD其中total = |S|^n，non_zero_total是B(x)^n的所有系数之和（即所有非零乘积的方案数）。注意，B(x)的构造只包含了 S 中的非零元素。
• 目标 c > 0：合法序列必须全部由非零元素构成。所以答案直接就是B(x)^n中对应指数t_c的系数。如果 S 中包含 0，那么任何包含 0 的序列都不合法，不会影响这个结果。

因此，在预处理构造B(x)时，我们只将 S 中的非零元素映射进去。has_zero标志位只用于最后计算总方案数total和答案合并。

4.5 快速幂中的多项式单位元

在计算多项式快速幂时，初始结果res应设为单位元。对于普通多项式乘法，单位元是[1]（即常数项为1，其他为0）。对于循环卷积mod (x^L - 1)，单位元同样是[1]。因为[1]与任何多项式P(x)做循环卷积，结果仍是P(x)。你可以验证：1 * P(x) mod (x^L - 1) = P(x)。

所以初始化很简单：vector<int> res(L, 0); res[0] = 1;。

5. 完整代码结构与调试要点

将以上所有部分组合起来，下面给出一个高层次的伪代码框架，并附上关键调试点。

#include <bits/stdc++.h> using namespace std; typedef long long ll; const int MOD = 1000000007; // 答案模数 const int NTT_MOD1 = 998244353, G1 = 3; const int NTT_MOD2 = 1004535809, G2 = 3; const int NTT_MOD3 = 469762049, G3 = 3; // 需要实现三套 NTT 函数，或者写一个模板函数 // 1. 辅助函数：快速幂，原根查找，离散对数表生成 ll qpow(ll a, ll b, ll p); int find_root(int p); void build_log_table(int m, int g, vector<int>& pow_g, vector<int>& log_g); // 2. 核心函数：三模数循环卷积 vector<int> multiply_circular_3mod(const vector<int>& a, const vector<int>& b, int L); // 3. 主函数 int main() { int n, m, c, s_size; // 序列长度，模数，目标值，集合大小 vector<int> S(s_size); // 读入 n, m, c, S ... // Step 1: 找到模 m 的原根 g，并建立映射表 int g = find_root(m); // 或硬编码已知原根 vector<int> pow_g(m), log_g(m, -1); build_log_table(m, g, pow_g, log_g); // Step 2: 构造多项式 B(x)，长度 L = m-1 int L = m - 1; vector<int> B(L, 0); bool has_zero = false; for (int x : S) { if (x % m == 0) { // 处理 0 has_zero = true; } else { int t = log_g[x % m]; B[t] = (B[t] + 1) % MOD; // 注意，如果S中同一个数出现多次？通常集合元素不重复。若可重复，此处是计数。 } } // Step 3: 计算 B(x)^n mod (x^L - 1)，使用快速幂 vector<int> res(L, 0); res[0] = 1; // 单位元 vector<int> base = B; ll exponent = n; // n 可能很大，用 ll while (exponent) { if (exponent & 1) { res = multiply_circular_3mod(res, base, L); } base = multiply_circular_3mod(base, base, L); exponent >>= 1; } // 现在 res 是 B(x)^n 的循环卷积结果，系数模 MOD // Step 4: 计算答案 ll total = qpow((int)S.size(), n, MOD); // 总方案数 |S|^n ll non_zero_total = 0; for (int coeff : res) non_zero_total = (non_zero_total + coeff) % MOD; ll ans = 0; if (c == 0) { // 方案数 = 总方案 - 全非零方案 ans = (total - non_zero_total) % MOD; if (ans < 0) ans += MOD; } else { // c > 0，必须全为非零元素 int t_c = log_g[c % m]; if (t_c == -1) { // c 不在模 m 的乘法群中？c=0 已处理，c>0且m为质数，则c一定在。 ans = 0; } else { ans = res[t_c]; } // 如果 has_zero 为真，那么包含零的序列已经自动被排除，因为 B(x) 不含零。 } cout << ans << endl; return 0; }

调试与验证要点：

1. 小数据暴力对拍：这是最有效的方法。写一个暴力 DP 程序，处理 n, m, |S| 都很小的情况（比如 n<=10, m<=10）。用随机生成的 S 和 n, c，对比你的 NTT 优化程序和暴力程序的输出。确保在数百组随机数据下完全一致。
2. 检查原根和映射表：对于给定的 m，验证你找到的 g 是否真的是原根（g^(m-1) ≡ 1，且对于 m-1 的所有真因子 d，g^d !≡ 1）。验证pow_g和log_g表是否正确，确保pow_g[log_g[x]] == x对所有非零 x 成立。
3. 验证循环卷积：单独测试multiply_circular函数。用两个小的多项式，手动计算循环卷积，与程序结果对比。
4. 验证快速幂：测试B(x)^1,B(x)^2，与暴力乘法的结果对比。
5. 关注模数：清晰区分NTT_MOD1/2/3和答案MOD。在 CRT 合并时，确保计算正确，最终系数取模MOD。
6. 零值处理：专门构造包含 0 和不包含 0 的 S，分别测试 c=0 和 c>0 的情况，验证计数逻辑是否正确。

6. 性能优化与扩展思考

对于竞赛场景，性能至关重要。这里有一些优化方向：

1. 预处理 NTT 旋转因子：在多次 NTT 调用中（快速幂的每次乘法都需要多次 NTT），可以预先计算好每个长度对应的旋转因子wn[]数组，避免重复计算。
2. 使用迭代 NTT：递归 NTT 简洁但常数大。迭代 NTT 通过位逆序置换实现，速度更快，是竞赛标准板子。
3. 减少拷贝：在multiply_circular函数中，避免不必要的向量拷贝。可以尝试原地操作或使用指针。
4. 选择适当的卷积长度：在快速幂过程中，res和base多项式的长度始终是 L。但线性卷积需要长度至少为2L-1。我们可以预先计算好需要的 NTT 长度N，并一次性初始化好对应的 rev 数组和旋转因子。

扩展思考：

1. 如果 m 不是质数：问题会变得复杂很多。因为模 m 的乘法不构成群，单位元 1 的逆元不一定存在。通常的解法是分解 m，然后利用中国剩余定理（CRT）将问题分解为模数为质数幂的子问题，再组合答案。这超出了本文范围，但知道这个方向很重要。
2. 如果操作不是乘法而是加法：问题就变成了更经典的“背包”计数，可以用生成函数和普通多项式快速幂解决，且不需要原根和下标映射，直接用 NTT/FFT 计算加法卷积的幂即可。
3. 动态集合 S：如果 S 不是固定的，而是可以动态添加删除元素，问题就变成了动态维护多项式幂。这可能需要用到更高级的数据结构，如线段树维护分治 NTT。

回过头看，“序列统计”这个问题是连接组合计数、数论（原根、离散对数）和多项式算法（NTT）的绝佳桥梁。它迫使你跳出 DP 的思维定式，用生成函数的视角看待序列构造，再用数论工具将乘法群同构到加法群，最终用 NTT 这把利刃高效解决。理解并实现一次，你对多项式优化和模运算计数的认识会上一个全新的台阶。在调试通过的那一刻，看着屏幕上对于超大 n 瞬间输出的正确结果，那种成就感，正是算法竞赛最迷人的地方之一。