『Python数値計算ノート』ではアフィリエイトプログラムを利用して商品を紹介しています。

双子素数

2019.10.21

双子素数

$1$ とその数自身以外に約数をもたない数のことを 素数 といいます：
 \[2,\ 3,\ 5,\ 7,\ 11,\ 13,\ 19,\ 23,\ 29,\ 31,\ 37,\ 41,\ ……\]
$2$ を除いたすべての素数は奇数です。隣り合う素数の差が $2$ となるような素数のペア $(p,\ p+2)$ を 双子素数 (twin prime) とよびます。$(2,\ 3)$ は、その差が $1$ なので双子素数には含まれません。双子素数を具体的に書き並べてみると
 \[(3,\ 5),\ (5,\ 7),\ (11,\ 13),\ (17,\ 19),\ (29,\ 31),\ (41,\ 43),\ \cdots\]
のようになります。

双子素数生成関数

Python で 100 未満の双子素数を生成してみましょう。
sympy.primerange() で 3 以上 100 未満の素数リストを生成して NumPy の配列に格納します。

# PYTHON_TWIN_PRIME

# In[1]

import numpy as np
from sympy import primerange

# 素数イテレータで素数リストを作成
pn = list(primerange(3, 100))

# 素数リストをNumPy配列に格納
pn = np.array(pn)

print(pn)
# [3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97]

位置をスライドさせた素数リストを生成します。

# In[2]

# 先頭～末尾の1つ手前までをスライス
x = pn[0:pn.shape[0]-1]

# 先頭の1つ後ろ～末尾までをスライス
y = pn[1:pn.shape[0]]

print("x:\n{}".format(x))
print("y:\n{}".format(y))

# x:
# [3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89]
# y:
# [5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97]

マスク操作によって、x と y の要素同士の差が 2 となるような配列を生成します。これは双子素数のペアの片側となります。

# In[3]

# xとyの要素同士の差が2の配列を生成
z=x[y-x==2]

print(z)
# [3  5 11 17 29 41 59 71]

この配列に 2 を加えると、ペアのもう片方が生成されます。

# In[4]

# zの各要素に2を加える
w = z + 2

print(w)
# [5  7 13 19 31 43 61 73]

numpy.vstack() で z と w を結合して転置すると、双子素数を格納した配列が完成します。

# In[5]

# zとwを結合して転置
tp = np.vstack((z, w)).T

print(tp)

'''
[[ 3  5]
 [ 5  7]
 [11 13]
 [17 19]
 [29 31]
 [41 43]
 [59 61]
 [71 73]]
'''

イテレータとして双子素数を格納することもできます。

# In[6]

# 双子素数のイテレータ
itp = iter(tp.tolist())

next(itp)
# [3, 5]

以上の処理を関数にまとめておきます。

# In[7]

# a以上b未満の双子素数と個数
def twin_prime(a, b, inv=False):
    pn = list(primerange(a, b))
    pn = np.array(pn)
    x = pn[0 : pn.shape[0]-1]
    y = pn[1 : pn.shape[0]]
    z = x[y - x == 2]
    w = z + 2
    tp = np.vstack((z, w)).T
    if inv == True:
        tp = 1/tp
    ct = tp.shape[0]
    return tp, ct

# twin_prime() をユニバーサル関数に変換
twin_prime = np.vectorize(twin_prime)

戻り値は双子素数と双子素数の個数です。オプション引数 inv を True に設定すると双子素数の逆数を返します。twin_prime() を使って 100 以上 200 未満の双子素数を探してみましょう。

# In[8]

# 100以上200未満の双子素数
tp = twin_prime(100, 200)

print(tp)
print(ct)

'''
[[101 103]
 [107 109]
 [137 139]
 [149 151]
 [179 181]
 [191 193]
 [197 199]]

7
'''

双子素数の予想

「双子素数は無数に存在するのか、あるいは有限個しかないのか」という問題 (双子素数の予想) は、いまだに数学上の未解決問題です。

双子素数の逆数和

双子素数の逆数和
 \[\left(\frac{1}{3}+\frac{1}{5}\right)+\left(\frac{1}{5}+\frac{1}{7}\right)+\left(\frac{1}{11}+\frac{1}{13}\right)+\cdots+\left(\frac{1}{p}+\frac{1}{p+2}\right)+\cdots\]
は収束することは証明されています。値はおよそ $1.9$ ぐらいで、ヴィーゴ・ブルンによって証明されたので ブルン定数 とよばれています。素数の逆数和は発散するので、双子素数は素数に比べて少ないことになりますが、無限級数が収束することもあるので、これだけで一概に双子素数が無限個あるとは断言できません。
 
双子素数の逆数和の収束は非常に遅いので、数値計算でブルン定数 $1.9$ 強に達するのは大変です。特に私の激安ノート PC だと本当にきついですが、ぎりぎり頑張って 1 億以下の双子素数の逆数和を計算してみました。

# In[9]

# 1億未満の双子素数と個数
tp = twin_prime(1, 1e8, inv=True)[0]

# 双子素数の逆数を全部足す
s = np.sum(tp)

print(s)
# 1.7588156210679744

… 全然足りていませんね。本気で計算してみたい人は、パワフルなデスクトップ PC を用意したほうがいいです。学生さんなら大学のコンピュータを使ってください。社会人のかたは会社の PC を …… もちろん冗談です (会社の PC を私用に使うのはやめましょう)。あと、twin_prime() の内部で使用している SymPy の primerange() があまり速くないような気がするので、NumPy でもっと効率の良い素数生成関数を作成したほうがいいかもしれません (良いプログラムが書けたらコメントで教えてください)。

双子素数の漸近公式

ハーディ (G.H.Hardy) とリトルウッド (J.E.Littlewood) は、$x$ 以下の双子素数の組の個数が漸近的に
 \[N(x)=2C\int_{2}^{x}\frac{x}{(\log x)^2}\tag{1}\]
で与えられると予想しています。$C$ は $3$ 以上の素数についての $1-1/(p-1)$ の総積
\[C=\prod_{p\gt 2}\left\{1-\frac{1}{(p-1)^2}\right\}\tag{2}\]
で計算される定数であり、ハーディ・リトルウッド定数 とよばれています。$C$ は mpmath.twinprime から取得できます。

# In[10]

from mpmath import mp, twinprime

# 精度を50桁に設定
mp.dps = 50; mp.pretty = True

# ハーディ・リトルウッド定数
c = twinprime

print(c)
# 0.66016181584686957392781211001455577843262336028473

漸近式 (1) に含まれる積分
 \[\int_{2}^{x}\frac{x}{(\log x)^2}\tag{3}\]
は普通の形で表すことのできない積分関数であり、一般に数値計算 (積分の近似計算) によって値を取得します。SciPy の integrate() を使えば、それほど難しくありません。

# In[11]

import numpy as np
from scipy import integrate

# 被積分関数を定義
f = lambda x : 1 / np.log(x)**2

# fを2からxまで積分
g = lambda x : integrate.quad(f, 2, x)[0]

# gをユニバーサル化
g = np.frompyfunc(g, 1, 1)

引数に配列を渡せるように、積分関数 g(x) をユニバーサル化しておきました。Matplotlib を使って、$100,000$ 未満の双子素数の個数の漸近曲線 (asymptotic curve) と真値 (true value) を確認しておきます。

# In[12]

import matplotlib.pyplot as plt

# プロットするデータを用意
t = np.arange(2, 100000, 1000)
y = 2*c*g(t)
z = twin_prime(1, t)[1]

# FigureとAxesの設定
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111)
ax.grid()
ax.set_title("Number of twin primes less than x", fontsize=16)
ax.set_xlabel("x", fontsize=14)
ax.set_ylabel("N(x)", fontsize=14)
ax.set_xlim(2, 100000)
ax.set_ylim(0, 1300)

# Axesに漸近曲線と真値をプロット
ax.plot(t, y, color = "blue", label="asymptotic curve")
ax.plot(t, z, color = "red", label=" true value")

# 凡例を表示
ax.legend()

双子素数 (twin prime) の個数とハーディ・リトルウッドの漸近線ハーディ・リトルウッドの漸近曲線を眺めると、双子素数は無限に存在するように思えますが、この漸近曲線自体があくまで予想なので、本当に正しい分布公式なのかどうかわかっていません。（参考：Wikipedia）

コメント