smectic_gの日記

2010-05-29

[]素数リストの計算

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兆は遠いので,確かに,面白いかも知れない.

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

2008-02-11

[]ほんの値段って一体

Matzにっき(2008-02-11)

去年この本をブックオフを2件まわって*1200円くらいで売り払った自分涙目.下手すると今ごろシュレッダーだと思うけど.

14800円ってどんだけ.

*1:最初の店では買い取りを拒否された

2007-12-10

[]Dynamic Creator続き

わかっている係数の組み合わせから乱数を作るgenmtrandを移植するのはMersenneTwisterとほぼ一緒ですぐに終わった.

ただ,係数の組み合わせを求めるプログラムの移植が非常に難儀.所々,Cだと当然だけど,Javaにそのまま移すのはどうだろう的なプログラミングをしてて,そこら辺をオブジェクト指向的に書き直していくのが何をやってるのかが場所によってはちんぷんかんぷんなので非常につらい.

でも,なんだかんだ言いつつ,eqdeg.c以外の移植は(テストはしてないけど,eclipseで警告が出ない程度には)完了.

eqdeg.cの作業で力尽きてきたので,C的な記述をそのまま移植するか,係数の計算はCでやってJavaでは乱数の生成だけをやらせればいいじゃないと割り切るかの二択で,後者を選択.現在テスト中.

ところで,イサコの4423ってメルセンヌ数(正確には2^4423-1がメルセンヌ数)なのね.意図的かな.

2007-12-06

[]Javaでunsignedがないのは不便?

Javaで符号無し整数がないと移植の時に困ると言うのを突如認識するはめに.

Mersenne Twister由来の並列対応エンジンであるDynamic CreatorをJavaに移植しようとソースコードを除いてみたら,当然のごとくunsigned intのあらし.元のMersenne Twisterからが基本unsigned intを使用してる.で,頭が痛いのはこれをどうやってJavaで実装するか.

参考にするために,Javaに移植されたMersenne Twisterのコードを眺めると(colt, 検索で引っかかるMersenneTwisterFast,奥村先生他著のJavaによるアルゴリズム辞典),どれも符号ありと無しの区別に関しては無視を決め込んでいる.これで問題がなければいいのだけど,少なくとも無視を決め込むと同じ種で初期化しても違う乱数が生成されることになる(実際に,mt19937ar.cとcern.jet.random.engine.MersenneTwister.javaで比較してみると全然違う).

これがどの程度の問題になるのかが見えてこないのだけど(世間的に許容されていると言うことは乱数としての統計には問題はないのだろう),Dynamic Creatorの場合にはどこかでmt19937.cで生成した乱数を使って種を決めているような気配があるので,やばそうだし,それ以上に複雑でかつ時間のかかる計算をしているのでunsigned intを無視すると手ひどいことになりそう.

追記

といっても,unsignedやsignedといっても,しょせん補数表示なわけだからビットとしては右シフトをきちんと論理シフトにすることさえ徹底すれば問題なさげなんだが,何で結果が違ってくるんだろ.

ちょっと計算し直し.

追記2

coltはビットシフトが全部算術シフトなのでおかしい.MersenneTwisterFastはちゃんと論理シフトになってるけど,初期化のアルゴリズムが古い.Javaによるアルゴリズム辞典はちゃんと論理シフトになってるし,初期化のアルゴリズムも最新版なのでOK.coltは他にもぞっとするようなコードが所々あって困る

2007-02-01

[]靴下を干す時

色違いの靴下を干そうとするとすごい間抜けなことになる.

ぬれた洋服の色って見づらいから,二つを並べないと違いがわからない.コード(Java)で書くと

class Socks {
  private final int color;
  private boolean isDamaged;

  public boolean isSameColor(Sock otherSocks) {
    return ( color == otherSocks.color);
  }

  public boolean isDamaged() {
   return isDamaged;
  }
}

しかインターフェースがないSocksクラスが入ったスタックを,

class Pair {
  private Socks a;
  private Socks b;

  Pair(Socks a, Socks b) {
    this.a = a;
    this.b = b;
  }
}

のリストにまとめるということになる.今どうやってるかと言うと,

public void hangSocks(Stack<Socks> items, List<Pair> hanger) {
  Stack temp = new Stack();

  while (!items.empty() ) {
    Socks socksA = items.pop();
    if (socksA.isDamaged()) {
      continue;  // 痛んでたら捨てて次のをとる
    }
    while (!items.empty()) {
       Socks socksB = items.pop();
       if (socksB.isDamaged() ) {
         continue;
       } else if (socksA.isSameColor(socksB)) {
         hanger.add(new Pair(socksA, socksB); // たまたま同じ色だったらハンガーにかける
         break;
       } else {
         temp.push(socksB); // 違う色だったらひとまず別のところに置く
       }
    }
    // このループを経てハンガーにかけられなかったら靴下は半端だから自動的に捨てる
    // 一時的においた靴下を元のスタックに戻す.
    for (Socks socks : temp) {
       items.push(socks);
    }
  }
}

というふうに,非常に場当たり的なかけ方をしているので死ぬほど効率が悪い.O(n^2)のはず.

最初に同じ色の靴下のリストを作ってしまう

public void hangSocks(Stack<Socks> items, List<Pair> hanger) {
  List<Socks> samples = new ArrayList<Socks>();
  List<Stack<Socks>> temp = new ArrayList<Stack<Socks>>();
  while (!items.empty()) {
    Socks socks = items.pop();
    if (socks.isDamaged()) continue;
    int i;
    for (i=0; i < samples.size(); i++) {
      if (socks.isSameColor(samples.get(i))) {
        temp.get(i).add(socks);
        break;
      }
    }
    if (i == samples.size() ) {
      samples.add(socks);
      Stack stack = new Stack();
      stack.push(socks);
      temp.add(socks);
    }
  }

  for (Stack<Socks> list : temp) {
    Socks stock;
    for (Socks socks : list) {
      if (stock == null) {
        stock = socks;
      } else {
        hanger.add(new Pair(stock, socks));
        stock = null;
      }
    }
  }
}

は色の数が少なければだいぶ効率的だが,メモリ消費量の点で採用していないのかもしれない.コードの長さ自体も長いし.

最初におもむろにPairを組んでhangerに突っ込んだあとにswapして揃える形式で効率の良いアルゴリズムがあればいいなあと思うが厳しいのかもしれない.