プログラミング最高！（数学夏祭り問4）

数学夏祭り

出典

攻略法

プログラムを組めばOK

ソースコード

from fractions import Fraction
MAX_N=100
isPrime=[1 for i in range(MAX_N+1)]
isPrime[0]=0;isPrime[1]=0
for i in range(MAX_N+1):
    if isPrime[i]==0:
        continue
    for j in range(2*i,MAX_N+1,i):
        isPrime[j]=0
nowmin=Fraction(1,1)
argmin=[]
for n in range(6,80+1,2):
    bunsi=0
    bunbo=0
    for p in range(1,n+1,2):
        bunbo+=1
        if isPrime[p] and isPrime[n-p]:
            bunsi+=1
    tmp=Fraction(bunsi,bunbo)
    if nowmin>tmp:
        argmin=[n]
        nowmin=tmp
    elif nowmin==tmp:
        argmin.append(n)
print(argmin,nowmin)

このプログラムの出力は

[68] 2/17

であるため、
$$n=68,P(n)=\frac{2}{17}$$
が答えである。

余談

この問題はゴールドバッハ予想が絡んでいるように見える。
$$f(n)=P(n)\cdot \frac{n}{2}$$
とおくと、素数の和で表すための場合の数になる。
Wikipediaによると、多少ガバガバな計算をすることでこれは{\frac{n}{\log^2{n}}}のオーダーっぽいということが分かるとのことである。
真偽はともかく、{n\leq 2^{20},n\leq 2^{25}}で実験すると以下のようになる。

f:id:shakayami:20200904002257p:plain
n≦2^20
f:id:shakayami:20200904002320p:plain
n≦2^25

これは
$$f(n)\cdot \frac{\log^2{n}}{n}$$
についての値をグラフにしている。つまり上記の「ガバガバな計算」が正しければ定数オーダーになるはずである。上限と下限が一定であるように見えるため、矛盾がないように思えるかもしれない。

ゴールドバッハ予想は未解決問題であるため、あまり変なことを言って間違いを書くことをやめたいので、ここで少々別の話題に変えてみる。

上記の図はpythonのプログラム(matplotlib)を使って書かれているが、プログラムを少々工夫しないと実行時間が年単位で掛かってしまう。この記事が今出せているということは、プログラムを工夫しているのである。

素数の和で表すための場合の数は以下のように書ける。

$$Q(x)=(x^2+x^3+x^5+x^7+x^{11}+x^{13}+\cdots)^2$$
としたときの{Q(x)}の{x^n}の係数

結局これは畳み込み積で書ける。するとFFT（高速フーリエ変換）を使うことで計算することができる。{1\leq k\leq 2^n}の範囲内での場合の数のリストを作成する場合、{n\cdot 2^n}程度の計算量で実行が終わる。一方、愚直にシミュレーションをすると{2^{2n}}程度の計算量がかかる。
問題は素数のリストを書くことだが、それはエラトステネスの篩を適用することで、{2^n\log{n}}程度の計算量となる。

最終的に図を出力するためのプログラムは以下のようになる。

import matplotlib.pyplot as plt
from fractions import Fraction
import math
import numpy as np
#MAX_N=1<<25
#にすると2^25の場合になる（実行時間がかかる）
MAX_N=1<<20
isPrime=[1 for i in range(2*MAX_N)]
isPrime[0]=0;isPrime[1]=0
for i in range(MAX_N):
    if isPrime[i]==0:
        continue
    for j in range(2*i,MAX_N,i):
        isPrime[j]=0
for i in range(MAX_N,2*MAX_N):
    isPrime[i]=0
PrimeSum=[int(i+0.5) for i in np.real(np.fft.ifft(np.fft.fft(np.array(isPrime))*np.fft.fft(np.array(isPrime))))]
result=np.array([PrimeSum[i]*math.log(i)**2/i for i in range(6,MAX_N+1,2)])
X=np.array([i for i in range(6,MAX_N+1,2)])
plt.plot(X,result,".",markersize=1)
plt.ylim(0,10)
plt.grid()
plt.show()

ちなみに{2^{25}}の場合の図を書く際、プログラムの実行時間は約2分であった。