数学をプログラミング

Python×Rustで数学体験!サービト・イブン=クッラの方法で友愛数を見つけよう

前回の記事「PythonとRustで!友愛数の組を見つけてみよう!」では、友愛数の定義からPythonとRustを使ってコードを実装し、友愛数の組を探しました。

今回は「友愛数入門! 数の面白シリーズ!」の記事の中で紹介したサービト・イブン=クッラの方法を使用して友愛数の組を探します。

友愛数とサービト・イブン=クッラの方法

まずは友愛数の定義を思い出します。

友愛数

2つの自然数ab友愛数であるとは、aの真の約数の和がbに等しく、かつbの真の約数の和がaに等しいことである。

以下は友愛数の組の例です。

  • (220,284)
  • (1184, 1210)
  • (2620, 2924)
  • (5020, 5564)
  • (6232, 6368)

友愛数を具体的に探そうとするとかなり大変です。
友愛数を探す手段として、次のサービト・イブン=クッラの方法があります。
ただし、この方法で探すことができる友愛数の組はかなり限定的です。
サービト・イブン=クッラの方法で探すことができない友愛数の組の例として、(6232, 6368) があります。

サービト・イブン=クッラの方法

整数n \geq 2に対して、次のようにp,q,rを定める:
p = 3 \times 2^{n-1} – 1
q = 3 \times 2^n -1
r = 9 \times 2^{2n-1} – 1

この時、p,q,rが素数であるならば、(2^npq,2^nr)は友愛数の組となる。

Pythonでサービト・イブン=クッラの方法を実装しよう

Pythonでサービト・イブン=クッラの方法を実装します。
ソースコードは以下になります。

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


# 友愛数を探す関数
# limit : サービト・イブン=クッラの方法に出てくる整数nの範囲の最大値
def thabit_ibn_qurra_method(limit):
    """友愛数のペアをサービト・イブン=クッラの方法で探す"""
    result = [] # 友愛数の組のリストを初期化

    for n in range(2, limit):
        # サービト・イブン=クッラの方法に従って計算
        p = 3 * (2 ** (n - 1)) - 1
        q = 3 * (2 ** n) - 1
        r = 9 * (2 ** (2 * n - 1)) - 1

        # p,q,rがすべて素数である必要がある
        if is_prime(p) and is_prime(q) and is_prime(r):
            # 友愛数の候補を計算
            a = 2 ** n * p * q
            b = 2 ** n * r
            # 友愛数の組をリストに追加
            result.append((a, b))

    return result


if __name__ == '__main__':
    amicable_pairs = thabit_ibn_qurra_method(20)
    for a, b in amicable_pairs:
        print(f"{a} と {b}")

実行すると次の結果が表示されます。

220 と 284
17296 と 18416
9363584 と 9437056

ソースコードについて、解説します。
is_prime関数は完全数のところでも出てきましたが、与えられた引数nが素数かどうかを判定するものです。
thabit_ibn_qurra_method関数がサービト・イブン=クッラの方法を使う関数です。
limitがサービト・イブン=クッラの方法に出てくる整数nの動く範囲の最大値になります。
この関数は次のような処理を行っています。

  1. 2からlimit-1までの各整数nをサービト・イブン=クッラの方法に出てくる3つの式に代入
  2. その3つの式により算出される値全てが素数であるかをis_prime関数でチェック
  3. 3つの値がすべて素数ならば、(2^npq,2^nr)は友愛数の組になるので、その2つの値を計算して友愛数の組のリストに追加
  4. limit-1までの全てのnを調べたら、友愛数の組のリストを返却

Rustでもサービト・イブン=クッラの方法を実装しよう

Rust版のソースコードは以下になります。

// 素数判定を行う関数
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 thabit_ibn_qurra_method(limit: u32) -> Vec<(u128, u128)> {
    let mut result = Vec::new(); //完全数を格納するベクタ

    for n in 2..limit {
        // サービト・イブン=クッラの方法に従って計算
        let p = 3 * 2u128.pow(n - 1) - 1;
        let q = 3 * 2u128.pow(n) - 1;
        let r = 9 * 2u128.pow(2 * n - 1) - 1;

        // p,q,rがすべて素数である必要がある
        if is_prime(p) && is_prime(q) && is_prime(r) {
            // 友愛数の候補を計算
            let a: u128 = 2u128.pow(n) * p * q;
            let b: u128 = 2u128.pow(n) * r;
            // 友愛数の組をリストに追加
            result.push((a, b));
        }
    }
    result
}

fn main() {
    let limit: u32 = 60;
    let amicable_pairs = thabit_ibn_qurra_method(limit);
    for (a, b) in amicable_pairs {
        println!("{} と {}", a, b);
    }
}

やっていることはPythonと同じです。
limitの値が60を超えるとサービト・イブン=クッラの方法に出てくる計算式がu128の値を超えるのでオーバーフローして実行途中でエラーが吐かれます。
こちらも結果は以下になります。

220 と 284
17296 と 18416
9363584 と 9437056   

まとめ

サービト・イブン=クッラの方法を使って、PythonとRustで友愛数の組を探しました。
この方法でn=2からn=60までで探すことができた組は

  • (220, 284)
  • (17296, 18416)
  • (9363584, 9437056)

でした。
友愛数の組の例を上で紹介していますが、(1184, 1210)(2620, 2924)など、この方法で探すことができていません。
そういう意味でサービト・イブン=クッラの方法は友愛数の組を探すという意味でかなり限定的な方法と言えます。

これまで

という3本の記事で、友愛数の定義や、友愛数の組を探す方法を解説して、プログラムで実際に定義とサービト・イブン=クッラの方法を使用して友愛数の組を探してみました。

次回からは、「トーラスの中を歩く!? 数学が教える“穴”の不思議な世界」という主題で位相幾何学という数学の一分野についてお話をします。
厳密性には欠けますが、位相幾何学という面白い数学の話に触れてみます!
お楽しみに!

ABOUT ME
MSK
数学・プログラミング・ゲーム制作が大好きな30代エンジニア。 趣味でUnityやC言語を使って数学を可視化したり、小さなアプリを作ったりしています。 教育・学びの楽しさにも関心があります。