ようこそ暇な方

過去のごちゃまぜ最適化ネタ(未整理&未完結)--
分岐/ 並列処理/ 同時実行/ 分岐の実験/ 分岐の実験少し/ 午後のコード

1. 分岐を減らす

パイプラインが進化したCPUでは分岐は禁物です. いかに分岐を減らすかが一つの大きなポイントです.

No. 001
まず簡単なところから

func( int i ){
    if( abs(i) <= 5){
        ....

こんなソースを考えます.
普通に考えると if() の中身は i <= 5 && i >= -5 ですから, 2回条件分岐することになります.
ところが

    if( (unsigned)( i + 5 ) <= 10 ){

とすると1回の分岐ですみます.

    mov     eax,i
    add     eax,5
    cmp     eax,10
    ja      .else
    本体
.else:
    その他

とできるのです.

No. 002
前回のはあまりに簡単だったかも知れません.
もう少し複雑な例を考えてみます.
unsigned int な x, y, z, wという変数にたいし, それらが全て0か1の時のみ処理したいという場合です.

    if( x <= 1 && y <= 1 && z <= 1 && w <= 1 ){
    ....

#ちなみに午後のこ〜だのソースの一部です.
さてどうしますか?
私は次のように変形しました.

    if ((x | y | z | w) <= 1) {

無論これが1ではなく6だったりするとうまくいきません.
こんな個別論的なものを考えても無駄だという意見もあるでしょうが, 逆にいえばいかにケースバイケースに対応するかが速度の別れ目になりえます.

No. 003
もう1つこんなのはいかが?

    if( ( c >= 0x81 && c <= 0x9F ) || ( c >= 0xE0 && c <= 0xFC ) ){
    ...

いわずと知れたSJISの第1byte判定コードです.
一般的には分岐を減らすためにテーブルを作成しておく方法があります.
無論これでも構いませんが意外とメモリアクセスに時間がかかることがあります.
ここでは分岐を1回ですむように変形してみます.

[ヒント1]下の表は`O'の部分がSJISの第1byteであることを示しています.

[00] ................
[10] ................
[20] ................
[30] ................
[40] ................
[50] ................
[60] ................
[70] ................
[80] .OOOOOOOOOOOOOOO
[90] OOOOOOOOOOOOOOOO
[A0] ................
[B0] ................
[C0] ................
[D0] ................
[E0] OOOOOOOOOOOOOOOO
[F0] OOOOOOOOOOOOO...

1回の分岐ではある範囲に入っているか否かしか判別出来ません.
上の様に2つに分かれていては無理です.
[ヒント2]xorを使って分かれている部分をくっつけてみましょう.

No. 004
入力値と0x20をxorすると0x80〜0x9Fの部分が0xA0〜0xBFに, 0xE0〜0xFFの部分が0xC0〜0xDFに移り1つの連続した部分になります (分かりますよね?).
後は絶対値のところのようにすればよいです.
結局

    if( ((unsigned char)((x^0x20)-0xA1)<=0x3B) )

とすると分岐1回ですむようになりました.

No. 005
adcとsbb
この2つの命令は CF (Carry Flag) を伴う add, sub です.
つまり

    adc eax,ebx

は

    if( CF ) inc eax;
    add eax,ebx

と言う意味です.
64bit演算のためにある命令ですが, これらの命令を使うと様々な場面で分岐を減らすことが出来ます.

    sbb eax,eax

という命令は

    eax = ( CF ) ? -1 : 0

に等価です.
例えば画面表示で入力パラメータの0クリッピングをする時

   eax = ( eax < 0 ) ? 0 : eax;

というコードが出てくるでしょう.
* ちなみに A = ( B ) ? C : D; は

  if( B ){
    A = C;
  }else{
    A = D;
  }

という意味です. これは以下の様に記述できます.

    cmp     eax,(1<<31)  ;CF = ( eax >= 0 ) ? 1 : 0;
    sbb     edx,edx
    and     eax,edx

No. 006
別のやり方もあります.

    mov     edx,eax
    sar     edx,31
    or      eax,edx
    sub     eax,edx

K6-2ではsbbはvector命令なのでこのようにshort命令に分解した方がやや速いようです.
ちなみにこの最初の2行はcdqと等価です. 昔のCPUだとその方がよいかもしれません.

No. 007
-x = ~x + 1 (`~'はnot演算子)を使うことで絶対値も可能です.

    mov     eax,[x]
    cdq             ; edx = ( x < 0 ) ? -1 : 0;
    xor     eax,edx ; eax = ( x < 0 ) ? ~x : x;
    sub     eax,edx ; eax = ( x < 0 ) ? -x : x;

No. 008
要はマスク(0xFFFFFFFF)をいかにうまく作るかが鍵となります.
BM386でも使ってるレイヤー合成部分の

    out = ( image1 ) ? image1 : image2;

を考えてみましょう. image1の`0'の領域を透過させてimage2を見せるということです.
image1が0か否かでマスクを作るのですね.

    mov     eax,[image1]
    mov     edi,[image2]
    cmp     eax,1
    sbb     edx,edx      ; edx = ( eax ) ?      0 :     -1;

これでマスクを作れました. 後はくっつけるだけです.

    and     edi,edx      ; edi = ( eax ) ?      0 : image2;
    not     edx
    and     eax,edx      ; eax = ( eax ) ? image1 :      0;
    or      eax,edi      ; eax = ( eax ) ? image1 : image2;

実は後半は3命令で実現できます. 分かりますか?

No. 009
答えはこうです.

    sub     edi,eax  ; edi = image2 - image1;
    and     edi,edx  ; edi = ( eax ) ?      0 : image2 - image1;
    add     eax,edi  ; eax = ( eax ) ? image1 : image2;

先に差分をとってマスク後にもう一度足すというこのテクニックは意外に使えます.

No. 010
実戦向け(?)として午後のこ〜だのHuffmanCodeL(in loop.c)から

    変数は全てunsigned int

    idx = 0;
    if( x ){
        cbits++;
        idx = x * 16;
    }
    if( y ){
        cbits++;
        idx += y;
    }
    ... ここ以後 x, yは破棄してよい

このルーチンは頻繁に呼ばれるので出来るだけ速くしたいのです.
まずidxはx, yが0であっても分岐する必要がありません.
つまり分岐をなくすためにはcbitsの処理が鍵となります.
私は以下のように書きました.
x = eax, y = edx, cbits = ebp, idx = eaxとして組んでいます.

    mov     eax,[x]
    mov     edx,[y]
    cmp     eax,1    ;   CF = (eax) ? 0 : 1
    sbb     ebp,-1   ; ebp += (eax) ? 1 : 0
    shl     eax,4
    add     eax,edx  ;eax=x * 16 + y
    neg     edx      ;    CF = (y) ? 1 : 0
    adc     ebp,0    ;cbits += (y) ? 1 : 0

xが0でないなら1を足すのですから, 0か否かでCFをたてなければなりません.
ここでは 1と比較しました.
CFによるincrement, decrementは4パターンあります.

	inc/decしたいCFの状態
	CF=1	CF=0
inc	adc eax,0	sbb eax,-1
dec	sbb eax,0	adc eax,-1

後半のyが0か否かを判定する部分ではnegを使ってみました.
前半同様cmpを使っても良いのですが, こちらの方が命令長が小さいのと yはもう不要だったからです.

2. まとめて処理する

MMX系命令が無いから並列処理なんて出来ない, わけではありません.
ここでは16bit color [R=G=B=5bit, 最上位bitは0]での簡単な画像処理を考えてみます.

No. 011
2つの画面を50%ずつ混ぜる処理を考えます.
C言語にはない構文ですが0b...で2進数を表すことにします.

#define maskB(x) ( (x) & 0b11111 )
#define maskG(x) ( ( (x) & ( 0b11111 << 5 ) ) >> 5 )
#define maskR(x) ( ( (x) & ( 0b11111 <<10 ) ) >>10 )

    col1:1枚目の画像の色データ, col2:2枚目の画像の色データ
    tmp1 = ( maskB(col1) + maskB(col2) ) / 2;
    tmp2 = ( maskG(col1) + maskG(col2) ) / 2;
    tmp3 = ( maskR(col1) + maskR(col2) ) / 2;
    output = ( tmp3 << 10 ) | ( tmp2 << 5 ) | tmp1;

かなり重い処理になりそうです.
この z = ( x + y ) / 2 をまとめて処理することを考えます.
z = x / 2 + y / 2と変形します.
2で割る処理は当然 shift で実現します. bitをずらした時に隣のデータと交じらないようにマスクをかけます.

    output = ( ( col1 & 0b111101111011110 ) >> 1 )
            |( ( col2 & 0b111101111011110 ) >> 1 )

ただこれは元の方法と結果が異なってしまいます.
整数を2で割った時の端数 0.5を無視してしまうからです.
さてどうしましょうか.

No. 012
1つ目の方法と2つ目の方法とで値が異なってしまうのは xとyが共に奇数の時です.
そのとき ( x + y ) / 2が x / 2 + y / 2よりも1大きくなります.
つまり

    adjust = ( x & 1 ) & ( y & 1)
    z = x / 2 + y / 2 + adjust;

とすれば完全に一致します.
この様にアルゴリズムを変更すれば 32bitの変数を使って隣接する2ドットの画像データをまとめて処理できます.

#define maskA ( 0b000010000100001 * ( 65536 + 1) )
#define maskB ( 0b111101111011110 * ( 65536 + 1) )

    adjust = col1 & col2 & maskA;
    output = ( ( col1 & maskB ) >> 1 ) + ( ( col2 & maskB ) >> 1 ) + adjust;

もとの方法では1pixel辺りbit / 整数演算が15回だったのに対し, この方法では 8/2 =4回と劇的に改良されました.

No. 013
Mamiyaさんから

    output = ( (col1 & maskB) + (col2 & maskB) >> 1 ) + adjust;

の方がいいねとの御指摘がありました. うかつでした.
更に,

    output = (col0 & col1) + ( (col0 ^ col1) & maskB ) / 2;

という方法がfj.comp.lang.c++にあったそうです.
いや, 参りました.
どうやって同じであることを納得するかですが, こんな感じかな.
* 簡単にするため最小単位のみの表記です.

   m = (-1)<<1;
   ( (x & m) + (y & m) ) / 2 + (x & y & 1)
  =( (x & m) + (y & m) ) / 2 - (x & y & m) + (x & y)   ; -1 = m | 1
  =( (x & m) + (y & m) ) / 2 - (x & y & m) / 2 - (x & y & m) / 2 + (x & y) ;偶数だから
  =( (x & m) - (x & y & m) + (y & m) - (x & y & m) ) / 2 + (x & y)
  =( (x & ~y & m) + (~x & y & m) ) / 2 + (x & y)       ; x = x | ~x
  =( (x ^ y) & m ) / 2 + (x & y)                       ; x^y = x & ~y | ~x & y

[補足]
小次郎さん, さ〜さんから x + y = (x & y) * 2 + (x ^ y) を使って説明したほうが分かりやすいのではとの御指摘がありました.

No. 014
今回はホワイトアウトです.
RGBからYCbCrに移行して明るさだけ変更する考え方もありますが, ここでは単に各RGBに1を足すことを繰り返す方法です.
各色それぞれについて31でないときのみ1を足すというコードです.
* やはり簡単のため最小単位のみの表記です.
31か否かは全てのbitが立っているかそうでないかで判定できそうです.

    x = 色の5bit成分=[x4:x3:x2:x1:x0]       ;x0〜x4は各bit
    y = x >> 1;
    x &= y;         =[*:x43:x32:x21:x10]   ;x43はx4 & x3の意味  *は値不明
    y = x >> 1;
    x &= y;         =[*:*:x432:x321:x210]  ;x432はx4 & x3 & x2の意味
    y = x >> 2;
    x &= y;         =[*:*:*:*:x43210]
    x &= 0b00001;

これで x = ( x == 31 ) ? 1 : 0 を作れました.
後は元の値に1を足してこのxを引けば望みのものが出来ます.
もっとうまい方法はないのか? 以下次回

No. 015
掲示板のではincrementよりも遙かに高度な話が出ていますが, ゆっくり行きましょう (笑)
どんな思考過程を経て作ったのかを述べます.
31か否かを判定するにはやはり1を足して繰り上がりを見るのがスマートです.
問題は繰り上がると隣のbitと干渉してしまう点にあります.
それを防ぐために A=0b111101111011110 でマスクすることを考えます.
偶数の時は1を足しても繰り上がりませんから, 1を足すのではなく x & B (B=0b000010000100001) を足せばよさそうです.
ただ x & A + x & B = x ですから, 結局
(x + (x & B))>>5 で繰り上がりを判定できます.
* x + x & 1は常に偶数になることに注意してください.
全部まとめて書けば

    x += B - (( (x + (x & B)) >> 5 ) & B);

となりました.
今回は忙しいのでこれだけです. m(__)m

No. 016
以前ビデオキャプチャ機能を使ってCGの中に人間を移り込ませる遊びをやったことがあります.
まず人がいない時の部屋の画像Aを取り込みます. 次にそこに人が通った時の画像Bとの差分を検出し, 同じ部分だけ別の画像Cを透過させるという単純な原理です.
例えばCにピアノの鍵盤を書いておき, 自分の指を透過させて空中で指を振って演奏するなんてことも一応できました (面白そうでしょ).
原理は簡単ですが, いざやってみると各ドットごとに完全一致するか否かの判定ではノイズのために実用にはならず, ある程度の誤差を許容することが必要でした.
つまり,

    if( (Aの青-Bの青) < D && (Aの赤-Bの赤) < D && (Aの緑-Bの緑) < D )...

というふうにしないといけないのです.
さて, 私が以前やったときは何分非力なCPUでしたので

    if( ( A - B + 0b110001100011) & 0b110001100011000 )

とやっておりました. つまり上のD=4の場合で, 但し各bit成分の差が -3以上4以下で桁ずれは無視(^^; というおおらかなものです.
まあ実用的にはそれほど問題がなかったのですが, 今なら厳密な方法を発見できるかもしれません.
どうでしょうか.

No. 017
かなり間が空いてしまいました. でいろいろアルゴリズムを考えてみましたがたいしたのはできなかったので保留にします.

[Back]

御意見は光成滋生<herumi@nifty.com>までお願いします