Pythonで素数定理｜Pythonで数学を学ぼう！第11回

はじめに

レポート

毎月開催中のZoomオンライン授業、桜井進の算数・数学教室

桜井進の魔法の算数教室
桜井進の数学浪漫紀行
桜井進のPython・UNIX教室（入門コース全３回）
桜井進のPython・UNIX・Math教室（応用コース）

桜井進のPython・UNIX教室を除く3クラスの7月共通が対数でした。タイトルは順に

小学生からはじめよう！対数
小数点「.」誕生402年記念　対数
Pythonで素数　素数定理

でした。

桜井進のPython・UNIX・Math教室（応用コース）の「Pythonで素数　素数定理」には対数がありませんが、素数定理が対数そのものです。

9歳の小学生から60歳過ぎまでの方に参加いただき、Pythonのコーディングとともに素数定理を自分の手と目と頭で検証しました。

次の結果は、自然数を10の10乗までの素数の個数と素数定理を比較したものです。

2つの比が1に近づくことが分かります。

素数定理

素数定理は今から250年以上前、1760年にオイラーによって発見されました。

オイラー以降もガウス、ルジャンドル、ディリクレ、チェビシェフといった数学者たちの手によって発見されています。

中でもガウスの発見物語はよく知られています。15歳のガウスは300万までの中にある素数リストを作成し、寝る前に素数リストを眺めることを習慣にしていたといいます。

ガウスが素数リストとともに眺めていたのが自然対数表です。神がかった計算力により、ランダムに出現する素数の中に隠れた法則を見つけだすことに成功しました。

x以下の自然数の中にある素数の個数をπ(x)と表すとき、π(x)はxを大きくしていくとx/log x に限りなく近づいていくというのが素数定理です。

ガウスから100年後の1896年に、プーサンとアダマールが独立に証明しました。

したがって、x以下の素数の個数π(x)と素数定理によるx/log x を計算すれば素数定理が検証できることになります。

log xはネイピア数e（=2.718281828459045 鮒一鉢二鉢一鉢二鉢至極美味しい）を底とする自然対数です。

プログラム「primetheorem1.py」

>>> 　# 素数定理
>>> 　# エラトステネスの篩（ふるい）により素数をリストアップ
>>> 　import math
>>> 　n =int(input("10のn乗以下の素数の個数π(10^n)を表示。n >> "))
>>> 　p = [i for i in range(10**n + 1)]
>>>
>>> 　# 素数リストp[]の生成
>>> 　for i in p[3:]:
>>> 　　　if p[i] % 2 == 0:
>>> 　　　　　p[i] = 0
>>>
>>> 　root = (10**n) ** 0.5
>>> 　for i in range(3, 10**n):
>>> 　　　if i > root:
>>> 　　　　　break
>>> 　　　if p[i] != 0:
>>> 　　　　　for j in range(i, 10**n + 1, 2):
>>> 　　　　　　　if i * j >= 10**n + 1 :
>>> 　　　　　　　　　break
>>> 　　　　　　　p[i * j] = 0
>>>
>>> 　plist = list(set(p))[2:]
>>> 　plistn = str(len(plist))
>>>
>>> 　print(' x=',10**n)
>>> 　print(' π(x)=',plistn)
>>>
>>> 　# 素数定理の計算
>>> 　print(f'x/log x= {(10**n)/math.log(10**n):.0f}')
>>> 　print(f'2つの比= {((10**n)/math.log(10**n))/int(plistn):.4f}')

ダウンロードはコチラから

第7回で紹介したエラトステネスの篩のコードを使って素数リストを生成します。するとlen()関数を用いて取得したリストの要素数が素数の個数となります。

大きなxに対して計算できるようにxには10のn乗を与えます。素数定理の計算は、自然対数をmath.log()関数を用いて行います。
>>> 　print(f'x/log x= {(10**n)/math.log(10**n):.0f}')

最後に、2つの比（素数定理の値/本当の値）を表示して素数定理と本当の値を比較します。

結果は次のように出力されます。10万までに素数は9592個ありますが、素数定理では8686個（本当の値の約90%）と算出されています。

プログラム「primetheorem2.py」

>>> 　# 素数定理
>>> 　# 代数ライブラリSymPyの関数primepi(n) n以下の素数の個数
>>> 　import math
>>> 　import numpy as np
>>> 　from sympy import primepi
>>>
>>> 　n = int(input("10のn乗以下の素数の個数π(10^n)を表示（高速）。n >> "))
>>>
>>> 　for i in range(1, n+1):
>>> 　　print(' x=',10**i)
>>> 　　print(' π(x)=',primepi(10**i))
>>>
>>> 　　# 素数定理の計算
>>> 　　print(f'x/log x= {(10**i)/math.log(10**i):.0f}')
>>>
>>> 　　print(f'2つの比= {((10**i)/math.log(10**i))/primepi(10**i):.4f}')
>>> print()

ダウンロードはコチラから

エラトステネスの篩プログラムで素数リストを生成するには限界があります。

そこで代数ライブラリSymPyの関数primepi()を用いてみます。

primepi(n)によりまさにπ(n)＝n以下の素数の個数を返します。本質的にprimetheorem1.pyと異なるのはここだけです。

結果は次のように出力されます。

nに5と与えれば、10、100、1000、10000、100000に対するπ(n) =primepi(n)と素数定理の値、および2つの比を出力します。

プログラム「primetheorem3.py」

>>> 　# 素数定理
>>> 　# 代数ライブラリSymPyの関数prime(n) n番目の素数
>>> 　import math
>>> 　from sympy import prime
>>>
>>> 　n = eval(input("n番目の素数p(n)を計算しよう。n >> "))
>>>
>>> 　p1 = int(n*math.log(n)) # n番目の素数の近似式
>>> 　p2 = prime(n) # n番目の素数
>>> 　e = ((p2-p1)/p2)*100 # 誤差
>>> 　print(f'p({n})≒', p1, '素数定理の近似式')
>>> 　print(f'p({n})=', p2, '本当の値')
>>> 　print(f'誤差{e:.2f}%')

ダウンロードはコチラから

代数ライブラリSymPyにはn番目の素数を与える関数prime(n)もあります。

そこでそれを用して、素数定理の式から導かれるn番目の素数p(n)の近似式との比較も行ってみます。