遅延評価によるエラトステネスの篩のお話

Programming

http://d.hatena.ne.jp/laciryl/20100712#1278918861 のお話。

対象となるのは以下の有名なHaskellコード。

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

これをScalaで書く。

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

// def sieve(x:Stream[Int]):Stream[Int] = Stream.cons(x.head,
//     sieve(x.tail.filter(i => i%x.head!=0)))
// lazy val primes = sieve(Stream.from(2))

// With debug output
def sieve(x:Stream[Int]):Stream[Int] = Stream.cons(x.head,
    sieve(x.tail.filter(i => {
        println(i + " % " + x.head)
        i%x.head!=0
    })))
lazy val primes = sieve(Stream.from(2))

primes.foreach{ i =>
    if(i>=100) exit(0)
    println(i)
}

Scalaは副作用を許容するので、剰余計算がいつ行われているかを出力するようにしている。

2
3 % 2
3
4 % 2
5 % 2
5 % 3
5
6 % 2
7 % 2
7 % 3
7 % 5
7
8 % 2
9 % 2
9 % 3
10 % 2
11 % 2
11 % 3
11 % 5
11 % 7
11
12 % 2
13 % 2
13 % 3
13 % 5
13 % 7
13 % 11
13

まあ、見てどのような計算が行われているかを推測してほしいが、これと同じ手順を踏む計算はC++で以下のように書ける。

// This Haskell Code in C++
// primes = sieve [2..]
// sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p /= 0]
#include <cstdio>
#include <vector>
using namespace std;

int main(int argc, char **argv) {
    vector<int> primes;
    for(int i = 2; ; i++) {
        bool f = true;
        for(int j = 0; j < (int)primes.size(); j++) {
            printf("%d %% %d\n", i, primes[j]);
            if(i%primes[j]==0) {
                f = false; break;
            }
        }
        if(f) {
            if(i>=100)break;
            printf("%d\n", i);
            primes.push_back(i);
        }
    }
    return 0;
}

これは全ての値を試し割りしていて遅いらしい。素数のみについて考えてみると、素数定理より素数はn/log(n)個あると考えられるが、これらの素数はそれぞれそれより小さい素数全てで割られるので、厳密な計算ではないが、(n/log(n))^2くらいにはなると思われる。だからオーダーはO(n^2)付近じゃないかと予想される。

さて、これに対する解決方法として提案されている
http://haskell.g.hatena.ne.jp/nobsun/20060618
のコードをC++に直すと以下のようになると考えられる。

// This Haskell Code in C++
// primes = sieve [2..]
// sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p /= 0]
#include <cstdio>
#include <vector>
using namespace std;

int main(int argc, char **argv) {
    vector<int> sievers;
    int siever_index = 0;
    sievers.push_back(2);
    int i = 2;
    int print_i = 0;
    while(true) {
        int siever = sievers[siever_index];
        for(; i < siever*siever; i++) {
            bool f = true;
            for(int j = 0; j < siever_index; j++) {
                printf("%d %% %d\n", i, sievers[j]);
                if(i%sievers[j]==0) {
                    f = false; break;
                }
            }
            if(f) {
                sievers.push_back(i);
            }
        }
        for(; print_i < (int)sievers.size(); print_i++) {
            if(sievers[print_i] >= 100) return 0;
            printf("%d\n", sievers[print_i]);
        }
        siever_index++;
    }
    return 0;
}

これの出力は以下のようである。

2
2
3
4 % 2
5 % 2
5 % 2
6 % 2
7 % 2
7 % 2
8 % 2
5
7
9 % 2
9 % 2
9 % 3
10 % 2
11 % 2
11 % 2

これは確かに速そうだが、はっきり言って、これも一般的にRAMマシンと等価なマシン上で用いられるエラトステネスの篩とは違うものだと思う。

C++など配列を保有する言語では、配列のランダムアクセス性を生かした篩を行うし、剰余計算は現れない。

#include <cstdio>
// common sieve

int main(int argc, char **argv) {
    int siever[100] = {0};
    for(int i = 2; i < 100; i++) {
        if(siever[i])continue;
        printf("%d\n", i);
        for(int j = i*i; j < 100; j+=i) {
            siever[j]=1;
        }
    }
    return 0;
}

これと可能な限り等価なものをHaskellで書くのであれば、永続配列などのデータ構造を用いて書くしかないのではないかと思う。
Haskellのことはよくわからないので、僕の仕事はここまで。