たまたま最速の素数判定プログラム C# Java C++ という記事を見まして、条件によってはそれよりも更に速い方法もあるよと思ったので記事を書きます。

そのコードはこちら。GMP使ってます。

# define MAX_ARRAY   (100000000)
# define MAX_ARRAY_SQRT (10000)

char    array[MAX_ARRAY];

void initarray( void )
{
    int n, i;

    /*  配列を初期化する    */
    for( i = 0; i < MAX_ARRAY; i++ )
        array[i] = 1;

    /*  配列をふるいにかける    */
    for( n = 2; n <= MAX_ARRAY_SQRT; n++ ) {
        if( array[n] == 0 ) {
            for( i = n * n; i < MAX_ARRAY; i += n )
                array[i] = 0;
        }
    }
}

int isprime( mpz_t n )
{
    mpz_t i, n2, tmp;
    int _n, ret;

    if( mpz_cmp_ui( n, MAX_ARRAY ) < 0 ) {
        _n = mpz_get_ui( n );
        return( array[_n] );
    }

    /*  2割り切れたら合成数    */
    if( mpz_even_p( n ) )
        return( 0 );

    mpz_init( i );
    mpz_init( n2 );
    mpz_init( tmp );

    /*  sqrt(n)を求める */
    mpz_sqrt( n2, n );

    mpz_set_ui( i, MAX_ARRAY + 1 );
    ret = 1;
    while( mpz_cmp( i, n2 ) < 0 ) {
        if( mpz_divisible_p( n, i ) ) {
            ret = 0;
            break;
        }
        mpz_add_ui( tmp, i, 2 );
        mpz_set( i, tmp );
    }

    mpz_clear( i );
    mpz_clear( n2 );
    mpz_clear( tmp );

    return( ret );
}

要するに、一定以下の整数については予め素数判定を行って結果を配列に保持しておき、その範囲の素数判定は配列をひくだけで行うということです。

「はー？」と思った方もいらっしゃるかもしれませんが、実際のところとしてゴールドバッハ予想を単純な方法で検証した際には、24時間かかっていた処理を22秒で終わらせるという成果をあげたりもしてます。

一つの整数を一度だけ素数判定するなら元記事の方法がまあ速いでしょうけれど、多数の整数を繰り返し素数判定するならこんな方法もありますよというお話でした。