前回のあらすじ
前回までに、素数のデータベースの取り込み、エラトステネスの篩を使った素数生成、両者の比較検証、などが済み、いよいよ問題に直接チャレンジできるところまできた。今回でけりを付ける。
コードを分ける
エラトステネスの篩の部分(q6_eratos.py)、関数の定義(q6_func.py)、そして数学の問題を解く部分(q6_main.py)、の3つにプログラムファイルを分けることにする。
まずは数学の問題を解く部分(q6_main.py)。
import q6_eratos as q6e import q6_func as q6f Nmax = 21000000 prime_table = [ True for i in range(Nmax)] prime_table[0:2] = [False, False, True] prime_table = q6e.eratos(prime_table,Nmax) # Negative-integer check IMAX=101 print("\n f max",q6f.f(-IMAX)) for i in range(IMAX): im = -i fm = q6f.f(im) if fm > 0: print(im,fm,prime_table[fm]) # Positive-integer check IMAX=271 print("\n f max",q6f.f(IMAX)) for i in range(IMAX): fi = q6f.f(i) if prime_table[fi] : print(i,fi,prime_table[fi])
「篩」の表の中に含む整数は21,000,000まで(2100万個）とした。このときの最大の素数は20,999,999と算出されたが、これがちゃんと素数であることはThe PrimePagesによって個別に確認した。この設定に対して、篩が見つけた素数の数は1,329,943個。表の中に含まれる整数の1/10程度のオーダーになっているようだ（これは色々と試しに計算してみて、経験的に会得した推測である）。
負の整数の場合は$n=-101$までを試しているが、実は問題で与えられた関数は$x\rightarrow\infty$で$f(x)\rightarrow x^3$のように振る舞うので、関数$f(x)$自体が素数、つまり正数となるのは限られた$x<0$の値だけであることが（やってみると）わかる。数値計算の結果がその様子を抉り出すようにプログラミングされている。
正の整数については$n=271$まで試している。これは$f(n)<2100万$となるように決めた値である（もちろん、試行錯誤で）。
次は、エラトステネスの篩の部分のコード（q6_eratos.py)である。
def eratos(prime_table, Nmax): Pmax = 1 MaxPrime = 2 PrimeMax = 2 while MaxPrime < Nmax: n = MaxPrime while n <= Nmax/MaxPrime: prime_table[MaxPrime * n] = False n += 1 while True: MaxPrime += 1 if MaxPrime > Nmax -1: break if prime_table[MaxPrime] == True: Pmax += 1 PrimeMax = MaxPrime break print("Number of primes found: ",Pmax) print("Largest prime found: ",PrimeMax) return prime_table
global変数とlocal変数の違いをきちんと理解しておかないと、意外に関数というのは使うのが難しい。今回は、素数表とその最大整数値をglobal変数とし、外から（引数として）読み取ることにした。返すのは、関数内で処理した素数表であるが、これ自体はlocalなので、q6_main.pyの中で「値渡し」している。関数内ではglobal変数を参照することはできるが、「いじる（変化させる）」ことはできない。
最後に、問題文で与えられた関数のコード(q6_func.py)である。
def f(x): return x*(x**2 + 10*x + 20)
多項式の計算は、本来はホーナーのアルゴリズムというのを使ってやるべきなのだが、計算するのは因数分解できる3次式なので、直感的にわかりやすい表式で済ませた。
これで以上である。実行してみると次のようになった。
% python3 q6_main.py Number of primes found: 1329943 Largest prime found: 20999999 f max 1134331 -3 3 True -4 16 False -5 25 False -6 24 False -7 7 True f max 20642341 1 31 True
前半部分が負の整数$n<0$を試した結果である。抜き出された結果は、関数の値が正値をとる場合のみである。予想通り、$n=-3,-4,-5,-6,-7$の５つの場合のみに限定されている。$n\rightarrow -\infty$で$f(n)\rightarrow n^3 < 0$であるから、どんなに大きな$n$を試しても無駄なことになる。計算では$n=-101$までを試しているが、$n=-10$まで調べればもう十分である。ちなみに、$n=-101$のとき、$f(-101)=-930,311$であった。
$n<0$に対し関数値$f(n)$が正となる、これら５つの場合のうち、$f(n)$が素数となるのは$n=-3,-7$の2つだけであった。Trueというのは、エラトステネスの篩で見つかった素数の表の中に含まれている素数である、という意味である。
後半部分は$n>0$の場合であり、こちらは$n=1$の場合のみが適合した。試した整数の最大値は$n=271$で、このとき$f(271)=20,642,341$であった。これは見つけた最大の素数20,999,999に肉薄する値である。$n$のオーダーが100程度であっても、３次関数の値自体は非常に大きくなってしまうのは、数値計算的には頭の痛い内容である。
とはいえ、この計算に要した時間はわずか5.8秒弱である。もうちょっと頑張って、$n=1000$くらいまで試してみてもいいかもしれない。
ということで、問題6の前半はおしまいである。後半に関しては丁寧に分析していくだけであろうから、数値計算の観点からはあまり面白そうに感じられない。来年の受験シーズンまではまだ時間あるので、ゆっくりやっていくことにする（つまり、気が向いたら戻ってくるが、しばらくはやらない予定、ということ）。
