フェルマー数を求める

フェルマー数

フェルマー数とはフェルマーが考案した素数を生成する式で、2²ⁿ+1(n≧0)というものです。

プログラムは例によってフェルマー数の定義そのままに作成しています。べき乗をループを回して計算してしまっていますが、これはGMPにmpz_pow()が存在しないためです(あっておかしくないと思うのですが、なぜかありません)。

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

int isprime( const mpz_t n );

int main( int ac, char *av[] )
{
    int n;
    mpz_t f, p, i;

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

    /*  フェルマー数を求める    */
    mpz_init_set_ui( p, 2 );
    mpz_pow_ui( p, p, n );

    mpz_init_set_ui( i, 1 );
    mpz_init_set_ui( f, 2 );
    while( mpz_cmp( i, p ) < 0 ) {
        mpz_mul_ui( f, f, 2 );
        mpz_add_ui( i, i, 1 );
    }
    mpz_add_ui( f, f, 1 );

    if( isprime( f ) )
        printf( "OK\n" );
    else
        printf( "NG\n" );

    mpz_clear( f );
    mpz_clear( p );
    mpz_clear( i );

    return( 0 );
}

int isprime( const mpz_t n )
{
    mpz_t i, tmp;
    int d, ret;

    /*  mpz_probab_prime_p() で確定するか確認 */
    switch( mpz_probab_prime_p( n, 25 ) ) {
        case  0:
            return( 0 );
        case  2:
            return( 1 );
    }

    /*  自乗がn以下の2,3の倍数以外での剰余が0かどうか調べる   */
    d = 2;
    mpz_init_set_ui( i, 5 );
    mpz_init( tmp );
    ret = 1;
    while( 1 ) {
        mpz_mul( tmp, i, i );
        if( mpz_cmp( tmp, n ) > 0 )
            break;

        if( mpz_divisible_p( n, i ) ) {
            ret = 0;
            break;
        }
        mpz_add_ui( i, i, d );
        d = ( d == 2 ? 4 : 2 );
    }

    mpz_clear( i );
    mpz_clear( tmp );

    return( ret );
}

nが0,1,2,3,4の場合、フェルマー数は素数となります。しかし、nが5の場合はフェルマー数は4294967297となりこれは合成数でした。この先はどうかというと、現在のところ5以上のフェルマー数については合成数しか見つかっていないそうです。

(2012/2/14追記)

メルセンヌ素数と同様に、GMPを用いてプログラムを書き直しました。

(2020/4/25追記)

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

(2026/3/7追記)

素数判定で mpz_probab_prime_p() を使用するように修正するとともに、ソースコードを少し整理しました。素数生成式を検証すると同様、mpz_probab_prime_p() は確定時のみ採用しています。