エラトステネスの篩 in Java - 似非学問的な手記

「落ちつけ！素数を数えるんだ！」

というセリフが思い浮かんだのは俺だけでいい。
あれって元ネタは何なんだろう。

Wikipedia「エラトステネスの篩」
http://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%A9%E3%83%88%E3%82%B9%E3%83%86%E3%83%8D%E3%82%B9%E3%81%AE%E7%AF%A9

おそらく最古のアルゴリズムなのでは？と思う。
別に歴史を詳しく知ってるわけじゃないし、なんの裏付けも用意してないけど。

さて、この「篩」のコードは、確かC++の勉強をしているときに練習で書いたような気がする。このブログ始めるよりずっと前のことだけど。

その時は効率のいい実装を自分で調べて書いたんですが、今回そのソースのページが見つからんかった。
おぼろげな記憶を頼りに書いたのが次のコード。
1.１〜ｎまでのビットフラグを用意して、２から検索していく。
2.ビットが立っているのが素数である。
3.素数を見つけたら以降のその倍数のビットを全て降ろす。

例）（o:ビットが立っている,x:ビットが降りている）
2o,3o,4o,5o,6o,7o,8o,9o,10o,...
→2o,3o,4x,5o,6x,7o,8x,9o,10x,...(primes:2)
→2o,3o,4x,5o,6x,7o,8x,9x,10x,...(primes:2,3)
→2o,3o,4x,5o,6x,7o,8x,9x,10x,...(primes:2,3,5)
→...

コードは以下の通り。
出力は、引数としてコマンドラインで与えられた自然数以下の素数の個数を表示。素数全部を表示してもよかったけど、引数がデカいときの結果を知りたかったんでそのような仕様にしました。

あと、せっかくなので、素数の個数の近似値を与えるx/ln(x)を計算して表示するようにした。

sieve.java

import java.util.*;
class sieve{
	public static void main(String args[  ]){
		try{
			//サイズ保存

			int size=Integer.parseInt(args[ 0 ]);
			
			//初期化を2ビットづつ行う為、奇数が来るとArrayIndexOutOfBoundsException例外が発生する。
			//その防止に、1ビット多めにメモリを確保する。
			boolean num[  ]=new boolean[ size+2 ];
			
			//例外的に小さい引数を与えられた時対策。

			if(size<2){
				System.out.println("0 primes found.");
			}
			if(size==2){
				System.out.println("1 primes found.");
				System.out.println("π(2)="+(2/Math.log(2)));
			}
			if(size>2){
				//タイムカウントスタート。
				double start = System.currentTimeMillis();
				
				//素数カウンタ。

				int count=1;
				
				//初期化のターン。
				//高速化のため、篩にかける前に初期化の段階で偶数番目のビットを降ろしておく
				num[ 0 ]=num[ 1 ]=false;
				for(int i=2;i<=size;i+=2){
					num[ i ]=false;
					num[ i+1 ]=true;
				}
				
				//篩部分。

				//i=3から開始。
				//num[ n ]がtrue => nが素数
				//後で使うので、for文用変数nはfor文の外で宣言。

				int n;
				int sqrt_size=(int)(Math.sqrt(size)+1);
				for(n=3;n<sqrt_size;n+=2){
					//nが素数なら、nの倍数のビットを降ろす
					if(num[ n ]){
						if(size>=n*n){
							for(int j=n*n;j<size;j+=(n*2) ){
								num[ j ]=false;
							}
						}
						count++;
					}
				}
				//nが偶数で前半のforループを抜けたとき、強制的に奇数にする。

				if( (n%2)==0){
					if(num[ n ])count++;
					n++;
				}
				for(;n<=size;n+=2){
					if(num[ n ])count++;
				}
				//タイムカウントストップ。

				double stop = System.currentTimeMillis();
				System.out.println(count+" primes found.");
				System.out.println("π("+size+")="+(size/Math.log(size)));
				System.out.println("Time:"+(stop-start)+"[ msec ]");
			}
		}catch(Exception e){e.printStackTrace();}
	}
}

引数をnとすると、
(1)倍数消去は√n以下までの素数すべてに対して行うことが必要十分。
(2)初期化の段階で偶数は全て消しておく。
(3)素数pが見つかったら、次に現れる最初のpの倍数はpの二乗で、次から可能性があるものとしてp*(p+2),p*(p+4),...と続く。

などなど、高速化する手段は色々盛り込みました。
特に、偶数を徹底的に無視することでかなり速くなります。

ですがこのコードだと、勿論「篩の大きさ」は、ブール値配列が確保できる大きさまでが限界なので、当然そこまで大きい数には適用できません。

一応、目標として「一億」までを引数に取れるようにしたかったのですが、実験したところ、

java sieve 1000
168 primes found.
π(1000)=144.764...
Time:0.0[msec]

java sieve 10000
1229 primes found.
π(10000)=1085.736...
Time:3.0[msec]

java sieve 100000
9592 primes found.
π(100000)=8685.889...
Time:6.0[msec]

java sieve 1000000
78498 primes found.
π(1000000)=72382.413...
Time:18.0[msec]

java sieve 10000000
664579 primes found.
π(10000000)=620420.688...
Time:258.0[msec]

java sieve 100000000
Exception in thread "main" java.lang.OutOfMemoryError...

o...rz

どうやら、60000000程度が限界のようです。

ですが、「あと２倍評価できるようになれば……」と考えたところで、どうせ偶数番目ビットは使ってないんだから最初からn番目のビットに奇数2n+1を対応させるようにしておけば、２分の１のメモリで済むじゃんってことに気付いた。

で、改良したのが次のコード。
アルゴリズムは全く同じだが、今述べた手法でメモリ量を節約している。

sieve2.java

import java.util.*;
class sieve2{
	public static void main(String args[  ]){
		try{
			//サイズ保存

			int size=Integer.parseInt(args[ 0 ]);
			
			//例外的に小さい引数を与えられた時対策。
			if(size<2){
				System.out.println("0 primes found.");
			}
			if(size==2){
				System.out.println("1 primes found.");
				System.out.println("π(2)="+(2/Math.log(2)));
			}
			
			if(size>2){
				int half_size=( (size-1)/2)+1;
				//num[ n ]には命題「2n+1が素数である」の真偽が入る

				boolean num[  ]=new boolean[ half_size ];
				
				//タイムカウントスタート。
				double start = System.currentTimeMillis();
				
				//素数カウンタ。

				int count=1;
				
				//初期化のターン。
				num[ 0 ]=false;
				for(int i=1;i<half_size;i++){
					num[ i ]=true;
				}
				
				//篩部分。

				//n=1(2n+1=3)から開始。
				//num[ n ]がtrue => 2n+1が素数
				//後で使うので、for文用変数nはfor文の外で宣言。

				int n;
				int sqrt_size=(int)(Math.sqrt(size)+1);
				for(n=1;n<sqrt_size;n++){
					//2n+1が素数なら、2n+1の倍数のビットを降ろす
					if(num[ n ]){
						if(half_size>=(2*n*(n+1))){
							for(int j=(2*n*(n+1));j<half_size;j+=(2*n+1)){
								num[ j ]=false;
							}
						}
						count++;
					}
				}
				for(;n<half_size;n++){
					if(num[ n ])count++;
				}
				//タイムカウントストップ。

				double stop = System.currentTimeMillis();
				System.out.println(count+" primes found.");
				System.out.println("π("+size+")="+(size/Math.log(size)));
				System.out.println("Time:"+(stop-start)+"[ msec ]");
			}
		}catch(Exception e){e.printStackTrace();}
	}
}

で、実行。

java sieve2 10000000
664579 primes found.
π(10000000)=620420.688...
Time:176.0[msec]

java sieve2 100000000
5761455 primes found.
π(100000000)=5428681.023...
Time:2657.0[msec]

成功。
ちょっぴり速くもなりました。

このあと、３の倍数も初期化の時点で省く方法を試して高速化したコードも書いたんですが、そちらはイマイチ効果が上がらなかったので省略。

こうなると１０億までに挑戦したくもなりますが、さすがにそれには大改造がいると思うんでまたの機会ということに。

それでは（・ω・）ﾉｼ