素数リストの計算

no title

この記事に触発されて,大学一年の時に計算機実習が暇だったからやってた素数プログラムをまた作ってた.なんか,自分のプログラムの作り方が悪いのか2週間で終わるなんて考えられない結果になってるけど,その状況を適当に記事に.

一番楽に短く書けるからrubyのコードで書いてるけど,実際の計算はJavaかCで書いてます.

#探索範囲
max = 100
primes = Array.new()
i = 2
primes<<i
i+=1
while (i<=max)
   is_prime = true
   primes.each do |p|
      if (i % p == 0)
         # 素数ではない
         is_prime = false
         break
      end
   end
   if (is_prime)
      puts i
      primes << i
   end
   #2より大きい素数は全部奇数
   i+=2
end

というのが,ありがちなコード.元記事に載っているのもこれ(と思ったら,ルートで打ちきってたので違う).ただ,これは致命的に遅い.

問題は,素数を判定するときに,その時まで見つかっている素数全部で検査するというところ.素数の密度¥pi¥pi = n/¥ln(n)と,数が大きくなってもほとんど減らない.なので,このコードの計算量はO(n^2/¥ln(n))¥sim O(n^2)となって,致命的に遅い.というか,実は篩は意外と速くなりそうに見えてO記法で言うとln(n)で割る効果しかない.確かにたいしたコード量がいらない割に定数分確実に速くなるので悪い取引ではないけど.

ルート

素因数分解の基本だけど,何かの数で割れるか確かめるときに平方根以下の数で試せば十分である.これを入れるとかなり改善する.

#範囲
max = 100
primes = Array.new()
i = 2
primes<<i
i+=1
while (i<=max)
   is_prime = true
   primes.each do |p|
      if (p*p > i)
         break
      end
      if (i % p == 0)
         # 素数ではない
         is_prime = false
         break
      end
   end
   if (is_prime)
      puts i
      if (i*i <= max)
         primes << i
      end
   end
   #2より大きい素数は全部奇数
   i+=2
end

かけ算して平方根を超えるとループを打ちきるというコードを入れるだけで,かなり速度が改善する.つまり,O¥bigl(n^{1.5}/ln(n)¥bigr).ちなみに,素数判定に使う素数リストはメモリに載ってないと問題外だが,出力する素数を全部突っ込むとメモリ使用量が爆発するので,計算する最大値の平方根より小さい素数のみリストに追加している.

昔は,素数のリストを作るならこれで十分に速いと思ってたけど,実際に10兆まで計算しようと思うとこれではかなり非現実的になる.手元で計算して見た結果だと,

f:id:smectic_g:20100529115829p:image:h300

となって,一番速かったコードでもt = 6.8e-11¥text{s}¥times n^{1.487} = 6.8e-11 ¥times 1e13^{1.487} = 46.2¥text{year}と,なんつーか,スーパーコンピュータープリーズという結果に.

まっとうな方法

篩なんて原始時代の方法に頼らずに現代的な方法ってないの?というと.

wikipedia:ミラー-ラビン素数判定法なり,wikipedia:AKS素数判定法なりがある.これらは,多項式時間である数nを素数かどうか判定できる.つまり,これを単純に繰り返すだけでn¥cdot¥log(n)^p素数リストが作れることになる.定数値はともかくO記法的には圧倒的に速い.

で,wikipediaのミラーラビン法のコードをそのままCに移植してやってみた結果が次.ミラーラビンは非決定アルゴリズムなので,この目的にはあんまよろしくないけどまあ感じを掴むくらいなら問題ないと思う.

f:id:smectic_g:20100529160403p:image:h300

うーん,篩よりははるかに速く次数強しという感じはある(1e13までの探索で87.8日)けど,まだ絶望的に遅いなあ.どうやったら1e9までの探索が3分とか20秒とかで終わるんだろう.才能の欠如を感じる.

ということで,せいぜいがO(n^{1.5})にしかならない試し割りでは,10兆までの素数リストは埒が明かない.たぶん,現代的な多項式時間での方法でやっているんだろうけど,かなりチューニングしないと10兆は遠いので,確かに,面白いかも知れない.

というか,麻薬のように時間をとられるのでそろそろ自分はやめとこうと思う.才能無いし.