素數(shù)分布的Python可視化分析

王德貴

素數(shù)分布是數(shù)論中研究素數(shù)性質(zhì)的重要課題。素數(shù)，也稱為質(zhì)數(shù)，是指一個大于1的整數(shù)，除1和它本身外，不能被其他的正整數(shù)所整除。研究各種各樣的素數(shù)分布狀況，一直是數(shù)論中最重要和最有吸引力的中心問題之一。

最初，一些數(shù)學(xué)家想找到一個函數(shù)，可以表示任意一個素數(shù)，但這類嘗試都失敗了，后來數(shù)學(xué)家們又開始研究素數(shù)整體，于是便出現(xiàn)了很多相關(guān)理論和猜想。比如，歐拉公式、高斯公式、孿生素數(shù)猜想、梅林?jǐn)?shù)、黎曼猜想等等，至今黎曼猜想還只是猜想，沒有得到證明（圖1）。

本文嘗試用Python簡單分析素數(shù)的分布規(guī)律，涉及中國電子學(xué)會等級考試Python六級“數(shù)據(jù)可視化”的相關(guān)內(nèi)容。

關(guān)于素數(shù)個數(shù)的研究是素數(shù)分布中最重要的問題之一。在國內(nèi)國際，素數(shù)的研究一直沒有中斷，前人也總結(jié)了很多經(jīng)驗和公式，也有很多猜想，不少理論屬于高等數(shù)學(xué)知識，在此不多贅述。

筆者一直想從中小學(xué)生的角度來介紹素數(shù)分布規(guī)律及相關(guān)知識，這篇文章也早就想寫，但確實有難度，一直沒有一個完美的思路，今天整理出來，與大家共享，歡迎斧正。

程序設(shè)計思路是統(tǒng)計一定區(qū)間內(nèi)的素數(shù)個數(shù)，找出與區(qū)間值的關(guān)系，然后通過Python可視化以圖表的形式顯示出來，也就是中學(xué)的函數(shù)圖像，這樣就能直觀地看到素數(shù)分布的大致規(guī)律。

在[2，n]區(qū)間上（高中數(shù)學(xué)集合表示法——區(qū)間表示法）素數(shù)和孿生素數(shù)（孿生素數(shù)就是指相差2的素數(shù)對，例如3和5，5和7，11和13）的個數(shù)統(tǒng)計，直觀來看，區(qū)間越大，個數(shù)也越多（圖2）。

當(dāng)區(qū)間變大時，素數(shù)個數(shù)也在增加，那么這個遞增是線性的嗎（圖3）？

素數(shù)分布理論的中心定理，是關(guān)于素數(shù)個數(shù)問題的一個命題：設(shè)x≥1，以π（x）表示不超過x的素數(shù)的個數(shù)，當(dāng)x→∞時，π（x）～Li（x）或π（x）～x/ln（x），ln（x）為對數(shù)，Li（x）為對數(shù)積分。

這里我們只研究π（x）與x的關(guān)系，即y=π（x）函數(shù)。

以π（x）表示不大于x的素數(shù)個數(shù)，例如，π（10）=4，π（100）=25，π（1000）=168。

編寫程序，統(tǒng)計每一個確定范圍內(nèi)的素數(shù)個數(shù)，然后用可視化形式顯示出來。我們以100個為一組，依次增加100個，共100組的情況，即2-100，2-200，2-300，……，2-10000這100組中素數(shù)個數(shù)（圖4）。

其中要說明的是，本區(qū)間的第一個數(shù)與上一個區(qū)間的最后一個素數(shù)，如果是孿生素數(shù)，則計算在本區(qū)間上。因為統(tǒng)計方法是按區(qū)間統(tǒng)計的，以避免重復(fù)計算。

從輸出結(jié)果看，隨著區(qū)間的增大，素數(shù)個數(shù)（藍(lán)色線）和孿生素數(shù)個數(shù)（綠色線）也隨之增大，但我們看到它不是直線，即不是標(biāo)準(zhǔn)的線性關(guān)系。

從這些數(shù)據(jù)，我們可以求出線性回歸方程（紅色直線），與y=π（x） 圖像很接近（圖5）。

這里簡單介紹線性回歸方程的求法。直線的斜截式方程為：Y=kX+b，我們需要求出斜率k和截距b，這是初中數(shù)學(xué)知識點，但求回歸直線方程是高中數(shù)學(xué)知識點。

對于一組離散數(shù)據(jù)對，（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

為了大家能看明白，筆者使用了在程序中變量運算表達(dá)式的形式。（px，py）是線性回歸方程的解，也是回歸直線必然經(jīng)過的點。表達(dá)式也可以書寫成下列形式，意義與上述形式完全相同。

在同一圖像上顯示三條曲線，需要用第67行的設(shè)置，一行一列的第一部分，共用x軸，Y1、Y2、Y3分別表示素數(shù)個數(shù)、孿生素數(shù)個數(shù)和線性回歸方程，同時設(shè)置Y2和Y3的顏色（74、75行）分別是綠色（g）和紅色（r）。

前面討論的是在一定區(qū)間內(nèi)所有素數(shù)的個數(shù)，與區(qū)間值的關(guān)系，即y=π（x）中，小于x的所有素數(shù)個數(shù)與x之間的關(guān)系。如果我們把一個區(qū)間，分成n個自然數(shù)個數(shù)相同（m個整數(shù)）或不同的區(qū)間，再統(tǒng)計素數(shù)個數(shù)，還有規(guī)律嗎？

我們可以理解為將10000個整數(shù)，100個為一組，依次統(tǒng)計100組的素數(shù)個數(shù)。即第一組是2-100，第二組是101-200……

程序如圖6、圖7：

與前一個程序相比，變化的地方是統(tǒng)計的數(shù)組su和ls無需求和，直接存儲每個區(qū)段的素數(shù)個數(shù)，再一個變化是輸出時，每段的范圍不同了。

運行結(jié)果如下。其中藍(lán)色線表示素數(shù)個數(shù)，綠色線表示孿生素數(shù)個數(shù)，紅色直線表示素數(shù)個數(shù)的線性回歸方程。可以看出變化趨勢是數(shù)值變大，素數(shù)個數(shù)在波動中變少，而且變化越來越慢（圖8）。

n=72，m=10000的情況下，變化曲線如圖9。

n=1000，m=1000的情況下，素數(shù)個數(shù)變化趨勢如下圖（由于數(shù)據(jù)太大，生成時間太長，就只計算素數(shù)的個數(shù)了，前面程序稍作修改即可，這里不再給出程序），可以看出，區(qū)間值很大時，素數(shù)個數(shù)變化不大，有一定的波動（圖10）。

如果用每個區(qū)間的素數(shù)個數(shù)的總體占比，1000×1000的圖像如下。前面程序稍作修改，每個區(qū)間統(tǒng)計的個數(shù)除以總數(shù)即可，這里不再給出程序。可以看出，區(qū)間值很大時，占比很小（圖11）。

這種分析，我們無法找到規(guī)律，區(qū)間值很大時，素數(shù)個數(shù)較少，無法進(jìn)行對比，于是數(shù)學(xué)家們想到了另外的統(tǒng)計方法。

前面介紹的是均勻分布情況，在數(shù)值變大時，無法找到規(guī)律，于是我們可以用非均勻的分布來統(tǒng)計素數(shù)個數(shù)。

將自然數(shù)劃分成36n（n+1）為界的一個個區(qū)間，顯然每個區(qū)間的自然數(shù)個數(shù)是不一樣的，但可以看到素數(shù)的分布規(guī)律：各區(qū)間的素數(shù)，以波浪形式漸漸增多，只有個別的區(qū)間比前面的少，造成這種現(xiàn)象的原因是，有些合數(shù)的因子多少和素數(shù)對區(qū)間無法整除之故。

對比其他人的相關(guān)資料，與筆者計算的稍有差別，經(jīng)過驗證，筆者用程序計算出來的數(shù)據(jù)，總體素數(shù)個數(shù)和孿生素數(shù)個數(shù)均沒有錯誤。

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

統(tǒng)計100個區(qū)間的素數(shù)個數(shù)，存入列表中。第39-45行，由于區(qū)間變化，用公式來表示，然后分別調(diào)用自定義函數(shù)，計算出素數(shù)個數(shù)存入列表。

第47-54行計算線性回歸方程的參數(shù)，第56-60行設(shè)置三條曲線的函數(shù)關(guān)系，第62-65行，設(shè)定三條曲線及顯示方式，再顯示在屏幕上。

輸出結(jié)果如圖14，從這些數(shù)據(jù)基本可以看出，雖然有一些波動，但有一定的線性關(guān)系， 與y=π（x）曲線一樣，都有線性相關(guān)。

為節(jié)省篇幅，生成數(shù)據(jù)不再輸出（圖14）。

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