1. 从一道“省队互测”题说起当数论遇上组合数学最近在整理一些经典的数论题目时又翻到了这道P6271。题目名字很有意思“一个人的数论”听起来有点孤独但内容却一点也不孤单它把数论里几个重量级的概念——莫比乌斯反演和伯努利数——给串了起来。很多朋友第一次看到这个组合可能会有点懵莫比乌斯反演不是用来处理数论函数求和、化简gcd、lcm相关式子的吗伯努利数不是出现在幂和公式或者泰勒展开里的吗它们俩怎么会扯上关系这道题的精妙之处就在于此。它表面上可能只是让你求一个带幂次的因子和函数但内核却是一个经典的“数论函数与多项式”的桥梁问题。简单来说题目通常会给你一个关于正整数n的函数比如所有与n互质的数的k次幂之和记作 $S_k(n) \sum_{\substack{1 \le m \le n \ \gcd(m, n)1}} m^k$。我们的任务就是求出这个 $S_k(n)$ 的一个简洁表达式特别是当n以质因数分解形式 $n \prod_{i1}^w p_i^{\alpha_i}$ 给出时如何快速计算。直接枚举1到n的所有数并判断互质复杂度是O(n)当n很大时比如是多个大质数的高次幂乘积完全不可行。这时莫比乌斯反演就登场了它能帮我们把“互质”这个条件巧妙地转化掉。而伯努利数则为我们提供了计算自然数幂和公式的系数让我们能把一个求和式写成一个关于n的多项式。两者一结合就能把原问题转化为一个关于n的各质因子幂次的多项式求值问题复杂度瞬间降低到只与质因子个数w和幂次k有关实现了从“一个人的暴力枚举”到“优雅公式推导”的飞跃。这篇文章我就想结合这道题把莫比乌斯反演和伯努利数在这个场景下的应用逻辑彻底拆解清楚。我们不止看公式更要看每一步“为什么这么做”以及在实际推导和计算中会遇到哪些坑怎么绕过去。无论你是正在备战算法竞赛还是单纯对数学的美感有兴趣相信这个“数论与组合数学的联姻”故事都会让你有所收获。2. 核心武器拆解莫比乌斯反演与伯努利数在进入具体解题之前我们必须把两件核心武器的原理和使用方法摸透。很多资料只给结论但不知道背后的动机用起来就会束手束脚。2.1 莫比乌斯反演化“互质”为“整除”莫比乌斯函数 $\mu(n)$ 大家应该不陌生当n有平方因子时为0否则若n有偶数个质因子则为1奇数个则为-1。它的核心价值在于这个性质$\sum_{d|n} \mu(d) [n1]$其中 $[P]$ 是艾弗森括号当命题P为真时值为1否则为0。这个看似简单的式子正是解决“互质”条件的关键。考虑我们想求所有满足 $\gcd(m, n)1$ 的m的和。注意到 $[\gcd(m, n)1]$ 这个条件利用上述性质我们可以把它改写 $$[\gcd(m, n)1] \sum_{d|\gcd(m, n)} \mu(d)$$ 而 $d|\gcd(m, n)$ 等价于 $d|m$ 且 $d|n$。于是我们的目标求和式就变成了 $$S_k(n) \sum_{m1}^{n} m^k \cdot [\gcd(m, n)1] \sum_{m1}^{n} m^k \sum_{\substack{d|m \ d|n}} \mu(d)$$接下来是关键一步交换求和顺序。先对d求和d是n的因子再对m求和m是d的倍数且不超过n $$S_k(n) \sum_{d|n} \mu(d) \sum_{\substack{m1 \ d|m}}^{n} m^k \sum_{d|n} \mu(d) \sum_{t1}^{\lfloor n/d \rfloor} (dt)^k \sum_{d|n} \mu(d) d^k \sum_{t1}^{N} t^k$$ 其中 $N \lfloor n/d \rfloor$。看我们成功地把一个带互质条件的求和转化成了一个关于n的因子d的求和以及一个标准的自然数幂和 $\sum_{t1}^{N} t^k$。这就是莫比乌斯反演在此类问题中的经典应用它像一把筛子把不符合互质条件的项通过带正负号的 $\mu(d)$ 给“筛”掉了。注意这里我们使用的是莫比乌斯反演的“直接形式”或者说“容斥原理形式”而不是教科书上常见的“倍数形式”或“因数形式”公式。在实际解题中根据求和指标的不同灵活运用 $\sum_{d|n}\mu(d)[n1]$ 这个核心等式进行变换往往比生搬硬套公式更直接有效。2.2 伯努利数给自然数幂和一个“公式”现在问题转化成了求 $\sum_{t1}^{N} t^k$即前N个自然数的k次幂之和。我们知道这个和可以写成一个关于N的(k1)次多项式这就是著名的Faulhaber公式 $$\sum_{t1}^{N} t^k \frac{1}{k1} \sum_{j0}^{k} \binom{k1}{j} B_j N^{k1-j}$$ 其中 $B_j$ 就是伯努利数这里采用 $B_1 1/2$ 的约定也有些约定 $B_1 -1/2$需特别注意。伯努利数可以通过生成函数来定义$\frac{x}{e^x - 1} \sum_{n0}^{\infty} B_n \frac{x^n}{n!}$。前几项是$B_01, B_1\frac{1}{2}, B_2\frac{1}{6}, B_30, B_4-\frac{1}{30}, B_50, B_6\frac{1}{42}, \ldots$。一个重要的性质是所有奇数项除了$B_1$都是0。为什么这个公式重要因为它把求和变成了多项式计算。对于固定的k我们可以预先算出多项式 $\sum_{t1}^{N} t^k a_{k1} N^{k1} a_k N^k \cdots a_1 N a_0$ 的各个系数 $a_i$这些系数由伯努利数组合而成。那么之前式子里的 $\sum_{t1}^{N} t^k$ 就可以用这个多项式 $P_k(N)$ 来代替。于是我们的 $S_k(n)$ 表达式进一步变为 $$S_k(n) \sum_{d|n} \mu(d) d^k P_k(\lfloor n/d \rfloor)$$ 由于 $n/d$ 是整数所以 $\lfloor n/d \rfloor n/d$。最终 $$S_k(n) \sum_{d|n} \mu(d) d^k P_k\left(\frac{n}{d}\right)$$到这里我们已经完成了理论上的化简。$S_k(n)$ 被表达为对n的所有因子d求和每一项是 $\mu(d) d^k$ 乘以一个关于 $(n/d)$ 的多项式 $P_k$。如果n的质因数分解已知我们可以枚举其所有因子通常不会太多因为因子个数是各质因子指数加1的乘积然后快速计算每一项的值。这就是整个算法的理论核心。3. 从理论到实现算法步骤与细节剖析理论很美但要把代码写对、写快中间还有不少细节需要抠。我们一步步来。3.1 伯努利数与幂和多项式的预处理首先我们需要得到 $P_k(N) \sum_{t1}^{N} t^k$ 这个多项式的系数。根据Faulhaber公式系数为 $$a_{k1-j} \frac{1}{k1} \binom{k1}{j} B_j, \quad j0,1,\ldots,k$$ 其中 $a_{k1-j}$ 是 $N^{k1-j}$ 项的系数。因此第一步是计算前 $k1$ 个伯努利数 $B_0, B_1, ..., B_k$。伯努利数可以通过递归公式计算 $$B_m 1 - \sum_{j0}^{m-1} \binom{m}{j} \frac{B_j}{m - j 1}$$ 或者用生成函数方法。在算法竞赛中k通常不会太大比如 $k \le 1000$所以直接用递归或动态规划计算是可行的。注意这里计算的是有理数分数我们需要在模意义下进行题目通常要求对一个大质数取模比如 $10^97$。这意味着我们需要处理模意义下的分数即计算分母的模逆元。假设模数为 $MOD$计算 $B_m \mod MOD$ 的步骤初始化 $B_0 1$。对于 $m$ 从 1 到 $K$计算 $sum 0$对于 $j$ 从 0 到 $m-1$$C \binom{m}{j} \mod MOD$ 组合数可以预处理阶乘和逆元$inv (m - j 1)$ 在模 $MOD$ 下的逆元$sum (sum C * B_j % MOD * inv % MOD) % MOD$$B_m (1 - sum MOD) % MOD$ 注意保证结果非负得到伯努利数后就可以计算多项式 $P_k$ 的系数了。我们通常用一个数组coef[0..k1]来存储其中coef[i]对应 $N^i$ 的系数。注意根据公式$P_k(N)$ 是 $k1$ 次多项式所以coef[k1]是最高次项系数。 对于 $j$ 从 0 到 $k$$i k1 - j$ $N^i$ 的指数$C \binom{k1}{j} \mod MOD$$inv (k1)$ 在模 $MOD$ 下的逆元coef[i] C * B_j \% MOD * inv \% MOD实操心得这里有一个极其容易出错的地方——伯努利数的约定。一定要确认你使用的公式和计算的 $B_1$ 值是否匹配。如果你用 $B_11/2$ 的约定那么上面的Faulhaber公式就是对的。如果你用的库或参考资料是 $B_1-1/2$那么公式可能需要调整。最稳妥的方法是用小的k值如k1,2,3手动验算一下$\sum_{t1}^N t^1 N(N1)/2$$\sum_{t1}^N t^2 N(N1)(2N1)/6$看看你计算出的多项式系数是否匹配。这个验证步骤能帮你提前发现大部分伯努利数相关的错误。3.2 利用质因数分解高效枚举因子现在我们有 $S_k(n) \sum_{d|n} \mu(d) d^k P_k(n/d)$。n以质因数分解形式给出$n \prod_{i1}^{w} p_i^{\alpha_i}$。直接枚举n的所有因子如果n很大质因子指数高因子个数可能爆炸。但注意$\mu(d)$ 有一个关键性质如果d有平方因子则 $\mu(d)0$。这意味着在求和中我们只需要枚举那些每个质因子至多取1次的因子即从每个质因子 $p_i$ 中选择“不取”指数0或“取一次”指数1。这样的因子共有 $2^w$ 个w是不同质因子的个数。题目中w通常不会太大比如几十以内所以 $2^w$ 的枚举是完全可行的。具体枚举可以通过DFS深度优先搜索或二进制位掩码实现。对于每个枚举出的因子d计算 $d \prod_{j \in S} p_j$其中S是选中的质因子下标集合。计算 $\mu(d) (-1)^{|S|}$。计算 $d^k \mod MOD$。由于k可能很大需要用快速幂计算。计算 $n/d$。因为n是质因子幂的乘积d是其中一些质因子的乘积所以 $n/d \prod_{i1}^{w} p_i^{\alpha_i - [i \in S]}$。这里 $[i \in S]$ 是指示函数如果质因子 $p_i$ 被选中即在d中则指数减1否则不变。计算多项式 $P_k(n/d)$。即把 $x n/d$ 代入多项式 $\sum_{i0}^{k1} coef[i] * x^i$ 中求值。这里 $x n/d$ 也是一个可能很大的数需要边计算边取模。计算 $x^i$ 同样可以用快速幂但更高效的方法是使用秦九韶算法Horners method $$P_k(x) coef[0] x*(coef[1] x*(coef[2] ... x*coef[k1])...))$$ 这样只需要O(k)次乘法和加法避免了重复计算幂次。将 $\mu(d) * (d^k) % MOD * (P_k(n/d)) % MOD$ 累加到答案中注意处理负号加MOD再取模。3.3 模运算下的完整计算流程假设模数 $MOD$ 为质数如 $10^97$。输入幂次k质因子个数w以及w对 $(p_i, \alpha_i)$。预处理计算阶乘fact[0..K1]和阶乘逆元invfact[0..K1]用于快速计算组合数 $\binom{n}{m}$。计算伯努利数 $B_0$ 到 $B_k$模 $MOD$。利用伯努利数计算多项式 $P_k$ 的系数数组coef[0..k1]。枚举因子并求和初始化答案ans 0。使用DFS或位运算枚举所有 $2^w$ 种质因子选择方案对应因子d。对于每种方案选中集合S计算mu 1如果 |S| 为奇数则mu MOD-1即-1模MOD。计算d_val 1。遍历所有质因子如果 $p_i$ 在S中则d_val d_val * p_i % MOD。计算d_pow_k pow_mod(d_val, k, MOD)快速幂。计算nd_val 1。遍历所有质因子计算 $p_i$ 的指数若 $p_i$ 在S中指数为 $\alpha_i - 1$否则为 $\alpha_i$。则nd_val nd_val * pow_mod(p_i, exp, MOD) % MOD。这里nd_val就是 $n/d$。用秦九韶算法计算poly_val evaluate_poly(coef, nd_val, k1)即 $P_k(n/d)$。累加ans (ans mu * d_pow_k % MOD * poly_val % MOD) % MOD。输出ans即为 $S_k(n) \mod MOD$。避坑指南模运算下处理负号。mu可能是1或-1。在模 $MOD$ 下-1等价于 $MOD-1$。所以当mu为-1时我们实际存储为MOD-1。在累加时mu * d_pow_k % MOD这部分如果mu是MOD-1相乘后可能会得到一个很大的数接近 $MOD^2$所以一定要先取模一次(mu * d_pow_k) % MOD再与poly_val相乘。同样最终累加后也要取模。使用long long类型或语言的64位整数来避免中间结果溢出。4. 深度思考为什么是多项式算法的本质与边界我们已经有了完整的算法但理解其本质能让我们更灵活地应用它。为什么最终 $S_k(n)$ 能表示成关于n的各质因子幂次的多项式这背后是数论函数的一个重要性质积性。4.1 积性函数的视角回顾 $S_k(n) \sum_{d|n} \mu(d) d^k P_k(n/d)$。让我们把它稍微改写一下。注意到 $P_k(x)$ 是一个关于x的多项式设 $P_k(x) \sum_{i0}^{k1} c_i x^i$。那么 $$S_k(n) \sum_{d|n} \mu(d) d^k \sum_{i0}^{k1} c_i \left(\frac{n}{d}\right)^i \sum_{i0}^{k1} c_i \sum_{d|n} \mu(d) d^k \left(\frac{n}{d}\right)^i \sum_{i0}^{k1} c_i n^i \sum_{d|n} \mu(d) d^{k-i}$$令 $g_i(n) \sum_{d|n} \mu(d) d^{k-i}$。这是一个狄利克雷卷积的形式$g_i \mu \cdot \text{id}_{k-i}$其中 $\text{id}_a(n)n^a$。我们知道莫比乌斯函数 $\mu$ 和幂函数 $\text{id}_a$ 都是积性函数而积性函数的狄利克雷卷积依然是积性函数。所以 $g_i(n)$ 是积性函数。对于积性函数如果 $n \prod p^{\alpha}$那么 $g_i(n) \prod g_i(p^{\alpha})$。计算 $g_i(p^{\alpha})$ 非常简单因为对于质数幂 $p^{\alpha}$其因子只有 $1, p, p^2, ..., p^{\alpha}$。但 $\mu(d)$ 在d有平方因子时为0所以只有当 $d1$ 或 $dp$ 时假设 $\alpha \ge 1$$\mu(d)$ 才非零。$d1$: $\mu(1)1$, $d^{k-i}1^{k-i}1$贡献为1。$dp$: $\mu(p)-1$, $d^{k-i}p^{k-i}$贡献为 $-p^{k-i}$。$dp^t, t\ge2$: $\mu(p^t)0$贡献为0。因此$g_i(p^{\alpha}) 1 - p^{k-i}$。注意这个结果与 $\alpha$ 无关只要 $\alpha \ge 1$$g_i(p^{\alpha})$ 的值只取决于质因子p本身和指数差 $k-i$。于是$g_i(n) \prod_{p|n} (1 - p^{k-i})$。代回原式 $$S_k(n) \sum_{i0}^{k1} c_i n^i \prod_{p|n} (1 - p^{k-i})$$这个形式更清晰地揭示了 $S_k(n)$ 的结构它是关于 $n^i$ 的线性组合而每个系数 $c_i$ 又被一个与n的质因子相关的乘积 $\prod_{p|n} (1 - p^{k-i})$ 所调制。当n的质因数分解给出时计算这个表达式是非常快的。4.2 算法复杂度与优化空间我们有两种计算方式枚举因子法复杂度 $O(2^w \cdot (k \log MOD))$。其中 $2^w$ 是因子枚举k是多项式求值秦九韶算法$\log MOD$ 是快速幂的代价。当w较小时比如w 20这种方法非常高效。积性函数直接计算法根据最终公式 $S_k(n) \sum_{i0}^{k1} c_i n^i \prod_{p|n} (1 - p^{k-i})$。我们需要计算所有 $c_i$ (k2个)。对于每个i (0到k1)计算 $n^i \mod MOD$可以用递推 $n^{i} n^{i-1} * n \mod MOD$。对于每个i计算 $\prod_{p|n} (1 - p^{k-i})$这需要对每个质因子p计算 $p^{k-i} \mod MOD$然后连乘。计算 $p^{k-i}$ 需要快速幂复杂度 $O(w \cdot \log MOD)$。总复杂度约为 $O(k \cdot w \cdot \log MOD)$。两种方法各有优劣。枚举因子法在w很小时是首选代码直观。而直接计算法在w相对较大但k不太大时可能更有优势因为它避免了 $2^w$ 的枚举。在实际题目中通常会对w或k的范围有所限制需要根据具体情况选择。经验之谈在竞赛中如果题目明确给出了n的质因数分解形式并且w不同质因子数不大那么枚举因子法通常是更稳妥和易于实现的选择。它的逻辑清晰每一步都可以单独验证。而直接使用积性函数公式虽然数学上更优美但实现时需要注意 $k-i$ 可能为负数的情况当i k时。在模运算中$p^{k-i}$ 需要计算 $p^{k-i} \mod MOD$如果 $k-i$ 是负数则意味着计算模逆元 $p^{i-k} \mod MOD$这要求p与MOD互质通常MOD是大质数p是n的质因子只要p不等于MOD就成立。处理指数为负的情况会稍微增加代码复杂度。4.3 边界情况与测试实现完成后必须用一些边界情况测试。k0$S_0(n)$ 是与n互质的数的个数即欧拉函数 $\varphi(n)$。我们的算法应该能正确计算出 $\varphi(n) n \prod_{p|n}(1 - 1/p)$。验证时注意多项式 $P_0(N) N$伯努利数 $B_01, B_11/2$公式给出 $P_0(N) \frac{1}{1}\binom{1}{0}B_0 N^1 \frac{1}{1}\binom{1}{1}B_1 N^0 N 1/2$等等这里出问题了。实际上对于k0自然数0次幂和 $\sum_{t1}^N t^0 \sum_{t1}^N 1 N$。但Faulhaber公式在k0时$P_0(N) \frac{1}{1}\sum_{j0}^0 \binom{1}{j}B_j N^{1-j} \binom{1}{0}B_0 N^1 N$。因为公式中j从0到k当k0时j只能取0。所以 $P_0(N)N$。代入我们的推导$S_0(n) \sum_{d|n} \mu(d) * d^0 * (n/d) n \sum_{d|n} \mu(d)/d$。利用狄利克雷卷积或者莫比乌斯反演的性质可以证明 $\sum_{d|n} \mu(d)/d \varphi(n)/n$。所以结果是正确的。但这也提醒我们在处理k很小的时候要小心公式的边界。n1只有一个数1且与1互质所以 $S_k(1)1^k1$。我们的算法中n1时质因子列表为空枚举因子只有d1$\mu(1)1, d^k1, P_k(n/d)P_k(1)\sum_{t1}^1 t^k 1$结果为1正确。n为质数p$S_k(p) \sum_{m1}^{p-1} m^k$因为1到p中只有p与p不互质。另一方面根据我们的公式因子只有d1和dp。$\mu(1)1, 1^k1, P_k(p/1)P_k(p)$$\mu(p)-1, p^k, P_k(p/p)P_k(1)1$。所以 $S_k(p) 11P_k(p) (-1)p^k1 P_k(p) - p^k$。而 $P_k(p) \sum_{t1}^p t^k$所以 $P_k(p) - p^k \sum_{t1}^{p-1} t^k$吻合。通过这些测试可以增强对代码正确性的信心。5. 代码实现的关键片段与解析光说不练假把式这里给出一些核心代码片段的思路和实现细节以C为例模数MOD1e97。typedef long long ll; const int MOD 1e9 7; // 快速幂 ll pow_mod(ll a, ll b) { ll res 1; while (b) { if (b 1) res res * a % MOD; a a * a % MOD; b 1; } return res; } // 预处理阶乘和逆元用于组合数 vectorll fact, inv_fact; void init_fact(int n) { fact.resize(n1); inv_fact.resize(n1); fact[0] 1; for (int i1; in; i) fact[i] fact[i-1] * i % MOD; inv_fact[n] pow_mod(fact[n], MOD-2); for (int in-1; i0; i--) inv_fact[i] inv_fact[i1] * (i1) % MOD; } ll C(int n, int m) { if (m0 || mn) return 0; return fact[n] * inv_fact[m] % MOD * inv_fact[n-m] % MOD; } // 计算伯努利数 B[0..k] vectorll get_bernoulli(int k) { vectorll B(k1); B[0] 1; for (int m1; mk; m) { ll sum 0; for (int j0; jm; j) { ll comb C(m1, j); // 注意这里公式是 C(m1, j)但更常见的是 C(m, j)需要根据递推公式调整 // 实际上更标准的递推是B_m 1 - sum_{j0}^{m-1} C(m, j) * B_j / (m - j 1) // 我们使用这个公式 comb C(m, j); ll inv pow_mod(m - j 1, MOD-2); sum (sum comb * B[j] % MOD * inv) % MOD; } B[m] (1 - sum MOD) % MOD; } return B; } // 计算多项式 P_k(x) 的系数 coef[0..k1]coef[i] 对应 x^i vectorll get_poly_coef(int k) { vectorll B get_bernoulli(k); vectorll coef(k2, 0); // 最高次为 k1 ll inv_k1 pow_mod(k1, MOD-2); for (int j0; jk; j) { ll comb C(k1, j); int exp k1 - j; // x^{exp} 的系数 coef[exp] comb * B[j] % MOD * inv_k1 % MOD; } // 注意j从0到k对应exp从k1到1。coef[0]项常数项呢 // 当jk时exp1。当j0时expk1。所以常数项coef[0]没有被赋值应为0。 // 验证k1时P_1(N)N(N1)/2 0.5*N^2 0.5*N常数项为0。 return coef; } // 秦九韶算法计算多项式值 ll eval_poly(const vectorll coef, ll x, int max_deg) { ll res 0; for (int imax_deg; i0; i--) { res (res * x % MOD coef[i]) % MOD; } return res; } // 主求解函数枚举因子法 ll solve(int k, const vectorpairll, int factors) { // factors: (p_i, alpha_i) vectorll coef get_poly_coef(k); int w factors.size(); ll n 1; for (auto [p, a] : factors) { // 实际n可能很大超出long long我们只需要n mod MOD但后面需要n/d所以这里先计算n模MOD // 更安全的做法是在计算过程中需要n/d时直接用质因子指数相减来构造。 } ll ans 0; // 使用DFS枚举所有因子d每个质因子选或不选 functionvoid(int, ll, ll, int) dfs [](int idx, ll d, ll nd, int mu) { if (idx w) { // 枚举到一个因子d ll d_pow_k pow_mod(d, k); ll poly_val eval_poly(coef, nd, k1); // nd n/d ans (ans mu * d_pow_k % MOD * poly_val % MOD MOD) % MOD; return; } ll p factors[idx].first; int a factors[idx].second; // 不选当前质因子 ll nd_mul pow_mod(p, a); // n/d 中这个质因子的贡献是 p^a dfs(idx1, d, nd * nd_mul % MOD, mu); // 选当前质因子指数最多为1因为mu(d)在平方因子时为0 // d 乘以 p // n/d 中这个质因子的贡献是 p^{a-1} ll nd_mul_sel (a 1) ? pow_mod(p, a-1) : 1; // 如果a1则选完后指数为0贡献为1 dfs(idx1, d * p % MOD, nd * nd_mul_sel % MOD, (mu 1) ? MOD-1 : 1); // mu取反 }; dfs(0, 1, 1, 1); // 初始d1, n/d1, mu1 return ans; }代码细节提醒伯努利数递推代码中使用的递推公式是 $B_m 1 - \sum_{j0}^{m-1} \binom{m}{j} \frac{B_j}{m-j1}$。注意分母是m-j1不要写成m-j。多项式系数coef数组下标i代表 $x^i$ 的系数。在get_poly_coef中我们只填充了exp从1到k1的系数coef[0]默认为0这与理论一致自然数幂和多项式没有常数项。枚举因子DFS函数中我们同时维护了d当前已选质因子的乘积和nd当前已计算的 n/d 的值。这样避免了最后再去计算n/d提高了效率。注意nd的计算如果不选该质因子则n/d中包含 $p^a$如果选则包含 $p^{a-1}$如果a1则 $p^01$。模运算所有乘法后都要立即取模。mu用1和MOD-1表示正负。在累加时(mu * d_pow_k % MOD * poly_val % MOD MOD) % MOD可以确保结果非负但更安全的写法是分步取模。n的计算在枚举过程中我们并不需要显式地计算出完整的n它可能巨大只需要n/d。我们通过质因子指数直接构造nd这是处理大数的关键。6. 总结与延伸思考回顾整个解题过程我们从一道要求互质数幂和的题目出发引入了莫比乌斯反演来化解互质条件又借助伯努利数得到的自然数幂和公式将问题转化为多项式求值。最终通过枚举n的平方自由因子即 $\mu(d) \ne 0$ 的因子我们能在 $O(2^w \cdot k)$ 的时间复杂度内解决问题其中w是n的不同质因子数。这道题的价值不仅在于其解法更在于它展示了一种典型的“组合数学数论”的解题范式通过莫比乌斯反演将数论条件转化为求和再通过已知的求和公式如伯努利数公式、斯特林数公式等将求和转化为多项式或其他封闭形式最后利用输入数据的结构如质因数分解进行高效计算。在这个过程中有几个点值得反复品味数学工具的选择为什么用莫比乌斯反演因为[gcd(m,n)1]这个条件天生适合用 $\sum_{d|gcd(m,n)} \mu(d)$ 来替换。为什么用伯努利数因为我们需要一个关于N的幂和 $\sum t^k$ 的精确多项式表达式。计算与优化理论公式到可执行代码之间有大量细节需要打磨。模运算下的分数处理、伯努利数的递推、多项式的求值秦九韶算法、利用质因数分解枚举因子并同时计算n/d这些都是将数学转化为高效算法的关键步骤。验证与调试对于这类推导复杂的题目用小的、可手算的案例如k0,1, n为小质数进行验证至关重要。它不仅能帮助发现代码错误更能加深对公式本身的理解。最后这个解法还可以进一步扩展。例如如果题目不是求互质数的k次幂和而是求gcd为某个定值g的数的k次幂和我们依然可以通过类似的莫比乌斯反演技巧将其转化为本题的形式。伯努利数也可以被斯特林数替代它们之间也有联系都能用于表达幂和。这些延伸就留给有兴趣的朋友去探索了。数论的世界就是这样一个精巧的问题往往能串联起多个美丽的数学概念而解决它的过程本身就是一次深刻的思维训练。
莫比乌斯反演与伯努利数：高效求解互质数幂和问题
1. 从一道“省队互测”题说起当数论遇上组合数学最近在整理一些经典的数论题目时又翻到了这道P6271。题目名字很有意思“一个人的数论”听起来有点孤独但内容却一点也不孤单它把数论里几个重量级的概念——莫比乌斯反演和伯努利数——给串了起来。很多朋友第一次看到这个组合可能会有点懵莫比乌斯反演不是用来处理数论函数求和、化简gcd、lcm相关式子的吗伯努利数不是出现在幂和公式或者泰勒展开里的吗它们俩怎么会扯上关系这道题的精妙之处就在于此。它表面上可能只是让你求一个带幂次的因子和函数但内核却是一个经典的“数论函数与多项式”的桥梁问题。简单来说题目通常会给你一个关于正整数n的函数比如所有与n互质的数的k次幂之和记作 $S_k(n) \sum_{\substack{1 \le m \le n \ \gcd(m, n)1}} m^k$。我们的任务就是求出这个 $S_k(n)$ 的一个简洁表达式特别是当n以质因数分解形式 $n \prod_{i1}^w p_i^{\alpha_i}$ 给出时如何快速计算。直接枚举1到n的所有数并判断互质复杂度是O(n)当n很大时比如是多个大质数的高次幂乘积完全不可行。这时莫比乌斯反演就登场了它能帮我们把“互质”这个条件巧妙地转化掉。而伯努利数则为我们提供了计算自然数幂和公式的系数让我们能把一个求和式写成一个关于n的多项式。两者一结合就能把原问题转化为一个关于n的各质因子幂次的多项式求值问题复杂度瞬间降低到只与质因子个数w和幂次k有关实现了从“一个人的暴力枚举”到“优雅公式推导”的飞跃。这篇文章我就想结合这道题把莫比乌斯反演和伯努利数在这个场景下的应用逻辑彻底拆解清楚。我们不止看公式更要看每一步“为什么这么做”以及在实际推导和计算中会遇到哪些坑怎么绕过去。无论你是正在备战算法竞赛还是单纯对数学的美感有兴趣相信这个“数论与组合数学的联姻”故事都会让你有所收获。2. 核心武器拆解莫比乌斯反演与伯努利数在进入具体解题之前我们必须把两件核心武器的原理和使用方法摸透。很多资料只给结论但不知道背后的动机用起来就会束手束脚。2.1 莫比乌斯反演化“互质”为“整除”莫比乌斯函数 $\mu(n)$ 大家应该不陌生当n有平方因子时为0否则若n有偶数个质因子则为1奇数个则为-1。它的核心价值在于这个性质$\sum_{d|n} \mu(d) [n1]$其中 $[P]$ 是艾弗森括号当命题P为真时值为1否则为0。这个看似简单的式子正是解决“互质”条件的关键。考虑我们想求所有满足 $\gcd(m, n)1$ 的m的和。注意到 $[\gcd(m, n)1]$ 这个条件利用上述性质我们可以把它改写 $$[\gcd(m, n)1] \sum_{d|\gcd(m, n)} \mu(d)$$ 而 $d|\gcd(m, n)$ 等价于 $d|m$ 且 $d|n$。于是我们的目标求和式就变成了 $$S_k(n) \sum_{m1}^{n} m^k \cdot [\gcd(m, n)1] \sum_{m1}^{n} m^k \sum_{\substack{d|m \ d|n}} \mu(d)$$接下来是关键一步交换求和顺序。先对d求和d是n的因子再对m求和m是d的倍数且不超过n $$S_k(n) \sum_{d|n} \mu(d) \sum_{\substack{m1 \ d|m}}^{n} m^k \sum_{d|n} \mu(d) \sum_{t1}^{\lfloor n/d \rfloor} (dt)^k \sum_{d|n} \mu(d) d^k \sum_{t1}^{N} t^k$$ 其中 $N \lfloor n/d \rfloor$。看我们成功地把一个带互质条件的求和转化成了一个关于n的因子d的求和以及一个标准的自然数幂和 $\sum_{t1}^{N} t^k$。这就是莫比乌斯反演在此类问题中的经典应用它像一把筛子把不符合互质条件的项通过带正负号的 $\mu(d)$ 给“筛”掉了。注意这里我们使用的是莫比乌斯反演的“直接形式”或者说“容斥原理形式”而不是教科书上常见的“倍数形式”或“因数形式”公式。在实际解题中根据求和指标的不同灵活运用 $\sum_{d|n}\mu(d)[n1]$ 这个核心等式进行变换往往比生搬硬套公式更直接有效。2.2 伯努利数给自然数幂和一个“公式”现在问题转化成了求 $\sum_{t1}^{N} t^k$即前N个自然数的k次幂之和。我们知道这个和可以写成一个关于N的(k1)次多项式这就是著名的Faulhaber公式 $$\sum_{t1}^{N} t^k \frac{1}{k1} \sum_{j0}^{k} \binom{k1}{j} B_j N^{k1-j}$$ 其中 $B_j$ 就是伯努利数这里采用 $B_1 1/2$ 的约定也有些约定 $B_1 -1/2$需特别注意。伯努利数可以通过生成函数来定义$\frac{x}{e^x - 1} \sum_{n0}^{\infty} B_n \frac{x^n}{n!}$。前几项是$B_01, B_1\frac{1}{2}, B_2\frac{1}{6}, B_30, B_4-\frac{1}{30}, B_50, B_6\frac{1}{42}, \ldots$。一个重要的性质是所有奇数项除了$B_1$都是0。为什么这个公式重要因为它把求和变成了多项式计算。对于固定的k我们可以预先算出多项式 $\sum_{t1}^{N} t^k a_{k1} N^{k1} a_k N^k \cdots a_1 N a_0$ 的各个系数 $a_i$这些系数由伯努利数组合而成。那么之前式子里的 $\sum_{t1}^{N} t^k$ 就可以用这个多项式 $P_k(N)$ 来代替。于是我们的 $S_k(n)$ 表达式进一步变为 $$S_k(n) \sum_{d|n} \mu(d) d^k P_k(\lfloor n/d \rfloor)$$ 由于 $n/d$ 是整数所以 $\lfloor n/d \rfloor n/d$。最终 $$S_k(n) \sum_{d|n} \mu(d) d^k P_k\left(\frac{n}{d}\right)$$到这里我们已经完成了理论上的化简。$S_k(n)$ 被表达为对n的所有因子d求和每一项是 $\mu(d) d^k$ 乘以一个关于 $(n/d)$ 的多项式 $P_k$。如果n的质因数分解已知我们可以枚举其所有因子通常不会太多因为因子个数是各质因子指数加1的乘积然后快速计算每一项的值。这就是整个算法的理论核心。3. 从理论到实现算法步骤与细节剖析理论很美但要把代码写对、写快中间还有不少细节需要抠。我们一步步来。3.1 伯努利数与幂和多项式的预处理首先我们需要得到 $P_k(N) \sum_{t1}^{N} t^k$ 这个多项式的系数。根据Faulhaber公式系数为 $$a_{k1-j} \frac{1}{k1} \binom{k1}{j} B_j, \quad j0,1,\ldots,k$$ 其中 $a_{k1-j}$ 是 $N^{k1-j}$ 项的系数。因此第一步是计算前 $k1$ 个伯努利数 $B_0, B_1, ..., B_k$。伯努利数可以通过递归公式计算 $$B_m 1 - \sum_{j0}^{m-1} \binom{m}{j} \frac{B_j}{m - j 1}$$ 或者用生成函数方法。在算法竞赛中k通常不会太大比如 $k \le 1000$所以直接用递归或动态规划计算是可行的。注意这里计算的是有理数分数我们需要在模意义下进行题目通常要求对一个大质数取模比如 $10^97$。这意味着我们需要处理模意义下的分数即计算分母的模逆元。假设模数为 $MOD$计算 $B_m \mod MOD$ 的步骤初始化 $B_0 1$。对于 $m$ 从 1 到 $K$计算 $sum 0$对于 $j$ 从 0 到 $m-1$$C \binom{m}{j} \mod MOD$ 组合数可以预处理阶乘和逆元$inv (m - j 1)$ 在模 $MOD$ 下的逆元$sum (sum C * B_j % MOD * inv % MOD) % MOD$$B_m (1 - sum MOD) % MOD$ 注意保证结果非负得到伯努利数后就可以计算多项式 $P_k$ 的系数了。我们通常用一个数组coef[0..k1]来存储其中coef[i]对应 $N^i$ 的系数。注意根据公式$P_k(N)$ 是 $k1$ 次多项式所以coef[k1]是最高次项系数。 对于 $j$ 从 0 到 $k$$i k1 - j$ $N^i$ 的指数$C \binom{k1}{j} \mod MOD$$inv (k1)$ 在模 $MOD$ 下的逆元coef[i] C * B_j \% MOD * inv \% MOD实操心得这里有一个极其容易出错的地方——伯努利数的约定。一定要确认你使用的公式和计算的 $B_1$ 值是否匹配。如果你用 $B_11/2$ 的约定那么上面的Faulhaber公式就是对的。如果你用的库或参考资料是 $B_1-1/2$那么公式可能需要调整。最稳妥的方法是用小的k值如k1,2,3手动验算一下$\sum_{t1}^N t^1 N(N1)/2$$\sum_{t1}^N t^2 N(N1)(2N1)/6$看看你计算出的多项式系数是否匹配。这个验证步骤能帮你提前发现大部分伯努利数相关的错误。3.2 利用质因数分解高效枚举因子现在我们有 $S_k(n) \sum_{d|n} \mu(d) d^k P_k(n/d)$。n以质因数分解形式给出$n \prod_{i1}^{w} p_i^{\alpha_i}$。直接枚举n的所有因子如果n很大质因子指数高因子个数可能爆炸。但注意$\mu(d)$ 有一个关键性质如果d有平方因子则 $\mu(d)0$。这意味着在求和中我们只需要枚举那些每个质因子至多取1次的因子即从每个质因子 $p_i$ 中选择“不取”指数0或“取一次”指数1。这样的因子共有 $2^w$ 个w是不同质因子的个数。题目中w通常不会太大比如几十以内所以 $2^w$ 的枚举是完全可行的。具体枚举可以通过DFS深度优先搜索或二进制位掩码实现。对于每个枚举出的因子d计算 $d \prod_{j \in S} p_j$其中S是选中的质因子下标集合。计算 $\mu(d) (-1)^{|S|}$。计算 $d^k \mod MOD$。由于k可能很大需要用快速幂计算。计算 $n/d$。因为n是质因子幂的乘积d是其中一些质因子的乘积所以 $n/d \prod_{i1}^{w} p_i^{\alpha_i - [i \in S]}$。这里 $[i \in S]$ 是指示函数如果质因子 $p_i$ 被选中即在d中则指数减1否则不变。计算多项式 $P_k(n/d)$。即把 $x n/d$ 代入多项式 $\sum_{i0}^{k1} coef[i] * x^i$ 中求值。这里 $x n/d$ 也是一个可能很大的数需要边计算边取模。计算 $x^i$ 同样可以用快速幂但更高效的方法是使用秦九韶算法Horners method $$P_k(x) coef[0] x*(coef[1] x*(coef[2] ... x*coef[k1])...))$$ 这样只需要O(k)次乘法和加法避免了重复计算幂次。将 $\mu(d) * (d^k) % MOD * (P_k(n/d)) % MOD$ 累加到答案中注意处理负号加MOD再取模。3.3 模运算下的完整计算流程假设模数 $MOD$ 为质数如 $10^97$。输入幂次k质因子个数w以及w对 $(p_i, \alpha_i)$。预处理计算阶乘fact[0..K1]和阶乘逆元invfact[0..K1]用于快速计算组合数 $\binom{n}{m}$。计算伯努利数 $B_0$ 到 $B_k$模 $MOD$。利用伯努利数计算多项式 $P_k$ 的系数数组coef[0..k1]。枚举因子并求和初始化答案ans 0。使用DFS或位运算枚举所有 $2^w$ 种质因子选择方案对应因子d。对于每种方案选中集合S计算mu 1如果 |S| 为奇数则mu MOD-1即-1模MOD。计算d_val 1。遍历所有质因子如果 $p_i$ 在S中则d_val d_val * p_i % MOD。计算d_pow_k pow_mod(d_val, k, MOD)快速幂。计算nd_val 1。遍历所有质因子计算 $p_i$ 的指数若 $p_i$ 在S中指数为 $\alpha_i - 1$否则为 $\alpha_i$。则nd_val nd_val * pow_mod(p_i, exp, MOD) % MOD。这里nd_val就是 $n/d$。用秦九韶算法计算poly_val evaluate_poly(coef, nd_val, k1)即 $P_k(n/d)$。累加ans (ans mu * d_pow_k % MOD * poly_val % MOD) % MOD。输出ans即为 $S_k(n) \mod MOD$。避坑指南模运算下处理负号。mu可能是1或-1。在模 $MOD$ 下-1等价于 $MOD-1$。所以当mu为-1时我们实际存储为MOD-1。在累加时mu * d_pow_k % MOD这部分如果mu是MOD-1相乘后可能会得到一个很大的数接近 $MOD^2$所以一定要先取模一次(mu * d_pow_k) % MOD再与poly_val相乘。同样最终累加后也要取模。使用long long类型或语言的64位整数来避免中间结果溢出。4. 深度思考为什么是多项式算法的本质与边界我们已经有了完整的算法但理解其本质能让我们更灵活地应用它。为什么最终 $S_k(n)$ 能表示成关于n的各质因子幂次的多项式这背后是数论函数的一个重要性质积性。4.1 积性函数的视角回顾 $S_k(n) \sum_{d|n} \mu(d) d^k P_k(n/d)$。让我们把它稍微改写一下。注意到 $P_k(x)$ 是一个关于x的多项式设 $P_k(x) \sum_{i0}^{k1} c_i x^i$。那么 $$S_k(n) \sum_{d|n} \mu(d) d^k \sum_{i0}^{k1} c_i \left(\frac{n}{d}\right)^i \sum_{i0}^{k1} c_i \sum_{d|n} \mu(d) d^k \left(\frac{n}{d}\right)^i \sum_{i0}^{k1} c_i n^i \sum_{d|n} \mu(d) d^{k-i}$$令 $g_i(n) \sum_{d|n} \mu(d) d^{k-i}$。这是一个狄利克雷卷积的形式$g_i \mu \cdot \text{id}_{k-i}$其中 $\text{id}_a(n)n^a$。我们知道莫比乌斯函数 $\mu$ 和幂函数 $\text{id}_a$ 都是积性函数而积性函数的狄利克雷卷积依然是积性函数。所以 $g_i(n)$ 是积性函数。对于积性函数如果 $n \prod p^{\alpha}$那么 $g_i(n) \prod g_i(p^{\alpha})$。计算 $g_i(p^{\alpha})$ 非常简单因为对于质数幂 $p^{\alpha}$其因子只有 $1, p, p^2, ..., p^{\alpha}$。但 $\mu(d)$ 在d有平方因子时为0所以只有当 $d1$ 或 $dp$ 时假设 $\alpha \ge 1$$\mu(d)$ 才非零。$d1$: $\mu(1)1$, $d^{k-i}1^{k-i}1$贡献为1。$dp$: $\mu(p)-1$, $d^{k-i}p^{k-i}$贡献为 $-p^{k-i}$。$dp^t, t\ge2$: $\mu(p^t)0$贡献为0。因此$g_i(p^{\alpha}) 1 - p^{k-i}$。注意这个结果与 $\alpha$ 无关只要 $\alpha \ge 1$$g_i(p^{\alpha})$ 的值只取决于质因子p本身和指数差 $k-i$。于是$g_i(n) \prod_{p|n} (1 - p^{k-i})$。代回原式 $$S_k(n) \sum_{i0}^{k1} c_i n^i \prod_{p|n} (1 - p^{k-i})$$这个形式更清晰地揭示了 $S_k(n)$ 的结构它是关于 $n^i$ 的线性组合而每个系数 $c_i$ 又被一个与n的质因子相关的乘积 $\prod_{p|n} (1 - p^{k-i})$ 所调制。当n的质因数分解给出时计算这个表达式是非常快的。4.2 算法复杂度与优化空间我们有两种计算方式枚举因子法复杂度 $O(2^w \cdot (k \log MOD))$。其中 $2^w$ 是因子枚举k是多项式求值秦九韶算法$\log MOD$ 是快速幂的代价。当w较小时比如w 20这种方法非常高效。积性函数直接计算法根据最终公式 $S_k(n) \sum_{i0}^{k1} c_i n^i \prod_{p|n} (1 - p^{k-i})$。我们需要计算所有 $c_i$ (k2个)。对于每个i (0到k1)计算 $n^i \mod MOD$可以用递推 $n^{i} n^{i-1} * n \mod MOD$。对于每个i计算 $\prod_{p|n} (1 - p^{k-i})$这需要对每个质因子p计算 $p^{k-i} \mod MOD$然后连乘。计算 $p^{k-i}$ 需要快速幂复杂度 $O(w \cdot \log MOD)$。总复杂度约为 $O(k \cdot w \cdot \log MOD)$。两种方法各有优劣。枚举因子法在w很小时是首选代码直观。而直接计算法在w相对较大但k不太大时可能更有优势因为它避免了 $2^w$ 的枚举。在实际题目中通常会对w或k的范围有所限制需要根据具体情况选择。经验之谈在竞赛中如果题目明确给出了n的质因数分解形式并且w不同质因子数不大那么枚举因子法通常是更稳妥和易于实现的选择。它的逻辑清晰每一步都可以单独验证。而直接使用积性函数公式虽然数学上更优美但实现时需要注意 $k-i$ 可能为负数的情况当i k时。在模运算中$p^{k-i}$ 需要计算 $p^{k-i} \mod MOD$如果 $k-i$ 是负数则意味着计算模逆元 $p^{i-k} \mod MOD$这要求p与MOD互质通常MOD是大质数p是n的质因子只要p不等于MOD就成立。处理指数为负的情况会稍微增加代码复杂度。4.3 边界情况与测试实现完成后必须用一些边界情况测试。k0$S_0(n)$ 是与n互质的数的个数即欧拉函数 $\varphi(n)$。我们的算法应该能正确计算出 $\varphi(n) n \prod_{p|n}(1 - 1/p)$。验证时注意多项式 $P_0(N) N$伯努利数 $B_01, B_11/2$公式给出 $P_0(N) \frac{1}{1}\binom{1}{0}B_0 N^1 \frac{1}{1}\binom{1}{1}B_1 N^0 N 1/2$等等这里出问题了。实际上对于k0自然数0次幂和 $\sum_{t1}^N t^0 \sum_{t1}^N 1 N$。但Faulhaber公式在k0时$P_0(N) \frac{1}{1}\sum_{j0}^0 \binom{1}{j}B_j N^{1-j} \binom{1}{0}B_0 N^1 N$。因为公式中j从0到k当k0时j只能取0。所以 $P_0(N)N$。代入我们的推导$S_0(n) \sum_{d|n} \mu(d) * d^0 * (n/d) n \sum_{d|n} \mu(d)/d$。利用狄利克雷卷积或者莫比乌斯反演的性质可以证明 $\sum_{d|n} \mu(d)/d \varphi(n)/n$。所以结果是正确的。但这也提醒我们在处理k很小的时候要小心公式的边界。n1只有一个数1且与1互质所以 $S_k(1)1^k1$。我们的算法中n1时质因子列表为空枚举因子只有d1$\mu(1)1, d^k1, P_k(n/d)P_k(1)\sum_{t1}^1 t^k 1$结果为1正确。n为质数p$S_k(p) \sum_{m1}^{p-1} m^k$因为1到p中只有p与p不互质。另一方面根据我们的公式因子只有d1和dp。$\mu(1)1, 1^k1, P_k(p/1)P_k(p)$$\mu(p)-1, p^k, P_k(p/p)P_k(1)1$。所以 $S_k(p) 11P_k(p) (-1)p^k1 P_k(p) - p^k$。而 $P_k(p) \sum_{t1}^p t^k$所以 $P_k(p) - p^k \sum_{t1}^{p-1} t^k$吻合。通过这些测试可以增强对代码正确性的信心。5. 代码实现的关键片段与解析光说不练假把式这里给出一些核心代码片段的思路和实现细节以C为例模数MOD1e97。typedef long long ll; const int MOD 1e9 7; // 快速幂 ll pow_mod(ll a, ll b) { ll res 1; while (b) { if (b 1) res res * a % MOD; a a * a % MOD; b 1; } return res; } // 预处理阶乘和逆元用于组合数 vectorll fact, inv_fact; void init_fact(int n) { fact.resize(n1); inv_fact.resize(n1); fact[0] 1; for (int i1; in; i) fact[i] fact[i-1] * i % MOD; inv_fact[n] pow_mod(fact[n], MOD-2); for (int in-1; i0; i--) inv_fact[i] inv_fact[i1] * (i1) % MOD; } ll C(int n, int m) { if (m0 || mn) return 0; return fact[n] * inv_fact[m] % MOD * inv_fact[n-m] % MOD; } // 计算伯努利数 B[0..k] vectorll get_bernoulli(int k) { vectorll B(k1); B[0] 1; for (int m1; mk; m) { ll sum 0; for (int j0; jm; j) { ll comb C(m1, j); // 注意这里公式是 C(m1, j)但更常见的是 C(m, j)需要根据递推公式调整 // 实际上更标准的递推是B_m 1 - sum_{j0}^{m-1} C(m, j) * B_j / (m - j 1) // 我们使用这个公式 comb C(m, j); ll inv pow_mod(m - j 1, MOD-2); sum (sum comb * B[j] % MOD * inv) % MOD; } B[m] (1 - sum MOD) % MOD; } return B; } // 计算多项式 P_k(x) 的系数 coef[0..k1]coef[i] 对应 x^i vectorll get_poly_coef(int k) { vectorll B get_bernoulli(k); vectorll coef(k2, 0); // 最高次为 k1 ll inv_k1 pow_mod(k1, MOD-2); for (int j0; jk; j) { ll comb C(k1, j); int exp k1 - j; // x^{exp} 的系数 coef[exp] comb * B[j] % MOD * inv_k1 % MOD; } // 注意j从0到k对应exp从k1到1。coef[0]项常数项呢 // 当jk时exp1。当j0时expk1。所以常数项coef[0]没有被赋值应为0。 // 验证k1时P_1(N)N(N1)/2 0.5*N^2 0.5*N常数项为0。 return coef; } // 秦九韶算法计算多项式值 ll eval_poly(const vectorll coef, ll x, int max_deg) { ll res 0; for (int imax_deg; i0; i--) { res (res * x % MOD coef[i]) % MOD; } return res; } // 主求解函数枚举因子法 ll solve(int k, const vectorpairll, int factors) { // factors: (p_i, alpha_i) vectorll coef get_poly_coef(k); int w factors.size(); ll n 1; for (auto [p, a] : factors) { // 实际n可能很大超出long long我们只需要n mod MOD但后面需要n/d所以这里先计算n模MOD // 更安全的做法是在计算过程中需要n/d时直接用质因子指数相减来构造。 } ll ans 0; // 使用DFS枚举所有因子d每个质因子选或不选 functionvoid(int, ll, ll, int) dfs [](int idx, ll d, ll nd, int mu) { if (idx w) { // 枚举到一个因子d ll d_pow_k pow_mod(d, k); ll poly_val eval_poly(coef, nd, k1); // nd n/d ans (ans mu * d_pow_k % MOD * poly_val % MOD MOD) % MOD; return; } ll p factors[idx].first; int a factors[idx].second; // 不选当前质因子 ll nd_mul pow_mod(p, a); // n/d 中这个质因子的贡献是 p^a dfs(idx1, d, nd * nd_mul % MOD, mu); // 选当前质因子指数最多为1因为mu(d)在平方因子时为0 // d 乘以 p // n/d 中这个质因子的贡献是 p^{a-1} ll nd_mul_sel (a 1) ? pow_mod(p, a-1) : 1; // 如果a1则选完后指数为0贡献为1 dfs(idx1, d * p % MOD, nd * nd_mul_sel % MOD, (mu 1) ? MOD-1 : 1); // mu取反 }; dfs(0, 1, 1, 1); // 初始d1, n/d1, mu1 return ans; }代码细节提醒伯努利数递推代码中使用的递推公式是 $B_m 1 - \sum_{j0}^{m-1} \binom{m}{j} \frac{B_j}{m-j1}$。注意分母是m-j1不要写成m-j。多项式系数coef数组下标i代表 $x^i$ 的系数。在get_poly_coef中我们只填充了exp从1到k1的系数coef[0]默认为0这与理论一致自然数幂和多项式没有常数项。枚举因子DFS函数中我们同时维护了d当前已选质因子的乘积和nd当前已计算的 n/d 的值。这样避免了最后再去计算n/d提高了效率。注意nd的计算如果不选该质因子则n/d中包含 $p^a$如果选则包含 $p^{a-1}$如果a1则 $p^01$。模运算所有乘法后都要立即取模。mu用1和MOD-1表示正负。在累加时(mu * d_pow_k % MOD * poly_val % MOD MOD) % MOD可以确保结果非负但更安全的写法是分步取模。n的计算在枚举过程中我们并不需要显式地计算出完整的n它可能巨大只需要n/d。我们通过质因子指数直接构造nd这是处理大数的关键。6. 总结与延伸思考回顾整个解题过程我们从一道要求互质数幂和的题目出发引入了莫比乌斯反演来化解互质条件又借助伯努利数得到的自然数幂和公式将问题转化为多项式求值。最终通过枚举n的平方自由因子即 $\mu(d) \ne 0$ 的因子我们能在 $O(2^w \cdot k)$ 的时间复杂度内解决问题其中w是n的不同质因子数。这道题的价值不仅在于其解法更在于它展示了一种典型的“组合数学数论”的解题范式通过莫比乌斯反演将数论条件转化为求和再通过已知的求和公式如伯努利数公式、斯特林数公式等将求和转化为多项式或其他封闭形式最后利用输入数据的结构如质因数分解进行高效计算。在这个过程中有几个点值得反复品味数学工具的选择为什么用莫比乌斯反演因为[gcd(m,n)1]这个条件天生适合用 $\sum_{d|gcd(m,n)} \mu(d)$ 来替换。为什么用伯努利数因为我们需要一个关于N的幂和 $\sum t^k$ 的精确多项式表达式。计算与优化理论公式到可执行代码之间有大量细节需要打磨。模运算下的分数处理、伯努利数的递推、多项式的求值秦九韶算法、利用质因数分解枚举因子并同时计算n/d这些都是将数学转化为高效算法的关键步骤。验证与调试对于这类推导复杂的题目用小的、可手算的案例如k0,1, n为小质数进行验证至关重要。它不仅能帮助发现代码错误更能加深对公式本身的理解。最后这个解法还可以进一步扩展。例如如果题目不是求互质数的k次幂和而是求gcd为某个定值g的数的k次幂和我们依然可以通过类似的莫比乌斯反演技巧将其转化为本题的形式。伯努利数也可以被斯特林数替代它们之间也有联系都能用于表达幂和。这些延伸就留给有兴趣的朋友去探索了。数论的世界就是这样一个精巧的问题往往能串联起多个美丽的数学概念而解决它的过程本身就是一次深刻的思维训练。