完全数を見つけるプログラムをもっと速く！メルセンヌ素数で改善してみよう

前回の記事「完全数をPythonで見つけよう！数論の不思議をプログラムで体感!」では、完全数の定義に従って、PythonとRustを使ってコードを実装し、1〜10000の範囲で完全数を探しました。

しかし、定義に忠実に約数をすべてリストアップしていく方法では、処理効率が非常に悪くなってしまいます。
そこで今回は、「完全数入門! 数の面白シリーズ!」で紹介した完全数の性質（メルセンヌ素数との関係）を活用して、より効率的に完全数を求める方法を試してみたいと思います。

Contents

完全数と完全数の性質（復習）
メルセンヌ数を使って完全数を探そう!
- （おまけ）Rustでもメルセンヌ数を使って完全数を探そう!
最後に

完全数と完全数の性質（復習）

まず、完全数の定義を思い出します。

完全数 (Perfect Number)

自身以外の正の約数の和が、ちょうど自分自身に等しい数のこと

完全数の例として

$6 = (1 + 2 + 3)$
$28 = (1 + 2 + 4 + 7 + 14)$
$496$

8128[katex]</li> <!-- /wp:list-item --></ul> <!-- /wp:list --> <!-- wp:paragraph --> <p>などがあります。 <strong>メルセンヌ<strong>素数</strong></strong>([katex]2^n-1

で表される素数)と2のべき乗の積は完全数になります。

命題1

$2^n-1$ が素数ならば、 $2^{n-1}(2^n-1)$ は完全数である。

この事実を使って完全数を探していきます。

メルセンヌ数を使って完全数を探そう!

まずはPythonで実装します。

# 素数判定を行う関数
def is_prime(n):
    # 2未満の数は素数ではない
    if n < 2:
        return False
    # 2から√nまでの数で割り切れるか調べる
    for i in range(2,int(n**0.5) + 1):
        if n % i == 0:
            # 割り切れたら素数ではない
            return False
    # 割り切れる数がなければ素数
    return True

# 偶数完全数を生成する関数
def generate_even_perfect_numbers(n):
    perfect_numbers = []  # 完全数を格納するリスト
    p = 2  # 素数指数の初期値 
    while True:
        # メルセンヌ数を計算
        mersenne_number = 2**p - 1
        # メルセンヌ素数かどうか調べる
        if is_prime(mersenne_number):
            # メルセンヌ素数であれば完全数を算出
            perfect_number = 2**(p-1) * mersenne_number
            # 調べる数の最大数まで調べたら終了
            if perfect_number > n:
                break
            # そうでなければ作成した完全数をリストに追加
            perfect_numbers.append(perfect_number)
        p += 1
    return perfect_numbers

if __name__ == '__main__':
    # 10000以下の偶数完全数を表示
    for pn in generate_even_perfect_numbers(10000):
        print(f"{pn}は完全数です")

プログラムでやっていることは、generate_even_perfect_numbers関数に与えられた数までメルセンヌ数を計算して、その数が素数であれば $2^{n-1}(2^n-1)$ を計算してリストに追加して、最後にそのリストの数を表示するということです。

実行したら分かるのですが、前回のプログラムに比べてかなり処理速度が向上しています。
試しにgenerate_even_perfect_numbers関数に与える数を1000000000000000にしてみましたが、かかった時間は0.003秒です！
以前のプログラムで100000までの完全数を調べたら、176.91041秒かかりました。
(1000000000000000で調べようと思いましたが、1時間以上まっても終わらなかったので断念しました。)

（おまけ）Rustでもメルセンヌ数を使って完全数を探そう!

上で作成したコードのRust版を作成します。

use std::time::Instant; // 時間計測のためのモジュールをインポート

// 素数判定を行う関数
fn is_prime(n: u128) -> bool {
    // 2未満の数は素数ではない
    if n < 2 {
        return false;
    }
    // 2から√nまでの数で割り切れるか調べる
    let sqrt_n = (n as f64).sqrt() as u128;
    for i in 2..sqrt_n {
        if n % i == 0 {
            // 割り切れたら素数ではない
            return false;
        }
    }
    // 割り切れる数がなければ素数
    true
}

// 偶数完全数を生成する関数
fn generate_even_perfect_numbers(n: u128) -> Vec<u128> {
    let mut perfect_numbers = Vec::new(); //完全数を格納するベクタ
    let mut p = 2;

    loop {
        //メルセンヌ数を計算
        let mersenne_number = 2u128.pow(p) - 1;
        // ルセンヌ素数かどうか調べる
        if is_prime(mersenne_number) {
            // メルセンヌ素数であれば完全数を算出
            let perfect_number = 2u128.pow(p - 1) * mersenne_number;
            // 調べる数の最大数まで調べたら終了
            if perfect_number > n {
                break;
            }
            // そうでなければ作成した完全数をリストに追加
            perfect_numbers.push(perfect_number);
        }
        p += 1;
    }
    perfect_numbers
}

fn main() {
    // 開始時間の記録
    let start = Instant::now();
    for pn in generate_even_perfect_numbers(1000000000000000000000000000000000000) {
        println!("{} は完全数です", pn);
    }
    // 経過時間を計算
    let duration = start.elapsed();
    println!("かかった時間: {:.6}秒", duration.as_secs_f64());
}

上のプログラムを実行すると、

6 は完全数です
28 は完全数です
120 は完全数です
496 は完全数です
8128 は完全数です
33550336 は完全数です
8589869056 は完全数です
137438691328 は完全数です
かかった時間: 0.001483秒

となりました。
PythonよりもRustが処理速度が速いので、このような結果になっています。
試しに、1000000000000000000000000000000000000をgenerate_even_perfect_numbersの引数に与えてみたのですが、以下のようになりました。

6 は完全数です
28 は完全数です
120 は完全数です
496 は完全数です
8128 は完全数です
33550336 は完全数です
8589869056 は完全数です
137438691328 は完全数です
2305843008139952128 は完全数です
かかった時間: 21.303184秒

かなり大きい数で調べてみたのですが、20秒程度でこの数までの完全数を探すことができています。