読者です 読者をやめる 読者になる 読者になる

どせいたんさき。

ナスダヨー

Octave でエラトステネスの篩

引数として与えられた数 n よりも小さな素数の列を返す関数を考える.Project Euler の問題を解くために何度も書いたのでさすがにソラで書けるようになった.

エラトステネスの篩

function lst = eratosthenessieve(n)
  sqn = sqrt(n);
  idx = ones(1,n); idx(1) = 0;
  for k =2:sqn; idx(k*k:k:n) = 0; endfor
  lst = (1:n)(idx==1);
endfunction

実行速度はこんなもん.

octave> tic(); eratosthenessieve(5e6); toc();
Elapsed time is 1.1 seconds.

偶数を省いて高速化

octave, matlab, IDL なんかではループを使わずにベクトルのまま計算したほうが速いのでとりあえず偶数を例外として処理してみる.

function lst = eratosthenessieve2(n)
  sqn = sqrt(n);
  idx = ones(1,n); idx(1) = 0;
  idx(4:2:n) = 0;
  for k =2:2:sqn; idx(k*k:k:n) = 0; endfor
  lst = (1:n)(idx==1);
endfunction

実行速度がちょっと上がった.

octave> tic(); eratosthenessieve2(5e6); toc();
Elapsed time is 0.74 seconds.

ループを分けて高速化

でかいループを回すよりもループを小分けにしたほうが効率的かもしれない.

function lst = eratosthenessieve3(n)
  sqn = sqrt(n);
  idx = ones(1,n); idx(1) = 0;
  idx(4:2:n) = 0; idx(9:3:n) = 0;
  for k =5:6:sqn; idx(k*k:k:n) = 0; endfor
  for k =7:6:sqn; idx(k*k:k:n) = 0; endfor
  lst = (1:n)(idx==1);
endfunction

実行速度がまた少し上がった.

octave> tic(); eratosthenessieve3(5e6); toc();
Elapsed time is 0.59 seconds.

おまけ:組み込みの関数 primes

octave にはまったく同じ働きをする組み込みの関数 primes がある.

octave> tic(); primes(5e6); toc();
Elapsed time is 0.15 seconds.

さすがに速い.