ゴールドバッハ予想

2001/11/28作成

ゴールドバッハ予想

ゴールドバッハの予想とは「全ての偶数の合成数は、2個の素数の和で表すことができる」というものです。コラッツ予想と同じく、予想ですからまだ証明はされていません。

偶数の合成数ですから、2は含まれず、4以上の偶数が対象になります。4は2と2の和で表される事は簡単に分かります。偶数の素数は2だけですから、4以外で偶数素数の和で表される偶数は無いことになります。と言うことは、ゴールドバッハ予想は「全ての6以上の偶数は、2個の奇素数の和で表すことが出来る」と等価である事になります。

以下のプログラムは、ゴールドバッハ予想の定義のままに作成しています。

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

#define START_R (0x4000)

int isprime( int n );

int main( int ac, char *av[] )
{
    int n, max_n, m1, m2;

    /*  コマンドラインから探索範囲を決定する    */
    if( ac < 2 )
        return( 1 );
    max_n = strtol( av[1], NULL, 10 );
    max_n = max_n % 2 == 1 ? max_n + 1 : max_n;

    /*  探索範囲の数を調べる    */
    for( n = 6; n <= max_n; n += 2 ) {
        /*  足してnになる整数組(m1,m2)が共に素数か調べる    */
        for( m1 = 3, m2 = n - 3; m1 <= m2; m1 += 2, m2 -= 2 ) {
            if( isprime( m1 ) && isprime( m2 ) )
                break;
        }
        if( m1 > m2 )
            printf( "%d NG\n", n );
    }

    return( 0 );
}

int isprime( int n )
{
    int i, n2, r, d;

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

    /*  2,3は素数 */
    if( n == 2 || n == 3 )
        return( 1 );

    /*  2,3で割り切れたら合成数    */
    if( n % 2 == 0 || n % 3 == 0 )
        return( 0 );

    /*  sqrt(n)を求める */
    r = START_R;
    n2 = 0;
    while( r ) {
        if( n >= ( n2 + r ) * ( n2 + r ) )
            n2 += r;
        r /= 2;
    }

    /*  n2以下の2,3の倍数以外での剰余が0かどうか調べる   */
    d = 2;
    for( i = 5; i <= n2; i += d, d = ( d == 2 ? 4 : 2 ) ) {
        if( n % i == 0 )
            return( 0 );
    }

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

素数の和で表されるまで単純に全ての数の組み合わせについて調べているので、処理速度は大変遅いです。素数の判定も単純な方法ですから2重に遅いことになります。それでも簡易なプログラムで安定してゴールドバッハ予想の探索が出来ると言うメリットがあります。

実際に1億以下の数について調べてみましたが、全ての数においてゴールドバッハ予想が成り立ちました。しかし、1億以下の範囲を探索するのに約10時間掛かっているため、これ以上の範囲について探索を行うのは困難です。

(2012/2/12追記)

素数と同様に、プログラムに冗長なところがあったので修正しました。それとともにAthlon 64 X2 5600+(2.9GHz)で1億以下を再計算してみたところ、1時間半程度で終了しました。

(2012/2/19追記)

素数と同様に素数判定の処理を高速化しました。Athlon 64 X2 5600+(2.9GHz)で1億以下を再計算したところ、約50分で終了しました。更に10億以下を検証したところ約24時間で全てゴールドバッハ予想が成り立ちました。

(2012/3/18追記)

素数と同様に素数判定において平方根を求める処理をきちんと実装しました。

素数テーブル

(2012/2/22追記)

ゴールドバッハ予想の検証を行うためには素数判定を何度も行う必要があります。上記の単純なプログラムでは、同じ数に対して素数判定を何度も行っているため、処理時間が無駄に掛かってしまっていました。そこで、事前にある程度の範囲までは素数一覧テーブルを作成し、その範囲の数については素数判定はテーブルを参照することで処理を高速化してみました。

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

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

char    array[MAX_ARRAY];

void initarray( void );
int isprime( mpz_t n );

int main( int ac, char *av[] )
{
    mpz_t n, max_n, m1, m2, tmp;

    /*  コマンドラインから探索範囲を決定する    */
    if( ac < 2 )
        return( 1 );
    mpz_init_set_str( max_n, av[1], 10 );
    mpz_init( n );
    mpz_init( m1 );
    mpz_init( m2 );
    mpz_init( tmp );

    initarray();

    /*  探索範囲の数を調べる    */
    mpz_set_ui( n, 6 );
    while( mpz_cmp( n, max_n ) <= 0 ) {
        /*  足してnになる整数組(m1,m2)が共に素数か調べる    */
        mpz_set_ui( m1, 3 );
        mpz_sub_ui( m2, n, 3 );
        while( mpz_cmp( m1, m2 ) <= 0 ) {
            if( isprime( m1 ) && isprime( m2 ) )
                break;
            mpz_add_ui( tmp, m1, 2 );
            mpz_set( m1, tmp );
            mpz_sub_ui( tmp, m2, 2 );
            mpz_set( m2, tmp );
        }
        if( mpz_cmp( m1, m2 ) > 0 )
            gmp_printf( "%Zd NG\n", n );

        mpz_add_ui( tmp, n, 2 );
        mpz_set( n, tmp );
    }

    mpz_clear( max_n );
    mpz_clear( n );
    mpz_clear( m1 );
    mpz_clear( m2 );
    mpz_clear( tmp );

    return( 0 );
}

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 );
}

この修正の効果は顕著で、Athlon 64 X2 5600+(2.9GHz)での10億以下のゴールドバッハ予想の検証が、わずか22秒で完了しました。

(2012/3/7追記)

高速化に成功したので、探索範囲を広げるためにGMP化しました。10億以下の検証が約6分、100億以下が約50分、1000億以下が約9時間、1兆以下が約94時間で完了しました。それぞれゴールドバッハ予想は成り立っています。

(2012/3/10追記)

initarray()における処理を素数と同様に修正しました。ただ、このプログラムにおいてinitarray()の処理時間の占める割合は微小なので、計算時間には影響はほとんどありません。


あおやぎのさいと2.0初歩の整数論プログラミング