莫比乌斯反演与伯努利数：解决互质幂和问题的数论利器

1. 项目概述当数论难题遇上组合数学的“瑞士军刀”看到这个标题很多搞算法竞赛或者对数论感兴趣的朋友可能会心头一紧。[湖北省队互测2014]一个人的数论这个名头听起来就很有分量再加上“莫比乌斯反演”和“伯努利数”这两个组合妥妥的一道硬核数论题。我第一次接触这类题目时感觉就像面对一个结构精密的密码锁你知道它肯定能被解开但手里只有几把形状各异的钥匙数学工具需要找到正确的组合方式和开锁顺序。这道题的核心简单来说是要求解一个关于正整数n的幂和函数S_d(n) 1^d 2^d ... n^d在n的质因数分解形式已知但n本身巨大甚至可以是形式表达式时的值或其某种算术性质。这里的d是一个给定的非负整数。如果n是一个具体的数字我们或许可以暴力计算或者用已知的幂和公式。但题目往往给出n p1^{a1} * p2^{a2} * ... * pk^{ak}这样的形式要求你给出一个关于这些质因子pi及其指数ai的表达式。这就把问题从“计算”提升到了“推导”和“表示”的层面。为什么需要莫比乌斯反演和伯努利数这就像是解决这个问题的两把关键钥匙。莫比乌斯反演擅长处理“整除”、“因子”关系下的求和转换它能把一个看似复杂的、带有限制条件的求和转化为另一个更容易处理的求和式。而伯努利数则为我们提供了将幂和S_d(n)表示为一个关于n的d1次多项式的系统工具。这个多项式的系数就由伯努利数决定。单独使用任何一个工具可能都无法优雅地解决这个“一个人的数论”问题但将它们串联起来就能构建出一条清晰的解题路径先用伯努利数公式将幂和多项式化再利用题目中n的质因数分解形式结合莫比乌斯反演处理多项式求值中涉及的因子和问题。接下来我会把自己在理解、推导和实现这道题思路过程中的所有细节、技巧和踩过的坑毫无保留地拆解开来。无论你是正在备战竞赛还是单纯想深入理解这两个优美的数学工具如何协同工作这篇内容都将带你从问题本质出发一步步走到最终那个简洁的表达式。2. 核心数学工具莫比乌斯反演与伯努利数深度解析要攻克这道题必须对这两件“武器”了如指掌。我们不能仅仅停留在会背公式的层面更要理解它们为何在此处能发挥作用。2.1 莫比乌斯反演穿透整除关系的“X光”莫比乌斯函数μ(n)定义很简单μ(1) 1如果n含有平方因子即存在质数p使得p^2 | n那么μ(n) 0否则n是k个不同质数的乘积则μ(n) (-1)^k它的核心价值体现在莫比乌斯反演公式上。最常见的形式有两种因数形式若F(n) Σ_{d|n} f(d)则f(n) Σ_{d|n} μ(d) F(n/d)。倍数形式若F(n) Σ_{n|d} f(d)则f(n) Σ_{n|d} μ(d/n) F(d)。在本题相关的推导中倍数形式更为常用。它为什么是“X光”想象一下F(n)是一个关于整数n的“总体测量”它包含了所有其因子或倍数d的贡献f(d)。这个求和关系可能很复杂。莫比乌斯反演就像一束X光通过μ函数的正负抵消作用这正是其定义的精妙之处能从F(n)中逆向透视并分离出我们真正关心的单个f(n)。这种“求和”与“反演”的转换在数论中处理包含gcd,lcm, 整除条件的计数问题时威力巨大。实操心得记忆和使用反演公式时务必厘清求和下标。d|n表示d整除n求和遍历n的所有正因子。n|d则表示n整除d求和遍历n的所有倍数。写代码或推导时下标写错一步满盘皆输。2.2 伯努利数为幂和穿上多项式的外衣计算1^d 2^d ... n^d是数学中一个经典问题。伯努利数B_k是一列有理数它们通过生成函数定义t/(e^t - 1) Σ_{k0}^{∞} B_k * t^k / k!。对我们而言最关键的是伯努利数将幂和表示为多项式的公式Faulhaber公式S_d(n) 1^d 2^d ... n^d 1/(d1) * Σ_{k0}^{d} C(d1, k) * B_k * n^{d1-k}其中C(a,b)是组合数。令d1-k j也可以写成S_d(n) 1/(d1) * Σ_{j1}^{d1} C(d1, j) * B_{d1-j} * n^j这个公式的颠覆性在于它告诉我们尽管S_d(n)是d次幂的求和但它本身却是n的一个d1次多项式。系数由d和伯努利数共同决定。对于固定的d这些系数是常数。例如d1:S_1(n) 12...n n(n1)/2 (1/2)n^2 (1/2)n对应B_01, B_1-1/2。d2:S_2(n) 1^22^2...n^2 n(n1)(2n1)/6 (1/3)n^3 (1/2)n^2 (1/6)n。在本题中的应用逻辑题目通常不会直接求S_d(n)而是求S_d(n)模某个数或者求Σ_{i1}^n [gcd(i, n)1] * i^d这类与n互质的数的幂和。这时我们可以先利用伯努利公式将i^d的求和转化为n的多项式再利用n的质因数分解形式结合包含gcd条件的莫比乌斯反演将“互质”的条件转化为对n的因子的求和从而最终得到一个关于pi和ai的表达式。注意事项伯努利数有B_1 -1/2和B_1 1/2两种约定取决于生成函数是t/(e^t-1)还是t/(1-e^{-t})。在算法竞赛和大多数数论文献中通常采用B_1 -1/2的约定即上面使用的。使用任何数学库或参考资料前务必确认其约定否则会导致多项式系数符号错误。3. 问题拆解与核心推导流程我们以一道典型的衍生问题为例进行完整推导计算F_d(n) Σ_{i1}^{n} i^d * [gcd(i, n) 1]其中[gcd(i, n)1]是艾弗森括号当gcd(i,n)1时为1否则为0。已知n的质因数分解n Π_{i1}^{k} pi^{ai}。3.1 第一步引入莫比乌斯函数转化互质条件这是最关键的一步。数论中有一个非常实用的等式[gcd(i, n) 1] Σ_{d | gcd(i, n)} μ(d)因为μ(d)的和在gcd(i,n)1时为1只有d1一项贡献1在gcd(i,n)1时根据莫比乌斯函数定义其因子和的净结果为0。因此我们可以重写求和式F_d(n) Σ_{i1}^{n} i^d * [gcd(i, n) 1] Σ_{i1}^{n} i^d * (Σ_{d | gcd(i, n)} μ(d))交换求和顺序这是处理此类问题的标准操作F_d(n) Σ_{d1}^{n} μ(d) * (Σ_{i1, d|gcd(i,n)}^{n} i^d)条件d | gcd(i, n)等价于d | i且d | n。由于d要能整除n所以外层求和d实际上只需遍历n的所有正因子F_d(n) Σ_{d | n} μ(d) * (Σ_{i1, d|i}^{n} i^d)3.2 第二步处理内层求和联系伯努利公式内层求和Σ_{i1, d|i}^{n} i^d表示对所有n以内d的倍数i求i^d的和。令i d * j则当i从1遍历到n且d|i时j从1遍历到floor(n/d)。代入Σ_{i1, d|i}^{n} i^d Σ_{j1}^{n/d} (d*j)^d d^d * Σ_{j1}^{n/d} j^d d^d * S_d(n/d)看这里出现了我们熟悉的幂和S_d(m)其中m n/d。3.3 第三步代入伯努利多项式公式根据伯努利公式S_d(m) 1/(d1) * Σ_{j1}^{d1} C(d1, j) * B_{d1-j} * m^j。 将其代入F_d(n) Σ_{d | n} μ(d) * d^d * [1/(d1) * Σ_{j1}^{d1} C(d1, j) * B_{d1-j} * (n/d)^j] 1/(d1) * Σ_{j1}^{d1} C(d1, j) * B_{d1-j} * n^j * [ Σ_{d | n} μ(d) * d^{d-j} ]注意d^d * (1/d)^j d^{d-j}。3.4 第四步利用积性函数性质化简因子求和现在我们面对的是形如G_j(n) Σ_{d|n} μ(d) * d^{d-j}的求和。μ(d)是积性函数d^{d-j}对于固定的d-j注意这里指数是d-jd是求和变量不是常数所以d^{d-j}本身不是积性函数这造成了主要困难。这里需要一点技巧。我们回到F_d(n)的另一个等价的、更简洁的表达式。由莫比乌斯反演的倍数形式有F_d(n) Σ_{i1}^n i^d * Σ_{d|i, d|n} μ(d) Σ_{d|n} μ(d) * d^d * S_d(n/d)这和之前一样。关键的观察点对于S_d(n/d)我们可以直接将其视为关于(n/d)的多项式记为P_d(x) 1/(d1) * Σ_{j0}^{d} C(d1, j) B_j x^{d1-j}这里调整了索引与之前等价。那么F_d(n) Σ_{d|n} μ(d) * d^d * P_d(n/d)由于n的质因数分解已知n Π p_i^{a_i}且μ(d)只在d为无平方因子数即d是某些不同质因子的乘积时非零。设d Π_{i in S} p_i其中S是{1,2,...,k}的一个子集。则n/d Π p_i^{a_i} / Π_{i in S} p_i Π_{i in S} p_i^{a_i -1} * Π_{i not in S} p_i^{a_i}。此时P_d(n/d)是一个关于n/d的d1次多项式。将n/d的表达式代入整个F_d(n)就可以展开为一系列关于各p_i的幂次的项的和。由于μ(d)的取值(-1)^{|S|}和d^d都可以明确写出理论上可以通过枚举所有子集S来计算。但当k质因子个数很大时枚举子集不可行。更优雅的做法注意到F_d(n)对于n是积性函数吗我们来检验。若gcd(a,b)1则F_d(ab) Σ_{i1}^{ab} i^d [gcd(i, ab)1] Σ_{i1}^{ab} i^d [gcd(i,a)1 且 gcd(i,b)1]。 这并不直接等于F_d(a)*F_d(b)因为i的取值范围和分解方式有问题。但是考虑中国剩余定理模ab的互质剩余类可以唯一分解为模a和模b的互质剩余类的笛卡尔积。实际上可以证明F_d(n)是积性函数。一种证明思路是利用狄利克雷卷积F_d(n) (id_d * μ)(n)其中id_d(n)n^d。而id_d和μ都是积性函数它们的狄利克雷卷积也是积性的。既然F_d(n)是积性函数那么对于n Π p_i^{a_i}有F_d(n) Π F_d(p_i^{a_i})。问题瞬间简化我们只需要计算F_d(p^a)对于单个质数幂的表达式。3.5 第五步计算单个质数幂情形F_d(p^a)F_d(p^a) Σ_{i1}^{p^a} i^d [gcd(i, p^a)1]。 在1到p^a中与p^a不互质的数就是p的倍数共有p^{a-1}个。因此我们可以用总数减去p的倍数的贡献F_d(p^a) S_d(p^a) - Σ_{i1, p|i}^{p^a} i^d S_d(p^a) - Σ_{j1}^{p^{a-1}} (p*j)^d S_d(p^a) - p^d * S_d(p^{a-1})看这里又出现了S_d(m)。将伯努利公式S_d(x) (1/(d1)) Σ_{j0}^{d} C(d1, j) B_j x^{d1-j}代入F_d(p^a) (1/(d1)) Σ_{j0}^{d} C(d1, j) B_j [ p^{a(d1-j)} - p^d * p^{(a-1)(d1-j)} ] (1/(d1)) Σ_{j0}^{d} C(d1, j) B_j * p^{(a-1)(d1-j) d} * [ p^{d1-j} - 1 ] (1/(d1)) Σ_{j0}^{d} C(d1, j) B_j * p^{a(d1-j) - (d1-j) d} * (p^{d1-j} - 1) (1/(d1)) Σ_{j0}^{d} C(d1, j) B_j * p^{(a-1)(d1-j) d} * (p^{d1-j} - 1)这个表达式虽然看起来复杂但对于给定的d和p^a它是一个关于p的明确表达式。因为d通常很小题目常限制d 100或类似所以这个求和项数有限可以快速计算。3.6 第六步整合最终答案由于F_d(n)是积性函数且n Π_{i1}^{k} p_i^{a_i}那么最终答案就是Ans F_d(n) Π_{i1}^{k} F_d(p_i^{a_i})其中每个F_d(p_i^{a_i})由上面的公式计算。至此我们完成了一个典型的、基于莫比乌斯反演和伯努利数求解“与n互质的数的幂和”问题的完整推导。原始题目[湖北省队互测2014]一个人的数论可能在此基础上有一些变化例如模一个特定数或者求的是Σ_{i1}^n [gcd(i, n)1] * i^d对某个大质数取模的值但核心的数学推导框架完全一致。踩坑实录在推导F_d(p^a)时最容易出错的地方是S_d(p^a) - p^d * S_d(p^{a-1})这一步的指数运算。(p*j)^d p^d * j^d所以第二项是p^d * S_d(p^{a-1})而不是p * S_d(p^{a-1})或S_d(p^{a-1})^d。务必仔细处理代数运算。4. 算法实现与关键代码剖析理论推导完成后我们需要将其转化为可运行的代码。这里以计算F_d(n) mod MM为大质数为例给出实现步骤和代码细节。假设输入为d幂次k质因子个数以及k对(p_i, a_i)。4.1 预处理伯努利数与组合数由于d较小我们可以预处理出前d2个伯努利数因为公式中用到B_j,j从0到d。伯努利数可以通过递归关系或生成函数求逆得到。递归关系Σ_{k0}^{m} C(m1, k) * B_k 0对于m 1。由此可递推B_m -1/(m1) * Σ_{k0}^{m-1} C(m1, k) * B_k同时需要预处理组合数C(n, m) mod M可以用递推公式C(n, m) C(n-1, m-1) C(n-1, m)或者预先计算阶乘和逆元。// 假设 MOD 是大质数 const int MOD 1e9 7; const int MAXD 105; // 根据d的最大值设定 long long bernoulli[MAXD]; // 伯努利数 B[0], B[1], ... B[d] long long comb[MAXD][MAXD]; // 组合数 long long inv[MAXD]; // 逆元 inv[i] i^(-1) mod MOD // 预处理逆元 void initInv() { inv[1] 1; for (int i 2; i MAXD; i) { inv[i] (MOD - MOD / i) * inv[MOD % i] % MOD; } } // 预处理组合数 void initComb() { comb[0][0] 1; for (int i 1; i MAXD; i) { comb[i][0] comb[i][i] 1; for (int j 1; j i; j) { comb[i][j] (comb[i-1][j-1] comb[i-1][j]) % MOD; } } } // 预处理伯努利数使用递归公式结果模 MOD // 注意这里计算的是模 MOD 下的伯努利数对于有理数 B_k p/q我们存储 p * q^(-1) mod MOD void initBernoulli(int d) { bernoulli[0] 1; for (int m 1; m d; m) { long long sum 0; for (int k 0; k m; k) { sum (sum comb[m1][k] * bernoulli[k]) % MOD; } bernoulli[m] (MOD - sum) * inv[m1] % MOD; } // 根据约定B_1 -1/2需要特别调整。 // 我们的递推公式基于 Σ C(m1,k)B_k0该定义下 B1 -1/2。 // 验证当 m1 时公式为 C(2,0)B0 C(2,1)B1 0 1*1 2*B1 0 B1 -1/2。 // 所以递推得到的结果已经是正确的。 }4.2 实现单点计算F_d(p^a) mod MOD根据公式F_d(p^a) (1/(d1)) Σ_{j0}^{d} C(d1, j) B_j * [ p^{a(d1-j)} - p^d * p^{(a-1)(d1-j)} ]我们实现一个函数。// 快速幂 long long pow_mod(long long x, long long n) { long long res 1; x % MOD; while (n) { if (n 1) res res * x % MOD; x x * x % MOD; n 1; } return res; } // 计算 F_d(p^a) % MOD long long calc_F(long long p, long long a, int d) { long long res 0; long long inv_d1 inv[d1]; // 1/(d1) for (int j 0; j d; j) { long long exponent1 a * (d 1 - j); long long exponent2 d (a - 1) * (d 1 - j); long long term1 pow_mod(p, exponent1); long long term2 pow_mod(p, exponent2); long long diff (term1 - term2 MOD) % MOD; // p^{a(d1-j)} - p^{d(a-1)(d1-j)} long long coeff comb[d1][j] * bernoulli[j] % MOD; res (res coeff * diff) % MOD; } res res * inv_d1 % MOD; return res; }4.3 主逻辑与最终计算主函数读取d和k然后读取k个(p_i, a_i)利用积性计算总答案。int main() { int d, k; cin d k; initInv(); initComb(); initBernoulli(d); // 预处理到 d 即可 long long ans 1; for (int i 0; i k; i) { long long p, a; cin p a; ans ans * calc_F(p, a, d) % MOD; } cout ans endl; return 0; }代码细节与优化指数运算exponent2的计算d (a-1)*(d1-j)来源于公式推导(a-1)(d1-j) d。务必确保与推导一致。负数处理term1 - term2可能为负数需要 MOD再取模。逆元预处理频繁使用inv[d1]和组合数中的除法预处理逆元能大幅提升效率。模运算一致性所有运算包括伯努利数的递推都应在模MOD意义下进行。因为我们最终只要模MOD的结果。伯努利数为有理数在模质数MOD下我们可以将分数p/q表示为p * q^(-1) mod MOD。这正是我们在initBernoulli中使用逆元inv[m1]的原因。5. 常见问题、调试技巧与扩展思考即使理解了原理和算法实现时依然会遇到各种问题。下面是我在多次实现和调试中积累的一些经验。5.1 典型错误与排查清单问题现象可能原因排查方法答案输出为01. 模运算中乘法溢出未取模。2.calc_F中diff计算错误导致结果为0。3. 伯努利数B_1的符号约定错误。1. 检查所有乘法和加法后是否及时% MOD。2. 输出中间变量term1,term2,diff检查指数计算是否正确。3. 验证B[1]的值应为(MOD - inv[2]) % MOD即-1/2。答案与暴力对拍小数据不符1. 组合数C(n,m)预处理错误或越界。2. 积性函数性质使用错误误用于非积性函数。3.F_d(p^a)公式推导错误。1. 编写暴力程序计算小n和d的F_d(n)进行对拍。2. 单独测试calc_F函数用几个小的p^a和d手动计算验证。3. 重新检查从S_d(n)到F_d(p^a)的每一步代数变形。程序运行超时1. 快速幂pow_mod未使用long long或效率低。2. 对每个质因子重复计算了高次幂未利用预处理。1. 确保pow_mod参数和返回值是long long。2. 本题中d很小指数a*(d1-j)可能很大但快速幂是O(log n)的可以接受。主要检查是否有不必要的重复计算。5.2 调试技巧构建暴力验证程序对于数论问题编写一个针对小范围的暴力程序进行对拍是黄金法则。// 暴力计算 F_d(n) sum_{i1..n} (i^d) * [gcd(i,n)1] long long brute_force(int n, int d, int mod) { long long res 0; for (int i 1; i n; i) { if (std::gcd(i, n) 1) { long long term 1; for (int j 0; j d; j) term (term * i) % mod; // 计算 i^d % mod res (res term) % mod; } } return res; }用这个函数去验证calc_F以及最终积性乘法的结果。可以从n1到30d0到5进行测试。一旦暴力与优化算法结果一致信心就大大增强了。5.3 扩展思考题目可能的变化形式“一个人的数论”这类题目本质是考察对积性函数、狄利克雷卷积、莫比乌斯反演和伯努利数的综合运用。除了上述标准形式还可能有以下变体求和条件变化求Σ_{i1}^n [gcd(i, n)g] * i^d其中g是n的某个因子。可以通过变量替换i g * j转化为gcd(j, n/g)1的问题。多次询问给定固定的d和不同的n以质因数分解形式给出要求快速计算多个F_d(n)。这时我们的算法已经是O(k*d)的非常高效预处理伯努利数和组合数后每次询问只需遍历每个质因子进行计算。模数非质数如果模数M不是质数那么求逆元会变得复杂。伯努利数是有理数需要分别处理分子分母模M。一种方法是使用Python的大整数直接计算有理数最后取模或者使用中国剩余定理将M分解为质数幂的乘积分别计算后合并。与欧拉函数结合当d0时i^d 1F_0(n)就是与n互质的数的个数即欧拉函数φ(n)。我们的公式应能退化到φ(n) n * Π_{p|n} (1 - 1/p)。验证这一点是检查公式正确性的好方法。5.4 关于伯努利数计算的进一步优化当d较大如几百时递推求伯努利数的O(d^2)复杂度可能成为瓶颈。可以使用生成函数t/(e^t-1)结合多项式求逆在O(d log d)的时间内计算出前d项伯努利数模质数。这需要用到快速傅里叶变换FFT或数论变换NTT。这在算法竞赛的高阶题目中有时会涉及。// 伪代码思路伯努利数 B(x) 的指数生成函数 EGF 是 x/(e^x-1)。 // 令 F(x) (e^x-1)/x Σ_{k0} x^k/(k1)!。 // 那么 B(x) 1/F(x)。因此可以通过多项式求逆得到 B(x) 的系数。 // 注意这是 EGF最终需要乘以 k! 得到通常的伯努利数 B_k。不过对于本题及大多数竞赛题d 100左右O(d^2)的递推完全足够。最后一点个人体会数论问题的推导过程往往比最终代码复杂得多。像“一个人的数论”这样的题目其价值在于训练我们将一个复杂的求和问题通过莫比乌斯反演进行转化再利用伯努利数等工具进行化简最终归结到积性函数的性质上从而得到高效的算法。这个过程锻炼的是一种“分解”和“转化”的数学思维能力。在实现时耐心和细心至关重要尤其是处理模运算和指数运算时一个符号的错误就可能让整个程序输出错误的结果。多设置中间变量输出多写暴力程序对拍是调试这类问题最有效的方法。