DSAS開発者の部屋

まくおさんをバージョンアップしましたのでお知らせします。
主な変更点は以下のとおりです。

【Ver1.1.3→Ver1.2.0】

• 応答速度と安定性の向上
• エラーログに"[error]"という文字列を追加
• メンバがタイムアウトしたときの状況を詳細に出力するようにした
• msync --statusで、実行中のコマンドを表示するようにした
• msyncを複数同時実行したときに挙動がおかしくなる不具合を修正
• ファイルクローズに失敗したときに異常終了する不具合を修正
• msyncにdeleteオプションを指定したときの結果表示がおかしくなる不具合を修正
• その他細かい不具合を修正

【ダウンロード】

• ダウンロードページ
• makuosan-1.2.0.tar.gz

「ところで、まくおさんってだれ？」という方は、こちらやこちらを見てやってください(^^;

前回 インタプリタ型言語を高速化する computed goto で紹介したcomputed gotoを、「phpでも使えないの？」という声が社内であったので、 php のソースコードを見たところ、特定のビルド手順で php でも computed goto が使えることが判りました。そのビルド手法とベンチマーク結果を紹介しておきます。

php の VM のソースコードは、phpソースパッケージ中の Zend/ というディレクトリの中にあります。zend_vm_で始まる幾つかのファイルのうち、 zend_vm_execute.h というファイルが命令ディスパッチが実装されているファイルで、このファイルは zend_vm_gen.php というスクリプトで生成されています。そして、 zend_vm_gen.php のオプションで、命令ディスパッチの方法を選択できます。(phpのビルドにphpが必要です)

computed goto を使うビルド手順は以下の通りになります。

$ tar xzf php-5.2.x
$ cd xzf php-5.2.x/Zend
$ php zend_vm_gen.php --with-vm-kind=GOTO
$ cd ..
$ ./configure --with-zend-vm=GOTO [その他オプション]
$ make

Zend/ ディレクトリの中に bench.php というファイルがあったので、これを使って通常の関数ポインタ配列版とcomputed goto版の速度を比較してみました。

# 関数ポインタ配列版.
$ ~/local/stow/php-5.2.9-func/bin/php bench.php
simple             0.319
simplecall         0.412
simpleucall        0.639
simpleudcall       0.693
mandel             0.991
mandel2            1.311
ackermann(7)       0.767
ary(50000)         0.047
ary2(50000)        0.037
ary3(2000)         0.534
fibo(30)           1.739
hash1(50000)       0.087
hash2(500)         0.098
heapsort(20000)    0.312
matrix(20)         0.257
nestedloop(12)     0.539
sieve(30)          0.232
strcat(200000)     0.034
------------------------
Total              9.047

# computed goto版
$ ~/local/stow/php-5.2.9-goto/bin/php bench.php
simple             0.283
simplecall         0.360
simpleucall        0.584
simpleudcall       0.645
mandel             1.022
mandel2            1.213
ackermann(7)       0.788
ary(50000)         0.046
ary2(50000)        0.036
ary3(2000)         0.485
fibo(30)           1.716
hash1(50000)       0.079
hash2(500)         0.090
heapsort(20000)    0.276
matrix(20)         0.237
nestedloop(12)     0.479
sieve(30)          0.234
strcat(200000)     0.035
------------------------
Total              8.610

5%くらい速くなっていますね。 Python 3.1 だけじゃなくて php をビルドするときも computed goto を試してみて下さい。

先日Python 3.1a1 がリリースされました。 Python 3.0 は Python 2.6 に比べてパフォーマンスが悪かったのですが、Python3.1はPython2.6よりも速くなるかもしれません。

Python3.1のパフォーマンス向上は、主に次の2点が影響しています。

• ioモジュールがC言語で書き直された
• computed goto の採用 (--with-computed-gotos というconfigureオプションで有効)

computed goto という名前を聞き慣れなかったのですが、調べてみると Ruby 1.9 の VM (YARV) や、 Perl6 の VM として開発されとうとうリリースされた Parrot にも採用されている手法でした。今回は簡単に computed goto の紹介をしてみます。

    goto label;
    puts("foo"); // 表示されない.
label: //ここにジャンプする.
    puts("bar");

    void *jump_target = &&label;
    goto *jump_target;
    puts("foo");
label:
    puts("bar");

for (;;) {
    code = *code_pc++;
    switch (code) {
    case CODE_ADD:
        // 足し算処理.
        break;
    case CODE_SUB:
        // 引き算処理.
        break;
    // ...
    case CODE_END:
        goto LOOP_END;
    }
}
LOOP_END;

    static const void *JUMPTABLE[] = {&&CODE_ADD, &&CODE_SUB, ... , &&CODE_END};

    code = *code_pc++;
    goto *JUMPTABLE[code];

CODE_ADD:
    //足し算処理.
    code = *code_pc++;
    goto *JUMPTABLE[code];

CODE_SUB:
    // 引き算処理
    code = *code_pc++;
    goto *JUMPTABLE[code];
    ...
CODE_END:

bitslice とは

Hack the Cell '09 に参加して知った、Cellに限らず一般的に使えるビット演算の高速化手法について紹介します。

Bitslice と呼ばれる手法では、ビット順を90度回転します。言葉で説明するよりもコードを見たほうが早いので、回転させるコードの例を見てください

int   x[32], y[32]; // x が元のデータ、y が回転後のデータ.
for (int i = 0; i < 32; ++i) {
    int t = 0;
    for (int j = 0; j < 32; ++j)
        t |= ((x[j] >> i) & 1) << j; // x[j] の i ビット目を
    y[i] = t;                        // y[i] の j ビット目にする
}

この変換をすることで、y[0] には x[0] - x[31] の最下位ビットが、 y[1] には 2番目のビットが格納されることになります。

回転することでビット演算が高速になる様を見てみます。

//--
// and 演算

// 通常
for (int i = 0; i < 32; ++i) {
    x[i] &= 3; // load, and, store, 各32回
}

// bitslice
for (int i = 2; i < 32; ++i) {
    y[i] = 0; // store 30回
}

//--
// xor 演算

// 通常
for (int i = 0; i < 32; ++i) {
    x[i] ^= 0xff; // load, xor, store, 各32回
}

// bitslice
for (int i = 0; i < 8; ++i) {
    y[i] = ~y[i]; // load, not, store 各8回
}

//--
// シフト演算
for (int i = 0; i < 32; ++i) {
    x[i] >>= 16; // load, 右シフト, store, 各32回
}

// bitslice
for (int i = 0; i < 16; ++i) {
    y[i] = y[i+16]; // load, store, 各16回
}

bitslice はなぜ速いか

無駄なビットに対する演算を省略できる

たとえば、 x ^= 1 (xは32bit整数) という演算は、最下位ビットを反転させているだけで、最下位以外の31bitにはまったく影響がありません。せっかくの32bit演算なのに、実質は1bit演算になってしまっているわけです。

bitslice なら最下位ビットだけを集めた32bitに対して演算するので、32bit演算の能力を無駄にすることがありません。

マスクがいらない

マスクは操作したくないビットを保護するためのものです。操作したいビットが集まっていれば、マスクが不要です。

マスク付きの通常ビット演算とbitsliceでの演算の対応は次のようになります。

• and演算は、マスクビットが0の場所への 0 代入に
• or演算は、マスクビットが1の場所への 1 代入に
• xor演算は、マスクビットが1の場所への not 演算に

load が減る

普通の and や or 演算においてロードが必要なのは、操作しない(マスクされた)ビットを元の値にしておくためです。bitsliceでは全bitに1や0を代入できるので、不要な load を減らすことができます。

load が減るということは、単純に命令数が削減する以上の意味があります。一般的に言って、store命令は遅延して次の命令と並列に実行できるのに対して、load命令は完了するまでloadした値に依存する命令を実行できません。load命令が減ることで、パイプラインストールを減らすことができます。

工夫次第でさらに速く

先ほどの例では、通常のシフト演算を配列内での移動に置き換えていました。しかし、 y[i] = y[i+16] の代わりに y += 16 が許されるのであれば、シフト演算は完全に省略できることになります。

同じように or 演算でも、 y[i] = 0 をする代わりに、次に y[i] を参照する場所を 0 に置き換えることができるかもしれません。状況によって、いろいろな工夫が考えられます。

bitslice が苦手な計算

bitslice では各ビットがバラバラになっているので、ビット単位の演算は得意な代わりに、ビットをまたがる演算が苦手です。たとえば加算や減算は繰り上がりや繰り下がりがビットをまたぐので、bitsliceでは遅くなります。試しに加算を実装してみます

// 通常の演算
for (int i = 0; i < 32; ++i) {
    x[i] += z;
}

// bitslice
int carry = 0;
for (int i = 0; i < 32; ++i) {
    int c;
    c = (carry & y[i]) | (y[i] & z) | (z & carry);
    y[i] = carry ^ y[i] ^ z;
    carry = c;
}

ただし、特定の条件では通常の演算に近い速度が出せることもあります。次の例では、CPUに bitcnt() という bit が 1 になっている数を数える命令があると仮定して、配列の合計値を計算しています

// 通常の演算
int sum = x[0]; // load 1回
for (int i = 1; i < 32; ++i) {
    sum += x[i];  // load, add, 各31回
}

// bitslice
int sum = bitcnt(y[0]); // load, bitcnt, 各1回
for (int i = 1; i < 32; ++i) {
    sum += bitcnt(y[i]) << i; // load, bitcnt, shift, add 各31回
}

Hack the Cell '09 に参戦していました。Fixstars社からの結果発表を待ってから成績報告をしようと思っていたのですが、他の参加者の皆さんがどんどんとスコアや素晴らしいテクニックを披露されているので、予定を前倒しして私のスコアとソースコードを公開します。（中身に関してはまた別の記事に書きます）

提出物とスコア

ソースコード (試行錯誤の跡が整理されていません)

スコア

ORIGNAL:         sum=3c927c56, 294030647 ticks
MINE:            sum=3c927c56, 4464381 ticks
ORIGNAL:         sum=2e987a4d, 424155603 ticks
MINE:            sum=2e987a4d, 6440068 ticks
ORIGNAL:         sum=ef1b6aef, 312102737 ticks
MINE:            sum=ef1b6aef, 4738888 ticks
ORIGNAL:         sum=eedd2516, 290055058 ticks
MINE:            sum=eedd2516, 4403913 ticks
ORIGNAL:         sum=f7e967a8, 14366822 ticks
MINE:            sum=f7e967a8, 218363 ticks
ORIGNAL:         sum=1f37a7db, 214216185 ticks
MINE:            sum=1f37a7db, 3252580 ticks
ORIGNAL:         sum=c7d41f36, 294964206 ticks
MINE:            sum=c7d41f36, 4478494 ticks
ORIGNAL:         sum=aa9d2e9f, 259565045 ticks
MINE:            sum=aa9d2e9f, 3941212 ticks
ORIGNAL:         sum=8abd398a, 250844221 ticks
MINE:            sum=8abd398a, 3808729 ticks
ORIGNAL:         sum=a374bd58, 6110290 ticks
MINE:            sum=a374bd58, 92978 ticks

倍率としては、約66倍といったところです。ただし、ORIGINALの速度が安定しない(プログラムのアライメントによって実行速度が変化する)ので、倍率は参考程度にとどめておいてください。

今回のコンテストの内容について

課題は、メルセンヌツイスターという乱数生成機の高速化でした。

最適化する対象となる関数はこちらです(一部省略)

#define N 624
#define M 397
#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
#define LOWER_MASK 0x7fffffffUL /* least significant r bits */

static unsigned long mt[N]; /* the array for the state vector  */
static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */

/* ... */

unsigned long genrand_int32(void)
{
    unsigned long y;
    static unsigned long mag01[2]={0x0UL, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */

    if (mti >= N) { /* generate N words at one time */
        int kk;

        for (kk=0;kk<N-M;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        for (;kk<N-1;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
        mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];

        mti = 0;
    }

    y = mt[mti++];

    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return y;
}

しかし、最適化したコードを実装する側のスケルトンとなる関数は次のような定義になっていました

unsigned int
genrand_mine(int num_rand)
{
  int r = 0;
  /* r = num_rand 個の乱数のチェックサム */
  return r;
}

この関数は戻り値として、 num_rand 個の生成された乱数のチェックサムを返す必要があります。

もともとの課題では、「※疑似乱数列は、メルセンヌ・ツイスタの実装 mt19937ar.c と同じ乱数列を生成してください。」となっていたのですが、参加者によって以下のように異なった解釈がされていました。

• 同じ32bit整数としての乱数を求められている (私はこの解釈をしていました)
• 乱数のビット順がどうなっていても良い
• そもそも配布されているプログラムがチェックサムしか見ていないのだから、チェックサムだけ合ってれば何をしてもOK

この解釈の違いが問題になっていたのですが、最終的には現在課題のページにあるように、一番ゆるい「チェックサムが合っていれば何をしてもOK」というルールで優勝・準優勝を決め、きちんと乱数を生成していた人にはFixstars社の判断で「フィックススターズ賞」が送られることになりました。

結果として、今回のコンテストにおいて最適化の方針として3つが生まれました。

1. 元のコードをそのままSIMD化する

   基本は元のプログラムのまま、mt[] の要素を4つひとまとめに処理する。 その後、命令数を削減したりパイプラインを埋めて速度を稼いでいく。 この方法だと、60倍〜70倍程度までを目指せます。

   元の乱数生成器と同じ32bit整数で乱数を生成するので、フィックススターズ賞ねらいになります。

2. ビット順を90度入れ替えた bitsliece と呼ばれる構造で乱数を計算する

   SPUの128bitレジスタを、32bit整数4つとしてではなくて、1bitが128個あると考えて最適化していきます。命令の組み合わせだけでなくデータ構造をどうするのかにも工夫が必要になります。この方法で100倍を超えた参加者もおられるようです。

   この方法では乱数が32bit整数の形では出現しないので、優勝・準優勝ねらいになります。

3. 乱数を生成しないで直接チェックサムを計算する

   mt[] の値は 32bit整数を乱数の種として初期化されるので、種と num_rand と組み合わせた 64bitの情報から結果のチェックサム32bitを求めるということが究極的な目標になります。 極端な例を言えば 32bit * 2^64 分のテーブルを用意すれば、O(1) でチェックサムの計算ができる事になります。

   私には全く思いつかなかったのですが、 この計算をメモリ制約の中で O(n) 未満で実行する方法が存在するかもしれません。

私は2月は忙しくて全く高速化できなかったので、ルール確定前までに書いていた 1 の方法のプログラムのまま提出しました。 提出締め切り後、同じ1の方法でもあと5%以上高速化する方法が他の参加者のブログで紹介されていて悔しい思いをしています。 次回があればまた参加したいと思います。

本日発売される日経Linux 4月号からKLab の社員で執筆した「2000円ルーターを改造しよう」という連載が始まります。

この連載では 2000円〜4000円で購入できる無線LANルーター La Fonera に DD-Wrt や OpenWrt などの代替ファームウェアに置き換えて活用する方法を紹介していきます。

4月号の第1回では、前準備としてFONルーターのブートローダーにシリアル接続する方法や、telnet でログインする方法、などを紹介しています。

これから安価で自由度の高い無線LANルーターを作ってみたいという方はぜひご覧にな って下さい。

すっかり花粉が飛び交う季節になってしまいました。 今や花粉症は5人に1人疾患する症状とも言われており、マスクや薬などによる花粉対策の準備を始めている方が多いのではないかと思います。

そこで、今回はプログラマに出来る花粉対策を紹介したいと思います。 花粉症対策にはとにかく花粉を吸わないことが重要です。 大量の花粉にばく露されることで、症状が悪化したり、今まで花粉症で無かった人も花粉症が発症する場合があります。 つまり花粉が多く飛んでいる時間帯をなるべく避けることが有効な花粉対策なのですが、 天気予報の花粉飛散情報は1日単位のデータだったりして、参考にならないことが多いです。

環境省が運用する花粉観測システム(はなこさん)では１時間毎の花粉飛散データが公開されています。 1時間毎に花粉飛散データが更新されるため外出のタイミングを決定するのに役立ちます。 プログラマであれば、このデータを cron で定期的に監視して警告メールを送ったり、Shell のプロンプトやemacs のミニバッファに表示させたい、と思うはずです。

そこで、この花粉観測システムはなこさんから花粉飛散データを取得する perl module を作成しました。

WWW::Hanako - Perl interface for Hanako(Pollen observation system at Japan) - search.cpan.org

この perlモジュールを使用する前に、まずお近くの観測所コードを調べる必要があります。 東京の観測所コードは以下の通りです。

例えば、日本医科大学付属病院での観測情報を得るコードは以下の様になります。

use WWW::Hanako;
my $hanako = WWW::Hanako->new(area=>3, mst=>51300200);
print $hanako->now()->{pollen} . "\n";

表示される、飛散量の数値は個／m^3中に観測された花粉数の１時間平均、となります。 この数値が1000を越えた場合、かなり危険なので十分注意して外出しましょう。

また、この花粉観測システム(はなこさん)のデータは転用や営利目的とした利用を禁止していますので、あくまで個人利用に限るよう注意してください。

参考サイト: 花粉症環境保健マニュアル−２００９年２月改定版−

WSGIとは

PythonでWebアプリを作るときに必ず出てくる単語にWSGIがあります。 WSGIとは、Web Server Gateway Interface の略で、WebサーバーとPython製Webアプリを つなげる標準インタフェースです。

WSGIの上で動くようにアプリケーションを作ると、そのアプリケーションは修正無しに Apache+mod_wsgi, Apache+mod_python, fastcgi, scgi, cgi, 等の環境で動かせるように なります。

他にもミドルウェアという考え方があります。例えばOpenID認証機能をWebフレームワークの プラグインとして開発した場合では他のWebフレームワークでは利用できないのですが、 WSGIミドルウェアとして開発すればWebフレームワークを問わずに利用できるようになります。

標準ライブラリのwsgirefモジュール

Python2.5から、標準ライブラリにwsgirefモジュールが追加されました。その中にはWSGIサーバーも 入っているので、自分で書いたWSGIアプリを手軽に立ち上げたいときに重宝します。

次のコードは、簡単なWSGIサーバー+アプリの例です。

from wsgiref.simple_server import *

def myapp(environ, start_response):
    start_response('200 OK', [('Content-Type', 'text/plain')])
    return ['hello']

server = make_server('0.0.0.0', 8000, myapp)
server.serve_forever()

WSGIServerのマルチスレッド化

上記のコードで立ち上がるWebサーバーは、シングルスレッド&シングルプロセスで動作しているので、 複数のクライアントからアクセスされた場合に問題があります。接続が遅かったり不安定だったりする クライアントが一つでも合った場合、その接続が完了するか切断するまで他の接続ができないので、 接続に失敗したり数十秒かかったりすることになります。

普通はWebアプリを公開する場合には標準ライブラリのWSGIServerではなくて別のサーバーを使う ものですが、お手軽な方法として3行追加と1行修正のみで、他のライブラリ無しに マルチスレッド化したりマルチプロセス化することが出来ます。 次のコードが上記のコードをマルチスレッド化したものになります。

from wsgiref.simple_server import *
from SocketServer import *
class ThreadingWsgiServer(ThreadingMixIn, WSGIServer):
    pass

def myapp(environ, start_response):
    start_response('200 OK', [('Content-Type', 'text/plain')])
    return ['hello']

server = make_server('0.0.0.0', 8000, myapp, ThreadingWsgiServer)
server.serve_forever()

このコードについて簡単に説明します。 wsgiref.simple_server.make_server は第三引数でサーバークラスを指定することができて、 デフォルト値は wsgiref.simple_server.WSGIServer になっています。 WSGIServerの継承関係は次のようになっています。

SocketServer.BaseServer
  ↑
SocketServer.TCPServer
  ↑
BaseHTTPServer.HTTPServer
  ↑
wsgiref.simple_server.WSGIServer

SocketServerモジュールには、ThreadingMixIn, ForkingMixIn というクラスがあり、 SocketServer.BaseServer を継承したサーバークラスであれば簡単にマルチスレッド化 したりfork化したりできるようになっています。 WSGIServer も BaseServer を継承しているので、ThreadingMixIn/ForkingMixIn クラスを Mix-in してやるだけでマルチスレッド化やマルチプロセス化が可能です。

小さなテクニックですが、Webアプリが完成するまではサーバーの設定の問題とアプリの 問題が混ざらないようにしたいときとか、Webサーバーの設定ファイルを書くのに慣れてなくて 後回しにしたいときに重宝します。

ちなみに、 flup というライブラリには、 thread pool や prefork を使ったWSGIサーバーがあって十分実用になる速度が出るので、 標準ライブラリのサーバーだと遅いという時に利用できます。