アルゴリズムの素晴らしさに気付かせてくれたのはエラトステネスの篩だった

Haskellによるエラトステネスの篩(sieve of Eratosthenes)の美しい実装を見て、初めてアルゴリズムの素晴らしさに気付かせてくれたのがエラトステネスの篩だったことを思い出した。たしか中学生の頃だ。そこで、当時を懐かしみながら簡単な実装を書き留めておくことにした。とりあえず、Haskell, C++, Pythonの実装を以下に示す。コードは比較的短いが、実行効率を優先させているわけではない。

Haskell

これはHaskellのサンプルコードでよく出てくる無限リストと遅延評価による実装だが、とても分かりやすいし、美しいコードだと思う。

primes = sieve [2..] sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p /= 0]

以下のように関数takeを使って必要な分だけ素数を取り出すことができる。ただし、実行効率はあまり良くない。

take 10 primes [2,3,5,7,11,13,17,19,23,29]

C++

次にC++を使って実装してみた。ここではC++0xを使っている。そのほうが簡潔で分かりやすく書けるし、これからはC++0xがより使われ、普及して欲しいという意味もある。因みにgccではバージョン4.5以降、Intelコンパイラではバージョン11.0以降、Visual C++では2010(16.0)以降でコンパイルできる。

std::vector<int> primes; for (int i = 3; i < 100; i += 2) primes.push_back(i); auto end = primes.end(); for (auto x = primes.begin(); *x * *x <= *(end-1); ++x) end = std::remove_if(x + 1, end, [&x](int p){ return p % *x == 0; }); primes.erase(end, primes.end()); for (auto p = primes.begin(); p != primes.end(); ++p) std::cout << *p << " "; std::cout << std::endl;

上記のコードでは3から100までの素数を表示する。ノートPC (Core 2 Duo T9800 2.93 GHz / Visual Studio 2010)を使って1000万までの素数を生成してファイルに書き出すのにかかった時間は約2秒だった。

そういえば以前にC++のテンプレートで素数を求めたことがあったな。エラトステネスの篩ではなかったけど。

Python

最後にPythonによる実装を示す。効率はあまり良くないかも。

primes = range(3, 100, 2) for i in range(len(primes)): if primes[i]**2 > primes[-1]: break primes[i+1:] = filter(lambda p: p % primes[i] != 0, primes[i+1:]) print primes

これをそのままワンライナーに。

python -c "import sys;globals().__setitem__('P',range(3,100,2));[P.__setslice__(i+1,len(P),filter(lambda p:p%P[i]!=0,P[i+1:]))for i in range(len(P))if i<len(P)and P[i]**2<=P[-1]];[sys.stdout.write('%d '%p)for p in P]"

上記のコードではそれぞれ、赤字の数値が求める素数の上限値となっている。