メルセンヌ素数を求める

2002/12/16作成

メルセンヌ素数

メルセンヌ数とは2p-1(p>0)で表される数です。これが素数であるとメルセンヌ素数と呼びます。

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

int LucasLehmerTest( mpz_t n, int p );

int main( int ac, char *av[] )
{
    int p;
    mpz_t m;

    /*  コマンドラインから指数を決定する    */
    if( ac < 2 ) {
        printf( "usage : mersenne p\n" );
        return( 1 );
    }
    p = strtol( av[1], NULL, 10 );

    /*  メルセンヌ数を求める    */
    mpz_init_set_ui( m, 2 );
    mpz_pow_ui( m, m, p );
    mpz_sub_ui( m, m, 1 );

    if( LucasLehmerTest( m, p ) )
        printf( "OK\n" );
    else
        printf( "NG\n" );

    mpz_clear( m );

    return( 0 );
}

int LucasLehmerTest( mpz_t m, int p )
{
    int i, ret;
    mpz_t s;

    mpz_init_set_ui( s, 4 );
    for( i = 2; i <= p - 1; i++ ) {
        mpz_mul( s, s, s );
        mpz_sub_ui( s, s, 2 );
        mpz_mod( s, s, m );
    }
    ret = ( mpz_cmp_ui( s, 0 ) == 0 );

    mpz_clear( s );

    return( ret );
}

1〜63までの値について調べてみた結果は、p=2, 3, 5, 7, 13, 17, 19, 31, 61の場合に素数になりました。メルセンヌはpが素数ならメルセンヌ素数になると予想を立てたそうですが、確かにそう考えたくなる数の並びです。でも現実にはこの予想は外れていたわけですが。(2020/3/4訂正)メルセンヌが予想したのはpが2, 3, 5, 7, 13, 17, 19, 31, 67, 127, 257の場合だけでした。(2020/3/4訂正終わり)

(2012/2/12追記)

GMPを用いてプログラムを書き直しました。処理手順は以前と変わりませんが、これにより探索範囲をメモリが許す限り広げることが出来るようになりました。しかし、素数判定を単純な方法で行っているため、計算時間が現実的ではないという問題が出てきました。これについては、メルセンヌ数が素数かどうかを判定するリュカ-レーマー法(Lucas-Lehmer Test)というものがありますので、これを取り入れることで解決できます。

(2012/3/9追記)

前述した通り、素数判定にリュカ-レーマー法(Lucas-Lehmer Test)を用いて書き直しました。当然ながら、この方法で高速に素数判定が行えるようになり、p<256で実行した結果pが3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127でメルセンヌ素数になることが検証できました。

ところで、リュカ-レーマー法(Lucas-Lehmer Test)を実装していて気が付いたのですが、これってp=2には使用できないのでしょうか。実用上問題はないのですが(2だけ手計算をすればいい)、そのことに触れている資料がなかったのが気になりました。

(2016/4/12追記)

素数判定が高速化できたので更に探索範囲を広げてみました。無計画にp<1000000でスタートしたのですが、あまりにも時間がかかってしまうのでp<100000が完了したところで一旦ストップ。しかしここまで2か月以上かかってしまいました。うーん、さすがに手ごわいですな。

p<100000の範囲でメルセンヌ素数であったのは次の通り。

3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243

当たり前ですが、全て既知のメルセンヌ素数です。ちなみにですが、p=100000を検証するのにTurion II NEO N54L(2.2GHz)で139.25秒掛かっています。p=100000に達するのに約2か月掛かりましたが、p=1000000まで完了するのに同じペースでも1500日程度かかることになります。実際にはpが大きくなるほど計算時間も長くなりますので、その何倍もの時間がかかることでしょう。なんか工夫が必要ですね。

(2020/4/25追記)

ソースコードを少し整理しました。

(2020/4/25追記)

pは素数だけを調べればいいということに今更ながら気が付きました。ということで改めて p < 100000 の範囲で調べたところ、4日程度で完了しました。上述の2か月から大幅に短縮されましたが、では p < 1000000 を探索するとすると100倍以上かかるでしょうから1年を超えてしまいます。これではちょっと現実的ではないですね。