基于容斥公式的Mertens函数M(n)界估计方法的可行性探讨
首先得说，你用容斥原理推导Mertens函数$M(n)$的思路相当巧妙，直接把$M(n)$和素数计数函数$\pi(n)$建立了直观联系，这是个很自然的切入角度。
先梳理下你推导的核心内容：
原始容斥表达式：
$$M(n)=1-\pi\left(n\right)+\sum_{p_{i}\leq\frac{n}{p_i}}\left(\pi\left(\frac{n}{p_{i}}\right)-i\right)-\sum_{p_{i}<p_{j}\leq\frac{n}{p_ip_j}}\left(\pi\left(\frac{n}{p_{i}p_{j}}\right)-j\right)+\sum_{p_{i}<p_{j}<p_{k}\leq\frac{n}{p_ip_jp_k}}\left(\pi\left(\frac{n}{p_{i}p_{j}p_{k}}\right)-k\right)-\dots$$
你后来将其拆分为两个大项：
$$M(n)=\left(-\pi\left(n\right)+\sum_{p_{i}\leq\frac{n}{p_i}}\pi\left(\frac{n}{p_{i}}\right)-\sum_{p_{i}<p_{j}\leq\frac{n}{p_ip_j}}\pi\left(\frac{n}{p_{i}p_{j}}\right)+\sum_{p_{i}<p_{j}<p_{k}\leq\frac{n}{p_ip_jp_k}}\pi\left(\frac{n}{p_{i}p_{j}p_{k}}\right)-\dots\right)-\left(\sum_{p_{i}\leq\frac{n}{p_i}}i-\sum_{p_{i}<p_{j}\leq\frac{n}{p_ip_j}}j+\sum_{p_{i}<p_{j}<p_{k}\leq\frac{n}{p_ip_jp_k}}k-\dots\right)$$
代入素数定理近似后得到了表达式(1)，而且你通过数值计算发现，拆分后的两个表达式在x轴附近呈现准对称振荡的特点——这个观察非常有价值，这种抵消效应或许能帮你简化后续的界估计工作。
关于方法可行性的分析
你的思路理论上完全可行：
- 像Dussart给出的$\pi(x)$紧界$\frac{x}{\log x -1} < \pi(x) < \frac{x}{\log x -1.1}$，代入后确实能给每一项套上上下界，进而得到$M(n)$的区间估计；
- 对于小$n$的场景，这种方法能给出精确的计算结果，数值观察到的抵消效应还能进一步优化计算效率。
但也需要注意几个核心难点：
- 截断误差控制：容斥表达式是无穷级数，实际计算只能取有限项，你需要严格估计截断后剩余项的贡献，否则得到的界会不够严谨。而且随着$n$增大，需要考虑的素数组合项数会急剧增长，误差控制的难度也会直线上升。
- 计算复杂度：如果要得到严格的界，每一项代入$\pi(x)$的界后，需要处理多层求和的上下界，计算量会非常大，尤其是大$n$场景下，这种方法的效率远低于成熟的解析数论方法。
- 现有方法的对比：目前学界对$M(n)$的非平凡界（比如$M(n)=O(n^{1/2+\epsilon})$，或者更紧的指数型界）都是通过解析数论方法得到的——利用黎曼ζ函数的零点分布、复积分技巧等，这些方法能更高效地得到大$n$下的渐近界，而容斥方法更适合小$n$的精确计算或者启发式分析。
后续建议
如果你坚持走容斥这条路，可以试试这几个方向：
- 先聚焦于估计截断误差，比如证明当素数乘积超过某个阈值后，剩余项的绝对值可以被某个小量控制；
- 利用你观察到的两个拆分表达式的准对称特性，尝试证明它们的绝对值相近，从而把$M(n)$的界转化为其中一个表达式的界，大幅简化计算；
- 先在小$n$范围内验证方法的有效性，再逐步推广到大$n$场景。
如果你的目标是得到大$n$下的非平凡界，那解析数论的方法会更成熟高效，建议参考相关经典文献（比如关于Mertens函数的渐近分析、ζ函数零点的应用）。
你用于数值实验的代码
@author: juanmoreno import math import matplotlib.pyplot as plt def sieve_of_eratosthenes(n): is_prime = [True] * (n + 1) is_prime[0] = is_prime[1] = False primes = [] for p in range(2, int(math.sqrt(n)) + 1): if is_prime[p]: for i in range(p * p, n + 1, p): is_prime[i] = False for i in range(2, n + 1): if is_prime[i]: primes.append(i) return primes def calculate_expression_1(n): return (-n / math.log(n)) def calculate_expression_2(n, primes): result = 0 result_2 = 0 for pi in primes: if pi < n / pi: result += (n / (pi * (math.log(n) - math.log(pi)))) result_2 -= primes.index(pi) return result, result_2 def calculate_expression_3(n, primes): result = 0 result_3 = 0 for i in range(len(primes)): for j in range(i + 1, len(primes)): pij = primes[i] * primes[j] if primes[j] < n / pij: result -= (n / (pij * (math.log(n) - math.log(primes[i]) - math.log(primes[j])))) result_3 += j return result, result_3 def calculate_expression_4(n, primes): result = 0 result_4 = 0 for i in range(len(primes)): for j in range(i + 1, len(primes)): pij = primes[i] * primes[j] if pij < n / pij: for k in range(j + 1, len(primes)): pijk = pij * primes[k] if primes[k] <= n / pijk: result += (n / (pijk * (math.log(n) - math.log(primes[i]) - math.log(primes[j]) - math.log(primes[k])))) result_4 -= k return result, result_4 def calculate_expression_5(n, primes): result = 0 result_5 = 0 for i in range(len(primes)): for j in range(i + 1, len(primes)): pij = primes[i] * primes[j] if pij < n / pij: for k in range(j + 1, len(primes)): pijk = pij * primes[k] if pijk < n / pijk: for l in range(k + 1, len(primes)): pijkl = pijk * primes[l] if primes[l] <= n / pijkl: result -= (n / (pijkl * (math.log(n) - math.log(primes[i]) - math.log(primes[j]) - math.log(primes[k])-math.log(primes[l])))) result_5 += l return result, result_5 # Create lists to store results and square roots of n results_1 = [] results_2 = [] sqrt_n = [] x_axis = [] for n in range(2, 10001): primes = sieve_of_eratosthenes(n) result1 = calculate_expression_1(n) result2 = calculate_expression_2(n, primes)[0] result3 = calculate_expression_3(n, primes)[0] result4 = calculate_expression_4(n, primes)[0] result5 = calculate_expression_5(n, primes)[0] final_result_1 =1 + result1 + result2 + result3 + result4 + result5 result21 = calculate_expression_2(n, primes)[1] result31 = calculate_expression_3(n, primes)[1] result41 = calculate_expression_4(n, primes)[1] result51 = calculate_expression_5(n, primes)[1] final_result_2 = result21 + result31 + result41 + result51 results_1.append(final_result_1) results_2.append(final_result_2) sqrt_n.append(math.sqrt(n)) x_axis.append(0) # Create a plot to compare results and sqrt(n) plt.figure(figsize=(10, 6)) plt.plot(range(2, 10001), results_1, label='Expression (1)', linewidth=2) plt.plot(range(2, 10001), results_2, label='Expression (2)', linewidth=2) plt.plot(range(2, 10001), sqrt_n, label='square root of n', linestyle='--', linewidth=2) plt.plot(range(2, 10001), x_axis, label='x-axis', linestyle='--', linewidth=2) plt.xlabel('n') plt.ylabel('Value') plt.legend() plt.title('Comparison of Final Result and x-axis') plt.grid(True) plt.show()
备注：内容来源于stack exchange，提问作者Juan Moreno
