孪生素数分布的大规模计算实验:分段筛实现与经验统计分析
首页 文章 正文 2026年02月22日 19:58 17

孪生素数猜想断言存在无穷多对形如 (p,p+2)(p,p+2) 的素数,但该猜想至今仍未解决。Hardy–Littlewood 的第一猜想给出了孪生素数计数函数的经典启发式渐近公式,而张益唐、Maynard 与 Polymath8 的工作则证明了存在无穷多对素数,其间隔有统一上界,但这一上界尚未降至 2 [1–3]。

摘要

孪生素数猜想断言存在无穷多对形如 (p,p+2)(p,p+2) 的素数,但该猜想至今仍未解决。Hardy–Littlewood 的第一猜想给出了孪生素数计数函数的经典启发式渐近公式,而张益唐、Maynard 与 Polymath8 的工作则证明了存在无穷多对素数,其间隔有统一上界,但这一上界尚未降至 2 [1–3]。

本文从实验计算角度出发,设计并实现了一套基于 C++ 的分段筛程序,对区间 [0,108][0,108] 内的素数、孪生素数对以及若干自定义筛选统计量进行大规模精确或近精确计算。程序首先生成不超过 NN 的基础素数表,然后对 [0,N][0,N] 进行分段筛;在此基础上,引入 wheel 条件 gcd⁡(n,2310)=gcd⁡(n+2,2310)=1gcd(n,2310)=gcd(n+2,2310)=1,并计算 ΩΩ 函数奇偶性,由此定义并统计 wheel 候选集、odd-odd 类以及经验比例

r(x)=π2(x)O(x),r(x)=O(x)π2(x),

其中 π2(x)π2(x) 表示不超过 xx 的孪生素数对数,O(x)O(x) 表示满足 wheel 条件且 Ω(n)Ω(n) 与 Ω(n+2)Ω(n+2) 均为奇数的整数对数量。

实验结果表明,在 x=108x=108 时,程序统计得到素数个数为 5,761,455,孪生素数对数为 440,312,最大相邻素数间隔为 220。基于当前实现,得到 r(108)≈0.3001158717r(108)≈0.3001158717。在长度为 107107 的窗口划分下,局部比例 rlocalrlocal 从约 0.3978 下降到约 0.2670,并在当前样本范围内与 C/(log⁡x)βC/(logx)β 形式具有较高拟合度 (R2≈0.99R2≈0.99)。此外,本文还将观测到的孪生素数对数与 Hardy–Littlewood 启发式预测作了有限范围比较 [4]。

需要强调的是,本文全部结果均属于有限范围内的计算实验与经验统计,不构成孪生素数猜想或 Hardy–Littlewood 渐近公式的严格证明。本文工作的主要意义在于给出一套可复现的分段筛实验框架,并为奇偶障碍背景下的经验量化分析提供基础数据。

关键词:孪生素数;分段筛法;wheel 筛选;奇偶障碍;Hardy–Littlewood 猜想;实验计算

1 引言

1.1 研究背景

孪生素数是指形如 (p,p+2)(p,p+2) 的素数对,例如 (3,5)(3,5)、(5,7)(5,7)、(11,13)(11,13)。孪生素数猜想断言这样的素数对有无穷多组,但目前尚无严格证明 [5]。

与这一猜想密切相关的经典启发式是 Hardy–Littlewood 第一猜想。其在孪生素数情形下预测:

π2(x)∼2C2x(log⁡x)2,π2(x)∼2C2(logx)2x,

其中 C2C2 为孪生素数常数。该公式长期以来被视为理解孪生素数分布的核心经验模型,但仍属于猜想而非定理 [4]。

在严格理论方面,张益唐于 2014 年证明存在无穷多对素数,其间隔小于 7×1077×107 [1];随后 Polymath8 项目将该上界降至 246 [2];Maynard 则提出了新的小间隔方法,证明了更一般的有界素数间隔结论 [3]。尽管这些成果极大推进了该方向,但“间隔 2”这一孪生素数情形仍未解决。

1.2 问题动机

筛法在素数问题中极其重要,但其在孪生素数问题上受到所谓“奇偶障碍”的限制。粗略地说,很多筛法只能有效排除带有小素因子的数,却难以在“真正素数”和“具有奇数个素因子的合数”之间作出足够精细的区分。这使得“看起来像素数的候选对”与“真正的孪生素数对”之间存在系统性混杂。

本文并不试图证明孪生素数猜想,而是提出一个可操作的实验问题:

在排除了前几个小素数模类干扰之后,那些“在奇偶层面看起来像素数”的候选对中,真正成为孪生素数对的比例有多大?该比例在有限范围内如何变化?

为此,本文引入 wheel 候选集与 odd-odd 类,并定义经验比例 r(x)r(x) 及局部比例 rlocalrlocal,对其进行数值统计与拟合分析。

1.3 本文贡献

本文的主要工作包括:

设计并实现了一套适用于 [0,108][0,108] 规模的分段筛实验程序,支持素数、孪生素数对和间隔统计;

引入 wheel 条件 gcd⁡(n,2310)=gcd⁡(n+2,2310)=1gcd(n,2310)=gcd(n+2,2310)=1,构造候选集并精确计算 ΩΩ 函数奇偶性;

定义并统计 odd-odd 类与经验比例 r(x)=π2(x)/O(x)r(x)=π2(x)/O(x);

对局部窗口中的 rlocalrlocal 进行趋势观察与 log-log 拟合;

将有限范围内的孪生素数对计数与 Hardy–Littlewood 启发式预测进行对照。

2 相关背景与符号定义

2.1 基本符号

本文使用如下记号:

π(x)π(x):不超过 xx 的素数个数

π2(x)π2(x):不超过 xx 的孪生素数对个数,即满足 p≤xpx 且 p,p+2p,p+2 同为素数的起点 pp 的个数

Ω(n)Ω(n):整数 nn 的素因子总数(计重数)

C2C2:孪生素数常数,其近似值为 0.6601618158468695

W=2⋅3⋅5⋅7⋅11=2310W=2⋅3⋅5⋅7⋅11=2310

2.2 wheel 候选集

定义 wheel 候选集

W(x)={n≤x:gcd⁡(n,W)=gcd⁡(n+2,W)=1}.W(x)={nx:gcd(n,W)=gcd(n+2,W)=1}.

该集合排除了所有因被 2,3,5,7,11 的模类阻断而不可能成为孪生素数起点的整数。

2.3 odd-odd 类

在 wheel 候选集中,进一步定义 odd-odd 类

O(x)={n≤x:gcd⁡(n,W)=gcd⁡(n+2,W)=1, Ω(n)≡Ω(n+2)≡1(mod2)}.O(x)={nx:gcd(n,W)=gcd(n+2,W)=1, Ω(n)≡Ω(n+2)≡1(mod2)}.

它表示:在排除小素数模类后,那些在 ΩΩ 奇偶性意义下同时“像素数”的整数对。

基于此定义经验比例

r(x)=π2(x)O(x).r(x)=O(x)π2(x).

这里必须强调,r(x)r(x) 是本文引入的经验统计量,不是现有理论中的标准常数,也不能直接推出孪生素数猜想的真伪。

3 算法设计与实现

3.1 总体流程

设扫描上界为 NN,分段长度为 BB。程序整体流程如下:

生成不超过 NN 的基础素数表;

将区间 [0,N][0,N] 分成若干长度约为 BB 的分段;

对每个分段执行分段筛,得到局部素性数组;

同时计算该分段内每个整数的 ΩΩ 奇偶性;

统计该分段中的素数数目、孪生素数对数、wheel 候选数、odd-odd 数;

进行窗口汇总与检查点输出;

最终输出全局统计量与局部窗口结果。

3.2 基础素数表

程序使用普通埃拉托斯特尼筛法生成所有不超过 NN 的素数。对应实现如下:

static vector<int> simple_sieve(int limit) {

   vector<bool> isPrime(limit + 1, true);

   isPrime[0] = isPrime[1] = false;

   for (int p = 2; 1LL * p * p <= limit; ++p) {

       if (isPrime[p]) {

           for (long long x = 1LL * p * p; x <= limit; x += p)

               isPrime[(int)x] = false;

       }

   }

   vector<int> primes;

   for (int i = 2; i <= limit; ++i)

       if (isPrime[i]) primes.push_back(i);

   return primes;

}

对于 N=108N=108,只需筛到 104104 量级,因此该步骤开销很小。

3.3 分段筛

对每个分段 [L,R][L,R],初始化 isPrimeSeg 为真,再使用基础素数逐个划去合数。伪代码如下:

算法 1:分段筛主过程

输入:分段区间 [L,R][L,R],基础素数表 basePrimes输出:分段素性数组 isPrimeSeg

初始化 isPrimeSeg[i]=true

对每个基础素数 pp

若 p2>Rp2>R,停止;

计算当前分段中第一个不小于 LL 的 pp 倍数;

从该位置起每隔 pp 标记一次合数

修正 0 和 1 的素性

返回 isPrimeSeg

该方法的优点是将内存需求从 O(N)O(N) 降到 O(B)O(B)。

3.4 ΩΩ 奇偶性的计算

对每个分段,程序不仅要知道某数是否为素数,还要计算 Ω(n) mod 2Ω(n)mod2。为此,对每个分段元素维护一个剩余值 rem[i],并用基础素数反复整除。每成功整除一次,就翻转该位置的奇偶位:

while (rem[idx] % p == 0) {

   rem[idx] /= p;

   parity_out[idx] ^= 1;

}

若最终 rem[i] > 1,则说明尚有一个大素因子未被剥离,再翻转一次。这样得到的 parity_out[i] 即为 Ω(L+i) mod 2Ω(L+i)mod2。

3.5 wheel 候选判定

wheel 条件采用

gcd⁡(n,2310)=1且gcd⁡(n+2,2310)=1gcd(n,2310)=1且gcd(n+2,2310)=1

即要求 nn 和 n+2n+2 均不被 2,3,5,7,11 整除。实现中直接使用欧几里得算法完成。

3.6 孪生素数对统计

对每个分段内的 xx,若 xx 与 x+2x+2 都为素数,则记作一对孪生素数。考虑到孪生对可能跨越分段边界,程序额外保存上一段末尾附近的素性状态,并做跨段补偿。该设计的目的是保证 π2(x)π2(x) 统计尽量完整。

需要说明的是,程序当前对孪生素数跨段补偿做了处理,但对 odd-odd 跨段补偿尚未完全对称实现。因此 π2(x)π2(x) 更接近精确统计,而 O(x)O(x) 和由此得到的 r(x)r(x) 应理解为“当前实现下的统计值”。

3.7 检查点与窗口设计

程序在若干检查点处输出累计统计量,例如 106,107,5×107,108106,107,5×107,108。

此外,将 [0,108][0,108] 划分为长度 107107 的 10 个窗口,分别统计每个窗口内:

孪生素数对数 twtw

odd-odd 对数 oooo

wheel 候选数 candcand

并定义

rlocal=twoo.rlocal=ootw.

3.8 复杂度分析

设 NN 为扫描上界,BB 为分段大小。则:

基础筛复杂度约为 O(Nlog⁡log⁡N)O(NloglogN);

分段筛主过程总复杂度近似为 O(Nlog⁡log⁡N)O(NloglogN);

ΩΩ 奇偶性计算额外增加一个与分段筛同量级的因子,但在 N=108N=108 范围内仍可在个人计算机上快速完成。

在本文实验参数下,程序运行时间约为 8 秒,说明该实现对 108108 范围是可行的。

4 实验设置

4.1 软硬件环境

硬件:Apple M1 芯片,16 GB 内存

操作系统:macOS

编译器:支持 C++17 的本地编译环境

实现方式:单线程

4.2 参数设置

上界:N=108N=108

分段长度:B=106B=106

wheel 参数:W=2310W=2310

窗口长度:107107

5 实验结果

5.1 全局统计量

程序最终输出如下核心结果:


统计量数值
素数总数 π(108)π(108)5,761,455
孪生素数对数 π2(108)π2(108)440,312
最大相邻素数间隔220
wheel 候选数5,844,137
odd-odd 数1,467,140
r(108)r(108)0.300115871696

其中,素数总数与孪生素数对数均属于当前程序实现下的直接累计值;r(108)r(108) 为基于 odd-odd 统计值得到的经验比例。

5.2 与 Hardy–Littlewood 启发式的有限范围比较

程序使用

2C2x(log⁡x)22C2(logx)2x

作为 Hardy–Littlewood 启发式预测值,并与观测值作对比。对于较大的检查点,得到如下结果:


x观测孪生数HL 预测obs/pred
10610614,8716,9172.15
10710764,04050,8221.26
5×1075×107243,139210,0651.16
108108440,312389,1071.13

可以看到,在当前实验范围内,obs/pred 随 xx 增大而明显下降,并向 1 靠近。这说明 Hardy–Littlewood 公式在有限范围内已经具有较好的数量级预测能力。

但必须指出,这种现象仅是有限范围的经验相容性,不能据此推出渐近公式已被证明。

5.3 局部窗口中的 rlocalrlocal

将 [0,108][0,108] 按 107107 划分为 10 个窗口,得到:


窗口区间孪生数odd-odd 数rlocalrlocal
(0,107](0,107]58,977148,2430.3978400329
(107,2×107](107,2×107]43,367132,2240.3279813045
(2×107,3×107](2×107,3×107]40,787131,7230.3096422037
(3×107,4×107](3×107,4×107]39,398132,2220.2979685680
(4×107,5×107](4×107,5×107]38,045131,6150.2890627968
(5×107,6×107](5×107,6×107]37,419132,0680.2833313142
(6×107,7×107](6×107,7×107]36,786131,6820.2793548093
(7×107,8×107](7×107,8×107]35,999131,8220.2730879519
(8×107,9×107](8×107,9×107]35,713131,2630.2720720995
(9×107,108](9×107,108]35,274132,0920.2670411531

在当前这 10 个窗口样本上,rlocalrlocal 呈现持续下降趋势。

应当强调,这一结论仅针对本文给定的窗口划分与上界范围,不应直接外推为普遍渐近规律

5.4 拟合结果

以窗口中点为 xx,对 rlocalrlocal 做如下 log-log 拟合:

r(x)≈C(log⁡x)β.r(x)≈(logx)βC.

程序得到:

β≈2.205039β≈2.205039

C≈162.554059C≈162.554059

R2≈0.990019R2≈0.990019

这说明在当前样本范围内,上述形式对 rlocalrlocal 的拟合度较高。

但这里的 ββ 只应被理解为有限区间经验拟合参数,不应被解释为严格的渐近指数。

5.5 相邻素数间隔

程序还统计了最大相邻素数间隔,得到:

最大间隔为 220

出现在 47,326,693 与 47,326,913 之间

这里需要特别说明:这只是 108108 内的有限区间统计量。它与 Zhang–Maynard–Polymath 所证明的“无穷多次出现的有界间隔上界”不是同一类结论,二者不应直接比较为“谁更强”或“谁更紧”。

6 讨论

6.1 关于 r(x)r(x) 的解释边界

本文的 r(x)r(x) 是一个自定义经验比例,用来衡量:

在通过 wheel 条件与 ΩΩ 奇偶性筛选后的候选对中,真正成为孪生素数对的比例。

从结果看,r(x)r(x) 在 108108 范围内约为 0.30,而局部窗口中的 rlocalrlocal 则呈下降趋势。这说明 odd-odd 类中包含大量“看起来像素数但并非孪生素数”的噪声对象。

这可以作为奇偶障碍的一种经验呈现,但不能被理解为对奇偶障碍的严格定量定理。

6.2 关于 Hardy–Littlewood 对照的解释

观测值与 Hardy–Littlewood 启发式预测的比值从较大范围看确实在向 1 靠近,但这只说明:在当前范围内,数据与该预测相容。

不能从这一点推出误差项已经被某个特定渐近式严格控制。

6.3 关于程序局限

当前实现存在以下边界条件:

检查点输出:若以“当前分段累计值”近似代替“严格截止到 xx 的累计值”,则小检查点的数据会偏粗;

跨段补偿:孪生素数跨段补偿已实现,但 odd-odd 跨段补偿尚未完全对称处理,因此 O(x)O(x) 与 r(x)r(x) 更适合写成“当前实现下的统计值”;

无符号整数下溢:若对 twin 起点差进行额外统计,必须确保跨段顺序严格单调,否则会发生无符号整数下溢。

这些问题并不影响本文作为“有限范围实验统计”的基本定位,但要求在论文中明确结果的解释边界。

7 结论

本文设计并实现了一套基于分段筛的孪生素数统计实验程序,对 [0,108][0,108] 范围内的素数、孪生素数对以及 wheel–odd-odd 经验统计量进行了系统计算。

主要结果包括:

统计得到 π(108)=5,761,455π(108)=5,761,455,π2(108)=440,312π2(108)=440,312;

计算得到 r(108)≈0.3001158717r(108)≈0.3001158717;

在长度为 107107 的窗口上,rlocalrlocal 从约 0.3978 下降到约 0.2670;

在当前样本范围内,rlocalrlocal 与 C/(log⁡x)βC/(logx)β 型函数具有较高拟合度;

观测到的孪生素数对数与 Hardy–Littlewood 启发式预测在 106106 至 108108 范围内逐步接近。

本文结果属于有限范围内的计算实验与经验分析,不构成孪生素数猜想或 Hardy–Littlewood 猜想的证明。其主要价值在于提供一套可复现的实验框架,以及一组可用于进一步讨论奇偶障碍与候选筛选效率的基础数据。

未来工作可沿以下方向继续展开:

将上界扩展到 109109 甚至更高;

完善 odd-odd 的跨段补偿,使 r(x)r(x) 统计更严格;

引入并行化提升吞吐性能;

对不同 wheel 参数下的经验比例进行比较分析。

参考文献

[1] Zhang, Y. Bounded gaps between primes. Annals of Mathematics, 2014, 179(3): 1121–1174.

[2] Polymath, D. H. J. Variants of the Selberg sieve, and bounded intervals containing many primes. Research in the Mathematical Sciences, 2014, 1:12. 文中给出无条件上界 H1≤246H1≤246.

[3] Maynard, J. Small gaps between primes. Annals of Mathematics, 2015, 181(1): 383–413.

[4] Hardy, G. H., Littlewood, J. E. Some problems of "Partitio Numerorum", III: On the expression of a number as a sum of primes. Acta Mathematica, 1923, 44: 1–70.

[5] Weisstein, E. W. Twin Prime Conjecture. Wolfram MathWorld. 关于孪生素数猜想与相关常数的综述。



程序调试结果


/Volumes/陈恩华/c/Demo2/cmake-build-debug/Demo2

============================================================

陈恩华·马虎算法(孪生素数证实程序)

============================================================

MAX_N=100000000 | BLOCK=1000000 | WHEEL=2310 | WINDOW_H=10000000

[中文调试] 扫描到 999999 | 素数=78498 | 孪生对=8169 | wheelOddOdd=15226 | r_global=0.5365164850 | maxGap=114 | 用时=0.08s

[检查点] x=100 | obs_twins=8169 | HL_pred≈6 | obs/pred=1312.139073 | C_est=866.224113 | errScaled(|obs-pred|/(x/log^3x))=7972.139 | r_global=0.5365164850

[检查点] x=1000 | obs_twins=8169 | HL_pred≈28 | obs/pred=295.231291 | C_est=194.900425 | errScaled(|obs-pred|/(x/log^3x))=2683.528 | r_global=0.5365164850

[检查点] x=10000 | obs_twins=8169 | HL_pred≈156 | obs/pred=52.485563 | C_est=34.648965 | errScaled(|obs-pred|/(x/log^3x))=626.097 | r_global=0.5365164850

[检查点] x=100000 | obs_twins=8169 | HL_pred≈996 | obs/pred=8.200869 | C_est=5.413901 | errScaled(|obs-pred|/(x/log^3x))=109.459 | r_global=0.5365164850

[检查点] x=1000000 | obs_twins=14871 | HL_pred≈6917 | obs/pred=2.149778 | C_est=1.419201 | errScaled(|obs-pred|/(x/log^3x))=20.973 | r_global=0.4902904619

[检查点] x=10000000 | obs_twins=64040 | HL_pred≈50822 | obs/pred=1.260081 | C_est=0.831857 | errScaled(|obs-pred|/(x/log^3x))=5.535 | r_global=0.3928231866

[中文调试] 扫描到 50999999 | 素数=3057494 | 孪生对=243139 | wheelOddOdd=749540 | r_global=0.3243842890 | maxGap=220 | 用时=4.11s

[检查点] x=50000000 | obs_twins=243139 | HL_pred≈210065 | obs/pred=1.157446 | C_est=0.764102 | errScaled(|obs-pred|/(x/log^3x))=3.685 | r_global=0.3243842890

[中文调试] 扫描到 100000000 | 素数=5761455 | 孪生对=440312 | wheelOddOdd=1467140 | r_global=0.3001158717 | maxGap=220 | 用时=7.92s

[检查点] x=100000000 | obs_twins=440312 | HL_pred≈389107 | obs/pred=1.131596 | C_est=0.747037 | errScaled(|obs-pred|/(x/log^3x))=3.201 | r_global=0.3001158717

--- (尾部窗口) r_local 序列:r_local = twin_in_window / oddodd_in_window ---

窗口#0 区间=(0,10000000] | twin=58977 | odd-odd=148243 | cand=584413 | r_local=0.3978400329

窗口#1 区间=(10000000,20000000] | twin=43367 | odd-odd=132224 | cand=525974 | r_local=0.3279813045

窗口#2 区间=(20000000,30000000] | twin=40787 | odd-odd=131723 | cand=525974 | r_local=0.3096422037

窗口#3 区间=(30000000,40000000] | twin=39398 | odd-odd=132222 | cand=525974 | r_local=0.2979685680

窗口#4 区间=(40000000,50000000] | twin=38045 | odd-odd=131615 | cand=525974 | r_local=0.2890627968

窗口#5 区间=(50000000,60000000] | twin=37419 | odd-odd=132068 | cand=525974 | r_local=0.2833313142

窗口#6 区间=(60000000,70000000] | twin=36786 | odd-odd=131682 | cand=525973 | r_local=0.2793548093

窗口#7 区间=(70000000,80000000] | twin=35999 | odd-odd=131822 | cand=525973 | r_local=0.2730879519

窗口#8 区间=(80000000,90000000] | twin=35713 | odd-odd=131263 | cand=525973 | r_local=0.2720720995

窗口#9 区间=(90000000,100000000] | twin=35274 | odd-odd=132092 | cand=525972 | r_local=0.2670411531

min r_local = 0.267041153136 发生在区间=(90000000,100000000]

--- (趋势拟合审计) r_local ≈ C / log(x)^beta ---

beta=2.205039 | C=162.554059 | R^2=0.990019

【解释】beta>0 -> 在当前范围更支持 r_local 缓慢趋零(不等价于孪生素数有限)。

============================================================

陈恩华·马虎算法 结果汇总(实证审计)

============================================================

Prime count = 5761455

Twin count (gap=2) = 440312

Max prime gap = 220 (47326693 -> 47326913)

gap<=246 total = 5761454 (应≈primeCount-1=5761454)

Wheel candidates = 5844137

Wheel odd-odd = 1467140

r_global = twins/odd-odd = 0.300115871696

--- twin-free gap(按孪生起点差)---

maxTwinStartGap = 18446744073708551878 (30999737 -> 29999999)

【关键提醒】

1) 从 maxGap<=220@1e8 不能推出 liminf(p_{n+1}-p_n)=2(证明层面需要击穿奇偶障碍)。

2) 从 C_est@1e8 不能推出 对所有x的误差不等式(需要解析误差项控制)。

3) 从 r≈0.3@1e8 不能推出 r(x)≥c>0;相反,本程序也会显示 r_local 更支持缓慢趋零。

4) 但 r(x)->0 并不等价于“孪生素数有限”(密度趋零仍可数量无穷)。

============================================================

用时=7.92s

进程已结束,退出代码为 0