憂鬱な午後のひととき

『歯車のおはなし』(中里為成)

図書館でなんとなく目についたので読んでみた. 円と直線, 円と円と接して動くといえばサイクロイドを思い浮かべる. 実際, 機械式時計では使われている(いた?).

ところが数学的にはよくてもちょっとの誤差でかみ合いにくくなるそうな. というわけで普通はサイクロイドは使われず, インボリュート曲線が利用されるらしい.

現実的な誤差を許容するとか作りやすいとかの条件により求める解が異なるというのは, 数学と工学の違いを表していておもしろい.

やや残念なのは文章が堅苦しいというか妙に古い感じ. 卓上計算機って, これ1996年発行の本なんやけど. あと文章と図の対応が分かりにくい. たとえば全歯たけやモジュールという重要そうな用語が言葉だけで説明していて図に書かれていない. 矢印でも引いて図の中に書いてくれればずっと分かりやすいと思うのだが.

Python Tools for Visual Studio

数年前から公開され, 名前はみたことはあったが今日初めて使ってみた. もっと早くインストールしておけばよかった. PythonはWindowsなら秀丸, Linuxならviで書いてたのだけどこれは便利だ.

通常のC/C++などの開発と同様, マウスを変数の上に置くだけで中身が見える. リストも配列も自作クラスでもOK. 補完もできるしbreak pointもマウスでOK. 設定も何もしなくても考えなくても普通のVCと同じように使えてとてもやりやすい. ごく簡単な型推論(?)もしてるようで, 書いてる途中で自動的にチップ表示される関数の引数にintやstrが表示されていた.

LinuxやMac OSのリモートデバッグもできるとある. そのうち試してみる.

高速な逆数平方根の近似計算

Fast inverse square rootというのを知った.

わずか数個の演算で相対誤差数%を実現する. もともと1990年代に発見されていたが, 一般に知られるようになったのは2002年以降らしい.

1998年にはAMDのK6-2に搭載された3D Now!に逆数平方根の近似計算命令が入っていたし, SSEにもあるので今, 直接使う機会はあまりないかもしれないがエレガントな方法だ.

float rsqrt1(float x)
{
    uint32_t i;
    memcpy(&i, &x, sizeof(i)); // i = *(uint32_t*)&x;

    i = 0x5f3759df - (i >> 1); // (***)

    memcpy(&x, &i, sizeof(x));
    return x;
}

このコードはWikipediaやオリジナルはlongを使っていたものをuint32_tに直してある. それからstrict aliasing ruleを破ってるのでmemcpyに置き換えた. もう一つ, この計算の後のNewton法による近似を入れてるが今回のポイントではないのでそこは省いた.

floatのポインタをintにキャストして値を取り出すコードがうまく動かなくなるような今どきの最適化コンパイラなら, この程度のmemcpyは最適化されて消える. 実際gcc 4.7.3では上記コードは

rsqrt1:
  vmovd   edx, xmm0
  mov eax, 0x5f3759df
  shr edx, 1
  sub eax, edx
  vmovd   xmm0, eax
  ret

とmemcpyはmovd1命令に置き換えられていた.

さて, どうやってこのコードが導かれたのだろう. 上記Wikipediaやそこの参考文献に書いてあるが自分の理解のために書いてみる.

まずfloat xの値をx = (1 + M_x/L)2^{(E_x - B)}と書くことにする. ここでM_xは仮数部, E_xは指数部を表す. L = 2²³, B = 127である.

y = 1/√xとするとlog₂(y) = -(1/2) log₂(x). yの仮数部をM_y, 指数部をE_yとして代入すると

log₂(1 + M_y/L) + (E_y - B) = -(1/2) log₂(1 + M_x/L) - (1/2) (E_x - B).

ここで大胆な近似を行う. log₂の中は1以上2未満なのでlog₂(1 + x) ≒ x + s(sはある定数)としてしまう. 2とe(=2.718...)を同じものとみなしてしまうのだから凄い. するとこの式は

(M_y/L + s) + (E_y - B) = -(1/2)(M_x/L + s) - (1/2)(E_x - B)となる. 整理すると

E_y L + M_y = (3/2)(B - s)L - (E_x L + M_x)/2.

E_y L + M_yはE_yが指数部でM_yが仮数部でLをかける, つまり23ビットシフトしてorしたもの. これはyのfloatとしてのbit表現が表す整数(I_y)である. つまり

I_y = (3/2)(B - s)L - I_x/2.

(3/2)(B - s)Lは定数なので(***)が導かれた.

普通なら1bit右シフトすると指数部の最下位ビットが仮数部の最上位ビットに入るので別処理したくなるが, このコードにはそれがないのがポイント. あとは定数sを決めればよい.

sの決め方として, そのままでよい近似となるものや, Newton法を一度かけたときによくなるものなどいろいろ研究されているようだ. 自分でもざっくりと計算してみた.

log₂(1 + x)とx + sが区間[0, 1]で二乗誤差が最小となるようにsを決めてみる. Mapleだとこんな感じ.

f := int((log(x+1)/log(2)-(x+t))^2,x=0..1);
s := solve(diff(f,t)=0,t);
evalf(s,15);
> .573049591110395e-1

このsとB=127, S=2²³を使うと(3/2)(B - s)L = 0x5f34ff58となる. この数値を使ってテストしてみた[rsqrt.cpp].

なんか私がもとめた数値の方がよく見えるのだけど多分何か勘違いしてるのだろう.

『応用数理ハンドブック』(日本応用数理学会監修)

「数論アルゴリズムとその応用」研究部会の幹事の関係で私も4ページ書かせて頂いた. 項目を決めるときに, せっかくならとペアリング暗号を入れさせてもらう. 機械学習の項目もあるし折紙工学というのもある. 「こんな研究があるのか～」とぱらぱらめくると楽しそうだがいかんせん結構なお値段(かつ多分重たいので寝っころがっては見られない). どうしようかな.

サービスプロバイダからの通達

家で利用してるインターネット接続プロバイダから「重要」という封筒が来ていた. 何かと思ったらうちが著作権侵害のものを公開してる可能性があるという通達だった.

なんやそれと思ったがほっとくと接続を止める可能性があると書かれてるので無視するわけにもいかない. 通達によるとBitTorrent関連らしいのでportを見てみるけど怪しそうなものはない. そもそもBitTorrentのソフトは入れてないし.

可能性は低いけどウィルスが勝手にやってたら恐ろしいなとスキャンしてみる. もちろん何も見つからない(まあスキャンする側も感染してたらどうしようもないが).

だから「無関係だ」と連絡しようと思った. けど ここに連絡しろというフリーダイアルをネットで調べるとヒットしない. もう一つの普通電話(電話料金がかかるほう)の方はちゃんとプロバイダの電話番号だった. フリーだけが嘘の詐欺メールの可能性もあるかなあと思って普通の方に掛けて聞いてみる.

わりとすんなりと話は進んでルータのMACアドレスや割り当てられていたIPアドレスを教えたらうちじゃないと判明. 何か別のところだったらしい. しかし向こうが知ってる情報だろうに. 関係なかったので一安心だがいろいろ気持ち悪い.

『オンラインゲームを支える技術』(中嶋謙互)

私はオンラインゲームを遊んだことがないのだけども, この本は非常に面白く読めた. ゲーム開発者だけでなく, ネットワーク, インフラに関わる人でも楽しめて参考になると思う.

ネットワークが実は結構遅いこと, TCPのヘッダサイズを考慮した転送速度の見積もり, 1/60秒でCPUができること, すべきことなど具体的な数値による考察が素晴らしい.

通信レイテンシの大きい環境の中で, どうやって複数の入力によるゲームの状態を矛盾なく維持すべきか. それをプレイヤーの満足感を優先する原則の中で解こうとするのも, (ゲームにとっては当たり前なんだろうけど)新鮮な考え方.

4章以降の架空のゲーム開発における企画資料から始まり, 必要なサーバ台数やDBの性能, マシンやネットワークリソースのコストの見積もりは圧巻. できるだけ仕様を簡潔にするプロトコルの実装の原則も, 理由と共に具体的に説明されていて分かりやすい.

1章「オンラインゲームの歴史と進化」もよく調べてまとめたなあと感心する. 富士通のHabitatも紹介されていて驚いた.

整体

最近あまりに腰が痛いのでたまらず初めてのところの整体に行ってみた.

結構きつめだったせいか, 疲れがたまりすぎていたせいか, 終わってから体がだるくて2日間全然動けない.

自己紹介の練習

本社で営業のプレゼン技術向上のための勉強会があったで参加してみた. はっきりとした声を出す, 重要なところはゆっくり大きく, 最初の1分が大事など.

言葉にするとあたりまえなんだけど, ちゃんとできてないことも多いので参考になる. 実際に講義の初めと終わりにみなの前での練習があったのがよい(普段殆ど交流のない人ばかりだったので初対面に近い).

2回目は声を出そうとしたら, 何をしゃべるか分からなくなってしまった. 難しい. ただひたすら練習あるのみだそうだ.

同僚に「光成さんは大事なところでつまらなそうに話す」と言われたことがあるので改善したい.

左薬指の爪

が何かの拍子で半分に欠けてしまって, 薄皮でつながってる状態で痛い. 当たると余計に痛いのでセロハンテープや絆創膏や貼ったりするとキーボードーを打ちにくい.

打ちにくいと普段よく使ってるコマンドとかが出てこなくて入力できない. 検索しようにもなんだか分からないので調べられない. おもしろいもんだ.

結局コマンド履歴や, バイナリがあるフォルダを眺めて思い出した.

アセンブラ短歌

夕方コードを書く気分にならなくてつい詠んでしまった[nop.asm].

Windowsなら
nasm -fwin32(or -fwin64) nop.asm && link nop.obj msvcrt.lib
Linuxなら
nasm -felf32(or -felf64) nop.asm -o nop.o && gcc -m32(or -m64) nop.o

一応解説. 1byteから7byteまでのnopを使った. そうすると全部で1 + 2 + ... + 7 = 28(byte). 5 7 5 7 7の境界さえ揃ってればよいようなので適当に並べる. そうすると残り3byte. retは通常0xc3で1byteなので2byte nopを使ってもいいけど, 既に使ったのでなんとなく嫌. だから16bit即値つきretだと3byte命令になるのでret 0をそれでエンコードした.

0f 1f 44 00 00       nop    DWORD PTR [eax+eax*1+0x0]; 5byte nop

0f 1f 80 00 00 00 00 nop    DWORD PTR [eax+0x0]      ; 7byte nop

66 90                xchg   ax,ax                    ; 2byte nop
0f 1f 00             nop    DWORD PTR [eax]          ;+3byte nop

90                   nop                             ; 1byte nop
66 0f 1f 44 00 00    nop    WORD PTR [eax+eax*1+0x0] ;+6byte nop

0f 1f 40 00          nop    DWORD PTR [eax+0x0]      ; 4byte nop
c2 00 00             ret    0x0                      ;+3byte ret

何もしたくない気分をnopのれんちゃんで表現した(既出なら忘れて下さい).

デバッガでRTMを使う

x86/x64最適化勉強会6で@shinhさんがデバッガで使えないって仰ってたけど, Software Developer's Manual vol.1の14.3.7にRTM-Enabled Debugger Supportという記述があった. DR7やIA32_DEBUGCTL_MSRなどを設定するとどうのこうのとあるので, ツールなどの対応が必要だけど何もないわけじゃなさそう.

@natto_heavenさんのHaswellのTSXのテスト: 追試編では AVXとSSEをまぜるなとの指摘もあるし(上記丸あるでは14.3.8.1), ちゃんとドキュメント読まないといろいろはまりそうだ. FPUやMMXも駄目っぽい.

glibcのstrlenの実装が賢い感じ

glibcの実装ってstrlen.Sだろうか. historyをみるとSSE4.2を使っていたけどこっちの方が速いぜ, うりゃってSSE2に戻したのかな. このコードだけとってきて自分でも比較してみた(gcc-4.8.1).

% /lib/x86_64-linux-gnu//libc.so.6
GNU C Library (Ubuntu EGLIBC 2.17-0ubuntu5) stable release version 2.17, by Roland McGrath et al.

strlen(2.7.2)は標準のやつ. strlen(latest)がgitレポジトリからとってきたもの. 確かに速くなってる. 特にHaswellでは2.7.2のやつでも速い.

テクニック的にはpcmpeqbでbyte単位の比較をしてpmovmskbやってbsfで見つける方法なので方針は私のstrlenSSE2の実装と大体に近いようだ. もともと私の実装でも長い文字列のところではSSE4.2よりもstrlenSSE2の方が速い. glibcではループアンロールしてるのでより速くなってるようだ. 短いところでも速くなってるのは当時は非アライメントな領域のアクセスが遅かったので処理が面倒で手を抜いていたからかもしれない.

速くならないどころか かなり遅くなったなーと思った。はて

というのは上記の私のベンチマークだと短い文字列や長い文字列それぞれでベンチマークしてるので混ざってると分岐が沢山あるglibcより、いつも同じように処理するSSE4.2のほうが有利からかなあ(状況が分かってないので自信なし). ってそういう意味じゃなくてAVXを使ったのにという意味か.

ERMSB(Enhanced Rep MovSB)

最適化勉強会で紹介したのだけど, IvyやHaswellではrep movsb, rep stosbが高速化されている(rep stosbはHaswellだけ? 要確認).

もともとrep movsbって何だったかというと, 8086の頃からあるbyte単位でのメモリ転送命令. これはrsi(転送元), rdi(転送先), rcx(回数)を引数にとる. word単位で行うmovsw, dword単位のmovsdもある.

  rep movsb ≒ memcpy(rdi, rsi, rcx);

80386の以降, dword単位での転送が速いのでn byteコピーするならn/4回rep movsd, (n%4)回rep movsbするのがセオリーだった(x64ならrep movsqとか).

ただ初期コストが大きいとか, SSE登場以降, そちらを使った方が速かったこともあったりして, 私はあんまし使わなかった(memcpyなんてコンパイラに任せればよいようになったのもある). 特にmovsbなんて端数処理専用で, 自分で書くコードに出番は殆どなかった.

ただIntelはこれらの命令を地道に改良していて, (どの世代か忘れたけどしばらく前から)4未満のnに対してだけmovsbが高速に動作するようにしたり, 初期コストを小さくしたりするなどの改良が入っていた.

で, Ivy以降はその殆ど端数処理専用となっていたmovsbを高速化するという力業に出た. 同じバイトコードなのにERMSBと名前を変えるぐらいの力の入れようだ.

実際, eax = 7, ecx = 0でcpuidをしたecxの9bit目が立っていればERMSB, そうでなければ従来のmovsbという区別がある.

詳しい速度は下のリンク先の資料を見てもらったり自分でベンチマークをとってもらうとして, ざっくりとした感じでは4KBを越えるとAVX2によるメモリ転送よりも高速に動作する. 凄い.

何よりアライメントや端数処理をしなくてすむ(若干遅くなるが)のが素晴らしい.

注意点はドキュメントによるとDF(direction flag)=0のときだけ速いようだ. DFなんて普段意識しないので忘却の彼方だったりするけど念のため.

x86/x64最適化勉強会6レポート

会社のブログに書いてみました.

x86/x64最適化勉強会6の発表資料と録画

昨日はx86/x64最適化勉強会6でした. 今回も会場を提供してustream中継をしてくださった@kzmogiさん, (株)AXELL社のご好意に感謝します. 発表者の方々もありがとうございました. 質疑応答でも鋭い発言をしてくださる方が多くてとても勉強になりました.

とりあえず資料へのリンク.

1. @herumi : Haswellサーベイと有限体クラスの紹介 [ust-1]
2. @i_saint : Runtime C++ Editing [ust-2]
3. @tanakmura : Cortex-A15 vs Jaguar [ust-3]
4. @shinh : TSXを使うとTBDがはやくなりました！みたいな [ust-4]
5. @starpoz : Intel TSX HLEを触ってみた [ust-5]
6. @k_nitadori : Xeon PhiとN体計算コーディング [ust-6]