素数というのは、1と自分自身以外の約数を持たない自然数の事です。 まずは、この定義通りにプログラムを作り素数を調べだしてみます。
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
int isprime( unsigned long long n );
int main( int ac, char *av[] )
{
unsigned long long n, min_n, max_n;
/* コマンドラインから素数探索範囲を決定する */
if( ac < 3 )
return;
min_n = strtouq( av[1], NULL, 10 );
if( min_n < 2 )
min_n = 2;
max_n = strtouq( av[2], NULL, 10 );
/* 探索範囲の数を調べる */
for( n = min_n; n <= max_n; n++ ) {
if( isprime( n ) )
printf( "%qu\n", n );
}
}
int isprime( unsigned long long n )
{
unsigned long long i, n2, pn2;
/* 1は素数ではない */
if( n == 1 )
return( 0 );
/* 2以外で2で割り切れたら合成数 */
if( n != 2 && n % 2 == 0 )
return( 0 );
/* sqrt(n)以上の出来るだけ小さい値を求める */
pn2 = n2 = n;
while( n2 * n2 > n || n2 * n2 < n2 ) {
pn2 = n2;
n2 /= 2;
}
if( pn2 >= n )
n2 = n - 1;
else
n2 = pn2;
/* n2以下の奇数での剰余が0かどうか調べる */
for( i = 3; i <= n2; i += 2 ) {
if( n % i == 0 )
return( 0 );
}
/* 素数であった */
return( 1 );
}
素数の定義通り、2以上の値で次々に割ってみて余りが出るかどうかを調べています。
素数の判定を効率よく行うため、2以外は奇数のみで割っているのと、平方根より少し大きい値までしか調べておりません。平方根についてはsqrt()というANSI標準関数がありますが、double型で計算するため64bit整数に対しては誤差が大きすぎ、平方根以下の値を返すことがあるため使用していません。
このプログラムは比較的性能が良く、1億以下の素数を約1時間で全て求めることが出来ました。しかし、より大きな値まで探索範囲を広げるには、このプログラムだけでは不十分だと思われます。
前節の単純な方法では、素数を求める事はできるものの非常に時間が掛かってしまいます。そこで、効率よく素数を求める方法としてよく知られている「エラトステネスのふるい」という方法を使ってみます。これは、自然数に対応するフラグの配列を用意し、自然数の倍数に印を付けていくというものです。最後まで印が付かなかったものが素数であることになります。
#include <stdio.h>
#define MAX_ARRAY (10000000)
char array[MAX_ARRAY];
int main( void )
{
unsigned long long n, i;
/* 配列を初期化する */
for( i = 0; i < MAX_ARRAY; i++ )
array[i] = 0;
/* 配列をふるいにかける */
for( n = 2; n <= MAX_ARRAY; n++ ) {
if( array[n] == 0 ) {
for( i = 2; i * n < MAX_ARRAY; i++ )
array[i * n] = 1;
}
}
/* ふるいで残った数は素数である */
for( n = 2; n < MAX_ARRAY; n++ ) {
if( array[n] == 0 )
printf( "%qu\n", n );
}
}
期待通り、高速に素数探索を行うことが出来ました。1000万以下の素数を求めるのに、先ほどのプログラムは2分半掛かったのですが、エラトステネスのふるいを使うと24秒で済み、約5倍の速さで計算することが出来ました。正直言って、期待した程の性能は出なかったのですが、これは探索範囲が狭いために、アルゴリズムによる差が出にくいということだと思われます。
このプログラムの欠点は探索範囲分のフラグ配列を必要とするため、メモリの許す程度の範囲しか探せないという事です。
(2005/4/18追記)
何年も前に指摘を頂いていた、ふるいに掛ける時に既にふるい落とされた数の倍数も繰返し調べるという全く無駄な処理をしていた点を修正しました。これにより更に4倍ほど処理速度が速くなりました。
メモリの許す範囲しか素数を探せないというエラトステネスのふるいの欠点を、ふるいを繰り返し使用するという方法で解決してみました。今回のプログラムでは、区間の大きさを1000万に設定してみました。もちろん区間が大きいほど効率が良くなるのですが、その為には大量のメモリが必要になってきます。
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#define MAX_ARRAY (10000000)
char array[MAX_ARRAY];
int main( int ac, char *av[] )
{
unsigned long long base, arraysize, n, min_n, max_n, i;
/* コマンドラインから素数探索範囲を決定する */
if( ac < 3 )
return;
min_n = strtouq( av[1], NULL, 10 );
if( min_n < 2 )
min_n = 2;
max_n = strtouq( av[2], NULL, 10 );
/* MAX_ARRAY分ごとの整数区間をふるいにかける */
for( base = min_n; base < max_n; base += MAX_ARRAY ) {
/* 整数区間配列の大きさを決める */
if( max_n - base + 1 > MAX_ARRAY )
arraysize = MAX_ARRAY;
else
arraysize = max_n - base + 1;
/* 配列を初期化する */
for( i = 0; i < arraysize; i++ )
array[i] = 0;
/* 配列をふるいにかける */
for( n = 2; n < ( base + arraysize ) / 2; n++ ) {
if( base <= n )
i = n * 2;
else if( base % n == 0 )
i = base;
else
i = ( base / n + 1 ) * n;
for( ; i < base + arraysize; i += n )
array[i - base] = 1;
}
/* ふるいで残った数は素数である */
for( n = 0; n < arraysize; n++ ) {
if( array[n] == 0 )
printf( "%qu\n", base + n );
}
}
}
この方法では、1億以下の素数を全て調べるのに6分程度で済みました。一番最初の単純な方法では同じ範囲を約1時間掛かりましたから、やはり探索範囲を広げるほどエラトステネスのふるいの性能がよく表れるようになっています。
これまでに求めた9999万9999以下の素数一覧データです。