Javaで素数を100万個求める
LL Ringの素数の話が盛り上がっている。
http://ll.jus.or.jp/2006/blog/doukaku1
Javaで参戦しようかと思ったがネタにも走れず、かといって短さにも走れず、速さだけなら前回のたらいの2番煎じになってしまう。
芸が無いので今回は見てるだけにしようと思った。
しかし色々な関連する記事を読んでいたら既にJavaで非常に優れたプログラムを書いている人を発見できた。そしてすぐに100万個の素数を計算できてしまった。簡単に説明しようと思う。
まずRubyの界隈でメモリと実行速度の最適化が行われているのを見てその内容に魅かれてしまった。
まず最初に知ったのはMatzさんの日記。
Matzにっき(2006-06-20)
こちらで紹介されているid:sumimさんの以下の記事を読んだ。
エラトステネスのふるいをコンパクトにする方法 - Smalltalkのtは小文字です
SmallTalkは全く読めないので良くわからなかったのだが、Sumimさんが紹介されているJohn Moyerさんの記事を読んでみてその内容のレベルの高さに驚いてしまった。
Print Prime Numbers
Johnさんは綿密な数学の理論とビット演算を用いることによりふるいに必要なメモリを圧縮している。
Sumimさんの説明をそのまま引用すれば、
換言すれば、８ビット（８個のフラグ）、つまり１バイトあれば、３０個分の整数の素数判定状況を保持することができるわけです。すげー。
実際にJohnさんの以下のC言語のソースを読むと、上の理論を拡張し、1バイトあたり38.5個までに効率を向上している。
C Source Code for a Sieve Program
/* expand scheme from 30 to 2310 30 == 2*3*5 and 8 == (2-1)*(3-1)*(5-1) 2*3*5*7*11 == 2310 1*2*4*6*10 == 480 bits == 60 bytes ==> prime or not prime for 38.5 integers per byte In each 2310 integers, for N >= 1, the numbers that might be prime are: N*2310+1, N*2310+13, N*2310+17, . . . N*2310+13*13, N*2310+13*17, . . . N*2310+2309 */
さてJohnさんはこれだけでは終わらず、ある整数Nより大きい10個の素数を表示するJavaアプレットを作成している。
http://www.rsok.com/~jrm/source/java/
このプログラムが凄いのはアルゴリズムがふるいではなく、ある数Nの次の素数を直ぐに計算できることだ。
このプログラムでは整数nがn < 341,550,071,728,321であるならば完全に正確な素数を算出できる。今回の問題の条件には楽勝で合致する。
アルゴリズムはOlivier Langloisさんによる以下のものだそうだ。
http://www.utm.edu/research/primes/prove/prove2_3.html
さすがにこのレベルになると難しくてまだ理解しきれていない。
しかしJavaのプログラムを読んでみると、少し改造すればすぐに100万個の素数を求めることが確認できた。
以下が100万個の素数を求めたプログラムである。JohnさんのNext_10_Primesを開始０、個数100万個に変えただけだ。
public class NextMillionPrimes {
public static void main(String[] args) {
int i;
long x = 0;
long y = 0;
for (i = 0; i < 1000000; i++) {
y = LFindNextPrime.nextprime(x);
System.out.println(y);
x = y;
}
}
}ただし、Johnさんのプログラムは数学的に大きな数にも対応できるようにBigDecimalで記述されている。
このせいで100万個を計算するのに思ったよりも非常に長い時間がかかってしまう。*1
このため少しでも高速にするためJohnさんのプログラムを全てlongで書き直した。*2
書き直すときにソースのB_をL_に直しているので上のソースを使うときには注意が必要だ。
それでもJohnさんのプログラムは多量のループで使うことを前提としていない。FindNextPrimeは毎回nがどのような値であるか大量のif文で確認しているし、あらかじめnの大きさで場合分けをしていないのでこのへんは高速化できると思う。
そのへんを直した上で、さらにlongでなくintにしてより高速にして計りたかったがもう時間がないのであきらめた。
実行結果は以下のとおりだ。
$ time java NextMillionPrimes | tail 15485761 15485773 15485783 15485801 15485807 15485837 15485843 15485849 15485857 15485863 real 0m34.843s user 0m5.217s sys 0m17.202s
