素數分布的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）。

每個程序都是一邊寫，一邊測試，沒有錯誤再繼續寫下一段。當計算數據很大的時候，需要測試更長時間，所以文中例子的數據量都不是太大。本文只是在一定范圍內的簡單分析，數據太大時，未完成測試，如果文中有紕漏的地方，還請各位同仁、朋友斧正。