憂鬱な午後のひととき

6月30日_
「柳家三三で北村薫。２０１２」
前回に続き2回目. 雰囲気に慣れたので落ち着いて楽しめた. プロによるクラリネットの演奏もよかった. 落語の演目が鰍沢なのも期待通り. 作中の謎出題者の出題前後の会話自体が謎解きに関わる多重構造はやっぱり凄い. ミステリとして一級品だと思う.
6月29日_
『鹿男あをによし』(万城目学)
久しぶりに再読. やっぱり面白い. 初読はドラマを見た後で, ドラマの印象が強くて「かりんとう」先生は女性の方がいいなあ. ドラマも秀逸だったのでまた見たい. 特に剣道の試合. たくさんの鹿が走るシーンもよかったな. というか奈良を散歩したい.
6月27日_
「都の排ガス規制強化後 脳卒中死亡率8.5%減」(朝日新聞6/26)
2006年に東京都がディーゼル車の排ガスの規制を強化した. それで大気中のPM2.5という汚染物質の濃度が低下した. そして, 規制強化前後33カ月の脳卒中の死者数は20460人から19728人に減ったらしい.
脳卒中死者数自体の変化や, 気温など様々な外部要因を除去すると, 東京都の1年あたり632人の脳卒中死を防いだことになるそうな.
ごく最近の規制でそんなに効果があるとは驚きだ. 600人ってすごい人数だと思うが. 東京ってよっぽど空気が悪いのか.
6月26日_
シフト量に負の数指定して良かったでしたっけ?
みんなつっこんでそうですけど, exponent <= 0のときの1 << -exponentなので問題ないのでは.
6月25日_

rank()で遊ぶの続き.

前回は128bitあたり64bitの追加データであった.1bitあたり0.5bitの増加というわけでかなり大きい. もっと容量を減らすことを考える.

最後は64bitのpopcntを使うとして, 途中何個かをpopcntでつなげるか, テーブル引きするか. いずれにせよ素朴にやると何回足すかの予測困難な分岐が発生する. marisa-trieのbit-vectorはmod 8で分岐しているようだ.

今回は別の方法を考えてみよう. popcntの128bit版なんてあれば嬉しいがそんなものはない. SIMD命令の中で代用できそうなものが無いか考えてみる. で, 思いついたのがpsadbw. コーデック屋さんじゃなければほとんど使う機会のない命令じゃなかろうか.
# 私はコーデック屋さんではないので思い出すのに苦労した.
Packed Sum of Absolute Differences (Byte to Word)の略で(uint8_t a[8], uint8_t b[8])からΣi|a[i] - b[i]|を求める.

片方を0にしておけば8個のuint8_tを足してくれる. つまりbyte単位でのmaskをしてからこの命令を呼べば最大8個までの和が求められる.

128bitレジスタなら16個!…と実装し始めて二つの点ではまった.

mask(n) = (1ULL << n) - 1;相当のことをSSEでやろうとすると0xff...ffをmod 16の余りrで右シフトすればよい. が, byte単位の128bit右シフト命令psrldqはシフト数を即値でしか取らない!

他のたとえばpsllqは即値とSSEレジスタの二つを選べる. psrldqでは即値しか設定できない. ということはシフトする数が実行時まで決まらない場合は 

__m128i shift_right_byte(__m128i x, int n)
{
    switch (n) {
    case 0: return x;
    case 1: return psrldq(x, 1);
    case 2: return psrldq(x, 2);
    case 3: return psrldq(x, 3);
    ...
    case 15: return psrldq(x, 15);
    }
}

と書かないといけないということだ. 条件分岐を減らそうとしてシフトで分岐が必要とは, そんな馬鹿な.
# ちなみにIntelの最適化マニュアルにはこのギャグみたいなコードが載ってたりする(Example 10-19 in sse4i_atol).

もう一点はそもそもpsadbwはSSEレジスタに対しては16個のuint8_tに対して総和をとってくれるのではなく(8個のuint8_t) x 2を実行する命令ということ.
# こちらは完全に失念していた(この命令使うの10年振りぐらいだし).
って検索したら, まるもさんが10年前に文句を書いていてちょっと笑った.

8個でしか使えないならpsllqじゃなくてpsrlqを使えばよいので即値問題も解決する.

結局512bitデータ毎の累積和(32 or 64bit)とそこからの64bitのpopcntの結果8bit x 8個を持てばよいということになった. 512bitに対して128bitの増加なので1bitあたり1/4bitとちょっと減った. 実装コードは次の通り(rank.hpp). 

inline uint32_t rank1(size_t idx) const
{
    const Block& blk = blk_[idx / 512];
    size_t ret = blk.rank;
    uint64_t q = (idx / 64) % 8;
    V128 vmask;
    vmask = pcmpeqd(vmask, vmask); // all [-1]
    V128 shift((8 - q) * 8);
    vmask = psrlq(vmask, shift);
    V128 v = V128((uint32_t*)blk.s8);
    v = pand(v, vmask);
    v = psadbw(v, Zero());
    ret += movd(v);
    uint64_t mask = (uint64_t(2) << (idx & 63)) - 1;
    ret += popCount64(org_[idx / 64] & mask);
    return ret;
}

/* テーブル初期化部 */
 uint32_t r = 0;
 size_t pos = 0;
 for (size_t i = 0; i < tblNum; i++) {
     Block& b = blk_[i];
     b.rank = r;
     for (size_t j = 0; j < 8; j++) {
         uint64_t v = pos < blkNum ? blk[pos++] : 0;
         uint64_t s8 = popCount64(v);
         r += s8;
         b.s8[j] = static_cast<uint8_t>(s8);
     }
 }

ベンチマークはこちら. 

--my_rank
 Total Time(count:10000000): 131666379.000000
 Throughputs(query/sec): 0.075950
 Unit Speed(sec/query): 13.166638clk
my chk=168115228

@maropuさんのが18clkぐらいで, 前回のasm修正版が12clkだったからいい感じ. ビット演算やメモリレイアウトをまじめに頑張れば数clkぐらい減らせるのかもしれない.

命令数が多少増えるが1024bit毎の累積和と128bitのpopcntの結果8bit x 8というデータの持ち方もある. キャッシュ効率も上がるので検討する価値はあるだろう.

6月24日_

22日のソースコード, VC11の場合: #if 0のとき 

  and     dl, 64 ; m
  lea     rcx, QWORD PTR [rcx+rcx*2]
  mov     r8d, 2
  lea     r9, QWORD PTR [rax+rcx*8]
  movzx   ecx, bl
  and     cl, 63
  shl     r8, cl
  mov     rcx, QWORD PTR [r9+8]
  dec     r8  ; mask
  neg     dl
  sbb     rdx, rdx ; m = (idx & 64) ? (-1) : 0
  and     rcx, rdx
  mov     rax, rdx
  and     rcx, r8
  or      rax, r8
  and     rax, QWORD PTR [r9]
  popcnt  rcx, rcx
  popcnt  rbx, rax

#if 1のとき 

  and     r8d, 64  ; idx & 64
  cmove   rdi, rdx ; m0 ← mask
  test    r8, r8   ; !(idx & 64)
  cmovne  rbx, rdx ; m1 ← mask

ほぼ期待通りのコード生成. #if 1のときはtest r8, r8は不要なのに惜しい.

6月23日_

まるも製作所のB-CAS関連の記事がご本人によりすべて削除された.

先日のx86/x64最適化勉強会での彼の資料やustreamへのリンクも削除しました. ご本人がそうせざるをえなくなってしまったことはとても残念だ.

詳細は知らないので単なる個人的な印象になってしまうが, 技術検証や議論をしていたページがそうなる風潮はきな臭い. 彼の活動の場の素早い復帰を望む.

6月22日_

rank()で遊ぶ.

先日のx86/x64最適化勉強会で簡潔データ構造に使うrankの高速化の話があったので私も試してみた. コードはrank.hpp. なお, まだ作りかけだしメモリ食いなのであしからず.

簡潔データ構造自体の説明はリンク先の資料などを参照. rank1(idx)は与えられたbit列の[0..idx]の範囲に入っている1の個数を数える関数.

rank.hppでは128bitごとのビットの数の累積数をBlock構造体のrankメンバに, 元の128bitのデータをb0, b1メンバ保持した配列を作る. 8byte境界に揃えるためにpaddingが入ってる. 本当はキャッシュラインに揃えるべきかも. でもpaddingで大きくなりすぎなのでやめた. 

  +---------+---------+-----+-----+
  |    64bit|    64bit|32bit|32bit|
  |       b0|       b1|rank |  pad|
  +---------+---------+-----+-----+

与えられたidxに対して128bitで割ったところのblockのrankを取り出し, 残り(idx % 128)の分だけ[b1:b0]の下位ビットをマスクしてpopcountを使う. 

uint32_t rank1(size_t idx) const
{
    const Block& blk = blk_[idx / 128];
    size_t ret = blk.rank;
    uint64_t r = idx % 64;
    uint64_t b0 = blk.b0;
    uint64_t b1 = blk.b1;
    uint64_t mask = (uint64_t(2) << r) - 1;
#if 0
    uint64_t m0 = -1;
    uint64_t m1 = 0;
    if (!(idx & 64)) m0 = mask;
    if (idx & 64) m1 = mask;
#else
    uint64_t m = (idx & 64) ? (-1) : 0;
    uint64_t m0 = mask | m; // (idx & 64) ? (-1) : mask;
    uint64_t m1 = mask & m; // (idx & 64) ? mask : 0;
#endif
    ret += popCount64(b0 & m0);
    ret += popCount64(b1 & m1);
    return (uint32_t)ret;
}

実際にベンチマークをとってみるとなんだか遅い(https://github.com/herumi/dbitv). 比較対象の"--rank"は@maropuさんので"--my_rank"が今回のやつ. ちなみに先日の勉強会の発表によると@maropuさんのrankはmarisa-trieの数倍速いらしい. @maropuさんのは256bitずつとより省エネ型なのに, メモリ食ってるmy_rankがそれと変わらないというのはなんだろう. 

--rank
 Total Time(count:10000000): 196085440.000000
 Throughputs(query/sec): 0.050998
 Unit Speed(sec/query): 19.608544clk
   chk=168115228
--my_rank
 Total Time(count:10000000): 183576868.000000
 Throughputs(query/sec): 0.054473
 Unit Speed(sec/query): 18.357687clk
my chk=168115228

-Sでasmコードを眺めるとマスク生成のところで条件分岐コードが挿入されている(下記***の部分). ここはcmov命令一択だろうと思ってので意外. 

   and edx, 64
   je  .L3209     ; ***
   and r10, rcx
   mov rcx, -1
   popcnt  r15, r10
.L3209:
   and r9, rcx
   add r8, r15
   add rax, 1
   popcnt  r9, r9
   add r8, r9
   add ebx, r8d

rank1()のコードをいかにもcmovを使ってくれな感じの#if 1にしてもわざわざjmpしてくれる. 仕方がないので-Sオプションで出力した*.sを直接手で修正することにした.
# gccのインラインアセンブラにはまるよりは手書きの方がましだろうと思う人.

-O3 -fomit-frame-pointerで関数はループの中に埋め込まれているからレジスタを壊さないように, 余計にレジスタを使わないように修正するのがちょっとだけめんどい. けどまあ, これぐらいならたいしたことはない. 


   sub rcx, 1
   ; ---------------
   or r15, -1
   and edx, 64
   cmovz r15, rcx
   cmovz rcx, rdx
   and r15, r9
   popcnt r15, r15
   and rcx, r10
   popcnt rcx, rcx
   add r8, r15
   add r8, rcx
   add rbx, r8
   add rax, 1

でベンチマークを取り直したら速くなった. 

--rank
 Total Time(count:10000000): 197397900.000000
 Throughputs(query/sec): 0.050659
 Unit Speed(sec/query): 19.739790clk
   chk=168115228
--my_rank
 Total Time(count:10000000): 123700761.000000
 Throughputs(query/sec): 0.080840
 Unit Speed(sec/query): 12.370076clk
my chk=168115228

6clkぐらい速くなったということはgccのもとの生成コードは3回に1回分岐予測が外れてる感じか(適当).

でも一桁clkにはいかないんだなあ. 256bitにするときにどれぐらい遅くならないようにするかがポイントか. ただマスク生成がよりややこしくなるので難しそう.

6月21日_
神田・公開鍵暗号ワークショップ参加
set unionやMQの拡張のMultivariate Cubicとかsignature hijackingとかいろいろ面白かった. それとは別にしばらく前から思っていた自分のしょうもない疑問も答えを教えて頂いた. ありがとうございます.
6月20日_

絵で見るペアリング

登場人物

有限体Fpとは歯車のことで, Fp上の演算とは歯車をぐるぐる回す感じ.

gear1

一つ目の図は歯が7個あるのでp=7に相当する. 歯車を回して決められた場所にぴったり合わせることを離散対数問題(DLP)を解くという. pが大きいときはDLPを解くのが難しい. pが1024ビットぐらいの素数のとき1024ビットRSAを解くぐらい難しいと言われている(演算コストは280ぐらいで今のところ無理).

歯車を複数組み合わせてより複雑なものを作ることもできる.

gear2

二つ目の図はFpを6個組み合わせたもの. これをFpの6次拡大体という. それぞれの歯車の状態の組み合わせで一つの状態を表すので状態数はp6個となる.

楕円曲線とは浮輪みたいなもの.

torus

浮輪をぐるぐる回して決められたある場所にぴったり合わせることを楕円離散対数問題(ECDLP)を解くという. 浮輪に書かれた"ます目"にもFpが使われる(ここかなり説明のギャップがある). ます目の数pが大きいとやっぱりECDLPを解くのは難しい.

ECDLPの方がDLPより難易度が高くて, 同じぐらいの難しさの問題なら歯車の数よります目の数が少なくてすむ. 一般に1024ビットのDLPは140～160ビットのECDLPと同じぐらいと言われている. ECDLPが難しいことを利用した暗号を楕円曲線暗号という.

pが小さい方が計算量やメモリが少なくてすむので嬉しい. ICカードでは楕円曲線暗号がよく使われる.

ここで浮輪の2点を決めると歯車の組み合わせの1点が決まる対応があることが知られている. これをペアリングと言って, なかなかおもしろい性質がある. この対応を利用した暗号をペアリング暗号という.

pairing

この図では浮輪の2点から歯車6枚への対応を表している.

ペアリング暗号を作るにはDLPとECDLPの両方とも解けないようなパラメータを選ばないといけない. 逆に言うとどちらか一方さえ破れればペアリング暗号は破れる.

さて歯車の話に戻るとpが大きいと計算が大変なのでとても小さいp(たとえば3)を使って, その代わり歯車を97個も組み合わせたりする手法がある.

これを使って浮輪を作る. 浮輪の点の個数は397～2153で150ビットぐらいのECDLPになる. 点の行き先はその歯車6枚なので(397)6～2920ビットぐらいのDLPになる.

ECDLPとしては150ビットぐらいあるし, DLPとしても920ビットぐらいあるので1024ビットにちょっと足りないけど「まあそこそこ安全じゃね?」と思われていた.

が, DLPとしては実は253程度の演算コストでできるということをISPEC2012で発表されていて, それを実証してみたというのが今回の発表だそうな.

ちなみに今回の攻撃手法はpが大きい素数のときには使えないのでたとえば2年前に実装したやつはまだ大丈夫らしい.

6月19日_

ビットマスク生成時のよくある落とし穴

与えられた整数nだけならんだビットマスクを作ることを考えよう. たとえばn = 2なら2進数で0b11(=3), n = 5なら0b11111(=31). 2n - 1とすればよいので 

int mask(int n)
{
    return (1 << n) - 1;
}

となる. 実際にn = 0のとき0b1, n = 1のとき0b11, ..がほしいならmask(n + 1)をつかえばよい.

ただこれでは二つ問題点がある.

一つは符号付き整数の演算結果がその整数の扱える範囲を超えるとアウト. 上記で言えばsizeof(int) = 4だとするとn = 31のとき 1<< 31 = 0x80000000となってINT_MAXを超えてしまう. 

unsigned int mask(unsigned int n)
{
    return (1U << n) - 1;
}

とすべき. もう一つはn = 32のときである. シフトする数がその整数のビット長以上のシフトはアウト. 実際x86では32bitレジスタのシフト命令のシフトする数は下位5bitしか使われない. つまりn = 32ならn = 0と同じに扱われてしまう. ただしCのコードでどう処理されるかはわからない. 

>cat t.cpp
#include <iostream>

int main(int argc, char *[])
{
        int c = (argc == 1) ? 32 : 3;
        std::cout << "c=" << c << std::endl;
        std::cout << (1U << 32) << std::endl;
        std::cout << (1U << c) << std::endl;
}
>clang++ t.cpp && ./a.out
c=32
3381226160
1
>% clang++ t.cpp -O1 && ./a.out
c=32
3140643504
8
>% clang++ t.cpp -O2 && ./a.out
c=32
0
8
% g++ t.cpp -O3 && ./a.out
c=32
0
0
% g++ t.cpp  && ./a.out
c=32
0
1

見事にばらばら. n == 32のときは別処理するのが安全. あるいは mask(n + 1) for n = 0, ..., 31が欲しいなら最初から1bitシフトしておいて 

unsigned int mask2(unsigned int n)
{
    assert(0 <= n && n < 32);
    return (2U << n) - 1;
}

としてしまうのも手.

6月18日_

次世代暗号の解読で世界記録を達成
F397上の楕円曲線の離散対数問題が解けたらしい. まだまだ先だろうと思ってたけど凄いですね.
ただ今回の攻撃でペアリング暗号そのものが破れたわけじゃないので一応注意. (3p)6でp=97は実装評価でよく使われるものでは一番小さい(とは言え繰り返しになるけどp=97を解いたのはすごい功績). 安全性を担保するためにp=239やp=509とか, (2256)12とか, 最近はそれより大きいのも扱われる.
たとえるならRSAは最近は1024ビットよりも大きい2048ビットなどの長い鍵を使うべしとなっている. 768ビットの鍵長が解かれたので, それより小さい512ビットなんて使ってると痛い目に会うというようなもの.
ペアリング暗号はまだ歴史が浅いので安全性の見積もりはRSAよりは精度よくないだろう. こういう攻撃によって, 安全なパラメータをより正確に考察できるようになる.
もちろん, ある日ものすごい攻撃アルゴリズムの発見により, 完全に破られる可能性はゼロじゃないけどね.

6月16日_

x86/x64最適化勉強会4終了

今回も快適な会場のご提供とustreamをしてくださった@kzmogiさんと(株)AXELL社のご好意に感謝します. 本当にありがたいです.

  1. maropu : x86-optimized rank/select dictionary for bit sequences[ust-1]
  2. starpoz : Intel TSXについて[ust-2]
  3. Constellation : Context Threading JIT[ust-3]
  4. nothingscosmos : x86向けHotspotの気持ちになって考える、JITコンパイラの適応的最適化について[ust-4]
  5. umezawa_takeshi : 「AMDで使うと遅いんだけど」[ust-5]
  6. w_o : Bulldozer入門[ust-6]
  7. kzmogi : Intel OpenCL SDK
  8. herumi : AA-sort with SSE4.1[ust-8]

以下感想など. 今回も多岐に渡ってお話が聞けてとても楽しかった.

  1. キャッシュにヒットしやすいデータ構造でrankを数倍高速化. 中村実さんの256bitではpopcnt 4回の方がよいのではという意見も興味深い.
  2. トランザクション入門としてわかりやすい. ありがたいです.
  3. 着実に高速化されていっていつも感心してる. あのしゃべりの速さにも慣れてきたのだけど本当はもっと時間をとってゆっくり聞きたい.
  4. 文字列処理はC++よりJavaが速くなることも多いというのもうなずけそう.
  5. キャッシュ競合は忘れたころにぶつかるととても悩む.
  6. CPU性能を見るベンチマークコードの作り方がうまい.
  7. キャッチーな題材で笑った. GPGPUの本当のパワーは専門の人の意見が聞きたかったな.
  8. iwiwiさんとmaropuさんに後続の研究を教えていただいたのでちょっと見てみる.
6月15日_
・明日のx86/x64最適化勉強会4, いろいろおもしろそうな話題が並んでいて楽しみ.
・先月買ったマシン, なぜか時々勝手に電源が入ってるので気味が悪かったが単にwake on LANが有効になっていたからだった.
6月13日_
MySQLの認証がうっかり通ってしまうバグ, memcmpの戻り値をcharにcastしてるので例えば256だと0になって一致してしまうせいだったみたい. しかし修正箇所にtestを使うのか. てか相変わらずまだそんなマクロがmy_global.hなんてglobalなヘッダにあるのね.
6月11日_
風邪引いた. 鼻水が止まらない. 微熱もあるようだ. 最近眠くてどこでも寝てしまうせいだろうか.
6月10日_
ネジ穴がつぶれてどうしても抜けない. webでいろいろ探したらネジ穴に1, 2滴垂らして回すと開くという不思議な薬品があるらしい. 近くのホームセンターに行ってみると何種類か売られていた. ついでに先が切り込むような特殊な形になってるドライバーもあったので一緒に買ってみる. 合わせて1100円ぐらい.
家に帰って早速試してみる. 片方づつ試して余計にネジ穴がつぶれても嫌だったので両方同時に試す. あっさり開いた. これは凄い. ただ両方一緒だったのでどちらの効果が強かったのかは分からない. 薬の方は細かい砂利が入ってる感じだった.
6月4日_

kaorunさんに電化製品の消費電力を計れるワットチェッカーを貸していただいた. 早速パソコンやテレビに挟んでみる.

まずは先月買った超低電圧版のi3マシン. 33Wぐらい. CPUを使う計算させても40Wぐらいかなあ. 次, 6年前に購入したtype R(Pentium D)はなんと150～180W. うひゃー. もっと早く買い換えるべきだったか. 一緒に買ったDELLのモニタは60W弱. 結構食ってる.

あと自宅サーバは30～40W. i3マシンと同じぐらい. i3マシンに統合して引退させようかな. テレビは試聴時45Wぐらい, 消してしばらくは25W(録画時もそれぐらい?), 5分ぐらいしたら0W.

6月2日_
x86/x64最適化勉強会4を6/16(土)に開催します. 懇親会のatndも用意したのでよろしくお願いします.
昔
[Back]
御意見は光成滋生<herumi@nifty.com>までお願いします