mitsunari@cybozu labs

ここでは何回かに分けてx86(IA-32)，いわゆる普通のPentiumパソコンで使われている機械語の説明をする予定です．
アセンブラ，アセンブリ言語としては拙作のXbyakを使うことにしました．
理由は，

• 普通のgasやNASMによる解説はありふれていること
• Windows，Linux，Intel Macで同じソースが使え，C++だけで閉じているためアセンブラの設定に悩まなくてすむこと
• JITなどの特殊な最適化ができること

などがあります(半分は自己満足ですね)．

内容は基本的なところから始めますが，場合によってはマニアックネタに走るかもしれません．最初のうちはx86アセンブリ言語入門と重複することも多いと思います．

以下は単なる私の価値観ですが，機械語を絶対に知っておくべきものであるとは思いません．けれども何事もかじってみるのは悪くないと思います．
必要だから勉強する，不要だから勉強しないという合理的な判断は好きではありません．知りたい，興味があるから勉強してみるという素朴な動機が大事だと思います．

閑話休題

まずXbyakをインストールしましょう．XbyakはIntel，AMDなどのPentium互換CPUが搭載されたWindows Xp，Vista，Linux，Macの32bit OSで動作し，開発にはC++コンパイラを使います．Visual C++ 6 + STLport(必須)，Visual Studio 2005 Professional/Expression, Linux gcc 4.x, Intel Macなどで動作確認をしています(もしかしたらsecure OSではセキュリティレベルを下げないと動かないかもしれません)．
zipファイルを上記ページからダウンロードしてきて展開します．

VC系の場合はxbyak.dswを開いてtest0をコンパイル，Linux，Macの場合はmakeでサンプルプログラムができます．

>./test0
xbyak version=1.06
0 + ... + 0 = 0
0 + ... + 1 = 1
0 + ... + 2 = 3
0 + ... + 3 = 6
0 + ... + 4 = 10
0 + ... + 5 = 15
0 + ... + 6 = 21
0 + ... + 7 = 28
0 + ... + 8 = 36
0 + ... + 9 = 45
0 + ... + 10 = 55
0 + 0 = 0
1 + 1 = 2
2 + 2 = 4
3 + 3 = 6
4 + 4 = 8
5 + 5 = 10
6 + 6 = 12
7 + 7 = 14
8 + 8 = 16
9 + 9 = 18
call atoi("123") = 123
jmp atoi("456") = 456
OK

の様に表示されればOKです．もし，ビルドエラーになる，動かないなどがあれば環境とエラー内容を教えていただけるとありがたいです．

インストールするにはLinux，Macの場合はrootになって

make install

とすると/usr/local/include/xbyakにインストールされます．

WindowsのVisual Studio 6の場合は[オプション]→[ディレクトリ]，Visual Studio 2005 Expressionなどの場合は[ツール]→[オプション]→[プロジェクトおよびソリューション]→[VC++ディレクトリ]→[インクルードファイル]でincludeディレクトリに展開したxbyakディレクトリを追加してください．

サンプルプログラムについて

test0

1から10までの和や，足し算，atoi()の呼び出しテストです．

quantize

JPEGやMPEGで使われる量子化のサンプルです．起動して1～100の数値を入力してください．最初にC版，次にxbyakによる演算時間が表示されます．

calc

簡単な関数計算デモです．引数に"x*(x+1)"などを与えて実行してください．このデモのコンパイルにはboostのインストールが必要です．

bf

BrainfuckJIT環境です．./bf hello.bfとしてください．

今回はXbyakのインストール方法について説明しました．次回はサンプルプログラムの説明などをする予定です．

前回言い忘れましたが，このシリーズの目標はLL魂2007デモコードのchaos.cpp程度のものを読めて理解できることを考えています．
具体的には

• レジスタやスタックを理解する
• 関数を作れる
• SIMD命令の基礎を知る
• JITのメリットを理解する

あたりを目指します．なお，JITアセンブラは通常のアセンブラに比べてワンクッション概念が必要になります．もしかしたらアセンブリ言語の本当の入門としては不適切かもしれませんが，まあご了承ください．

それでは，最初に他の言語と同様，"Hello Xbyak!"を表示してみましょう．

#include "xbyak/xbyak.h"
#include <stdio.h>
 
struct HelloGenerator : public Xbyak::CodeGenerator { // (A)
    HelloGenerator() // (B)
    {
        push((int)"Hello Xbyak!"); // (C1)
        call((int)puts); // (C2)
        add(esp, 4); // (C3)
        ret(); // (C4)
    }
};
 
int main()
{
    HelloGenerator hello; // (D)
    void (*code)() = (void (*)())hello.getCode(); // (E)
    code(); // (F)
    return 0;
}

VC系ではコンソールアプリとしてコンパイルしてください．
gccでは-fno-operator-namesをつけるのを忘れないでください．これはandやorを演算子ではなく関数として扱いたいためです．

> g++ t.cpp -g -fno-operator-names
> ./a.out
> Hello Xbyak!

# C++的には

#include <stdio.h>
int main(int argc, char *argv[])
{
    argc == 1 and puts("need option");
}

ってできたのご存じでした?

それはともかくコードの説明に入ります．Xbyakでは実行時に命令をアセンブルしてバイトコード(以下マシン語と呼ぶことにします)を生成します．それを生成するためのクラスがXbyak::CodeGeneratorであり，これを継承する(A)ことでそのクラス内でx86の殆どの命令を記述できるようになります．

今後も繰り返すことになるかと思いますが，コンパイル時ではなく実行時であることに注意してください．ここではHelloGeneratorのコンストラクタにXbyak命令(面倒なので以後asmと呼びます)を記述している(B)ので，main()内のインスタンスhelloが生成されたとき(D)にアセンブルされます．

アセンブルされたマシン語はgetCode()メソッドで取得できます(E)．この返り値は生成されたマシン語へのポインタですので，適当な関数ポインタ(詳細は後述)にキャストし，その関数を呼び出すことでマシン語を実行します(F)．

以上がXbyakを使った場合の大まかな流れです．

続いて実行時の挙動を追いかけてみましょう．

VC系ではDebugモードでコンパイルし，(F)のところにブレークポイントを置いてデバッグ実行します．停止したところでVC6なら[表示]→[デバッグウィンドウ]→[混合モード]，VC8なら[デバッグ]→[ウィンドウ]→[逆アセンブル]で

19:       code();
00401696 8B F4                mov         esi,esp
00401698 FF 95 F4 FE FF FF    call        dword ptr [ebp-10Ch]
0040169E 3B F4                cmp         esi,esp
004016A0 E8 2B 14 01 00       call        __chkesp (00412ad0)

というようなウィンドウに入ってください．[F10]を一度押し，callのところで[F11]を押すと

003730A4 68 B0 01 44 00       push        offset string "Hello Xbyak!" (004401b0)
003730A9 E8 62 FA 09 00       call        puts (00412b10)
003730AE 83 C4 04             add         esp,4
003730B1 C3                   ret

が表示されると思います．左端の数値は違うかもしれません．意味も今は分からなくて結構です．ただ上記サンプルの(C1)～(C4)に対応していることを確認してください．

gccの場合はgdbを使っておっかけましょう．まず起動してgetCode()のところでブレークポイントを置きます．

>gdb ./a.out
>b HelloGenerator.getCode
Breakpoint 1 at 0x804b2b6: file xbyak/xbyak.h, line 1307.

実行した後，[n]と[Enter]でステップ実行し，code()の直前まで進みます．

(gdb) r
Starting program: /home/shigeo/Program/xbyak/a.out
Breakpoint 1, Xbyak::CodeGenerator::getCode (this=0xbffff2e0) at xbyak/xbyak.h:1307
1307                    assert(!hasUndefinedLabel());
(gdb) n
1309                    return top_;
(gdb)([Enter]を押す)
1310            } ([Enter]を押す)
(gdb) ([Enter]を押す)
main () at t.cpp:19
19              code();

ここでコードがどうなっているか見ます．その前にデバッグに見やすいaliasを作っておきます．

(gdb) define al
Type commands for definition of "al".
End with a line saying just "end".
>x/11i $pc
>end

と定義してみてください．

(gdb) al
0x8048904 <main+52>:    mov    eax,DWORD PTR [ebp-12]
0x8048907 <main+55>:    call   eax
0x8048909 <main+57>:    mov    ebx,0x0
0x804890e <main+62>:    lea    eax,[ebp-0x108]
0x8048914 <main+68>:    mov    DWORD PTR [esp],eax
0x8048917 <main+71>:    call   0x8049c2c <~HelloGenerator>
0x804891c <main+76>:    mov    DWORD PTR [ebp-0x10c],ebx
0x8048922 <main+82>:    jmp    0x8048952 <main+130>
0x8048924 <main+84>:    mov    DWORD PTR [ebp-0x110],eax
0x804892a <main+90>:    mov    ebx,DWORD PTR [ebp-0x110]
0x8048930 <main+96>:    lea    eax,[ebp-0x108]

siを使ってasmレベルでのステップ実行をします．

(gdb) si
0x08048907      19              code();
(gdb) ([Enter]を押す)
0x0804c190 in ?? ()

call eaxが実行された瞬間からXbyakが生成したマシン語に突入しています．alで確認しましょう．

(gdb) al
0x804c194:      push   0x804b44b
0x804c199:      call   0x8048734 <puts@plt>
0x804c19e:      add    esp,0x4
0x804c1a1:      ret

やはり上記サンプルの(C1)～(C4)に対応していることを確認してください．
うっかり行き過ぎたら，最初からやり直してマシン語を一行ずつ実行する感じをつかんでください．gdbを使うときはアセンブラで遊ぶ時に便利なgdb設定を参考にすると便利だと思います．上記alもこのリンク先の.gdbinitで定義されているのを拝借しました．

今回はXbyakの基本的な書き方と実行時に確認する方法の説明をしました．長くなりましたので一端区切ります．次はコードの説明に入ります．

前置きが長くなりましたが，x86(IA-32)について説明を始めます．他のCPUでも概ね似たようなものですが，アセンブリ言語(asm)を記述するときに最低限必要な知識は概ね次の四つです．

• アドレス
• レジスタ
• インストラクションポインタ
• スタック

順次説明します．

アドレス

レジスタ

インストラクションポインタ

スタック

要はCPUがやっているのは，メモリをレジスタに読みこんで演算し，結果をメモリ/レジスタに格納する，結果によってインストラクションポインタを変更する．これだけです．アセンブリ言語で開発するということはこれらの手続きを一つ一つ丁寧に記述するということです．難しいのではなく，面倒なのですね．

さて，"Hello Xbyak!"を解説します．

    push((int)"Hello Xbyak!"); // (C1)
    call((int)puts); // (C2)
    add(esp, 4); // (C3)
    ret(); // (C4)

このプログラムはCでいうところの

void hello()
{
    puts("Hello Xbyak!");
}

と同じです(関数名は関係ないですがとりあえずhelloとします)．

C1.

文字列"Hello Xbyak!"はメモリ上のstaticな領域に格納されています．これをputsの引数として与えるには，そのアドレスをスタック領域に保存します．なぜスタック領域に保存するのかについては今後説明します．とりあえずここでは関数の引数はスタックに格納する必要があることを覚えてください．

そして，データをスタック領域に保存するにはpush命令を使います．intにキャストしているのはXbyakの文法のせいです．つまりこの行は"Hello Xbyak!"が格納されたアドレスをスタック領域に確保することを意味します．スタックで述べたようにスタック領域は0に向かいますので，pushするとespの値が4(byte)減ります．デバッガでpushの前後でespの値が4減っていることを確認してください．

C2.

引数がスタック領域に確保にされたのでCのputs関数を呼び出します．Xbyakの文法上，intにキャストしています．この行が実行されるとコンソールに文字列が表示されていることをデバッガで確認してください．

C3.

C1でスタック領域に引数を渡したので，このままではスタック領域は4byte減ったままになります．スタックは貴重です．この関数ではもうこの領域は使いませんから開放する必要があります．espに4を足してスタック領域を基に戻しましょう．add esp, 4はC表記のesp += 4;と同じです．

C4.

hello()関数を抜けて元のmain()関数に戻ります．そのためにはret命令を使います．デバッガでこの命令を実行するとmainに戻ることを確認してください．

以下にgdbで追いかけたときのログを表示します．実際に自分でも実行してみると理解しやすいでしょう．

(gdb) al
0x804b190:      push   0x804a756
0x804b195:      call   0x8048710 <puts@plt>
0x804b19a:      add    esp,0x4
0x804b19d:      ret
(gdb) p $esp
$1 = (void *) 0xbffff2ac      // (*)
(gdb) si
0x0804b195 in ?? ()           // (C1)を実行した
(gdb) p $esp
$2 = (void *) 0xbffff2a8      // (*)に比べてespの値が4減っている
(gdb) ni                      // (C2)callを実行する
Hello Xbyak!                  // 文字列が表示された
0x0804b19a in ?? ()
(gdb) si
0x0804b19d in ?? ()           // (C3)を実行
(gdb) p $esp
$3 = (void *) 0xbffff2ac      // espの値が(*)に戻った
(gdb) ni
main () at t.cpp:19
19              return 0;

今回はアドレス，レジスタとスタック，puts("Hello Xbyak!")の中身を説明しました．次回はサンプルプログラムの説明を続けます．

Cからasm，逆にasmからCを呼び出すために必要な手続きを呼び出し規約といいます．たとえば前回puts()を使う際に，文字列のポインタをpushしましたが，そういう手続きのことです．
x86での規約は大きく分けて3種類あり，細かい部分についてはコンパイラ依存であることも多く，注意が必要です．しかし，ほぼすべてのコンパイラで共通である規約は簡単で覚えやすく，しかも移植性が高くなります．計算処理などOSに依存しにくい部分ではその規約のみに従って記述すれば，Windows/Linux/Macで同一のソースとすることができメンテナンスもしやすくなります．

まずは基本の規約(cdecl)をしっかりと使いこなし，必要になってから他のその他の規約を学べばよいでしょう．詳細は呼出規約(Wikipedia)などを参考にしてください．拙文でも多少触れています．

汎用性の高い規約

• C++での関数宣言にはextern "C"をつける．
• 引数の型はvoid, int(32bit)，またはポインタのみとする．構造体，浮動小数は扱わない．
• 返り値はvoid, int，またはポインタのみとする．

たとえば

double func(double a, double b);

とするよりは

void func(double *out, const double *a, const double *b);

とすれば移植性が高くなるということです(もちろんasmとのやりとりが必要な部分のみです)．原理的に後者の方がやや遅くなる可能性がありますが，高速化を目的とした処理では通常，関数内に大きなループが存在するため，この二つの差が問題になるケースはまずありません．逆に，問題になるぐらいならその前後の処理も含めてasmの対象したほうがよい可能性が高いです．

呼び出す側の規約

• funcの引数を右側からpushする．
• call funcする(Xbyakではcall(int(func));)．
• 引数の数 x 4だけespを大きくする．
• 関数の返り値はeaxに入っている．void型の場合はeaxの値は不定．

int func0(int a, char *b, double *c);

の場合は

push(ポインタcの値);
push(ポインタbの値);
push(aの値);
call(int(func0));
add(esp, 3 * 4); /* 引数が3個なので3 * 4 */

とします．引数がなければpushとaddは省略できます．たとえば

void func();

の場合は

call(int(func));

となります．

呼び出される側の規約

• 返り値はeaxに代入する(返り値がvoidならeaxの値は不定)．
• 汎用レジスタのうちebx, esi, edi, ebp, espの値は関数を抜けるときに呼び出されたときの状態に戻す．ecx, edxの値を保存する必要はない．
• ret();で関数から戻る．

関数の引数についてもう少し詳しく説明します．例としてfunc0()を考えましょう．

func0の呼び出し手続きに従ってfunc0が呼ばれたときのスタックの状況を考えます．c, b, aと順にpushされたのでスタックには大きい方から小さい方へc, b, aと並んでいるはずです．ということは

esp + 8 : c
esp + 4 : b
esp + 0 : a

なのでしょうか．実はちょっと違いまして，callを実行したとき，こっそりと実行していた命令の次の命令の先頭アドレスがスタックに格納されています．つまり正解は

esp + 12: c
esp + 8 : b
esp + 4 : a
esp + 0 : 次の命令のアドレス

となります．デバッガで確認してみましょう．

extern "C" void func(int a, int b, int c)
{
}
 
int main()
{
    func(1, 2, 3);
}

をfuncのところでステップ実行してみます．下記はcallを実行する直前でのレジスタとスタックの状況です．スタック領域に1, 2, 3が格納されていることを確認してください．

<asm>
00401B97   push        3
00401B99   push        2
00401B9B   push        1
00401B9D   call        @ILT+1400(func) (0040157d)
00401BA2   add         esp,0Ch
 
<レジスタ>
 EAX = CCCCCCCC EBX = 7FFDF000 ECX = 00000000 EDX = 00371538
 ESI = 00000000 EDI = 0012FF70
 EIP = 00401B9D ESP = 0012FADC EBP = 0012FF80 EFL = 00000212
 
<スタック>
0012FADC  01 00 00 00 02 00 00 00 03 00 00 00 00 00 00 00
0012FAEC  00 00 00 00 00 F0 FD 7F CC CC CC CC CC CC CC CC
0012FAFC  CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC

ここでcallを実行します．

<レジスタ>
 EAX = CCCCCCCC EBX = 7FFDF000 ECX = 00000000 EDX = 00371538
 ESI = 00000000 EDI = 0012FF70
 EIP = 00401B30 ESP = 0012FAD8 EBP = 0012FF80 EFL = 00000212
 
<スタック>
0012FAD8  A2 1B 40 00 01 00 00 00 02 00 00 00 03 00 00 00
0012FAE8  00 00 00 00 00 00 00 00 00 F0 FD 7F CC CC CC CC
0012FAF8  CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC

espの値が4減ってスタック領域にcallの次の命令のアドレス0x00401ba2が格納されたことがわかります．ret()でfunc()を抜けるときに，この値を参照して帰るべき場所を決定します．そのとき自動的にespも4増えています．

今回は関数の呼び出し規約について説明しました．納得するまで何度もデバッガで確認するのがよいと思います．

x86の解説をいざ始めてみると，どうもblogという媒体はやりにくいので別ページで進めることにしました．すいません．まとめ直すまでしばらくお待ちください．あと基本的なことばかり続いたので，ちょっとマニアックネタに走ってみます．

というわけで突然ですがクイズです．そこそこ高い難易度に設定したつもりですが，いかがでしょう．初心者の方は全然分からなくても大丈夫です．あえて曖昧な記述をしている部分もあります．後半の答えは凄いものがあるといいなあ．あと，難問奇問募集中．

以下は断りがない限り，

• 環境は32bit OS上のPentium4以降のx86 CPU
• 関数の呼び出し規約は__cdecl
• 配列は16byte alignmentされていて複数の配列はオーバーラップしていない
• ループは4の倍数と仮定してよい

ものとします．CPUに依存する場合は明記してください．

Q1(5点)

Q2(6点)

; input :  eax
; output:  eax
; destroy: edx
 
mov    edx, eax
sar    edx, 31
xor    eax, edx
sub    eax, edx

Q3(各4点)

/* 3.1 */
.label:
...
shr   eax, 2
jnz    .label

/* 3.2 */
.label:
...
sub   ecx, byte 1
jnz    .label

*すいません，jzはjnzのtypoでした．

Q4(8点)

Q5(8点)

magic dd 0x01800000
x dd 1.5 ; float
 
movss    xm0, [x]
paddd    xm0, [magic]

Q6(12点)

Q7(10点)

int getBit(unsigned int x)
{
    for (int i = 31; i >= 0; i--) {
        if (x & (1 << i)) return i;
    }
    return 0;
}

をbsrを使って最適化せよ．

Q8(8点)

sub    esp, 8192
mov    [esp], eax
add    esp, 8192

Q9(15点)

/*
    -1 <= in[i] <= 1とする．
*/
void func(int *out, const float *in, int n)
{
    for (int i = 0; i < n; i++) {
        out[i] = int(in[i] * 32767);
    }
}

Q10(20点)

void log2(int *dest, const int *src, int n)
{
    for (int i = 0; i < n; i++) {
#if 0
          /* わずかな誤差で意図しない値になることがあったので修正 */
//        dest[i] = (int)(log(src[i]) / log(2));
#else
        int j;
        int x = src[i];
        for (j = -1; x > 0; j++) {
            x >>= 1;
        }
        dest[i] = j;
#endif
    }
}

問題編はこちらのx86カルトクイズです．

Q1

lea    eax, [eax + eax*8]
lea    eax, [eax + eax*4]

45 = 9 * 5 = (8 + 1) * (4 + 1)を利用する．乗算命令が速くなったとはいえ，まだadd/sub/shift/leaの組み合わせを使った方が速いことが多い．

Q2

cdq
xor    eax, edx
sub    eax, edx

条件分岐は可能な限り使わない．次の二つの恒等式

-x = x ^ (-1) - (-1)
 x = x ^ 0 - 0

と

unsigned int m = (x < 0) ? -1 : 0

を組み合わせることで絶対値を

unsigned int m = x >> 31;
return x ^ m - m;

と実現する．ここでmを取得する部分は

mov    edx, eax
sar    edx, 31

であるが，これに限り32bit除算の準備のために用意されたcdqを使うことができる．

cdq ⇔ edx = eax < 0 ? -1 : 0

P4系ではどちらも同じμOPに分解されるがcdqの方がコードサイズは小さくてすむ．

Q3

/* 3.1 */
.label:
...
shr    eax, 1
shr    eax, 1
jnz    .label

Core系では二つの1bitシフトに分割したほうが速くなることがある．それまでのP4では分割しない方がよい．
P4以前では分割しない方がよい．また8倍までの左シフトに対してはshlよりもaddを使った方がよいことが多い．

3.2はdecに変えない方がよい．sub ecx, 1とdec ecxの違いはecx == 0のとき前者はCF(キャリーフラグ) = 1となるが，後者はCFが変化しない点にある．そのためdec/incはフラグに対してパーシャルレジスタストールを起こし遅くなる．バイト長が長くなるがsub/addを使うべきである．

Q4

処理内容は変わらないが，前者のコードが8B 0C 45 00000000と7byteなのに対して，後者は8B 0C 00と3byteですむ．アドレス生成回路としてもシフトが使われるか，加算が使われるかの差があるかもしれない．NASMでは前者のニーモニックを記述しても自動的に後者に変換される．XbyakでもReg32.optimize()において同様の処理を行っている．

Q5

整数として加算するのでxm0の値は0x3fc00000，floatとしては12.0である．

Q6

整数として0x01800000の加算は，xのfloatとしての指数部に3を足すことになる．したがって一般にこの操作はxのfloatとしての値を8倍することになる．これはFLT_MIN(= 0x00800000) <= |x| <= FLT_MAX(= 0x7f7fffff) / 8の範囲で正しい．

通常ならmulss xm0, [_8f]を使うが，padddの方がレイテンシが短い．ただし0は0ではなくFLT_MIN * 4になる．またデノーマル数(0 < |x| < FLT_MIN)，あるいは処理後がオーバーフローする(FLT_MAX/8 < |x|)ときも値が異なることに注意する．

Q7

mov    edx, [esp+4]
bsr    eax, edx
cmovz  eax, edx
ret

bsr命令は入力値edxが0のとき出力値eaxが不定となる仕様であるため単純には使えない．
そのため

    mov    eax, [esp+4]
    test   eax, eax
    jz     .skip
    bsr    eax, eax
    jmp    .exit
.skip:
    mov    eax, [0の時の特殊値]
.exit:

という形にしたくなる．が，bsrは入力値が0のときZF(ゼロフラグ)を1にする．これを利用すれば，ZF == 1のときのみ移動する条件付き移動命令cmovzを使って

    mov    edx, [esp+4]
    bsr    eax, edx
    cmovz  eax, [0の時の特殊値]

とできる．更に特殊値が0ならば答えのように短くできる．

Q8

Windowsではespの指す部分がコミットされていない場合，アクセスがあった時点で4KB単位でコミットされる．その際コミット領域が連続している必要がある．したがって確実にコミットさせるために

sub    esp, 8192
mov   [esp+4096], eax ; dummy write
mov   [esp], eax
add    esp, 8192

などのように4KBずつアクセス(readでも可)しておく必要がある．

Q9

floatからintへのキャストがボトルネックとなる．

通常使う四捨五入では0.5は常に切り上げとなる．そのため統計的にわずかであるが演算結果が大きくなってしまう．これを防ぐために浮動小数から整数へのデフォルトの変換命令は0.5を偶数方向に丸める四捨五入を採用している．たとえば1.5は2，2.5は2，3.5は4となる．

だがこの仕様はCの常に0方向に切り捨てる仕様と相容れない．そのためキャストごとにFPUの制御モードを変更する必要があり，そのコストが非常に大きい．

FPUの場合

cwOrg dw ; デフォルトのFPU制御ワード
cwTrunc dw ; cwOrg |= (3 << 10) ; truncate to 0
 
; ループ前
fldcw    [cwTrunc] ; 丸めモードを0への切り捨てに変更する
; ループ内
fistp    dword [...]
; ループ後
fldcw    [cwOrg] ; デフォルトモードに復元

SSEの場合

mxcsrOrg dd ; デフォルトのMXCSRレジスタ
mxcsrTrunc dd msxcrOrg |= (3 << 11) ; truncate to 0
; ループ前
ldmscsr    [mxcsrTrunc] ; 丸めモードを0への切り捨てに変更する
; ループ内
cvtps2pi [...] ; float x 2 → int x 2
; ループ後
ldmxcsr    [mxcsrOrg] ; デフォルトモードに復元

ただし，FPU制御ワードの変更に比べてmscsrの変更はコストが大きい．ループが小さい場合はSSEを使わずFPUで閉じた方がよいケースもある．

SSE2の場合

cvttps2piを使えばMSCSRレジスタを変更することなく切り捨てを行える

SSE3搭載CPUでFPUを使う場合

fisttpを使えばFPU制御ワードを変更することなく切り捨てを行える

Q10

xに対してxを超えない最大の2巾の指数を求めるにはQ7のbsrを使えばよい．
今回はx>0が保証されているため例外処理は不要である．ループアンロールなどはここではしない．

proc log2A
    mov   edx, [esp+4]
    mov   eax, [esp+8]
    mov   ecx, [esp+12]
    push  ebx
.lp:
    mov   ebx, [eax]
    add   eax, byte 4
    bsr   ebx, ebx
    mov   [edx], ebx
    add   edx, byte 4
    sub   ecx, byte 1
    jnz   .lp
    pop   ebx
    ret

またxを2^e * f(1 <= f < 2)という形で考えると，bsrはeを求めることに相当する．したがってxをdoubleに変換し，その指数部を取り出す方法もある．doubleの仮数部は11bit，指数部は52bitなので52(=32+20)bit右シフトした後，1023を引けばよい．bsrはレイテンシが大きいのでこちらの方が速い．

_1023   dd 1023, 1023, 1023, 1023
 
proc log2B
    mov        edx, [esp+4]
    mov        eax, [esp+8]
    mov        ecx, [esp+12]
.lp:
    cvtpi2pd   xm0, [eax]      ; [D1_H:D1_L:D0_H:D0_L]
    cvtpi2pd   xm1, [eax+8]    ; [D3_H:D3_L:D2_H:D2_L]
    add        eax, 16
    shufps     xm0, xm1, (3 * 64 + 1 * 16 + 3 * 4 + 1)
                               ; [D3_H:D2_H:D1_H:D0_H]
                               ; 32bit右シフト相当
    psrld      xm0, 20         ; 20bit右シフト
    psubd      xm0, [_1023]
    movaps     [edx], xm0
    add        edx, 16
    sub        ecx, byte 4
    jnz        .lp
    ret

なお，もし0 < x < (1 << 23)という条件があるならばdoubleではなくfloatに変換することで更に効率よく処理を行える．floatの仮数部は8bit，指数部は23bitなので23bit右シフトした後，127を引けばよい．

_127 dd 127, 127, 127, 127
 
proc    log2C
    mov        edx, [esp+4]
    mov        eax, [esp+8]
    mov        ecx, [esp+12]
.lp:
    cvtpi2ps   xm0, [eax]      ; [*:*:F1:F0]
    cvtpi2ps   xm1, [eax + 8]  ; [*:*:F3:F2]
    punpckldq  xm0, xm1        ; [F3:F2:F1:F0]
    add        eax, 16
    psrld      xm0, 23
    psubd      xm0, [_127]
    movaps     [edx], xm0
    add        edx, 16
    sub        ecx, byte 4
    jnz        .lp
    ret

よく「アセンブラを使っても最近の最適化コンパイラに勝てない」という意見が聞かれるが，それは単に「自分では確認せず，おうむ返しをしている」か，「コンパイラと同じ(or それ以下の)コードを書いている」かのどちらかであろう．もちろん「コストに見合わない」こともあるがそれはまた別の話である．

一般に最適化ではQ6, 7, 10のような，コンパイラで自動生成されることのない部分を考慮しつつ，さぼれるだけさぼれるロジックを考える．そしてそういうロジックを自在に記述するためにアセンブリ言語やイントリンジック関数を使う．実際，これらのテクニックは音声や動画のコーデックの開発でしばしば使われる．その際Q3などの知識がクリティカルとなることはまずない．瑣末な知識を覚えるよりも(無論覚えるに越したことはない --- 実践に暗記は必要である)，常に代替ロジックが無いかを考える方がよりよいコードに繋がると思う．

ペアリング暗号とは2000年頃に登場し，近年盛んに研究されている暗号で，ペアリングとはその暗号で使われる演算の名前です．
この暗号は現在まだ普及していませんが，既存のブロック暗号や，公開鍵暗号ではできないことも実現可能になるため，今後いろいろな分野で実用化されると想像されます．

ただ当初は演算処理が極めて重く，私も以前高速な実装を試みたりしていましたが，なかなか大変なものでした．
ところが数年前からかなり効率のよいアルゴリズムが考え出され，未踏ユースで挑戦してる人も現れ，今ではRSAと同じかそれ以上の速度が出るようにもなってきているようです．

とはいえ，ソースが公開されているものは殆ど見かけません．
というわけで，今回自分で一から開発した実装を公開しました．
速度優先のため，あまりきれいなものではありませんが，何かの参考になればと思います．
コードはESSE3を利用しているためCore 2 Duo専用ですが，Core 2 Duo 1.8GHzで288usecと世界最速レベルだと思います(今までの最速は400usec程度? 違ったらすいません)．

実装の詳細についてはSEA & FSIJ 合同フォーラムでも多少触れる予定です．

2008/4/16にSEA & FSIJ 合同フォーラムでビット演算による最適化の妙味とJITアセンブラというタイトルで発表したときの資料を公開します(pptx, pdf)．
鋭いご質問，ご指摘をいただきありがとうございました．
またいろいろ議論できたらと思います．

よろしくお願いします．

ビットリバースのところで，『ハッカーのたのしみ』(ヘンリー・S・ウォーレン，ジュニア)に短いアルゴリズムがあったよと教えていただいたのでチェックしました(昔にちらっと見ただけなので忘れてました)．
それは

typedef unsigned int uint32;
 
#ifdef WIN32
uint32 bsr(uint32 x)
{
    __asm {
        bsr eax, x
    }
}
#else
 
#define bsr(x) __builtin_clz(x)
 
#endif
 
void put(uint32 x)
{
    for (int i = 0; i < 32; i++) {
        putchar((x & (1 << (31 - i))) ? '1' : '0');
    }
    printf("\n");
}
 
int main()
{
    uint32 x = 0;
    int len = 31;
    for (int i = 0; i < 32; i++) {
        put(x);
        uint32 s = 31 - bsr(~x);
        x = ((x << s) + (1 << len)) >> s;
    }
}

という感じのコードで，これはこれで面白いです．ただ今回のケースに限っては資料に書いた方法の方が命令数が短くて速いですね．

SEA & FSIJ 合同フォーラムでビット演算による最適化の妙味とJITアセンブラの中でデモに使ったVMを紹介します．
JITの紹介のために前日に2時間ででっちあげたVMなので本当に小さい(200行程度)ですが，エッセンスは楽しめるかなと思います．

ソースはXbyak.zipです．この中のxbyak/sample/toyvm.cppが今回作ったVMです(Win, Linuxと多分Intel Macでも動きます)．
このサンプルはフィボナッチ数列を計算して表示するだけのものです．
ここではどのように作ったかの説明をします．一つ前のエントリの資料も参考にしてください．

話の流れ

1. toyVMのスペック，命令セットと命令フォーマットを決める
   
2. toyVMのアセンブラを作る
   
3. toyVMの実行部分を作る
   
4. toyVM用のフィボナッチ数列プログラムを作って実行する
   
5. toyVMのマシン語をx86に変換するリコンパイラを作る
   
6. パフォーマンスを見る
   
7. リコンパイラを改良する
   

1. toyVMのスペックを決める

• 32bitレジスタA, Bの二つ．あとPC(program counter)
  
• メモリは4byte単位のみでのアクセス．4byte x 65536
  
• すべての命令は4byte固定
  
• 即値は全て16bit
  

命令群は次のものを用意しました．

メモリとレジスタの間の転送命令と加算，減算命令にレジスタの内容を出力する命令と，分岐命令だけです．
スタックが無いのでcall/retもありません．興味があれば自分で作ってみるのもよいかもしれません．
命令は4byteのうち先頭1byteに命令種別(code)，次の1byteにレジスタ種別(r)，最後の2byteに即値(imm)を入れることにします．
使われない場合は全て0にします．(code, r, imm)のペアと4byteデータの変換方法は次のようにします．簡単ですね．

void decode(uint32& code, uint32& r, uint32& imm, uint32 x)
{
    code = x >> 24;
    r = (x >> 16) & 0xff;
    imm = x & 0xffff;
}
void encode(Code code, Reg r, uint16 imm = 0)
{
    uint32 x = (code << 24) | (r << 16) | imm;
    code_.push_back(x);
}

2. toyVMのアセンブラを作る

void vldi(Reg r, uint16 imm) { encode(LDI, r, imm); }
void vld(Reg r, uint16 idx) { encode(LD, r, idx); }
void vst(Reg r, uint16 idx) { encode(ST, r, idx); }
void vadd(Reg r, uint16 idx) { encode(ADD, r, idx); }
void vaddi(Reg r, uint16 imm) { encode(ADDI, r, imm); }
void vsub(Reg r, uint16 idx) { encode(SUB, r, idx); }
void vsubi(Reg r, uint16 imm) { encode(SUBI, r, imm); }
void vjnz(Reg r, int offset) { encode(JNZ, r, static_cast(offset)); }
void vput(Reg r) { encode(PUT, r); }

3. toyVMの実行部分を作る

void run()
{
    uint32 reg[2] = { 0, 0 }; // A, B
    const uint32 end = code_.size();
    uint32 pc = 0;
    for (;;) {
        uint32 code, r, imm;
        decode(code, r, imm, code_[pc]);
        switch (code) {
           ...
        }
        pc++;
        if (pc >= end) break;
    } // for (;;)
}

switch (code) {
case LDI:
    reg[r] = imm;
    break;
case LD:
    reg[r] = mem_[imm];
    break;
case ST:
    mem_[imm] = reg[r];
    break;
case ADD:
    reg[r] += mem_[imm];
    break;
case ADDI:
    reg[r] += imm;
    break;
case SUB:
    reg[r] -= mem_[imm];
    break;
case SUBI:
    reg[r] -= imm;
    break;
case PUT:
    printf("%c %8d(0x%08x)\n", 'A' + r, reg[r], reg[r]);
    break;
case JNZ:
    if (reg[r] != 0) pc += static_cast<signed short>(imm);
    break;
default:
    assert(0);
    break;
}

void fibC(uint32 n)
{
    uint32 p, c, t;
    p = 1;
    c = 1;
lp:
    t = c;
    c += p;
    p = t;
    n--;
    if (n != 0) goto lp;
    printf("c=%d(0x%08x)\n", c, c);
}

Fib(int n)
{
    vldi(A, 1); // c = 1
    vst(A, 0); // p = 1
    vldi(B, n);
    vst(B, 2); // n
// lp
    vst(A, 1); // t = c
    vadd(A, 0); // c += p
    vld(B, 1); // mem[1]の値をBを経由してmem[2]に移動する
    vst(B, 0); // p = t
    vld(B, 2);
    vsubi(B, 1);
    vst(B, 2); // n--
    vjnz(B, -8); // PCを8減らせばlpのところにもどる．
    vput(A);
}

Fib fib(10);
fib.run();
>A      144(0x00000090)

void recompile()
{
    push(ebx);
    push(esi);
    push(edi);
 
    const Reg32 reg[2] = { esi, edi };
    const Reg32 mem(ebx);
 
    xor(reg[0], reg[0]);
    xor(reg[1], reg[1]);
    mov(mem, (int)mem_);
    const uint32 end = code_.size();
    uint32 pc = 0;
    uint32 labelNum = 0;
    for (;;) {
        uint32 code, r, imm;
        decode(code, r, imm, code_[pc]);
    L(toStr(labelNum++));
        switch (code) {
        ...
        pc++;
        if (pc >= end) break;
    } // for (;;)
 
    pop(edi);
    pop(esi);
    pop(ebx);
    ret();

    switch (code) {
    case LDI:
        mov(reg[r], imm);
        break;
    case LD:
        mov(reg[r], ptr[mem + imm * 4]);
        break;
    case ST:
        mov(ptr [mem + imm * 4], reg[r]);
        break;
    case ADD:
        add(reg[r], ptr [mem + imm * 4]);
        break;
    case ADDI:
        add(reg[r], imm);
        break;
    case SUB:
        sub(reg[r], ptr [mem + imm * 4]);
        break;
    case SUBI:
        sub(reg[r], imm);
        break;

    L(toStr(labelNum++));
        switch (code) {
        ...
        case JNZ:
            test(reg[r], reg[r]);
            jnz(toStr(labelNum + static_cast<signed short>(imm)));
            break;

.lp:
  mov   dword ptr [ebx+4],esi 
  add   esi,dword ptr [ebx] 
  mov   edi,dword ptr [ebx+4] 
  mov   dword ptr [ebx],edi 
  mov   edi,dword ptr [ebx+8] 
  sub   edi,1 
  mov   dword ptr [ebx+8],edi 
  test  edi,edi 
  jne   .lp

const Reg32 memTbl[] = { eax, ecx, edx };
const size_t memTblNum = NUM_OF_ARRAY(memTbl);
for (size_t i = 0; i < memTblNum; i++) xor(memTbl[i], memTbl[i]);

case ADD:
    if (imm < memTblNum) {        //
        add(reg[r], memTbl[imm]);    // 追加部分
    } else {                         //
        add(reg[r], ptr [mem + imm * 4]);
    }
    break;

.lp:
   mov         ecx,esi 
   add         esi,eax 
   mov         edi,ecx 
   mov         eax,edi 
   mov         edi,edx 
   sub         edi,1 
   mov         edx,edi 
   test        edi,edi 
   jne         .lp

最近，Visual C++ のことを高機能なマクロアセンブラだと思っている光成です．
その考えを64bit Windows/Linuxにも押し進めるため，64bitに対応したJITアセンブラXbyakを公開しました．
64bit Visual Studioではインラインアセンブラが廃止されたため，何かと便利になるのではないかと思います．

ところでWikipediaのAMD64などには64bit Windowsに関して

64ビットアプリケーションではx87命令・MMX命令及び3DNow!命令をサポートしない（x87レジスタをコンテキストスイッチの際にセーブしない）

テスト方法
test_mmx.cppをコンパイルして(binary)
コマンドプロンプトを二つ開いて

test_mmx 1
test_mmx 3

と実行します．それぞれスレッドを起動して，1, 2, 3, 4という数値が表示され続けました．

もしMMXレジスタがコンテクストスイッチで保存されていなければ全て同じものになるはずなので，保存されていると推測されます．

Legacy Floating-Point Supportにも

The MMX and floating-point stack registers (MM0-MM7/ST0-ST7) are preserved across context switches

とあるので保存されているように思われます．

一つ気になるのは私がテストしたのは32bit WindowsのVMware上の64bit Xpだったという点です．素の64bit環境でテストした場合どうなるかどなたか試していただけますでしょうか．

追記
奥さんから英語版Wikipediaには使えると書いてますよという情報をいただきました．
#見るの忘れてました．迂闊だった…．

strlen()とmemchr()のSIMD版を作ってみました．

今回は最速よりもお手軽さを重視したのでアセンブリ言語ではなくintrinsic関数を使っています．そのためVisual Studio 2008, gcc 4.xの両方でコンパイルでき32-bit, 64-bit OS上で動作します．
WindowsとLinuxでのみ確認していますが恐らくIntel Mac OS X上でも動作するでしょう(sample source)．

ベンチマークはランダムな長さの文字列の平均長(average length)を変化させつつ取りました．数値は1byteあたりにかかった処理時間比で小さいほど速いことを表します．
strlenが3種類(ANSI, BLOG, SSE2)とmemchrが2種類(ANSI, SSE2)あります．BLOGというのは今回試してみようというきっかけになったCounting Characters in UTF-8 Strings Is Fastの中のCバージョンです．

結果は極端に短いところではSIMDが使えず，かつそのチェックのためのオーバーヘッドがありそれほど速くなりませんが，32文字以上ではは2～5倍高速であることが分かります．
なお，gccのstrlen()は今となっては遅い命令を使っているためよくありません(なぜこうなっているのかちょっと理解に苦しみます)．

一つ細かい注意点としては，高速化のために与えられた範囲外の空間を少しアクセスします(readはするがその値は捨てます)．16byteアライメントされたところを超えることはないのでバグで16byte超えることがあったのですが修正しましたのでページングエラーにはなりませんが，purifyなどのツールを使うと警告がでるかもしれません．

intrinsic関数はVCとgccとで共通化され，しかも32bit, 64bit両対応なため汎用性も高いです．strlenSSE2のように30行未満の関数でも十分は効果はでます．利用しない手はないと思われます．

size_t strlenSSE2(const char *p)
{
    const char *const top = p;
    __m128i c16 = _mm_set1_epi8(0);
    /* 16 byte alignment */
    size_t ip = reinterpret_cast(p);
    size_t n = ip & 15;
    if (n > 0) {
        ip &= ~15;
        __m128i x = *(const __m128i*)ip;
        __m128i a = _mm_cmpeq_epi8(x, c16);
        unsigned long mask = _mm_movemask_epi8(a);
        mask &= 0xffffffffUL << n;
        if (mask) {
            return bsf(mask) - n;
        }
        p += 16 - n;
    }
    /*
        thanks to egtra-san
    */
    assert((reinterpret_cast(p) & 15) == 0);
    if (reinterpret_cast(p) & 31) {
        __m128i x = *(const __m128i*)&p[0];
        __m128i a = _mm_cmpeq_epi8(x, c16);
        unsigned long mask = _mm_movemask_epi8(a);
        if (mask) {
            return p + bsf(mask) - top;
        }
        p += 16;
    }
    assert((reinterpret_cast(p) & 31) == 0);
    for (;;) {
        __m128i x = *(const __m128i*)&p[0];
        __m128i y = *(const __m128i*)&p[16];
        __m128i a = _mm_cmpeq_epi8(x, c16);
        __m128i b = _mm_cmpeq_epi8(y, c16);
        unsigned long mask = (_mm_movemask_epi8(b) << 16) | _mm_movemask_epi8(a);
        if (mask) {
            return p + bsf(mask) - top;
        }
        p += 32;
    }
}

追記(2008/6/13)
64bit Windowsで_BitScanForwardの引数の型を勘違いしていたので修正．
ソースにライセンスを修正BSDライセンスと明記．

追記(2008/7/16)
strlenSSE2でループアンロールのため32byte単位で読み込んでいたため，ページングの境界を超えてしまうことがあったバグを修正(thanks to egtraさん)．