素数分布的Python可视化分析

王德贵

素数分布是数论中研究素数性质的重要课题。素数，也称为质数，是指一个大于1的整数，除1和它本身外，不能被其他的正整数所整除。研究各种各样的素数分布状况，一直是数论中最重要和最有吸引力的中心问题之一。

最初，一些数学家想找到一个函数，可以表示任意一个素数，但这类尝试都失败了，后来数学家们又开始研究素数整体，于是便出现了很多相关理论和猜想。比如，欧拉公式、高斯公式、孪生素数猜想、梅林数、黎曼猜想等等，至今黎曼猜想还只是猜想，没有得到证明（图1）。

本文尝试用Python简单分析素数的分布规律，涉及中国电子学会等级考试Python六级“数据可视化”的相关内容。

关于素数个数的研究是素数分布中最重要的问题之一。在国内国际，素数的研究一直没有中断，前人也总结了很多经验和公式，也有很多猜想，不少理论属于高等数学知识，在此不多赘述。

笔者一直想从中小学生的角度来介绍素数分布规律及相关知识，这篇文章也早就想写，但确实有难度，一直没有一个完美的思路，今天整理出来，与大家共享，欢迎斧正。

程序设计思路是统计一定区间内的素数个数，找出与区间值的关系，然后通过Python可视化以图表的形式显示出来，也就是中学的函数图像，这样就能直观地看到素数分布的大致规律。

在[2，n]区间上（高中数学集合表示法——区间表示法）素数和孪生素数（孪生素数就是指相差2的素数对，例如3和5，5和7，11和13）的个数统计，直观来看，区间越大，个数也越多（图2）。

当区间变大时，素数个数也在增加，那么这个递增是线性的吗（图3）？

素数分布理论的中心定理，是关于素数个数问题的一个命题：设x≥1，以π（x）表示不超过x的素数的个数，当x→∞时，π（x）～Li（x）或π（x）～x/ln（x），ln（x）为对数，Li（x）为对数积分。

这里我们只研究π（x）与x的关系，即y=π（x）函数。

以π（x）表示不大于x的素数个数，例如，π（10）=4，π（100）=25，π（1000）=168。

编写程序，统计每一个确定范围内的素数个数，然后用可视化形式显示出来。我们以100个为一组，依次增加100个，共100组的情况，即2-100，2-200，2-300，……，2-10000这100组中素数个数（图4）。

其中要说明的是，本区间的第一个数与上一个区间的最后一个素数，如果是孪生素数，则计算在本区间上。因为统计方法是按区间统计的，以避免重复计算。

从输出结果看，随着区间的增大，素数个数（蓝色线）和孪生素数个数（绿色线）也随之增大，但我们看到它不是直线，即不是标准的线性关系。

从这些数据，我们可以求出线性回归方程（红色直线），与y=π（x） 图像很接近（图5）。

这里简单介绍线性回归方程的求法。直线的斜截式方程为：Y=kX+b，我们需要求出斜率k和截距b，这是初中数学知识点，但求回归直线方程是高中数学知识点。

对于一组离散数据对，（x1，y1），（x2，y2），……，（xn，yn），求得x，y的平均值px=（x1+x2+…+xn）/n，py=（y1+y2+…+yn）/n，则：

k=[（x1-px）*（y1-py）+（x2-px）*（y2-py）+……+（xn-px）*（yn-py）]/

[（x1-px）* （x1-px）+ （x2-px）* （x2-px）+……+（xn-px） （xn-px）]

b=py-k*px

为了大家能看明白，笔者使用了在程序中变量运算表达式的形式。（px，py）是线性回归方程的解，也是回归直线必然经过的点。表达式也可以书写成下列形式，意义与上述形式完全相同。

在同一图像上显示三条曲线，需要用第67行的设置，一行一列的第一部分，共用x轴，Y1、Y2、Y3分别表示素数个数、孪生素数个数和线性回归方程，同时设置Y2和Y3的颜色（74、75行）分别是绿色（g）和红色（r）。

前面讨论的是在一定区间内所有素数的个数，与区间值的关系，即y=π（x）中，小于x的所有素数个数与x之间的关系。如果我们把一个区间，分成n个自然数个数相同（m个整数）或不同的区间，再统计素数个数，还有规律吗？

我们可以理解为将10000个整数，100个为一组，依次统计100组的素數个数。即第一组是2-100，第二组是101-200……

程序如图6、图7：

与前一个程序相比，变化的地方是统计的数组su和ls无需求和，直接存储每个区段的素数个数，再一个变化是输出时，每段的范围不同了。

运行结果如下。其中蓝色线表示素数个数，绿色线表示孪生素数个数，红色直线表示素数个数的线性回归方程。可以看出变化趋势是数值变大，素数个数在波动中变少，而且变化越来越慢（图8）。

n=72，m=10000的情况下，变化曲线如图9。

n=1000，m=1000的情况下，素数个数变化趋势如下图（由于数据太大，生成时间太长，就只计算素数的个数了，前面程序稍作修改即可，这里不再给出程序），可以看出，区间值很大时，素数个数变化不大，有一定的波动（图10）。

如果用每个区间的素数个数的总体占比，1000×1000的图像如下。前面程序稍作修改，每个区间统计的个数除以总数即可，这里不再给出程序。可以看出，区间值很大时，占比很小（图11）。

这种分析，我们无法找到规律，区间值很大时，素数个数较少，无法进行对比，于是数学家们想到了另外的统计方法。

前面介绍的是均匀分布情况，在数值变大时，无法找到规律，于是我们可以用非均匀的分布来统计素数个数。

将自然数划分成36n（n+1）为界的一个个区间，显然每个区间的自然数个数是不一样的，但可以看到素数的分布规律：各区间的素数，以波浪形式渐渐增多，只有个别的区间比前面的少，造成这种现象的原因是，有些合数的因子多少和素数对区间无法整除之故。

对比其他人的相关资料，与笔者计算的稍有差别，经过验证，笔者用程序计算出来的数据，总体素数个数和孪生素数个数均没有错误。

先看程序（图12、图13）。

统计100个区间的素数个数，存入列表中。第39-45行，由于区间变化，用公式来表示，然后分别调用自定义函数，计算出素数个数存入列表。

第47-54行计算线性回归方程的参数，第56-60行设置三条曲线的函数关系，第62-65行，设定三条曲线及显示方式，再显示在屏幕上。

输出结果如图14，从这些数据基本可以看出，虽然有一些波动，但有一定的线性关系， 与y=π（x）曲线一样，都有线性相关。

为节省篇幅，生成数据不再输出（图14）。

每个程序都是一边写，一边测试，没有错误再继续写下一段。当计算数据很大的时候，需要测试更长时间，所以文中例子的数据量都不是太大。本文只是在一定范围内的简单分析，数据太大时，未完成测试，如果文中有纰漏的地方，还请各位同仁、朋友斧正。