[Python]素数をごにょごにょすると円周率になるやつを計算してみた
先週NHKでやってた 「素数」 - 笑わない数学 という番組でやってた式をPythonで計算してみました。
計算式
これが番組で出てた式。超有名な数学者オイラー先生が発見した法則です。
まとめると、以下のようになる。
この式はオイラー積という名前で、調べると以下の式3の方がよく出てきます。どちらも結果は一緒なので、式2の方をPythonで実装してみます。
コード
まず、素数のリストを作る必要があります。素数はSymPyというPythonのモジュールを使うと簡単に出せるのですが、今回はコードを書く練習のためSymPyモジュールは使わずに普通にコードを書いていきます。
SymPyについてはこちら。
素数リスト作成関数
def make_prime_number_list(number: int) -> list:
prime_number_list = []
n = 2
while len(prime_number_list) < number:
for p in prime_number_list:
if n % p == 0:
break
else:
prime_number_list.append(n)
n += 1
return(prime_number_list)
Pythonはforにelseをつけることができます。(知ってたけど実際に使ったのは初めて)
このelseは途中でbreakせず、forで回しているイテレータが最後まで到達したときのみ実行されます。フラグ変数を使わないですむというのがメリットですが、forやwhileのelseはリーダブルではないので使わないほうがよいと「EffectivePython」という本に書かれています。たしかに、このコードは初見だと直前のifに対応するelseじゃないの?って読み手が混乱するので使わないほうが良さそうです。
「EffectivePython」ではforのelseの代わりにヘルパー関数を使えと言っています。そうすると、以下のようなコードになります。
# 素数判定関数。引数が素数ならTrue、素数でなければFalseを返す
def is_prime_number(value: int) -> bool:
for i in range(2, value):
if value % i == 0:
return False
return True
def make_prime_number_list(number: int) -> list:
prime_number_list = []
n = 2
while len(prime_number_list) < number:
if is_prime_number(n):
prime_number_list.append(n)
n += 1
return(prime_number_list)
is_prime_numberという素数かどうかを判定する関数を追加しました。
こちらの方が読みやすいですね。ただ、このコードはis_prime_number関数内で引数より小さい全ての自然数との割り算をするforループを組んでいるので、求めたい素数の数が大きくなってくると指数関数的に実行時間が増えてしまいます。最初のコードは自分より小さい素数だけを割り算しているので、そこまで実行時間は増えません。今回は素数の個数を多くして計算したいので、最初のコードで実装しています。
オイラー積を計算する関数
素数のリストが作れるようになったので、式1の左辺を計算する関数を実装します。
def euler_product(number: int) -> float:
value = 1
for p in make_prime_number_list(number):
value = value * p**2 / (p**2 - 1)
return value
引数は素数の数です。先程作ったmake_prime_number_list関数で素数のリストを作成し、それをforで回して掛け算していきます。
計算結果
単体で実行
まずは素数10個でやってみます。素数10個だと以下のように2~29です。
euler_product関数に入れると、計算結果は約1.633になりました。
π^2/6を計算させると、約1.645くらいです。
いい感じなので、素数の数を増やしていきます。最終的には表にしたいので、ここからはpandasのデータフレームにしながら計算してみます。
素数の個数を振って傾向を見る
まずは計算する素数の個数を設定します。10個、100個、1000個、10000個のリストを設定して、euler_product_resultという名前のデータフレームを作成します。
次に、先程作ったeuler_product関数を使って、合計値を計算していきます。行毎に関数を実行するにはDataFrameのmapメソッドを使います。素数の個数をeuler_product関数の引数に代入した結果が新しい列(合計値)として追加されます。
最後に、π^2/6との比を右に追加してみます。pandasは行毎の演算結果を使って新しい列として追加できます。
個数が増えていく毎にπ^2/6に近づいていくのが分かります。素数10000個だと99.9999%となり、ほとんど一致ですね。
まとめ
というわけで、全く役に立たないコードですが、オイラー積の式の値を計算することができました。ついでにpandasの使い方も少し復習できました。
もっとスマートもしくは高速なコードが書けそうな気もしますが、沼が深そうなのでこれくらいにしておきます。
また面白いネタがあったらコード書いてみようと思います。
最後まで読んで頂き、ありがとうございました。