素数
素数というのは、1と自分自身以外の約数を持たない自然数の事です。まずは、この定義通りにプログラムを作り素数を調べだしてみます。
#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;
/* コマンドラインから素数探索範囲を決定する */
if( ac < 3 ) {
printf( "prime1 min_n max_n\n" );
return( 1 );
}
min_n = strtol( av[1], NULL, 10 );
max_n = strtol( av[2], NULL, 10 );
/* 探索範囲の数を調べる */
for( n = min_n; n <= max_n; n++ ) {
if( isprime( n ) )
printf( "%d\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以上の値で次々に割ってみて余りが出るかどうかを調べています。
素数の判定を効率よく行うため、2以外は奇数のみで割っているのと、平方根より少し大きい値までしか調べておりません。平方根についてはsqrt()というANSI標準関数がありますが、double型で計算するため大きな整数に対しては誤差が出るため使用していません。
このプログラムは比較的性能が良く、1億以下の素数を約1時間で全て求めることが出来ました。しかし、より大きな値まで探索範囲を広げるには、このプログラムだけでは不十分だと思われます。
(2012/2/12追記)
プログラムに冗長なところがあったので修正しました。Athlon 64 X2 5600+(2.9GHz)で同じく1億以下を求めてみたのですが、6分ほどで完了しました。以前に実行した時のCPUが何であったか覚えてないのですが、随分速くなったものです。
(2012/2/15追記)
探索範囲を10億以下まで広げてみました。同じくAthlon 64 X2 5600+(2.9GHz)で約4時間掛かりました。
(2012/2/16追記)
素数判定を行う際、従来は2の倍数を特別扱いして除数から弾いていたのですが、これに3の倍数も加えて弾くようにしました。これによって、1億以下の素数を求めるのが約6分、10億以下で約2時間半と多少高速化することが出来ました。ならば5の倍数、7の倍数、11の倍数と弾く数をどんどん増やしていけばいいようにも思いますが、そうすると除数の数列を求める計算がどんどん複雑になってしまって、トータルでは高速化にはならないものと予想されるため、3の倍数までの除外としています。
(2012/2/22追記)
素数判定で5以上の倍数を除外する理由として、除数を求める処理が複雑化することを挙げていましたが、よく考えるとこれは間違いでした。除数列の隣り合う数の差は循環しますので、これをテーブル化することで処理速度は問題なくなります。例えば、2、3、5の倍数を除外する処理は次のようになります。
int isprime( int n )
{
int i, n2, dn, d[] = { 4, 2, 4, 2, 4, 6, 2, 6 };
/* 1は素数ではない */
if( n == 1 )
return( 0 );
/* 2,3,5は素数 */
if( n == 2 || n == 3 || n == 5 )
return( 1 );
/* 2,3,5で割り切れたら合成数 */
if( n % 2 == 0 || n % 3 == 0 || n % 5 == 0 )
return( 0 );
/* sqrt(n)以上の出来るだけ小さい値を求める */
n2 = 2;
while( n2 * n2 < n ) {
n2 *= 2;
}
/* n2以下の2,3,5の倍数以外での剰余が0かどうか調べる */
dn = 0;
for( i = 7; i <= n2; i += d[dn], dn = ( dn + 1 ) % ( sizeof(d)/sizeof(int) ) ) {
if( n % i == 0 )
return( 0 );
}
/* 素数であった */
return( 1 );
}
しかし、この修正を行っても1億以下で約5分とさほど高速化はされません。それはなぜかというと、除外される数字の比率が低いからです。
2の倍数を除数から除外すると、処理量は1/2になりますので2倍高速化されます。では3の倍数を除外すると処理量は2/3になるかというと、そういうわけではありません。既に2の倍数は除外していますので、2と3の公倍数は除外済みなのです。ですので3の倍数で2の倍数でないものだけしか除外できません。この数の比率は1/6ですので17%程度の高速化にとどまります。これに更に5の倍数を除外すると1/30しか除外できませんので、3%程度しか高速化できません。7の倍数、11の倍数と除外していっても、その除外できる比率はどんどん少なくなるわけです。
(2012/3/18追記)
nの素数判定を行う際は√nまでの数で割り切れるか調べればよいのですが、これまで平方根を求める処理を簡略化して近似値を用いていました。これを改め、きちんと平方根を求めるように修正しました。余分な判定を行わなくてよくなった分多少高速化し、10億以下の素数を求めるのに約2時間となりました。
(2020/4/25追記)
ソースコードを少し整理しました。
エラトステネスのふるい
前節の単純な方法では、素数を求める事はできるものの非常に時間が掛かってしまいます。そこで、効率よく素数を求める方法としてよく知られている「エラトステネスのふるい」という方法を使ってみます。これは、自然数に対応するフラグの配列を用意し、自然数の倍数に印を付けていくというものです。最後まで印が付かなかったものが素数であることになります。
#include <stdio.h>
#define ARRAY_SIZE (100000000)
#define ARRAY_SIZE_SQRT (10000)
char array[ARRAY_SIZE];
int main( 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;
}
}
/* ふるいで残った数は素数である */
for( n = 2; n < ARRAY_SIZE; n++ ) {
if( array[n] == 1 )
printf( "%d\n", n );
}
return( 0 );
}
期待通り、高速に素数探索を行うことが出来ました。1000万以下の素数を求めるのに、先ほどのプログラムは2分半掛かったのですが、エラトステネスのふるいを使うと24秒で済み、約5倍の速さで計算することが出来ました。正直言って、期待した程の性能は出なかったのですが、これは探索範囲が狭いために、アルゴリズムによる差が出にくいということだと思われます。
このプログラムの欠点は探索範囲分のフラグ配列を必要とするため、メモリの許す程度の範囲しか探せないという事です。
(2005/4/18追記)
何年も前に指摘を頂いていた、ふるいに掛ける時に既にふるい落とされた数の倍数も繰返し調べるという全く無駄な処理をしていた点を修正しました。これにより更に4倍ほど処理速度が速くなりました。
(2012/2/12追記)
こちらもAthlon 64 X2 5600+(2.9GHz)で再実行してみたところ、0.72秒でした。
(2012/2/16追記)
メモリが潤沢に使用できるようになったので、配列のサイズを1億に増やしました。Athlon 64 X2 5600+(2.9GHz)での実行時間は8秒でした。
(2012/3/10追記)
ふるいにかけるのは√(MAX_ARRAY)の倍数までで十分なので、そのように修正しました。また、nの倍数をふるいにかけるときに2nから開始していたのですが、これは2の倍数で既にふるい落とし済みで不要な処理でした。ですのでn*nから始めるように修正しました。この2点の修正で、Athlon 64 X2 5600+(2.9GHz)で1億以下を求めるのに約5秒となりました。
(2020/4/25追記)
ソースコードを少し整理しました。
エラトステネスのふるいを繰り返し使う
メモリの許す範囲しか素数を探せないというエラトステネスのふるいの欠点を、ふるいを繰り返し使用するという方法で解決してみました。今回のプログラムでは、区間の大きさを1000万に設定してみました。もちろん区間が大きいほど効率が良くなるのですが、その為には大量のメモリが必要になってきます。
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#define ARRAY_SIZE_MAX (100000000)
char array[ARRAY_SIZE_MAX];
int primelistcount;
int *primelist;
int initprimelist( void );
int isqrt( int x );
int main( int ac, char *av[] )
{
int base, arraysize, n, min_n, max_n, i, j;
/* コマンドラインから素数探索範囲を決定する */
if( ac < 3 ) {
printf( "usage : prime3 min_n max_n\n" );
return( 1 );
}
min_n = strtol( av[1], NULL, 10 );
max_n = strtol( av[2], NULL, 10 );
if( !initprimelist() )
return( 0 );
/* ARRAY_SIZE_MAX分ごとの整数区間をふるいにかける */
for( base = min_n; base < max_n; base += ARRAY_SIZE_MAX ) {
/* 整数区間配列の大きさを決める */
if( max_n - base + 1 > ARRAY_SIZE_MAX )
arraysize = ARRAY_SIZE_MAX;
else
arraysize = max_n - base + 1;
/* 配列を初期化する */
for( i = 0; i < arraysize; i++ )
array[i] = 1;
/* 配列をふるいにかける */
for( i = 0; i < primelistcount; i++ ) {
n = primelist[i];
if( base <= n )
j = n;
else if( base % n == 0 )
j = base;
else
j = ( base / n + 1 ) * n;
j = j < n * n ? n * n : j;
for( ; j < base + arraysize; j += n )
array[j - base] = 0;
}
/* ふるいで残った数は素数である */
for( n = 0; n < arraysize; n++ ) {
if( array[n] == 1 )
printf( "%d\n", base + n );
}
}
return( 0 );
}
int initprimelist( void )
{
int n, i;
int sqrt_int_max;
/* sqrt(INT_MAX/2)を求める */
sqrt_int_max = isqrt( INT_MAX );
/* 配列を初期化する */
for( i = 0; i <= sqrt_int_max; i++ )
array[i] = 1;
/* 配列をふるいにかける */
for( n = 2; n <= sqrt_int_max; n++ ) {
if( array[n] == 1 ) {
for( i = n * n; i <= sqrt_int_max; i+=n ) {
array[i] = 0;
}
}
}
/* 素数一覧配列にコピー */
primelistcount = 0;
for( n = 2; n <= sqrt_int_max; n++ ) {
if( array[n] == 1 )
primelistcount++;
}
primelist = calloc( primelistcount, sizeof(int) );
if( primelist == NULL )
return( 0 );
for( n = 2, i = 0; n <= sqrt_int_max; n++ ) {
if( array[n] == 1 )
primelist[i++] = n;
}
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 );
}
この方法では、1億以下の素数を全て調べるのに6分程度で済みました。一番最初の単純な方法では同じ範囲を約1時間掛かりましたから、やはり探索範囲を広げるほどエラトステネスのふるいの性能がよく表れるようになっています。
(2012/2/12追記)
こちらもAthlon 64 X2 5600+(2.9GHz)で再実行してみたところ、1億以下で48秒でした。
(2012/2/15追記)
こちらも探索範囲を10億以下まで広げてみました。同じくAthlon 64 X2 5600+(2.9GHz)で約17分で終了しました。単純な方法に比べるとやはり高速ですね。
(2012/2/16追記)
こちらも配列のサイズを1億に増やしました。1億以下の探索は50秒、10億以下の探索は11分と探索範囲が広い場合には多少高速化しました。
(2012/3/11追記)
このプログラムではINT_MAXまでの素数しか探索しませんので、予め√(INT_MAX)までの素数一覧を作成しておき、ふるいにかけるのはこの素数一覧の倍数のみで行うように修正しました。この修正により、Athlon 64 X2 5600+(2.9GHz)で1億以下で4秒、10億以下の素数を求めるのに44秒と大変高速化できました。
(2012/3/18追記)
先頭のブログラムと同様に、平方根を求める処理をきちんと行うようにしました。ライブラリのsqrt()を用いても32bit整数の範囲内では問題ないのですが、64bit整数化する際に誤差が発生してしまうため、こちらも自前で平方根を求めています。
(2020/4/25追記)
ソースコードを少し整理しました。
素数一覧
ここにあった素数一覧データは素数一覧に移動しました。