ゴールドバッハ予想
ゴールドバッハの予想とは「全ての偶数の合成数は、2個の素数の和で表すことができる」というものです。コラッツ予想と同じく、予想ですからまだ証明はされていません。
偶数の合成数ですから、2は含まれず、4以上の偶数が対象になります。4は2と2の和で表される事は簡単に分かります。偶数の素数は2だけですから、4以外で偶数素数の和で表される偶数は無いことになります。と言うことは、ゴールドバッハ予想は「全ての6以上の偶数は、2個の奇素数の和で表すことが出来る」と等価である事になります。
以下のプログラムは、ゴールドバッハ予想の定義のままに作成しています。
#include <stdio.h>
#include <stdlib.h>
int isprime( int n );
int isqrt( int x );
int main( int ac, char *av[] )
{
int n, min_n, max_n, m1, m2;
/* コマンドラインから探索範囲を決定する */
if( ac < 3 ) {
printf( "usage : goldbach min_n max_n\n" );
return( 1 );
}
min_n = strtol( av[1], NULL, 10 );
min_n = min_n % 2 == 1 ? min_n + 1 : min_n;
max_n = strtol( av[2], NULL, 10 );
max_n = max_n % 2 == 1 ? max_n + 1 : max_n;
/* 探索範囲の数を調べる */
for( n = min_n; 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, 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)を求める */
n2 = isqrt( n );
/* 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 );
}
int isqrt( int x )
{
int s, t;
if( x == 0 ) return 0;
s = 1;
t = x;
while( s < t ) {
s <<= 1;
t >>= 1;
}
do {
t = s;
s = ( x / s + s ) >> 1;
} while( s < t );
return( t );
}
素数の和で表されるまで単純に全ての数の組み合わせについて調べているので、処理速度は大変遅いです。素数の判定も単純な方法ですから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追記)
素数と同様に素数判定において平方根を求める処理をきちんと実装しました。
(2020/4/25追記)
ソースコードを少し整理しました。
素数テーブル
(2012/2/22追記)
ゴールドバッハ予想の検証を行うためには素数判定を何度も行う必要があります。上記の単純なプログラムでは、同じ数に対して素数判定を何度も行っているため、処理時間が無駄に掛かってしまっていました。そこで、事前にある程度の範囲までは素数一覧テーブルを作成し、その範囲の数については素数判定はテーブルを参照することで処理を高速化してみました。
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#define ARRAY_SIZE (100000000)
#define ARRAY_SIZE_SQRT (10000)
char array[ARRAY_SIZE];
void initarray( void );
int isprime( int n );
int isqrt( int x );
int main( int ac, char *av[] )
{
int n, min_n, max_n, m1, m2;
/* コマンドラインから探索範囲を決定する */
if( ac < 3 ) {
printf( "usage : goldbach2 min_n max_n\n" );
return( 1 );
}
min_n = strtol( av[1], NULL, 10 );
min_n = min_n % 2 == 1 ? min_n + 1 : min_n;
max_n = strtol( av[2], NULL, 10 );
max_n = max_n % 2 == 1 ? max_n + 1 : max_n;
initarray();
/* 探索範囲の数を調べる */
for( n = min_n; 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 );
}
void initarray( void )
{
int n, i;
/* 配列を初期化する */
for( i = 0; i < ARRAY_SIZE; i++ )
array[i] = 1;
/* 配列をふるいにかける */
for( n = 2; n <= ARRAY_SIZE_SQRT; n++ ) {
if( array[n] == 1 ) {
for( i = n * n; i < ARRAY_SIZE; i += n )
array[i] = 0;
}
}
}
int isprime( int n )
{
int i, n2, d;
if( n < ARRAY_SIZE ) {
return( array[n] );
}
/* 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)を求める */
n2 = isqrt( n );
/* 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 );
}
int isqrt( int x )
{
int s, t;
if( x == 0 ) return 0;
s = 1;
t = x;
while( s < t ) {
s <<= 1;
t >>= 1;
}
do {
t = s;
s = ( x / s + s ) >> 1;
} while( s < t );
return( t );
}
この修正の効果は顕著で、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()の処理時間の占める割合は微小なので、計算時間には影響はほとんどありません。
(2020/4/23追記)
素数配列のフラグ処理を間違えていて、配列中の全ての値が素数であると判定してしまっているという、とてつもなく巨大なバグがあることに8年たってようやく気が付きました。大変申し訳ないです。
プログラムを修正し再度計測してみますと、1億以下の範囲では単純な方法では1時間程度かかっていたものが、素数配列を使ったものでは6秒と大幅に高速化しています。しかし、10億以下の範囲では単純な方法では30時間程度かかっていたものが、素数配列を使ったものでは29時間程度と、微小な高速化にとどまっています。10億以下の範囲の数値に対して1億以下の数値は素数配列で判定しても対象となるのは全体の1割程度です。でれば素数判定処理も1割程度の高速化にしかならないと見込まれるので、こちらの結果が妥当でしょう。
結果として、素数配列を用いた素数判定の高速化は、素数配列の範囲を超えた場合には効果は限定的になるということで、大きな効果があるとしたこれまでの結論は間違っていました。お詫びして訂正いたします。
合わせて、これまではより広範囲に探索が可能と期待してGMP化したプログラムを載せていましたが、あまり意味が無いので非GMPのプログラムに戻しています。