莫比乌斯反演与伯努利数：解决大规模自然数幂和问题的数论组合方法

1. 项目概述当数论难题遇上组合数学的“瑞士军刀”看到这个标题“P6271 [湖北省队互测2014]一个人的数论”很多搞算法竞赛的朋友可能会心一笑这味儿太对了。它完美融合了那种经典的、带着一丝孤独感的竞赛题名风格以及背后硬核的数学内核。这道题的核心是要求我们计算一个与幂和相关的数论函数在给定模数下的值。具体来说它通常形如求 $S_d(n) \sum_{i1}^{n} i^d$ 在 $n$ 以质因数分解形式给出且 $d$ 和模数 $M$ 给定的情况下的值。直接计算$n$ 可以巨大无比比如是一个大整数的质因数分解形式暴力求和是天文数字级的复杂度。这就是为什么我们需要请出标题里的两位“大神”莫比乌斯反演和伯努利数。你可以把这道题理解为一个“公式推导高效计算”的经典模板。莫比乌斯反演是处理数论求和、解决“计数问题”的利器它能将复杂的、带有整除条件的求和转化为相对简单的形式。而伯努利数则是一套系统处理自然数幂和 $1^k 2^k ... n^k$ 的强力工具它能将这个和表示为一个关于 $n$ 的 $k1$ 次多项式。这道题的巧妙之处就在于将这两个强大的工具串联起来先通过莫比乌斯反演化简问题形式再利用伯努利数提供的多项式公式进行高效计算。它非常适合有一定数论和组合数学基础想要深入理解如何将理论工具应用于解决具体、复杂计算问题的选手或爱好者。接下来我们就一层层剥开它的外壳看看里面的精妙构造。2. 核心思路拆解从暴力求和到公式化降维打击面对求 $S_d(n) \sum_{i1}^{n} i^d$ 这个问题最朴素的想法是循环累加。但当 $n$ 的规模达到 $10^9$ 甚至以质因数分解形式给出意味着其数值可能极大时$O(n)$ 的复杂度是完全不可接受的。我们需要寻找与 $n$ 大小无关或者只与 $n$ 的质因数个数相关的计算方法。2.1 为什么是莫比乌斯反演题目中“一个人的数论”这个背景往往暗示着求和可能带有某种“互质”或“因子相关”的条件。一个常见的变形是求 $\sum_{i1}^{n} [\gcd(i, n) 1] i^d$即求与 $n$ 互质的所有 $i$ 的 $d$ 次幂之和。这正是莫比乌斯反演的典型舞台。莫比乌斯函数 $\mu(n)$是一个定义在正整数上的函数其值为$\mu(1) 1$若 $n$ 含有平方因子即存在质数 $p$ 使得 $p^2 \mid n$则 $\mu(n) 0$若 $n$ 是 $k$ 个不同质数的乘积则 $\mu(n) (-1)^k$它的核心性质是$\sum_{d \mid n} \mu(d) [n1]$其中 $[n1]$ 是艾弗森括号当 $n1$ 时为1否则为0。这个性质是反演的基石。对于条件 $[\gcd(i, n) 1]$我们可以利用这个性质进行转换 $$[\gcd(i, n) 1] \sum_{d \mid \gcd(i, n)} \mu(d) \sum_{d \mid i, d \mid n} \mu(d)$$于是原求和式可以转化为 $$F(n) \sum_{i1}^{n} [\gcd(i, n) 1] i^d \sum_{i1}^{n} i^d \sum_{d \mid i, d \mid n} \mu(d)$$交换求和顺序这是关键一步 $$F(n) \sum_{d \mid n} \mu(d) \sum_{i1}^{n} i^d \cdot [d \mid i]$$令 $i d \cdot j$则当 $i$ 从1到 $n$ 且 $d \mid i$ 时$j$ 从1到 $\lfloor n/d \rfloor$。因此 $$F(n) \sum_{d \mid n} \mu(d) \sum_{j1}^{\lfloor n/d \rfloor} (d \cdot j)^d \sum_{d \mid n} \mu(d) d^d \sum_{j1}^{\lfloor n/d \rfloor} j^d$$看我们成功地将一个带有互质条件的复杂求和转化为了一个对 $n$ 的所有因子 $d$ 的求和而内层求和变成了标准的自然数幂和 $\sum j^d$。这就是莫比乌斯反演在此类问题中的威力化条件求和为因子枚举。2.2 伯努利数如何登场现在问题转化为计算形如 $S_d(m) \sum_{j1}^{m} j^d$ 的式子其中 $m \lfloor n/d \rfloor$。这正是伯努利数的用武之地。伯努利数 $B_k$是一个有理数序列它通过生成函数定义$\frac{x}{e^x - 1} \sum_{k0}^{\infty} B_k \frac{x^k}{k!}$。它最著名的应用之一就是给出自然数幂和的封闭表达式即Faulhaber公式 $$S_d(m) \sum_{j1}^{m} j^d \frac{1}{d1} \sum_{k0}^{d} \binom{d1}{k} B_k m^{d1-k}$$这个公式告诉我们$S_d(m)$ 是一个关于 $m$ 的 $d1$ 次多项式其系数由伯努利数和二项式系数共同决定。这意味着只要我们预先计算出前 $d$ 个伯努利数在模 $M$ 意义下那么对于任意给定的 $m$我们都可以在 $O(d)$ 的时间内计算出 $S_d(m)$这与 $m$ 的大小无关这实现了从 $O(m)$ 到 $O(d)$ 的复杂度飞跃。思路串联问题转化利用莫比乌斯反演将原问题可能与互质相关转化为对 $n$ 的所有因子 $d$ 求和内层是标准幂和 $S_d(\lfloor n/d \rfloor)$。核心计算利用伯努利数提供的 Faulhaber 公式将 $S_d(m)$ 的计算复杂度从 $O(m)$ 降至 $O(d)$。整体求解枚举 $n$ 的所有因子 $d$其个数是 $n$ 的因子数通常远小于 $n$ 本身对于每个 $d$用 $O(d)$ 时间计算内层和再乘以 $\mu(d) d^d$ 并累加最后模 $M$ 得到答案。这个框架清晰地将一个看似需要巨大计算量的问题分解成了因子枚举和多项式求值两个可高效解决的子问题。3. 关键细节与实操要点解析理论框架搭建好了但魔鬼藏在细节里。要实现这个方案有几个关键的技术点必须妥善处理否则极易出错或无法通过。3.1 伯努利数的计算与存储伯努利数是有理数而我们的计算是在模 $M$ 下进行的。这意味着我们需要计算伯努利数在模 $M$ 意义下的值。伯努利数可以通过递推公式计算 $$B_0 1, \quad \sum_{k0}^{m} \binom{m1}{k} B_k 0 \quad (m \ge 1)$$ 由此可以解出 $B_m -\frac{1}{m1} \sum_{k0}^{m-1} \binom{m1}{k} B_k$。实操要点模意义下的除法递推公式中存在除以 $m1$ 的操作。这要求 $m1$ 在模 $M$ 下有乘法逆元。如果 $M$ 不是质数或者 $d$ 较大时 $d1$ 可能与 $M$ 不互质这将导致致命错误这是本题一个非常关键的陷阱。通常竞赛题中 $M$ 会取一个质数如 $10^97$来规避这个问题。如果 $M$ 非质数则需要更复杂的处理如中国剩余定理分解模数。计算范围我们只需要前 $d1$ 个伯努利数$B_0$ 到 $B_d$。预处理复杂度为 $O(d^2)$在 $d$ 较小如 $d \le 1000$时是可接受的。存储格式直接存储模 $M$ 意义下的整数值即可。注意网上有些资料会提到用生成函数 $x/(e^x-1)$ 和FFT来 $O(d \log d)$ 求伯努利数但对于竞赛场景和通常的 $d$ 范围$d \le 1000$$O(d^2)$ 的递推实现更简单可靠且常数小。3.2 枚举因子与莫比乌斯函数值获取题目通常以 $n p_1^{a_1} p_2^{a_2} ... p_k^{a_k}$ 的形式给出 $n$ 的质因数分解。我们需要枚举 $n$ 的所有因子 $d$并同时知道每个因子 $d$ 对应的莫比乌斯函数值 $\mu(d)$。高效枚举方法采用DFS深度优先搜索遍历质因子的指数。初始因子 $d 1$, $\mu(d) 1$。对于每个质因子 $p_i$我们可以选择不乘 $p_i$因子 $d$ 不变$\mu$ 不变。乘 $p_i^1$因子 $d$ 变为 $d \times p_i$$\mu$ 变为 $-\mu$。乘 $p_i^t (t \ge 2)$此时 $d$ 含有平方因子根据定义 $\mu(d) 0$且该分支的所有后续乘积 $\mu$ 值均为0可以直接剪枝无需继续DFS。递归地进行上述选择即可生成所有 $\mu(d) \ne 0$ 的因子 $d$即无平方因子数并记录其 $\mu(d)$ 值。为什么只枚举无平方因子数因为对于含有平方因子的 $d$$\mu(d)0$它在求和式 $F(n) \sum_{d \mid n} \mu(d) ...$ 中的贡献为零所以根本不需要计算。这大大减少了需要枚举的因子数量。$n$ 的无平方因子数个数等于 $2^k$其中 $k$ 是 $n$ 的不同质因子个数。题目中 $k$ 通常不会太大。3.3 Faulhaber公式的高效求值公式 $S_d(m) \frac{1}{d1} \sum_{k0}^{d} \binom{d1}{k} B_k m^{d1-k}$ 需要计算一个 $d1$ 项的多项式求和。直接计算的时间复杂度是 $O(d^2)$如果每次独立计算组合数。我们可以优化到 $O(d)$。优化技巧预处理组合数可以预处理出 $C[d1][0...d]$但这需要 $O(d^2)$ 空间。更优的方法是在循环计算求和时利用组合数的递推关系 $\binom{d1}{k} \binom{d1}{k-1} \cdot \frac{d1-k1}{k}$ 来递推计算每一项的组合数这样只需 $O(1)$ 的额外空间和 $O(d)$ 的时间。快速幂计算 $m^{d1-k}$需要计算 $m$ 的从 $1$ 到 $d1$ 次幂。可以预处理 $m$ 的所有次幂或者从高次幂向低次幂计算时利用 $m^{t-1} m^t / m$ 的关系在模 $M$ 下做乘逆元。但更简单且高效的方法是在循环中维护一个变量pow_m 1从 $kd$ 循环到 $0$每次循环先使用当前的pow_m即 $m^{d1-k}$然后更新pow_m pow_m * m % MOD为下一次循环$k-1$做准备。这样也是 $O(d)$。综合起来计算一个 $S_d(m)$ 可以在 $O(d)$ 时间内完成。3.4 模运算的注意事项整个计算过程涉及大量模运算需要特别注意负数取模$\mu(d)$ 可能是 $-1$计算 $\mu(d) \cdot (\text{某项})$ 时可能产生负数。在加/减/乘运算后一定要(x % MOD MOD) % MOD确保结果在 $[0, MOD-1]$ 范围内。大数乘法$d^d$ 可能非常大即使取模前也需要小心溢出。在计算 $d^d \bmod MOD$ 时应使用快速幂算法并在乘法时使用(long long)中间类型或类似机制防止溢出。逆元的使用Faulhaber公式中的 $\frac{1}{d1}$ 需要计算 $d1$ 的模逆元。务必确保 $d1$ 与 $MOD$ 互质。4. 完整实现流程与代码框架下面我们以一个典型的题目设定为例勾勒出完整的实现步骤。假设求 $F(n) \sum_{i1}^{n} [\gcd(i, n)1] i^d \pmod{M}$$n$ 以质因数分解形式给出$n \prod_{i1}^{k} p_i^{a_i}$$d, M$ 给定$M$ 为质数如 $10^97$。4.1 步骤一预处理伯努利数const int MOD 1e97; const int MAXD 1005; // 假设d的最大范围 long long B[MAXD]; // 伯努利数 B[0]...B[d] long long inv[MAXD]; // 逆元 inv[i] i^(-1) mod MOD long long C[MAXD][MAXD]; // 组合数可选用于清晰表达 void initBernoulli(int d) { // 预处理逆元 inv[1] 1; for (int i 2; i d2; i) { inv[i] (MOD - MOD / i) * inv[MOD % i] % MOD; } // 计算伯努利数 B[0] 1; for (int m 1; m d; m) { long long sum 0; for (int k 0; k m; k) { long long comb 1; // 可以递推计算组合数 C(m, k)这里为清晰直接使用预计算或函数 // 假设有函数 getComb(m, k) sum (sum getComb(m1, k) * B[k]) % MOD; } B[m] (MOD - sum) * inv[m1] % MOD; } }4.2 步骤二DFS枚举因子并计算贡献long long p[105], a[105]; // 存储质因子p_i和指数a_i int k; // 不同质因子个数 long long d_val; // 题目中的d long long ans 0; // DFS函数枚举无平方因子数 // current: 当前因子值 // mu: 当前因子对应的莫比乌斯函数值 // idx: 当前处理到第几个质因子 void dfs(long long current, long long mu, int idx) { if (idx k) { // 枚举到一个因子 d current, 其莫比乌斯值为 mu if (current 1) return; // 因子1通常可能需要特殊处理视题目而定 long long m n / current; // 这里n是全局变量需要根据current计算 // 计算 S_d(m) 使用伯努利数公式 long long S calcPowerSum(d_val, m); // 计算贡献mu * (current^d_val) % MOD * S % MOD long long contrib mu * fastPow(current, d_val) % MOD; contrib contrib * S % MOD; ans (ans contrib) % MOD; return; } // 分支1不选当前质因子 dfs(current, mu, idx1); // 分支2选一次当前质因子指数为1 dfs(current * p[idx], -mu, idx1); // 指数2时mu0直接剪枝不继续递归 } // 计算 S_d(m) 的函数 long long calcPowerSum(long long d, long long m) { m % MOD; // 注意取模 long long res 0; long long pow_m fastPow(m, d1); // 计算 m^(d1) long long inv_m fastPow(m, MOD-2); // 计算 m 的逆元用于降幂 for (int k 0; k d; k) { long long comb getComb(d1, k); // 获取组合数 C(d1, k) long long term comb * B[k] % MOD; term term * pow_m % MOD; res (res term) % MOD; // 为下一轮准备pow_m m^(d1 - (k1)) m^(d-k) pow_m * inv_m % MOD pow_m pow_m * inv_m % MOD; } res res * inv[d1] % MOD; // 乘以 1/(d1) return res; }4.3 步骤三主逻辑与初始化int main() { // 读入 d_val, MOD, k, 以及 p[i], a[i] long long n 1; for (int i 0; i k; i) { n n * fastPow(p[i], a[i]) % MOD; // 这里计算n主要是为了后面求 m n/d注意n可能很大我们通常只在需要时用质因子表示计算 // 实际上在dfs中计算 m n / current 时我们并不直接计算巨大的n而是利用质因子分解 // n / current 相当于每个质因子指数相减。更稳妥的做法是在dfs时维护 current 的质因子指数状态从而直接得到 m 的质因子分解形式。 // 但为了代码简洁说明原理此处假设可以计算。 } // 初始化伯努利数需要用到 d_val initBernoulli(d_val); // 初始化组合数如果需要 initCombination(d_val5); // 执行DFS ans 0; dfs(1, 1, 0); // 从因子1mu1开始 // 注意dfs中可能没有处理因子d1的情况mu(1)1需要根据题目公式决定是否加上。 // 例如如果原公式是 F(n) sum_{d|n} mu(d) * d^d * S_d(n/d)那么d1的贡献是 1^d * S_d(n) S_d(n) // 我们可以在dfs外部加上这部分 long long S_n calcPowerSum(d_val, n); ans (ans S_n) % MOD; // 加上 d1 的贡献 // 输出答案 cout (ans MOD) % MOD endl; return 0; }实操心得在实际编码中n可能非常大直接计算n/current会导致溢出即使在模意义下除法也不直接。更标准的做法是在DFS过程中不仅维护因子current还维护当前因子计算后剩余的“商”的质因子表示即每个质因子指数a[i]减去当前因子选取的指数。这样我们可以直接得到m n/d的质因子分解形式进而可以用快速幂计算出m % MOD的值用于伯努利公式的计算。这是实现中的一个关键细节考验对数论和代码的结合能力。5. 常见问题、调试技巧与扩展思考即使理解了算法实现时也难免踩坑。下面记录一些常见问题和排查思路。5.1 结果错误或为负检查点1模运算与负数。这是最常见的问题。确保每一次加法、减法、乘法后都立即取模并且对于可能为负的结果使用(x % MOD MOD) % MOD修正。检查点2伯努利数计算。验证递推公式是否正确特别是模逆元的计算。可以手动计算前几个伯努利数$B_01, B_1-1/2, B_21/6, B_30, B_4-1/30$在模 $M$ 下的值进行对比。检查点3因子枚举与贡献公式。确认题目要求的原始公式是否与你实现的公式完全一致。是 $\sum_{d \mid n} \mu(d) d^d S_d(n/d)$ 还是 $\sum_{d \mid n} \mu(d) d^k S_d(n/d)$$k$ 可能为 $d$ 或其他仔细核对每一项的指数。检查点4$d1$ 的特例。确认因子 $1$ 是否被正确处理。在DFS中因子 $1$ 通常对应mu1但current1时current^d也是1S_d(n/1)S_d(n)。确保这部分贡献被计入。5.2 运行超时优化点1减少模运算。在内部循环如计算 Faulhaber 公式的求和中可以适当减少取模次数例如累积多次乘法后再取模但要注意防止溢出。优化点2预处理幂次。在calcPowerSum函数中我们通过乘以逆元来递推 $m$ 的幂。如果d很大频繁计算逆元可能稍慢。可以预处理出 $m^0, m^1, ..., m^{d1}$ 的所有值虽然空间复杂度 $O(d)$但可能更快。优化点3检查枚举范围。确认DFS正确剪枝了 $\mu(d)0$即含有平方因子的情况。如果 $k$质因子个数很大比如20$2^k$ 的枚举也可能压力大但这类题目通常会将 $k$ 控制得较小。5.3 模数 $M$ 非质数的处理这是一个进阶难点。如果 $M$ 不是质数那么计算逆元可能失败。此时需要分解模数将 $M$ 分解为若干质数幂的乘积 $M \prod q_i^{e_i}$。分别求解对于每个 $q_i^{e_i}$ 作为模数运行整个算法此时需要保证计算过程中的分母与 $q_i$ 互质否则需要特殊处理线性递推求序列。这要求伯努利数递推、Faulhaber公式中的除法都能在该模数下进行。中国剩余定理CRT合并得到每个模数下的答案后使用中国剩余定理合并得到模 $M$ 下的最终答案。这大大增加了代码复杂度通常只会在更高级别的比赛中出现。5.4 扩展与变式掌握了这个框架你可以解决一系列类似问题变式1求 $\sum_{i1}^{n} i^d [\gcd(i, n)g]$。可以通过变量代换 $i g \cdot j$转化为 $\sum_{j1}^{n/g} (gj)^d [\gcd(j, n/g)1]$即 $g^d \cdot F(n/g)$其中 $F$ 就是我们上面求解的函数。变式2$d$ 非常大如 $10^5$。此时 $O(d^2)$ 预处理伯努利数不可行。需要使用基于生成函数和FFT/NTT的 $O(d \log d)$ 方法求伯努利数或者寻找其他数学变换。变式3求和范围不是 $1$ 到 $n$而是 $a$ 到 $b$。利用 $S_d(b) - S_d(a-1)$ 即可。这道题的精髓在于它教会我们如何将复杂的数论条件求和通过莫比乌斯反演进行分解再借助伯努利数这类组合数学工具将大规模求和问题转化为小规模多项式求值和因子枚举问题。这种“转化与归约”的思想在算法竞赛和数学问题求解中至关重要。在实际敲代码时耐心调试模运算和边界情况理解每一个步骤的数学含义才能最终稳稳地拿下这类题目。