history

青木日記 RSS

<前の日 | この月 | 次の日>

2003-04-25

たらいまわしべんち最終章

しつこくたらいをまわしてみる。

うーむ、原先生ファイナルバージョンがかなり速いですねえ。 第三引数だけ遅延評価するとは凄いアイデアです。

ただ、やっぱりそれでも引数をデカくすると Hugs に負けるんですよー。

~/c/tarai % time ruby yield-tarai.rb                       aamine@harmony
768
ruby yield-tarai.rb  2.20s user 0.01s system 101% cpu 2.187 total
 
~/c/tarai % time runhugs tarai.hs                          aamine@harmony
768
runhugs tarai.hs  2.13s user 0.01s system 101% cpu 2.118 total

これは (768,384,0) で試した場合で、引数を倍にするともっと差が広がります。 笹田さん版も同じですね。 やはりどうやっても完全な遅延評価をしている Haskell 版には敵わないようです。

Ruby / memoise

ちなみに笹田さんの多段ハッシュ版は次のように書いたほうが短くて速いです。

$memo = {}
 
def tarai( x, y, z )
  (($memo[x] ||= {})[y] ||= {})[z] ||= tarai0(x,y,z)
end
 
def tarai0(x, y, z)
  if x <= y
  then y
  else tarai(tarai(x-1,y,z),
             tarai(y-1,z,x),
             tarai(z-1,x,y))
  end
end
 
puts tarai(768, 384, 0)
 
 
~/c/tarai % time ruby mhash-tarai.rb     # 笹田さんのオリジナル
768
ruby mhash-tarai.rb  15.21s user 0.25s system 100% cpu 15.452 total
 
~/c/tarai % time ruby mhash-tarai2.rb    # 青木改良版
768
ruby mhash-tarai2.rb  10.41s user 0.08s system 100% cpu 10.487 total

また負の整数がないと考えてよいのなら 配列を使うとさらに高速化できます。

$memo = []
 
# x だけは -1 になることがあるので +1 しておく。
def tarai( x, y, z )
  (($memo[x+1] ||= [])[y] ||= [])[z] ||= tarai0(x,y,z)
end
 
def tarai0(x, y, z)
  if x <= y
  then y
  else tarai(tarai(x-1,y,z),
             tarai(y-1,z,x),
             tarai(z-1,x,y))
  end
end
 
puts tarai(768, 384, 0)
 
 
~/c/tarai % time ruby mhash-tarai3.rb                      aamine@harmony
768
ruby mhash-tarai3.rb  8.15s user 0.09s system 100% cpu 8.230 total

まだ「第三引数のみyield」には敵わないものの、かなり速くなりました。

C / memoise (1)

C の力任せに計算するやつが Ruby memo 版にすら負けているのが気にくわなかったので、 C でも多次元配列に計算結果をキャッシュするようにしてみました。

#include <stdio.h>
#include <stdlib.h>
 
#define K1 192
#define K2 96
#define K3 0
 
static int tarai(int, int, int);
static inline int tarai0(int, int, int);
static void initialize_cache_table(int);
 
#define max3(a,b,c) max(max(a,b),c)
#define max(x,y) ((x) > (y) ? (x) : (y))
 
int
main(int argc, char **argv)
{
    initialize_cache_table(max3(K1, K2, K3) + 2);
    printf("%d\n", tarai(K1, K2, K3));
    exit(0);
}
 
static int *cache;
static int table_width;
 
static void
initialize_cache_table(int width)
{
    int i, j, k;
 
    table_width = width;
    cache = malloc(width * width * width * sizeof(int));
    if (!cache) {
        fprintf(stderr, "malloc failed\n");
        exit(3);
    }
    for (i = 0; i < width; i++) {
        for (j = 0; j < width; j++) {
            for (k = 0; k < width; k++) {
                cache[(i * width * width) + (j * width) + k] = -2;
            }
        }
    }
}
 
static int
tarai(int x, int y, int z)
{
    int result;
    int key;
    int w = table_width;
 
    key = ((x+1)*w*w) + ((y+1)*w) + (z+1);
    result = cache[key];
    if (result != -2)
        return result;
    result = tarai0(x, y, z);
    cache[key] = result;
    return result;
}
 
static inline int
tarai0(int x, int y, int z)
{
    if (x <= y) return y;
    return tarai(tarai(x-1, y, z),
                 tarai(y-1, z, x),
                 tarai(z-1, x, y));
}

これもかなり速かったんですが、キャッシュの密度が粗すぎるため、 (768,384,0) では malloc() が失敗してしまいました。がくり。

C / memoise (2)

やはり多次元配列では芸がなさすぎました。 イリッフィベクタにしてみます。

#include <stdio.h>
#include <stdlib.h>
 
#define K1 768
#define K2 384
#define K3 0
 
static void *new_cache(void*);
static int tarai(int, int, int);
static inline int tarai0(int, int, int);
 
#define max3(a,b,c) max(max(a,b),c)
#define max(x,y) ((x) > (y) ? (x) : (y))
 
void **cache;
int cache_size;
 
int
main(int argc, char **argv)
{
    cache_size = max3(K1, K2, K3) + 2;
    cache = new_cache(0);
    printf("%d\n", tarai(K1, K2, K3));
    exit(0);
}
 
static void*
new_cache(void *val)
{
    void **p;
    int i;
 
    p = malloc(sizeof(void**) * cache_size);
    if (!p) {
        fprintf(stderr, "malloc failed\n");
        exit(3);
    }
    for (i = 0; i < cache_size; i++) {
        p[i] = val;
    }
    return p;
}
 
static int
tarai(int x, int y, int z)
{
    void **p1, **p2;
    int result;
 
    if (!(p1 = cache[x+1])) {
        cache[x+1] = p1 = new_cache(0);
    }
    if (!(p2 = p1[y])) {
        p1[y] = p2 = new_cache(-2);
    }
    if ((result = p2[z]) == -2) {
        p2[z] = result = tarai0(x,y,z);
    }
    return result;
}
 
static inline int
tarai0(int x, int y, int z)
{
    if (x <= y) return y;
    return tarai(tarai(x-1, y, z),
                 tarai(y-1, z, x),
                 tarai(z-1, x, y));
}

ちっ、こっちのほうが多次元配列を使うより短かかったか。

~/c/tarai % gcc -O6 -finline-functions -omemo-tarai memo-tarai2.c
memo-tarai2.c: In function `tarai':
memo-tarai2.c:54: warning: passing arg 1 of `new_cache' makes pointer from integer without a cast
memo-tarai2.c:56: warning: assignment makes integer from pointer without a cast
memo-tarai2.c:57: warning: assignment makes pointer from integer without a cast

黙れコンパイラ。俺が正義だ!

~/c/tarai % time ./memo-tarai                              aamine@harmony
768
./memo-tarai  0.57s user 1.68s system 101% cpu 2.227 total
 
~/c/tarai % time runhugs tarai.hs                          aamine@harmony
768
runhugs tarai.hs  2.11s user 0.02s system 100% cpu 2.122 total

おしい。もうちょっとなんですけどね。 まだキャッシュの密度が薄いみたいなので、 データ構造を工夫すればまだ速くできそうです。

ただ、こうやって速くしたバージョンが本当に正しい値を計算してるか わからないんですよね。Haskell 版にしても Ruby yield 版にしても、 プログラムが正しそうであることはちょっと見れば確かめられます。 でもここで書いた C 版のプログラムが本当に正しいかどうかは ちょっと見ただけではわかりません。 少なくとも引数に負の数を与えると SEGV するわけだし。

C / memoise (3)

ええい、ここまでやったら行くとこまで逝け! st_table を利用した C / memoise version 3。

#include <stdio.h>
#include <stdlib.h>
#include "st.h"
 
#define K1 768
#define K2 384
#define K3 0
 
static int tarai(int, int, int);
static inline int tarai0(int, int, int);
 
st_table *cache;
 
int
main(int argc, char **argv)
{
    cache = st_init_numtable();
    printf("%d\n", tarai(K1, K2, K3));
    exit(0);
}
 
static int
tarai(int x, int y, int z)
{
    st_table *t1, *t2;
    int result;
 
    if (!st_lookup(cache, x, &t1)) {
        t1 = st_init_numtable();
        st_add_direct(cache, x, t1);
    }
    if (!st_lookup(t1, y, &t2)) {
        t2 = st_init_numtable();
        st_add_direct(t1, y, t2);
    }
    if (!st_lookup(t2, z, &result)) {
        result = tarai0(x,y,z);
        st_add_direct(t2, z, result);
    }
    return result;
}
 
static inline int
tarai0(int x, int y, int z)
{
    if (x <= y) return y;
    return tarai(tarai(x-1, y, z),
                 tarai(y-1, z, x),
                 tarai(z-1, x, y));
}
 
 
~/c/tarai % gcc -O6 -Wall -finline-functions memo-tarai3.c st.c alloc.c -omemo-tarai
memo-tarai3.c: In function `tarai':
memo-tarai3.c:28: warning: passing arg 3 of `st_lookup' from incompatible pointer type
memo-tarai3.c:30: warning: passing arg 3 of `st_add_direct' makes integer from pointer without a cast
memo-tarai3.c:32: warning: passing arg 3 of `st_lookup' from incompatible pointer type
memo-tarai3.c:34: warning: passing arg 3 of `st_add_direct' makes integer from pointer without a cast
memo-tarai3.c:36: warning: passing arg 3 of `st_lookup' from incompatible pointer type
 
~/c/tarai % time ./memo-tarai                              aamine@harmony
768
./memo-tarai  1.34s user 0.07s system 101% cpu 1.387 total

やたー

ちなみに

真の最速はコレみたいです。

~/c/tarai % ghc -O -o tarai tarai.hs                       aamine@harmony
~/c/tarai % time ./tarai                                   aamine@harmony
768
./tarai  0.01s user 0.00s system 579% cpu 0.002 total

to_int

見逃してましたが、4/13 の中田さんのやつが結構速いですね。 基本的には lambda を使うのと同じですけど、to_int を使って むりやりオブジェクト指向ぽくしてるのが特徴ですか。

# ちょっと修正したものを再掲
 
class Tarai
 
  include Comparable
 
  def initialize(x, y, z)
    @x, @y, @z = x, y, z
    @t = 0
  end
 
  def tarai_inspect(level = 3)
    return "?" if (level -= 1) < 0
    s = "Tarai(" <<
        [@x, @y, @z].map {|x| Tarai === x ? x.tarai_inspect(level) : x}.join(", ") <<
        ")"
    s << ("%#d" % @t) unless @t.zero?
    s
  end
 
  def +(val)
    @t += val
  end
 
  def -(val)
    @t -= val
  end
 
  def <=>(val)
    to_int() <=> val.to_int
  end
 
  def coerce(val)
    case val
    when Integer
      [val, to_int()]
    else
      super
    end
  end
 
  def to_int
    @t += if !@x
            0
          elsif @x <= @y
            @y
          else
            Tarai.new(Tarai.new(@x-1, @y, @z),
                      Tarai.new(@y-1, @z, @x),
                      Tarai.new(@z-1, @x, @y)).to_i
          end
    @x = nil
    @t
  end
 
  alias to_i to_int
 
  def to_s
    to_i.to_s
  end
 
end
 
 
t = Tarai.new(*ARGV.map {|n| n.to_i })
puts "#{t.inspect} = #{t}"

たらいまわしパック

全部まとめました。

思う存分たらいまわしてください。

文字列を整数に変換

めも

# Ruby
"15".to_i        # チェックが甘い。e.g. "a".to_i == 0
Integer("15")    # チェックが厳しい。 e.g. Integer("a") は例外
 
; Scheme
(string->number "15")
 
 -- Haskell
read "15" :: Int
 
 -- Haskell (2) HaXml から抜粋
atoi :: String -> Int
atoi = foldl (\x y-> x*10 + ord y - ord '0') 0

Scheme / delay force

遅ればせながら shiro さんとこの Wiki も見てきました。

おおっ、delay/force なんてので遅延評価できるんですね。 Scheme って凄いなー。

まだ続くらしい

うーむ、原先生最新バージョンは……これは…… デフォルト引数を使うことで第三引数の無駄な評価を さらに一段減らしたんですね。す、凄すぎ。 これで (768,384,0) までは Hugs に勝てるようになりました。

本日のツッコミ(全8件) [ツッコミを入れる]
ささだ (2003-04-25 08:16)

失礼ながら,笑ってしまいました.いやぁ,ようやる.つまり,Cで最速を目指すなら,GHC相当を書け,と.

ちなみに,mswin版だとスタック足りないっぽいです>768
インチキとして,tarai(z-1,x,y) は消してしまう,とかありますけど(卑怯).まぁ,そのための遅延評価なわけですが.

しかし,一般的にどういうシチュエーションでこの発想を利用するんだろう?

はら (2003-04-25 12:44)

ピンポイント遅延評価ってのはどうでしょう:

def tarai( x, y, z = nil )
  if x <= y
  then y
  else
    z ||= yield
    tarai(tarai(x-1, y, z), tarai(y-1, z, x)){ tarai(z-1, x, y) }
  end
end

puts tarai(ARGV[0].to_i, ARGV[1].to_i, ARGV[2].to_i)

はら (2003-04-25 12:54)

本来のHakellの話に戻すと(^^;、、、質問です:

(1)Haskell はmemoiseはしないのでしょうか。
(2)Haskell はローカルなmemoise(引数を2度評価しないこと)はしますよね?
(3)gosh の lambda 版が速くないのは、クロージャ生成のコストのせいでしょうか。

ruby-math に行った方がいいかなあ。

あおき (2003-04-25 13:13)

Haskell ML はどうでしょ。
http://www.sampou.org/haskell/
に案内があります。

shiro (2003-04-25 13:57)

(3)に関して:ちゃんと測ったことは無いのですが、delayも環境を保存しているので、生成のコストはクロージャとあまり変わらないと思います。lambda版が速くないのは単純に計算量の問題だと思います。
実装れべるで言えば、delayで作成されたpromiseは最初のforceの時にその内容が評価されて記憶されます。以降のforceは記憶された内容を返すだけです。クロージャの場合は呼び出す度に内容を評価してしまいます。Rubyでも結果をキャッシュするpromiseオブジェクトを作れば同じ様な計算量になると思います。
(キャッシュするというと副作用めいていて関数的でないですが、セマンティクスとしては、そもそも関数的な式は何度評価しても同じなので、その結果をキャッシュするのは処理系による最適化ということになると思います)。

はら (2003-04-25 15:38)

内容を評価するタイミングについてはlambda版もHaskellも同じでは。違いは(僕の言葉でいうところの)ローカルなメモを取るかどうかではないかと想像するのですが。例えば

(define (tarai x y z)
  (let ((x0 (x)) (y0 (y)))
    (if (<= x0 y0)
        y0
        (tarai (lambda () (tarai (lambda () (- x0 1)) y z))
         (lambda () (tarai (lambda () (- y0 1)) z x))
         (lambda () (tarai (lambda () (- (z) 1)) x y))))))

とすると劇的に速くなります。

shiro (2003-04-25 17:21)

確かに、ローカルで一回評価しておくと、計算量は遅延評価と同じ、線形になりますね。こちらで検証してみました:
http://www.shiro.dreamhost.com/scheme/wiliki/wiliki.cgi?Scheme%3a%a4%bf%a4%e9%a4%a4%a4%de%a4%ef%a4%b7%a4%d9%a4%f3%a4%c1

nobsun (2003-04-25 19:59)

うぇーっ。すごいことになっていますねぇ。
Haskellの簡約は、outermost graph reduction といって、同じ変数に束縛
した式は高々一度しか簡約しないというものです。要するに、変数は式を指す
ポインタでみたいなものになっているわけです。
Introduction to Functional Programming using Haskell 2ed. §7.1

名前
メールアドレス

<前の日 | この月 | 次の日>
2002|04|05|06|07|08|09|10|11|12|
2003|01|02|03|04|05|06|07|08|09|10|11|12|
2004|01|02|03|04|05|06|07|08|09|10|11|12|
2005|01|02|03|04|05|06|07|08|09|10|11|12|
2006|01|02|03|04|05|06|07|08|09|10|11|12|
2007|01|02|03|04|05|06|07|08|09|10|11|12|
2008|01|02|04|05|06|09