エラトステネスの篩

user-pic
0

アマゾンで何気なく本をみていたら、1000000までの素数を印刷した本があり、それなりに人気なのだという。素数好きな私は、まあそういうのも面白いと思うけど、やはりプログラムで素数を生成する方が楽しいと思った。

そこで、素数の生成アルゴリズムの定番といえば、なんといってもエラトステネスの篩である。そのためのプログラムは今まで何度も書いてきたが、ここで改めて書くからには多少なりとも工夫してみることにした。

工夫のその1は、ビット演算を行うことで、素数かどうかの情報を1ビットであらわすことにした。メモリの合理的な使用としては、これ以上の合理化はありえないはずである。

そのため、整数配列をメモリ・アロケートすることにした。すると、整数配列の1つの要素は4バイトだから、1つの要素は32ビット、すなわち32個の整数が素数かどうかを表すことができる。かりに10億までの素数をリストしようとすると、必要な配列数は10億/32 = 3.125×107程度なので、必要なメモリは、これに4を掛けて、1.25×108バイト程度である。

1.25億バイトはすごく大きいようだが、要するに125MBにすぎない。リソースモニタをみると、1500MBくらい使用可能とあるのだが、流石にそれを丸々アロケートできるとは思えない。それでも、1000MBくらいまでは使えそうで、すると320億くらいまでの素数が計算できそうであるが、それは可能性として保留しておいて、とりあえず10億までの素数をリストすることにした次第である。

計算を楽にするため、特定のビットに1を立てるためのマスクを32個用意することにした。これがあると、4バイトの要素に、特定のマスクを論理ORするだけでよい。

また、素数の情報を取り出すには、4バイトの要素を32回1ビットシフトさせて、端のビットが0かどうかその都度調べればよい。0だと、要するに素数ということである。

このような方針で作成したプログラムのソースコードは以下のとおりである。すると、10億までの素数の計算が、わずか38秒程度で完了した。それによって生成された素数は、50847535個なのだから、なかなかの生産性であるといえる。そのときのシステム仕様は、Windows 10 64ビット、Intel Core i5、2.7GHzなのだから、いまどき普通である。

ちなみに、もし300億までの素数を計算したときの時間を、O(n2)の計算量を想定して推定すると、38×30×30 = 34200秒 = 9.5時間なので、ほぼ一日のworking timeくらいコンピュータを走らせておけばよい、ということになる。

#include <stdio.h>
#include <string.h>
#include <math.h>

 unsigned int MASK[32];
 int IsPrime(int n);

int main(void)        {
 int i,j,k,n;
 int a,b;
 int N, LIM, AR;
 int *a_num;	

// MASKの用意
 n = 1;
 i = 0;
 do { MASK[i++] = n;
      n <<= 1;
 } while ( i < 32 );
 
N = 1000000000;
AR = N / 32 + 1;
LIM = int(sqrt((double) N)) + 1;

 a_num = new int[AR];
 
// 初期化

 for ( i = 0; i < AR; i++ )
  a_num[i] = 0;


//篩ビットのセット
 for ( i = 2; i < LIM; i++ )    {
   j = 2;
   if ( IsPrime(i) == 1 )  {
     while (( k = i*j) <= N )    {
     	a = k/32;
    	b = k%32;
	    a_num[a] = a_num[a] | MASK[b];
	    j++;
    }
   }
 }
 
 
k = 0;
 for ( i = 0; i < AR; i++ )                {
   if ( (a_num[i] & 1) == 1 ) k++;
   else { printf("%d\n",k++); }
   for ( j = 1; j < 32; j++ )  {
	  a_num[i] >>= 1;
	if ((a_num[i] & 1) == 1 )      {
	    k++;
	} else { printf("%d\n",k++);           }
	 if ( k > N-1 ) goto EXT;
   }
 }
EXT:

delete[] a_num; 
return 1;
}

int IsPrime(int n)     {
 int lim;
 int i = 2;
 
 lim = int(sqrt((double) n)) + 1;
 
 if ( n == 2 ) return 1;
 
 while ( i <= lim )     {
  if (( n%i) == 0 ) return 0;	
  i++;
 }
 
 return 1;

}

トラックバック(0)

トラックバックURL: http://commutative.world.coocan.jp/mt/mt-tb.cgi/4836

コメントする

月別 アーカイブ

このブログ記事について

このページは、あやたろうが2018年3月 7日 05:07に書いたブログ記事です。

ひとつ前のブログ記事は「プレーオフ、豊島一勝」です。

次のブログ記事は「ソープ・ケース」です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。