2006年06月17日 21:30 [Edit]

LLR2006 - 1,000,000(番目|まで)の素数

LLR2006
キミならどう書く 2.0 - ROUND 1 - ? Lightweight Language Ring
お題は「100までの整数から素数を列挙せよ」です.

に対して

mputの日記。 - キミならどう書く 2.0 - ROUND 1 -
100までではちょっと上が小さすぎる。「最初の1,000,000個」とかに変更すべき。ここまで大きければHaskellでも素朴なsieveでは表示できなくなる*1ので、腕の見せ所となる。

というツッコミが来たので、Haskell記事がunder constructionということもあってやってみた。

Javascriptを追記


ただし、他との整合性から、「最初の1,000,000個」ではなく、「1,000,000までの素数」にしてある。この程度であれば、特に工夫を凝らさなくても...

エラトステネスのふるい @ Ruby - 32nd diary (2006-06-16)
追記.ええと,「最初の1,000,000個」とかに変更すべき。という意見もあるようですが,少人数で腕を競いたいわけではなく,みんなで盛り上がろうという趣旨ですので100までということになりました.ご理解ください.

...皆で楽しむことが出来るので。また、以下のソースコードを「最初の1,000,000個」に変更するのはTrivialだ。

perl 5.8.8

use strict;
use warnings;
my @primes;
sub is_prime{
    my $n = shift;
    for my $d (@primes){
        last if $d*$d > $n;
        return unless $n % $d;
    }
    push @primes, $n;
    return $n;
}
my $max = shift || 100;
is_prime($_) for (2 .. $max);
print join(" ", @primes), "\n";
Mac OS X 10.4.6 @ MacBook Pro Intel Core Duo 2.0GHz
        8.76 real         8.08 user         0.02 sys
         0  maximum resident set size
         0  average shared memory size
         0  average unshared data size
         0  average unshared stack size
         0  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
         0  signals received
       521  voluntary context switches
         0  involuntary context switches

MacOS Xの/usr/bin/time -lは見てのとおりあまり参考にならないので、以下-lなしのみ。全部521 voluntary context switchesってのもなんだかなあ。

FreeBSD 6.1 @ Dual Xeon 2.8MHz

        9.96 real         9.84 user         0.03 sys
      6912  maximum resident set size
         7  average shared memory size
      1492  average unshared data size
       127  average unshared stack size
      1419  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
         0  signals received
       406  voluntary context switches
      1328  involuntary context switches

python 2.4.3

pythonのjoinの仕様は毒々しい。

import sys
primes = []
def is_prime(n):
    for d in primes:
        if d * d > n:
            break
        if n % d == 0:
            return 0 
    primes.append(n)
    return n
max = 100
if len(sys.argv) == 2:
    max = int(sys.argv[1])
for i in range(2,max):
    is_prime(i)
# not primes.join(" ") -- I love to hate python!
# not " ".join(primes)
print " ".join(map(str, primes))
Mac OS X 10.4.6 @ MacBook Pro Intel Core Duo 2.0GHz
       10.30 real        10.21 user         0.05 sys
FreeBSD 6.1 @ Dual Xeon 2.8MHz

        9.07 real         8.71 user         0.09 sys
     20156  maximum resident set size
       657  average shared memory size
     17582  average unshared data size
       128  average unshared stack size
      4843  page reclaims
        40  page faults
         0  swaps
        17  block input operations
         0  block output operations
         0  messages sent
         0  messages received
         0  signals received
       470  voluntary context switches
      1198  involuntary context switches

ruby 1.8.4

コードはきれいだが、ちょっと遅め。OS Xだと高速なのはなぜ?

@primes = []
def is_prime(n)
    @primes.each do |d|
        break if d * d > n
        return 0 if n % d == 0
    end
    @primes.push(n)
    return n
end
max = ARGV.size == 1 ? ARGV[0].to_i : 100
(2 .. max).each{|i| is_prime(i) }
puts @primes.join(" ")
Mac OS X 10.4.6 @ MacBook Pro Intel Core Duo 2.0GHz
       13.65 real        13.49 user         0.02 sys
FreeBSD 6.1 @ Dual Xeon 2.8MHz
       24.23 real        23.60 user         0.02 sys
      3036  maximum resident set size
         4  average shared memory size
       885  average unshared data size
       128  average unshared stack size
       551  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
         0  signals received
       415  voluntary context switches
      2691  involuntary context switches

ANSI C

やはりぶっちぎりの速さ!

他との整合性を保つためになんちゃってjoin()を実装して使っているが、モノホンが欲しいなあ。

#include <stdio.h>
#include <string.h>
#define PRIMEBUFSIZ 1024*1024
#define CHARBUFSIZ  1024*1024

static int nprimes = 0;
static int primes[PRIMEBUFSIZ];

int is_prime(int n){
    int i, d, p;
    for (i = 0; i < nprimes; i++){
        p = primes[i];
        if (p * p > n)  break;
        if (n % p == 0) return 0;
    }
    primes[nprimes++] = n;
    return n;
}

char *join(char *dst, int *array, int items){
    char dbuf[16];
    int i;
    for (i = 0; i < items; i++){
        sprintf(dbuf, "%d ", array[i]);
        strncpy(dst, dbuf, 16);
        dst += strlen(dbuf);
    }
    return dst;
}
int main(int argc, char **argv){
    int i;
    int n = 100;
    char obuf[CHARBUFSIZ], *tail;
    if (argc > 1) n = atoi(argv[1]);
    for (i = 2; i <= n; i++) is_prime(i);
    tail = join(obuf, primes, nprimes);
    *tail = '\n';
    fputs(obuf, stdout);
}
Mac OS X 10.4.6 @ MacBook Pro Intel Core Duo 2.0GHz
        0.34 real         0.29 user         0.01 sys
FreeBSD 6.1 @ Dual Xeon 2.8MHz
        0.33 real         0.24 user         0.02 sys
      1404  maximum resident set size
         4  average shared memory size
       559  average unshared data size
      1219  average unshared stack size
       267  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
         0  signals received
       407  voluntary context switches
        87  involuntary context switches

Haskell (GHC 6.4.1/6.5-Beta on OS X)

そして、問題のHaskell。まずNaiveな実装。

import System.Environment

main = do args <- getArgs
          let n =  if length args == 1 then read $ args !! 0  else 100
          putStrLn $ unwords $ map show $ takeWhile (\x -> x <=  n) primes

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

やってみればわかるが、1,000,000までの素数とかだと、遅くて話にならない。これは、n-1以下の全素数をない方表記、じゃない内包表記中で使っているため。ちゃんと「sqrt(n)以下の素数までチェックすればOK」という条件をなんとか実装したらどうなるか?

import System.Environment

main = do args <- getArgs
          let n =  if length args == 1 then read $ args !! 0  else 100
          putStrLn $ unwords $ map show $ takeWhile (\x -> x <= n) primes

primes = 2 :  filter sieve [3..]
    where
         sieve x = all primeTo $ takeWhile lessThanSquare primes
                where  primeTo y  = mod x y /= 0 
                       lessThanSquare z = z * z <= x
Mac OS X 10.4.6 @ MacBook Pro Intel Core Duo 2.0GHz
        2.94 real         2.91 user         0.02 sys
FreeBSD 6.1 @ Dual Xeon 2.8MHz
        6.27 real         6.26 user         0.00 sys
      5828  maximum resident set size
       301  average shared memory size
        36  average unshared data size
       128  average unshared stack size
      1213  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
       313  signals received
         1  voluntary context switches
        76  involuntary context switches

OS Xが高速なのは、6.5 Betaを使っているからだろうか?

1,000,000個目の素数は15,485,863

ちなみに、takeWhile (\x -> x <= n) primestake n primesに書き換えるだけで、このscriptは「素数1,000,000個」を表示するようになる。

FreeBSD 6.1 @ Dual Xeon 2.8MHz
      216.53 real       215.31 user         0.36 sys
     58772  maximum resident set size
       300  average shared memory size
        36  average unshared data size
       128  average unshared stack size
     14437  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
     10771  signals received
         1  voluntary context switches
      7929  involuntary context switches

そんな悪くないではないか。ちなみにCではこうなった。高速化のため、なんちゃってjoin()は使っていない。mainのみ示す。

int main(int argc, char **argv){
    int i;
    int n = 100;
    if (argc > 1) n = atoi(argv[1]);
    for (i = 2; nprimes < n; i++) is_prime(i);
    for (i = 0; i < n; i++){
        printf("%d ", primes[i]);
    }
    printf("\b\n"); /* \b is such a dirty hack */
}
FreeBSD 6.1 @ Dual Xeon 2.8MHz
       11.01 real        10.98 user         0.02 sys
      4504  maximum resident set size
         4  average shared memory size
      4100  average unshared data size
       128  average unshared stack size
      1036  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
         0  signals received
         1  voluntary context switches
       557  involuntary context switches

これならPerlでも行けそうだ。やはりmain相当の部分のみ示す。

for (my $n = 2; @primes < $max; $n++){
    print $n, " " if is_prime($n)
}
print "\b\n";
FreeBSD 6.1 @ Dual Xeon 2.8MHz
      345.64 real       344.65 user         0.13 sys
     26000  maximum resident set size
         7  average shared memory size
     14892  average unshared data size
       127  average unshared stack size
      6145  page reclaims
         0  page faults
         0  swaps
         0  block input operations
        62  block output operations
         0  messages sent
         0  messages received
         0  signals received
         6  voluntary context switches
     30162  involuntary context switches

ざっとこんなところかな。

Dan the Prime Lightweight Linguist

Javascript

追記した。やはりbrowserでdemoを見れるのはいい。

Javascriptなので、Prototypeとして作ってある。

function Prime(){
  this.primes = [];
  this.init   = function() {
    this.primes = [];
  };
  this.is_prime = function(n){
    for (var i = 0; i < this.primes.length; i++){
      var d = this.primes[i];
      if (d * d > n)  break;
      if (n % d == 0) return 0;
    }
    this.primes.push(n);
    return n;
  };
  this.find_primes = function(n){
    this.init();
    for (i = 2; i <= n; i++) this.is_prime(i);
    return this.primes
  };
}
までの素数:

2 3 5 7


この記事へのトラックバックURL

この記事へのトラックバック
キミならどう書く 2.0 - ROUND 1 - LL Ring の前哨戦として...
Brainf*ckで100までの素数を列挙してみるテスト【TAKESAKO @ Yet another Cybozu Labs】at 2006年06月26日 16:17
キミならどう書く 2.0 - ROUND 1 - LL Ring の前哨戦として...
Brainf*ckで100までの素数を列挙してみるテスト【TAKESAKO @ Yet another Cybozu Labs】at 2006年06月26日 16:15
エラトステネスの篩もよいけれど、別の問題もやろうよ。ということで「完全数」です。 def perfect(n = 100, &block) sum = Array.new(n + 1, 1) (2..n).each do |i| if i == sum[i] block.call(i) end k = i + i while k <= n sum[k] += i k += i end end sum end
10000までの完全数を列挙せよ【rubyco(るびこ)の日記】at 2006年06月20日 05:55
素数を列挙するにあたって、エラトステネスのふるいは比較的高速で、たとえば、404 Blog Not Found:LLR2006 - 1,000,000(番目|まで)の素数 で、いくつかの言語で例示されているのとほぼ同じことをする Squeak を使った Smalltalk の次のコード… [ | max primes isPrime | m
[Squeak][Smalltalk][OOPL] エラトステネスのふるいをコンパクトにする方法【sumim’s smalltalking-tos】at 2006年06月19日 17:13
この記事へのコメント
Watermelonさん、
直しました。この場合、primesでも動いちゃうのが問題ですね。join()の定義が別ファイルとかにあればすぐにわかる潜在バグではありますが。
Dan the Enbugger
Posted by at 2006年06月18日 02:58
細かい話で恐縮ですが、Cのコードのjoin関数内…

sprintf(dbuf, "%d ", primes[i]);
->
sprintf(dbuf, "%d ", array[i]);

では。「なんちゃって」でも気になっちゃいました。
Posted by watermelon at 2006年06月18日 01:05