エラトステネスの篩

http://haskell.g.hatena.ne.jp/route150/20101224/1293122299を見て気になったので、エラトステネスの篩を書いてみた。

n以下の素数の和を表示する。

import Control.Monad
import Data.Array.ST
import Data.Array.Unboxed
import Data.Array.Base(unsafeRead, unsafeWrite, unsafeAt)
import System.Environment

sieve :: Int -> UArray Int Bool
sieve n = runSTUArray $ do
  arr <- newArray (0, n) True
  let end = floor $ sqrt (fromIntegral n)
  writeArray arr 0 False
  writeArray arr 1 False
  loop arr end 2
  return arr
  where
    loop arr end k
      | k > end = return ()
      | otherwise = do
          v <- unsafeRead arr k
          when v $ remove_multiples arr k (2 * k)
          loop arr end (k + 1)
    remove_multiples arr p k
      | k > n = return ()
      | otherwise = do
          unsafeWrite arr k False
          remove_multiples arr p (k + p)

total :: Int -> UArray Int Bool -> Int
total n arr = loop 2 0
  where
    loop k acc
      | k > n = acc
      | arr `unsafeAt` k = loop (k + 1) (acc + k)
      | otherwise = loop (k + 1) acc

main = do
  [n] <- fmap (map read) getArgs
  print $ total n $ sieve n

# include <stdio.h>
# include <math.h>
# include <stdlib.h>

int main(int argc, char **argv)
{
  if(argc <= 1) return 1;
  int n = atoi(argv[1]);
  bool *arr = new bool[n+1];

  arr[0] = false;
  arr[1] = false;

  for(int i = 2; i <= n; i++)
    arr[i] = true;

  int end = (int)(sqrt(n));
  for(int i = 2; i <= end; i++)
    if(arr[i])
      for(int j = 2 * i; j <= n; j += i)
        arr[j] = false;

  long sum = 0;
  for(int i = 2; i <= n; i++)
    if(arr[i])
      sum += i;
  printf("%ld\n", sum);
}

どちらもコンパイルオプションは-O2。

n=1億で実行。

> time ./eratosthenes-c++ 100000000
279209790387276
./eratosthenes-c++ 100000000  2.73s user 0.05s system 99% cpu 2.788 total
> time ./eratosthenes 100000000
279209790387276
./eratosthenes 100000000  1.87s user 0.01s system 99% cpu 1.889 total

Haskellの方が30%ほど速い。たぶんSTUArray s Boolが一要素1ビットなのが効いている。

追記

bool[]の代わりにstd::vectorを使ってもまだHaskellに及ばない。謎。

> time ./eratosthenes-c++-vector 100000000
279209790387276
./eratosthenes-c++-vector 100000000  2.05s user 0.00s system 99% cpu 2.061 total

追追記

別に謎ではなかった。最初にtrueを代入するのをやめて「std::vector arr(n+1, true);」と初期化すれば良かった。

> time ./eratosthenes-c++-vector 100000000
279209790387276
./eratosthenes-c++-vector 100000000  1.85s user 0.01s system 99% cpu 1.857 total