[Python]素数をごにょごにょすると円周率になるやつを計算してみた

ホッタ

2022年7月18日 21:18

先週NHKでやってた 「素数」 - 笑わない数学 という番組でやってた式をPythonで計算してみました。

計算式

これが番組で出てた式。超有名な数学者オイラー先生が発見した法則です。

まとめると、以下のようになる。

この式はオイラー積という名前で、調べると以下の式３の方がよく出てきます。どちらも結果は一緒なので、式２の方を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ループを組んでいるので、求めたい素数の数が大きくなってくると指数関数的に実行時間が増えてしまいます。最初のコードは自分より小さい素数だけを割り算しているので、そこまで実行時間は増えません。今回は素数の個数を多くして計算したいので、最初のコードで実装しています。

オイラー積を計算する関数

素数のリストが作れるようになったので、式１の左辺を計算する関数を実装します。

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