憂鬱な午後のひととき

10月31日_
・SSEの効果が活かせるよう他の箇所も見てみたら, いじりたくなる部分が沢山ありすぎたので止めた.
・medianでおもしろかったのは最初手を抜いてminpsを使ってはまったこと. これはfloatとして比較してminを求めるものだが, 正のintで一定範囲内ならintにも使える(3D Now!のときは常套手段). ただintrinsic関数で使うには引数の型が__m128でないといけない. 今回は__m128iなので型が違って使えない. つまり 
  __m128i x, y, z;
  z = _mm_min_ps(x, y);
とはできず, 
  __m128i x, y, z;
  __m128 f;
  f = _mm_min_ps(*(__m128i)&x, *(__m128i)&y);
  z = *(__m128i*)&f;
としていた. なんてださい…. まあ, こう書いてもVCだとメモリアクセスせずにminps 1命令に変換されるので妥協していた. ただgccだとメモリアクセスが発生するのでどうしたもんかとヘッダを見ていたら, コンパイラに型変換を指示する(実際のコードは生成されない)便利な関数の存在を知った. つまりこう書ける. 
  x = _mm_castps_si128(
          _mm_min_ps(_mm_castsi128_ps(x), _mm_castsi128_ps(y)
      );
え, 常識? でも_mm_castps_si128や_mm_castsi128_psを googleで見ても ほとんどヒットしないからあまり知られてないんちゃうやろか.
とにかく長年のうっとうしい疑問が解けてよかった --- のはよいがここに大きな落とし穴があった.
median算出の中でminpsを使うとCore2duoにおいて強烈な(10倍以上の)速度低下を招くのだ. 確かにfloat系と整数系の命令をまぜるなというのはマニュアルにあるが実際ここまで顕著な差を見たのは初めて. intのビットパターンなのでfloatとしてはイリーガルなせいもあるかもしれないが少なくともPentiumDでは遅くならない(むしろ下記のコードより速い). 今まで書いてた自分の他のコードも見直してみよう.
ちなみにintのminはSSE4.1にしかないので(最初から入れとけ) 
  // mask = (x > y) ? -1 : 0
  // x = (y & mask) | (x & ~mask)
  __m128i mask = _mm_cmpgt_epi32(x, y);
  x = _mm_or_si128(_mm_and_si128(y, mask),
                   _mm_andnot_si128(mask, x));
とした.
10月30日_
・うーむ, まだSSEに落としただけで細かい最適化はしていないとは言え, 期待した速度(2倍程度)は出なかった(USE_SSEあり/なし).
PenDなら多少速く, Core2なら(特にgccだとVCに比べて分岐の最適化がよくてSSEの最適化が多少劣るので)逆に遅くなるかもしれない. なお, このコードはmedianのみに注目してるのでshinhさんとこのまともな最適化は何もしていない. 局所的なmedianのベンチでは1.5～2.5倍速いので何か別の部分が悪さをしてるのかな. ざっと見てるとSSEレジスタが無駄にメモリに保存されている? コスト比で考えるなら, SSEは使わずまっとうに速くするのがよさげ.
・ただ求めたいのは実は median じゃなくて
そうですね. で, 当初頭にあったのは 
  if (play[0] == median) out[0] += s;
  if (play[1] == median) out[1] += s;
  if (play[2] == median) out[2] += s;
をSIMDのままやってしまうというもの. こんな感じ. 
  mask = _mm_cmpeq_epi32(x, play.x_);
                 // [x == median ? -1 : 0 for x in play[0..2]]
  x = to_score(x);

  x = _mm_and_si128(x, mask);  // [x == median ? x : 0 for x]
  x = _mm_add_epi32(x_, x);    // [out + x for x]

  int result = 2;
  if (v[0] == median) {
      result = 4;
  }
  if (v[1] == median) {
      result -= 1;
  }
  if (v[2] == median) {
      result -= 1;
  }
のほうもv[0]なら2, v[1], v[2]なら-1というoffsetを加算してからあとでまとめることで同様に処理可能. ただし最後の水平加算が速くないのでとんとんかも.
・64bit VCだとcmov命令を生成していた. 結果, 32bit Xpよりも64bit Xp(VMware) on 32bit Xpの方が3割ほど速いというおもしろい結果に(レジスタ増えたせいもあるだろう).
10月29日_
次の表を見れば分かった.
状態min(a, b)min(b, c)min(c, a)median
a <= b <= cabab
a <= c <= bacac
b <= a <= cbbaa
b <= c <= abbcc
c <= a <= bacca
c <= b <= abccb
3個のminの値はmedianを除いたものが同じ. つまり 
  median(a, b, c) = min(a, b) ^ min(b, c) ^ min(c, a)
ただこれに変更しても遅くなった. VCだとcmov命令を使わないので論外(前から思ってたのだが何故VCはこの種の最適化をおこなわないんだろう. /arch:SSE2をつければcmovを使うとマニュアルにはあるのに. 何か条件が要るのか?) ここだけinline asmで書いても関数呼び出しコストが大きくて遅くなる. 
__declspec(naked) int median(int x0, int x1, int x2)
{
    __asm {
        mov       eax, dword ptr [esp + 4] // x0
        mov       edx, dword ptr [esp + 8] // x1
        cmp       eax, edx
        cmovnc    eax, edx // min(x0, x1)
        mov       ecx, dword ptr [esp + 12] // x2
        cmp       edx, ecx
        cmovnc    edx, ecx // min(x1, x2)
        xor       eax, edx
        mov       edx, dword ptr [esp + 4]
        cmp       ecx, edx
        cmovnc    ecx, edx
        xor       eax, ecx
        ret
    }
}

__declspec(naked) int medianSSE(const int *x)
{
    __asm {
        mov     ecx, dword ptr [esp + 4]
        movd    xmm0, dword ptr [ecx + 0]
        movd    xmm1, dword ptr [ecx + 4]
        movd    xmm2, dword ptr [ecx + 8]
        movaps  xmm3, xmm0
        minss   xmm0, xmm1 // min(x0, x1)
        minss   xmm1, xmm2 // min(x1, x2)
        minss   xmm2, xmm3 // min(x2, x0)
        pxor    xmm0, xmm1
        pxor    xmm0, xmm2
        movd    eax, xmm0
        ret
    }
}
gccではcmov命令を生成するのにも関わらず遅くなった. 分岐予測が当たってるんだろう. minpsを使えばなんとかなるかもしれないけどそのためには少なくとも付近もろともSSEを使うように変更しないと水平操作のコストが大きすぎる.
とある方からこのmedianのアルゴリズムはJack BenkualさんやIntel IPPチームが見つけたものというメールを頂いた. 
 median(a, b, c) = a + b + c - min(a, b, c) - max(a, b, c)
という別解も教えていただいた. ありがとうございます.
10月28日_
ロジックの修正やまっとうな高速化はshinhさんに任せてmedianを考えてみる.
scoresをSSEレジスタ一つで表現したいから, 問題は「a = x[0], b = x[1], c = x[2]が与えられたとき分岐無しでmedianを求めよ. ただしmin, max, (x > y) ? -1 : 0を使ってよい.」
α = a > bβ = b > cγ = c > a状態medianidx
000a = b = cany0
001a <= b <= c, a < cb1
010c <= a <= b, c < ba0
011a < c < bc2
100b <= c <= a, b < ac2
101b < a < ca0
110c < b < ab1
111not herex-
上記7パターンで分類できるから, α:β:γを3bitとしてtbl[] = { 0, 1, 0, 2, 2, 0, 1 };に対するテーブル引きでmedianのindexはとれる(3個の同時比較はSSE1命令). でも, 3bitにまとめて表引きするコストが大きいのでこのままでは遅い. もっとシンプルに, あるいは直接medianの値を求める方法を考える必要がある. [α^β:α^γ]として2bitにまではできるけどそれ以上まとめられない….
10月27日_
単に
scoresをboost::array<int, 3>に変えると10倍以上速くなりました. あとMWinGameのコンストラクタは参照渡しにしないとまずいような気が. ところでint x 3が与えられたときの真ん中の値を取るmedianの最適化ってどうするんだろう. 考えてみる.
10月23日_
家のLinuxサーバのkernelを2.4.35から2.6.27.3へ移行(大変だった…). gcc 4.xで2.4系をmakeするときはその場しのぎパッチが必要だったけど, これで何もしなくていい. どう考えてもクリーンインストールの方が楽だがPlamo-3.0の頃(7年前?)からいじってるので愛着があるというか依怙地になってる.
10月22日_
Pythonのraw文字列って\を最後におけないんだな. 変な仕様だ. r'\\'だと\二つ分になってしまうし.
10月21日_
『九月の恋と出会うまで』(松尾由美)
導入部の不動産屋さんの粋な計らい, アパートに住む怪しげなキャラクタたち, 意外な展開ときてこれは凄い期待大, とワクワクしたのに, そこからありきたりな話の単調現象で収束. 実にもったいない.
10月20日_
寝ても覚めても最適化な人たちと会う. えらいマニアックで(よく分からん)会話してるのを見てるのも面白かった.
10月17日_
That’s even faster than my C++ version
論点はそこじゃないので野暮なつっこみだが, C++のコードがちょっと酷いなあ. boost::shared_ptrは全く不要だ(なぜわざわざstd::listをnewする. 排他制御もついてまわるぞ). それ消すだけでblogのC++の2倍になった. ついでに(これはアンフェアやけども)template化すれば30倍になった. 
template<int N>
struct array {
    array() {}
    array(const int *x)
    {
        for (int i = 0; i < N; i++) x_[i] = x[i];
    }
    int x_[N];
};

template<int N, int aN, int bN>
array<N> merge(const array<aN>& a, const array<bN>& b)
{
    array<aN + bN> c;
    int ia = 0, ib = 0, ic = 0;
    while (ia != aN && ib != bN) {
        if (a.x_[ia] > b.x_[ib]) {
            c.x_[ic++] = b.x_[ib++];
        } else {
            c.x_[ic++] = a.x_[ia++];
        }
    }
    memcpy(&c.x_[ic], &a.x_[ia], (aN - ia) * sizeof(a.x_[ia]));
    memcpy(&c.x_[ic + aN - ia], &b.x_[ib],
           (bN - ib) * sizeof(b.x_[ib]));
    return c;
}

template<int N>
array<N> mergesort(const array<N>& a)
{
    if (N <= 1) return a;
    const int m = N / 2;
    return merge<N>(mergesort(array<m>(&a.x_[0])),
                    mergesort(array<N - m>(&a.x_[m])));
}
10月15日_
・暴走族を珍走団とお呼び
ページがごちゃごちゃしてて中身を読んでないのだけどタイトルを思いついただけで凄いと思った.
・下記のコードはgcc 4.3.0でコンパイルエラーになる. コピーコンストラクタをpublicにすればよいんだけど呼ばれてないのにエラーになるのはちょっとやだなあ. 
#include <stdio.h>
struct A {
    A(const char *p) : p_(p) {}
private:
    A(const A&) { puts("not call"); }
    const char *p_;
};

static const A tbl[] = {
    "abc",
    "def"
};

int main()
{
}
10月9日_
『七瀬ふたたび』(筒井康隆)
がドラマでやってる…. ちらっと見た限りでは原作と雰囲気が違うそうな感じだったけど, とりあえず撮っておく. 見られる時間取れるかな.
10月7日_
・ノーベル物理学賞、素粒子研究の日本人３氏に(via Yoichi's diary)
おお, ついに. 隅っこで講義を受けてた頃既に, 何時貰えるの?みたいな話があったけど素晴らしい.
授賞式が初の海外、筋金入り英語嫌い
当時から「自分の論文を日本語で世界の人に読んでもらいたい」というのが夢。そのせいもあり、英語嫌いは筋金入りで、これまでパスポートを取ったことがない。学会など海外からの招待はすべて断り、小林さんに任せてきた。
凄い.
・↓2倍ぐらいしか速くならないのもあるし, 当然のことながら遅くなるのもある. 簡単のに限ったとしても難しいなあ.
10月6日_
久しぶりに時間をとってへなちょこ正規表現ライブラリいじった. 結局単純なのしかサポートせず, しかもVMをすっとばしてダイレクトにx86を生成することにした(弱気). その代わり多少速くなった. 作為的なパターン 
const char *text = "a<div>jevifjaesdasdfasdf5asdsafasdfas\
                    dfivafasdfkaaaa3aaajasdfasdvslevrjase\
                    ivjfasiajevifjaesdasdfasdfas-d>safasd\
                    fasdfivafasdfk<aaaaaaabjasdfasdvslevr\
                    jaseivjfasijfajafveiabdevifjaesd<asdf\
                    a$sdfasd>safasdfasdfivafasdfkjasdfasd\
                    vslevrjaseivjfasijfajafveiabda<a5aab>\
                    sdfasdfXaw<html>ei</div>fvasvfas";
const char *reg = "<div>.*?</div>";
だとboost vs my regex = 3.250 : 0.270@Core Duo 1.8GHz.
10月5日_
久しぶりに川の側を歩いてみた. が, コンクリートで埋められてて嬉しくない. 京都の鴨川や比叡山を散歩したいなあ. 大台ヶ原へサイクリング再挑戦もしたい.
10月2日_
電車の終点で車掌さんに起こされてしまった. 眠い.
10月1日_
・うーーーん
最初に「ぱっと目に付く」のはnew char[]でサイズを指定していないことじゃないかしら. あとoperator=()は自己代入に対する考慮(配慮?)じゃなくて例外安全なswapを実装してから作るのがよいかと. とか書きながら昔自分が書いたコードを思い出して叫びたくなるやつ.
・「こんにゃく入りゼリー」よりものどに詰まって死亡した件数が多い危険な食べ物ベスト10
適度に皮肉も交えて面白い記事.
※こんにゃく入りゼリーの方が、ご飯などと比較して口にする割合が少ないので、「ご飯はこんにゃく入りゼリーよりもはるかに高い危険が食事のたびに付きまとっている」わけではないです。
当たり前なんだけどこういう注意を書かない(書けない?)記事がなんと多いことか. 件のすぐ辞めた大臣の発言を見ても, いかに(統計)教育が大事かがよく分かる.
昔
[Back]
御意見は光成滋生<herumi@nifty.com>までお願いします