メルセンヌ素数を探す

メルセンヌ数とは2^p-1(p>0)で表される数です。これが素数であるかどうか調べてみます。

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

int isprime( unsigned long long n );

int main( int ac, char *av[] )
{
    unsigned long long  n, m, i;

    /*  コマンドラインから指数を決定する    */
    if( ac < 2 )
        return;
    n = strtouq( av[1], NULL, 10 );

    /*  メルセンヌ数を求める    */
    m = 1;
    for( i = 0; i < n; i++ )
        m *= 2;
    m--;

    if( isprime( m ) )
        printf( "%qu OK\n", m );
    else
        printf( "%qu NG\n", m );
}

int isprime( unsigned long long n )
{
    unsigned long long  i, n2, pn2;

    /*  1は素数ではない */
    if( n == 1 )
        return( 0 );

    /*  2以外で2で割り切れたら合成数    */
    if( n != 2 && n % 2 == 0 )
        return( 0 );

    /*  sqrt(n)以上の出来るだけ小さい値を求める */
    pn2 = n2 = n;
    while( n2 * n2 > n || n2 * n2 < n2 ) {
        pn2 = n2;
        n2 /= 2;
    }
    if( pn2 >= n )
        n2 = n - 1;
    else
        n2 = pn2;

    /*  n2以下の奇数での剰余が0かどうか調べる   */
    for( i = 3; i <= n2; i += 2 ) {
        if( n % i == 0 )
            return( 0 );
    }

    /*  素数であった    */
    return( 1 );
}

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