Pythonで双子素数を計算する
差が2である2つの素数の組が構成する素数を双子素数といいます。例えば(3,5),(5,7),(11,13),(17,19)・・・が双子素数になります。双子というくらいなので本当はとなり同士の数の組にしたいとことですが、素数は2,3,5,7・・となり、2以外の偶数は素数にならず (2, 3)以外の組は考えられないので、差が2の組としたようです。
双子素数の生成とクレメントの双子素数判定法
双子素数の生成
双子素数を求めるため、まず素数のリストを作成し、そこから双子素数を見つけ出すプログラムを作成します。ここでは、ひとまず10,000までの素数を扱います。
10000までの素数を計算する
- import sympy
- prime_l=list(sympy.primerange(3,10000))
- len(prime_l)
1228
ここでは、いちいち表示せず、len関数で個数だけを調べることにします。
双子素数をリストにする
- prev=2
- twin_l=[]
- for i in prime_l:
- if i == prev + 2:
- twin_l.append((prev,i))
- prev=i
- twin_l
[(3, 5), (5, 7), (11, 13), (17, 19), (29, 31), (41, 43), (59, 61), (71, 73), （中略） (9719, 9721), (9767, 9769), (9857, 9859), (9929, 9931)]
全部で205組になります。素数が1228に対し、かなりの数になります。
クレメントの双子素数判定法
双子素数を生成する方法にクレメントの双子素数判定法があります。
クレメントの双子素数判定法
さっそく、プログラムにしてみます。
- fact = 1
- twin_c = []
- for n in range(2, 10000):
- fact *= n-1
- num = (fact+1)*4+n
- if num % (n*(n+2)) == 0:
- twin_c.append((n, n+2))
- twin_c == twin_l
True
前項で作成したtwin_lと比べると、ぴったりと一致します。
三つ子素数
三つ子素数は素数が3個の組で、素数pに対して(p,p+ 2,p+ 6)または(p,p+ 4,p+ 6)となるものをいいます。 (p,p+ 2,p+ 6)のように2つ飛びにならないのは、偶数を2つおきに3つ並べると、どれか1つは3の倍数となり素数とならないためです。
三つ子素数をリストにする
三つ子素数は、最大値と最小値の差が6であることから、prime_lからそのような組を抜き出します。
三つ子素数の生成
- prev1=3
- prev2=2
- tri_l=[]
- for i in prime_l:
- if i == prev2 + 6:
- tri_l.append((prev2,prev1,i))
- prev2=prev1
- prev1=i
- tri_l
[(5, 7, 11), (7, 11, 13), (11, 13, 17), (13, 17, 19), (17, 19, 23), (37, 41, 43), (41, 43, 47), (67, 71, 73), （中略） (9337, 9341, 9343), (9431, 9433, 9437), (9433, 9437, 9439), (9461, 9463, 9467)]
全部で112組になります。双子素数の半分くらいに減ります。
三つ子素数の分類
三つ子素数は、(p,p+ 2,p+ 6)型、(p,p+ 4,p+ 6)型の2つがあります。そこで、前節で求めたtri_lを分割してみます。
- tri_11=[]
- tri_13=[]
- for i in tri_l:
- if i[1]==i[0]+2:
- tri_11.append(i)
- else:
- tri_13.append(i)
- print(tri_11,len(tri_11))
- print(tri_13,len(tri_13))
[(5, 7, 11), (11, 13, 17), (17, 19, 23), (41, 43, 47), (101, 103, 107), (107, 109, 113), (191, 193, 197), (227, 229, 233), (311, 313, 317), (347, 349, 353), (461, 463, 467), (641, 643, 647), (821, 823, 827), (857, 859, 863), (881, 883, 887), 55 [(7, 11, 13), (13, 17, 19), (37, 41, 43), (67, 71, 73), (97, 101, 103), (103, 107, 109), (193, 197, 199), (223, 227, 229), (277, 281, 283), (307, 311, 313), (457, 461, 463), (613, 617, 619), (823, 827, 829), (853, 857, 859), (877, 881, 883) 57
(p,p+ 2,p+ 6)型が55個、(p,p+ 4,p+ 6)型が57個となり、概ねバランスしています。
ブルン定数の計算
双子素数(3, 5), (5, 7), (11, 13)・・・を3, 5, 5, 7,11, 13というように分解し、その逆数の合計をします。数式にすると次の通りです。
ブルン定数を計算するため、SymPyモジュールを使い素数の一覧を作成し、そこから双子素数を抜き出した上で、一次元のリストに変換します。
双子素数を求め、一次元のリストにする
求める素数の上限は10のmax乗とします。コンピュータの性能にもよりますがここでは10としておきます。
双子素数を求め、リストに分解する
- prime_l
- sympyモジュールで計算した素数を格納
- twin_l
- 素数のうち、前の素数との差が2のものについて双子素数の一次元のリストとして格納
- digits
- 10のべき乗のコントロールをする
- brun_l
- 素数が10のべき乗ごとに、計算されたブルン定数を格納
- import sympy
- max=4
- brun_l=[]
- digits=10
- prime_l=list(sympy.primerange(3,10**(max)))
- prev=2
- twin_l=[]
- for i in prime_l:
- if i == prev + 2:
- twin_l.append(prev)
- twin_l.append(i)
- if i >= digits:
- brun_l.append(prev)
- digits*=10
- prev=i
2. 求める双子素数の上限を指定します。ここでは、$10^6$まで求めることとし、max=6とします。
5. SymPyモジュールで、3（双子素数の最小値）から$10^max$までの素数をリストprime_lに格納します。
6. prevは仮に2を代入しておき、15.で処理した素数の1つ前の値を代入しておきます。
9.10. 読込んだ素数（i）がprev（1つ前の素数）との差が2の時に、リストtwin_lにprevおよびiを追加します。
14. digitsが10のn乗を超えたときに、その時の素数の値をリストburn_lに追加します。digitsが切りの良い10のべき乗を超えたら10倍しておきます。
ブルン定数を計算する
- brun=0
- for i in twin_l:
- if i in brun_l:
- print(i,sigma)
- brun+=1/i
- print(i,sigma)
11 0.8761904761904762 101 1.330990365719087 1019 1.5180324635595912 9931 1.6168935574321999
2. 前節で作成した双子素数のリストtwin_lを処理します。
5. 変数brunでブルン定数を計算するため、1/双子素数を足しこんでいきます。
3. brun_lに双子素数が$10^n$になった時にはそのときのbrun定数を出力します。
今回は$10^4$まで計算しましたが、ブルン定数は次の通りになり、1.902160583104に収束するといわれています。
$10^2$ 1.330990365719
$10^4$ 1.616893557432
$10^6$ 1.710776930804
$10^8$ 1.758815621067
$10^{10}$ 1.787478502719
$10^{12}$ 1.806592419175.
$10^{14}$ 1.820244968130
$10^{15}$ 1.825706013240
$10^{16}$ 1.830484424658
1996年、Thomas R. Nicelyという学者が $10^{14}以下の双子素数までの部分和を計算しました。このとき、CPUに関するバグであるPentium FDIV バグを発見しました。なぜ、ブルン定数の計算で見つかったかは想像もつきませんが、異常にゆっくりと収束するのが原因ではないかと邪推しています。