【C#】リュカ–レーマーテストでメルセンヌ素数を特定する

執筆の動機

中学校を卒業して、春休みの自由時間ができたため図書館を徘徊していると、偶然「世界は素数でできている」(著 小島寛之)に出会い、素数に興味を持ったため․

メルセンヌ素数

現在発見されている世界最大の素数は、

2^82589933-1 (10進法表記で24,862,048桁)

である․ これは、2018年12月7日にGIMPSでPatrick Laroche氏により発見された․

このように、10進法表記において2^M-1とすることができる素数を メルセンヌ素数 といい、 M82589933 のように表記する․

2^M-1で表記できる数は、Lucas–Lehmer primality test(リュカ–レーマーテスト)と呼ばれるテストにより少ない計算リソースで素数判定が行える․
ゆえに、コンピュータの登場以降に発見されている巨大の素数のほとんどがこのメルセンヌ素数である․

一般的な素数判定

現時点では、一部の例外を省いて、任意の自然数nが与えられたときにそれが素数であるかを判定するには n mod (2~sqrt(n)) === 0が一つでも存在するか確認する必要がある․

具体的には、以下のような判定により行える(C#)․
なお、ここでは自然数nがintやdecimalを超えるものである可能性を考慮して、BigIntegerを使用する․

/*!
 * 巨大な数字に対する素数判定プログラム
 * (c) 2023 ActiveTK.
 */

using System.Numerics;
using static System.Console;

static class Test4Qiita
{
    /*!
     * エントリーポイント
     */
    static void Main()
    {
        BigInteger num = new BigInteger(1729), sqrt = Sqrt(num);
        int i;
        for (i = 2; i < sqrt; i++)
            if (num % BigInteger.Parse(i.ToString()) == 0)
                break;
        WriteLine(num.ToString() + "は素数: " + (i == sqrt));
        ReadLine();
    }

    /*!
     * 平方根を求める
     * @param BigInteger 0以上の実数․
     * @return BigInteger 入力値の平方根․
     */
    static BigInteger Sqrt(this BigInteger n)
    {
        if (n > 0)
        {
            int bitLength = Convert.ToInt32(Math.Ceiling(BigInteger.Log(n, 2)));
            BigInteger root = BigInteger.One << (bitLength / 2);
            while (!(n >= root * root && n < (root + 1) * (root + 1)))
            {
                root += n / root;
                root /= 2;
            }
            return root * root == n ? root : root + 1;
        }
        else if (n == 0) return 0;
        else throw new ArithmeticException("NaN");
    }
}

出力例

1729は素数: False

1729はいわゆる タクシー数 である․
1729 = 7・13・19 と素因数分解できるから、これは合成数であり、上記結果は正しい․

Miller–Rabin(ミラーラビン)素数判定法

一般的な素数判定よりはかなり高速であるが、テストに通っても確実に素数であるとは限らない ․
ここではソースコードのみ掲載する(C#)․

/*!
 * 素数判定プログラム
 * (c) 2023 ActiveTK.
 */
using System.Numerics;
using static System.Console;

static class Test4Qiita
{
    /*!
     * エントリーポイント
     */
    static void Main()
    {
        BigInteger num = new BigInteger(1729);
        WriteLine(num.ToString() + "は素数: " + IsPrime(num));
        ReadLine();
    }

    /*!
     * Miller–Rabin(ミラーラビン)素数判定
     * @param BigInteger 0以上の実数․
     * @param int 試行回数․ 大きいほど判定を誤る確率が減るが、その分時間がかかる․ 一般的には20～30を指定する․
     * @return BigInteger 入力値の平方根․
     */
    static bool IsPrime(BigInteger n, int k = 20)
    {
        if (n < 2)
            return false;
        if (n != 2 && n % 2 == 0)
            return false;
        Random rand = new Random();
        for (int i2 = 0; i2 < k; i2++)
        {
            BigInteger a = new BigInteger();
            do
            {
                a = BigInteger.Remainder(BigInteger.Abs(new BigInteger(rand.Next())), n);
            } while (a == 0);
            BigInteger x = BigInteger.ModPow(a, n - 1, n);
            if (x != 1)
                return false;
        }
        return true;
    }
}

メルセンヌ数の素数判定

さて、Lucas–Lehmer primality test(リュカ–レーマーテスト)の実装をC#で行う․
例として、M31(2^31-1)を素数判定する․

/*!
 * メルセンヌ数の素数判定プログラム
 * (c) 2023 ActiveTK.
 */

using System.Numerics;
using static System.Console;

static class Test4Qiita
{
    /*!
     * エントリーポイント
     */
    static void Main()
    {
        // M31
        // 2^31-1(=2,147,483,647)が素数であるか判定する․
        int num = 31;
        WriteLine("2^" + num + "-1は素数: " + LucasLehmerTest(num));
        ReadLine();
    }

    /*!
     * Lucas–Lehmer primality testを行う
     * @param int Mの値․ 例えば、2^1000-1が素数であるか判定したければ、1000を入力する․
     * @return BigInteger 入力値の平方根․
     */
    static bool LucasLehmerTest(int p)
    {
        // ここで、mは 2^p-1 である․
        // 通常の Pow ではpの値が大きな時に時間がかかってしまうため、シフト演算を用いる․
        BigInteger m = (BigInteger.One << p) - BigInteger.One;
        BigInteger s = 4;
        for (int i = 0; i < p - 2; i++)
        {
            s = (s * s - 2) % m;
            if (s == 0)
                return true;
        }
        return s == 0;
    }
}

出力例

2^31-1は素数: True

2^31-1は、数値に直すと2,147,483,647である․
このような大きい数を通常の方法で素数判定するには、多くのリソースが必要となる․
しかし、このアルゴリズムにより、容易に巨大な素数を発見することができた․

また、プログラムでは実装しなかったが、メルセンヌ数が素数になるのはMの値が素数である時に限られるため、そもそも入力値が素数であるか判定するとリソースの無駄使いを防ぐことができる ことを把握してほしい․

これは、上記「一般的な素数判定」をMに対して行うことにより実装できる․

世界最大の素数を探索

現在発見されている世界最大の素数はM82589933 (2^82589933-1)であり、さらに発見から4年以上の時間が経過していることに鑑みれば、世界最大の素数を探し当てるには少なくともM100,000,000以上の値域を探さなければならない․
これは10進法表記で3000万桁以上になり、非常に難易度が高くなっている․

先のLucas–Lehmer primality testを用いても、 M100,000,000以上のメルセンヌ数が素数であるか判定するためには、3000万桁の数を法とするmod計算が1億回必要 であるため、たいへん多くのリソースを要求する․

例えば、1996年からメルセンヌ素数を探索しているプロジェクト「GIMPS」を用いても、私のPCでは一つのメルセンヌ数が素数であるのか判定するだけで24時間程度を要する(i9 12900Kで1スレッド/メモリ96GB)․

メルセンヌ素数が無数に存在するのかは証明されていない上、メルセンヌ数が素数である確率は定かではないが、一般的な素数定理に基づけば、

\pi (x)\sim \frac{x}{\log (x)}

よって、2^100Mから2^100M+100000の間には約3,520個の素数があると予測される․
濃度によって考えるべきであると数学者に怒られそうだが、ここでは確率を用いて考えてみよう․

3520/100000=0.0352より、この付近で素数を発見できる確率は0.0352%ほどであると分かった․
メルセンヌ数が素数である確率と一定の範囲内に存在する素数の確率が一定であるかは未知数だが、ここではそのように仮定する․

1回以上出る確率が50%になる試行回数nを求めると、

0.5 = 1 - (1 - 0.000352)^n
(1 - 0.000352)^n = 0.5
n × ln(1 - 0.000352) = ln(0.5)
n = ln(0.5) / ln(1 - 0.000352)
n ≈ 1965

ゆえに、メルセンヌ数が素数であるのかの判定を約1965回行えば、メルセンヌ素数を一つ以上発見できる確率が50%になると予想される․
ただし、これはあくまでも期待値に過ぎない․

一つの2^100,000,000(誤差＋0~100000)-1を判定するのに24 CPU時間必要であるとすれば、世界最大のメルセンヌ素数を発見できる50%の期待値を得るには、24・1965=47160 CPU時間(3100年分)が必要である․

3100年のCPU時間と聞くと現実味を失うが、GoogleとCWIがSHA-1を衝突させた際に使用したCPU時間は6500年であるから、その約半分である․

予算さえあれば、クラウドでリソースを購入して試しても面白いかもしれない․

GPUによる挑戦

CPUは複雑な命令を高速で処理することに特化しており、コア数が少ない․
それに対して、GPUは当初グラフィックス処理を目的に開発されていただけあって、簡単で膨大な計算をこなせるよう、コア数や並列処理などが最適化されている․

そこで、私はGPUでLucas–Lehmer primality testを実行できないかと考えた․
複雑な処理のように思われるが、テストの中身としては巨大な数を法としたmod計算に過ぎない․
仮にGPUでテストを実行できたとすれば、CPU時間の100分の1以下の時間で実行できるのではないかと思う․

C#は.NETの中間言語を介するため、ハードに対して複雑な命令を行うことができず、GPUの扱いが困難である․
そこで、私は現在C#で書いたコードをC++で書き直して、GPUで処理できないか試みている․
完成次第、AWSかGCPのGPUプランを契約して、世界最大の素数を探索してみたいと思う․