Raspberry Pi2
会社にあったのでちょっと試してみる. clangのバージョンが3.0だったのでjessieにupgradeして試す. と思ったが容量が無い. 16GBあるのに3GBぐらいしか見えてない. ディスクをフルに認識させるにはsudo raspi-configでExpand Filesystemらしい.
C++のコンパイルが遅い. Makefileは修正が必要だがコードはいじらなくてよかった. そうなるよう作ったつもりだったのが確認できてよかった. LLVMのおかげ.
CPUは900MHz quad-core ARM Cortex-A7だそうな. ペアリングは35msecだった. 遅い. 測定方法が間違ってるのか, コードが悪いのか.
DSLへの移行完了
ようやくC++版LLVM IRコード生成の移行が終わった. これからはこちらで最適化する. 多分大分作業がしやすくなる.
LLVM llcでPIC対応オプション
コードを書いていて急に共有ライブラリのbuildでエラーがでるようになった.
動く前と後を比べると今までinlineで呼ばれていた関数がcallで呼ばれるような出力になっていた. LLVM IRではinlineされるかされないかは最適化の状況次第. alwaysinlineというキーワードはあるけど今回はつけたくない.
-fPICはつけているからいいと勘違いしていたがcpp→sだけでなくll→sの問題. でllcで*.sを出力するときにPIC対応にしないといけないのだった. オプションは-relocation-model=pic .
こうすることで関数名に@PLTがつくようになる.
Windowsのselect
Linuxでselect(2)を使うときはfdが範囲外にあるとメモリアクセス違反を起こしてしまうかもしれないのでFD_SETSIZEより小さいことを確認しないといけない. でもWindowsではその引数は無視される. 逆に確認してはいけないのだった.
LLVM IR DSL
汎用的につくるといろいろ足りないところがあるが, とりあえず自分が今欲しい部分だけ追加していく.
Window 10 upgradeリトライ
2カ月ほど前にやって戻したけど, もう不具合は直ったかなとやってみる. 残念ながらまだリブートやサスペンドまわりがおかしい. 仕方がないのでまたWindows 7に戻す.
ISEC参加
今日から研究専門委員会というものに参加することになった.
DSLの続き
>mul(z, x, y)よりはz = mul(x, y)の方が見た目が近いので捨てがたい.
この方針で進めた.
%p = call i192 @makeNIST_P192()
%pe = zext i192 %p to i256
%c0 = sub i256 %z, %pe
%c1 = lshr i256 %c0, 192
%c2 = trunc i256 %c1 to i1
%c3 = select i1 %c2, i256 %z, i256 %c0
%c4 = trunc i256 %c3 to i192
store i192 %c4, i192* %out
ret void
が
Operand p = call(makeNIST_P192);
p = zext(p, 256);
Operand zp = sub(z, p);
Operand c = lshr(zp, 192);
c = trunc(c, 1);
z = select(c, z, zp);
z = trunc(z, 192);
store(z, out);
ret(Void);
と書けるようになって大分視認性も書きやすさもアップした. しかも値を返す形にしたので
Operand zp = sub(z, zext(p, 256)); Operand c = trunc(lshr(zp, 192), 1); store(trunc(select(c, z, zp), 192), out);
と入れ子にもできた. これが
%r5 = sub i256 %r2, %r4 %r6 = lshr i256 %r5, 192 %r7 = trunc i256 %r6 to i1 %r8 = select i1 %r7, i256 %r2, i256 %r5 %r9 = trunc i256 %r8 to i192 store i192 %r9, i192* %r3
と展開される. これは予想外の利点.
コード生成DSL続き
mul(z, x, y)よりはz = mul(x, y)の方が見た目が近いので捨てがたい. 次のようにやれば一応実現できる. mul(x, y)で文字列を出力するのではなく, 右辺を評価したものをEvalという型にして返す.
Operatorはoperator=(const Eval&)で自身のレジスタの情報を出してからEvalの結果を表示する.
Eval mul(const Operator& x, const Operator& y)
{
Eval e;
e.str = "mul(x, y)を文字列化したもの";
}
Operator::operator=(const Eval& e)
{
std::cout << "自身のレジスタ" << "=" << e.str << std::endl;
}
こんな感じ.
するとcallでEvalを返してもoperator=でアサインしない限り表示されない問題が起こる. こちらはEvalにbool used;を追加してoperator=が発動するとtrueにし, Evalのデストラクタでused = falseならstrを表示すればよい.
この程度なら難しくはないのでこちらを採用しよう. もっとややこしいことをしないといけなくなりそうだったらやめる.
コード生成のDSL
いい加減,
define void @mcl_fpDbl_mulPre$(bit)(i$(unit)* %pz, i$(unit)* %px, i$(unit)* %py) {
%y0 = load i$(unit)* %py
%xy0 = call i$(bu) @mulPv$(bit)x$(unit)(i$(unit) *%px, i$(unit) %y0)
%z0 = trunc i$(bu) %xy0 to i$(unit)
store i$(unit) %z0, i$(unit)* %pz
%t0 = lshr i$(bu) %xy0, $(unit)
こんな感じのコードを書くのはしんどくなってきたのでもうちょっと楽にできるようにC++で書き直す. どんなDSLにするか考えよう.
レジスタはサイズが分かっていれば明示する必要がないので
Operand x(Int, 64); Operand y(Int, 64); Operand z; z = add(x, y);
のように書きたい.
SSAであるLLVM IRは変数名を考えるのが大変なので
x = add(x, x);
とかけるようにもしたい. つまり出力は
%x1 = add i64 %x0, %x0
みたいになってほしい.
でもそもそもz = add(x, y)はちょっと考えないと実現が難しい. なぜならadd(x, y)をadd i64 %x0, %y0に変換するのは容易だがそれを出力するより前に左辺のレジスタを%z0を表示しないといけないから.
DSLの文法としてadd(z, x, y);を採用すればずっと簡単なのでそちらにするか.
ただcall(mulPv, x, y);は値を返すときと返さないときがある.
call(z, mulPv, x, y);とcall(mulPv, x, y);にするか. 悩ましい.
値返しの続き
場合によってはそちらのほうがよいというわけではないようだ. またnoalias(C99のstrictみたいなの)をつけた方がregister spillが発生してしまうこともある. なぜだろう.
それでもどうにかx64環境で256bit × 256bit → 512bit乗算でXbyakによる手書きアセンブラとほぼ同等なコードを生成できるようになった. mcl_fpDbl_mulPreがmulx16個, add7個, adc24個, mov16個, push/pop6個のコードになる.
汎用コードを書いてそれでうまく行く, というわけではなかった. まあそんなうまい話はなかなかない.
ちょっとその場凌ぎのフィルタコマンドでは辛くなってきたのでC++で書き直す.
『ラプラスの魔女』(東野圭吾)
1章最初のエピソードとタイトルでキーアイデアの見当がついた. で, それがたいしたひねりもなく使われずに終わった. 『ジェノサイド』(高野和明)もだけど, 「みえる」からといって「狙った通りにできる」わけではないと思うので, 落ち葉や風船のはなしを入れるなら「常人ばなれした器用さ」について何か言及が欲しい.
LLVM IRで値返しにするかポインタ経由で値を返すか
たとえば
uint256_t mul(uint192_t x, uint64_t y);
という関数を
void mul(uint256_t *z, uint192_t x, uint64_t y);
にしたほうがよいのかという話. いくつか実験したところ小さいときはどちらでも一緒に見えた. それで前者の方がstackに明示的に確保しなくてよいので便利だから使っていた. けれでも, サイズが大きいところいろいろ実験するとどうも後者の方が効率がよいコードが生成されることが多い感じ.
LLVMが頑張って不要なメモリ移動を極力減らすのだけど, xとyの乗算中にzを破壊してよいのかは判断できないので安全側に倒される. ここではrestrictを使えるわけでもないので差が出ることがあるようだ.
そのあたりに注意しながら書き直す.
DLLをPythonから使う
MS Office暗号化ツールのDLL化をしたのでPythonから使ってみる.
ctypesを使うと簡単にできる. Excelで適当にplain.xlsxを作り msoc.dllをダウンロードして
from ctypes import *
lib = cdll.LoadLibrary('msoc.dll')
outFile = c_char_p("enc.xlsx")
inFile = c_char_p("plain.xlsx")
ps = c_char_p("test")
lib.MSOC_encryptA(outFile, inFile, ps, None)
これだけでmsoc.dllが提供するMSOC_encryptAを呼び出して"test"というパスワードつきのenc.xlsxができる. Pythonの文字列をC文字列にするにはc_char_p, unicode文字列を渡したい場合はc_wchar_pを使う. 4番目のオプションは通常使わないのでNoneを渡す.
Windowsではmain関数に渡されたコマンドライン引数は(日本語に設定していれば)cp932(Shift JIS)に変換されている. したがってコマンドラインでShift JISで表現できない文字, たとえばハート♡を渡すと文字化けがおこる.
Cならint wmain(int argc, wchar_t **argv)という関数を使えばよいがPythonにはそういうのはなさそうだ. 仕方がないのでWin APIのGetCommandLineとCommandLineToArgvWを使って作る.
from ctypes import *
from ctypes.wintypes import LPWSTR
def getWargv():
cmd = cdll.kernel32.GetCommandLineW()
argc = c_int(0)
CommandLineToArgvW = windll.shell32.CommandLineToArgvW
CommandLineToArgvW.restype = POINTER(LPWSTR)
argv = CommandLineToArgvW(cmd, byref(argc))
ret = []
for i in xrange(argc.value):
ret.append(argv[i])
return ret
GetCommandLineWはkernel32.dllにあって引数はないのでそのまま呼び出せる. CommandLineToArgvWは返り値の型がLPWSTR*なのでrestypeでPOINTER(LPWSTR)と指定してから呼び出す. int変数argcのポインタを渡すにはbyrefを使う.
改良続き
64bit環境で256bit x 64bit乗算も改良された. 今まで
mulxq (%rsi), %r8, %r9
mulxq 8(%rsi), %rcx, %rax
mulxq 16(%rsi), %r11, %r10
mulxq 24(%rsi), %rdx, %rsi
addq %r9, %rcx
adcq $0, %rax
adcq $0, %rdx
adcq $0, %rsi
addq %r11, %rax
adcq %r10, %rdx
adcq $0, %rsi
movq %r8, (%rdi)
movq %rcx, 8(%rdi)
movq %rax, 16(%rdi)
movq %rdx, 24(%rdi)
movq %rsi, 32(%rdi)
retq
だったのが
mulx r8, rax, qword ptr [rsi]
mov qword ptr [rdi], rax
mulx rcx, r9, qword ptr [rsi + 8]
mov qword ptr [rdi + 8], r9
mulx rax, r10, qword ptr [rsi + 16]
mov qword ptr [rdi + 16], r10
mulx rsi, rdx, qword ptr [rsi + 24]
mov qword ptr [rdi + 24], rdx
add r8, r9
adc rcx, r10
adc rax, rdx
adc rsi, 0
mov qword ptr [rdi + 24], rax
mov qword ptr [rdi + 16], rcx
mov qword ptr [rdi + 8], r8
mov qword ptr [rdi + 32], rsi
ret
と出力されるようになった. LLVM IRの型をいちいち書くのが面倒. Xbyakライクなラッパを用意すると便利かもしれない.
LLVM最適化
6日のレジスタspillの件, 128bit x 32bitの乗算ルーチンで1行1行何が起こってるかを観察する. これぐらいのサイズだとspillは起こっていない. しかし(***)で囲った部分に無駄がある. なぜこんなコードが生成されてしまうのだろう. LLVMのcarryの扱いはどうもよくわからない.
/*
[z4:z3:z2:z1:z0] = [x3:x2:x1:x0] * y
esp 4 20 16 12 8 24
x3 x2 x1 x0
y
-------------------
eax:ecx
edi:esi
ebx:ebp
edx:eax
*/
mul128x32
push ebp
push ebx
push edi
push esi
mov edx, dword ptr [esp + 40] ; y
mulx eax, ecx, dword ptr [esp + 24] ; [eax:ecx] = x0 * y
mulx edi, esi, dword ptr [esp + 28] ; [edi:esi] = x1 * y
mulx ebx, ebp, dword ptr [esp + 32] ; [ebx:ebp] = x2 * y
add esi, eax
mulx edx, eax, dword ptr [esp + 36] ; [edx:eax] = x3 * y
(***)-----------------------
adc edi, 0 ; 不要
adc eax, 0 ; 不要
adc edx, 0 ; 不要
add edi, ebp ; adc edi, ebpにする
(***)-----------------------
adc eax, ebx
mov ebx, dword ptr [esp + 20] ; ptr to z
mov dword ptr [ebx], ecx
mov dword ptr [ebx + 4], esi
mov dword ptr [ebx + 8], edi
mov dword ptr [ebx + 12], eax
adc edx, 0
mov dword ptr [ebx + 16], edx
pop esi
pop edi
pop ebx
pop ebp
ret
いろいろがんばってLLVM IRを改良したら望みのコードがでるようになった.
mul128x32:
push ebp
push ebx
push edi
push esi
mov eax, dword ptr [esp + 20] ; &z
mov ecx, dword ptr [esp + 24] ; &x
mov edx, dword ptr [esp + 28] ; y
mulx edi, esi, dword ptr [ecx] ; [H0:L0] = [edi:esi] = x[0] * y
mov dword ptr [eax], esi ; z[0] = esi
mulx ebx, esi, dword ptr [ecx + 4] ; [H1:L1] = [ebx:esi] = x[1] * y
mov dword ptr [eax + 4], esi ; z[1] = L1
add edi, esi ; H0 += L1
mulx ebp, esi, dword ptr [ecx + 8] ; [H2:L2] = [ebp:esi] = x[2] * y
mov dword ptr [eax + 8], esi ; z[2] = L2
adc ebx, esi ; H1 += L2
mulx edx, ecx, dword ptr [ecx + 12] ; [H3:L3] = [edx:ecx] = x[3] * y
mov dword ptr [eax + 12], ecx ; z[3] = L3
adc ebp, ecx ; H2 += L3
mov dword ptr [eax + 12], ebp ; z[3] = H2
mov dword ptr [eax + 8], ebx ; z[2] = H1
mov dword ptr [eax + 4], edi ; z[1] = H0
adc edx, 0 ; H3
mov dword ptr [eax + 16], edx ; z[4] = H3
pop esi
pop edi
pop ebx
pop ebp
ret
次はこのLLVM IRを生成するようにコードを修正しよう.
MS Office暗号化ツールのDLL化
Excelで食ってるとおぼしき方から欲しいという要望があったので作った(msoffice).
const wchar_t *inFile, *outFile, *pass; MSOC_encrypt(outFile, inFile, pass, NULL); MSOC_decrypt(outFile, inFile, pass, NULL);
で暗号化や復号ができる. ファイル名やパスワードの型がconst char*ではなくconst wchar_t*なのに注意. 久しぶりにCでサンプルコードを書いたので辛かった.
VCの/MTと/MD
今までライブラリを/MDでbuildすることが多かったのだけどもろもろの事情で/MTに変更することにした. なんか/MDをつかても_MDが定義されないのだけどそういうものだっけ(MSDN).
#include <stdio.h>
int main()
{
#ifdef _MT
puts("MT");
#endif
#ifdef _MD
puts("MD");
#endif
}
>cl t.cpp /MD && t
>MT
『新世界』(柳広司)
原爆を開発したロスアラモスが舞台. どこまでが史実でどこが創作なのかわかりにくいが力作. 天才科学者ばかりの中で狂ってるのはいったいだれなのか. 黙示録の章の描写が圧巻.
『違和感から始まる社会学』(好井裕明)
社会学は自分の身の回りのことの常識を観察したり疑ってみたりすることから始まると紹介する. スマホを見てる学生の例がよく取り上げられるが, 何をいいたいのかよくわからなかった.
ただわりと新しめの題材も多くてよみやすい. 大学のときに受講した社会学の講義は『命題コレクション』を順番に読んでいただけ(? すまん. あまり覚えてない)でそんなに面白いと思わなかったが, この人の講義なら面白そう. 巻末に参考文献もいろいろ紹介されてるので見てみよう.
「素数」の会だった. お二人で交代で合計約6時間. お疲れさまでした. グリーン・タオの定理の証明の流れが聞けただけでもありがたい.
ただ時間が足りないのだから, クイズなどの演出は省いてもらって本題を少しでも長くできるようにしてもらったほうがありがたかったかな.
2歳ぐらいのお子さんを連れてきている方もいらして驚いた. 少し前はうちの会社に子ども連れで就職活動に来ていた方もいたようで, だんだんそういうのが普通になっていくのかもしれない.
mclの32bit対応
mclは最初から32bitでも動かすことを念頭に入れていたので1箇所修正するだけで動いた. つまり多分ARMやその他の環境でも動くようになったはず. ただx86で測定するとそんなに速くない. LLVMの出力するコードを見るとregister spillがめちゃくちゃ多くてメモリ移動が大量発生していた.
うーん, これはうまく書かないとよきにしてくれないのか. 高速なコードを出力させるにはちょっと試行錯誤が必要そうだ.
AVX-512の検出
マニュアルを見ながらテストする. AVX512はたくさんのグループに分かれてるようだ.
Fが基本で512bit対応しているか. CDはConflict Detectionというものらしい. あと逆数用のERIというのもある. 将来はIFMA, VBMIなどもでてくるようだ(よく見てない). それからbyte, wordに対応しているかのBW, 32/64bit対応かのDQ, 128/256に対応しているかのVLというフラグもある. 128bitのCD関係の命令を使いたかったらFとCDとVLのフラグが立ってるかを確認しないといけないようだ. かなり複雑.
フラグの検出にはAVX以降はOSXSAVEが使えるかを確認してからxgetbvで該当ビットが立っているかを見る. OKなら更にAVX512Fが対応してるかを確認してから個別のフラグを見る(xbyak_util.h).
if (type_ & tOSXSAVE) {
// check XFEATURE_ENABLED_MASK[2:1] = '11b'
uint64 bv = getXfeature();
if ((bv & 6) == 6) {
if (data[2] & (1U << 28)) type_ |= tAVX;
if (data[2] & (1U << 12)) type_ |= tFMA;
if (((bv >> 5) & 7) == 7) {
getCpuid(7, data);
if (data[1] & (1U << 16)) type_ |= tAVX512F;
if (type_ & tAVX512F) {
getCpuidEx(7, 0, data);
if (data[1] & (1U << 17)) type_ |= tAVX512DQ;
if (data[1] & (1U << 21)) type_ |= tAVX512IFMA;
if (data[1] & (1U << 26)) type_ |= tAVX512PF;
if (data[1] & (1U << 27)) type_ |= tAVX512ER;
if (data[1] & (1U << 28)) type_ |= tAVX512CD;
if (data[1] & (1U << 30)) type_ |= tAVX512BW;
if (data[1] & (1U << 31)) type_ |= tAVX512VL;
if (data[2] & (1U << 1)) type_ |= tAVX512VBMI;
}
}
}
}
Intelのエミュレータsdeで試す.
// Skylakeのクライアント版 xbyak>..\sde\sde -skl -- sample\test_util.exe 64bit vendor intel mmx mmx2 cmov sse sse2 sse3 ssse3 sse41 sse42 popcnt aesni rdtscp osxsave(xgetvb) pclmulqdq avx fma avx2 bmi1 bmi2 lzcnt enh_rep rdrand adx rdseed smap f16c movbe popcnt ok family=6, model=E, stepping=0, extFamily=0, extModel=4 display:family=6, model=4E // Skylakeのサーバー版 xbyak>..\sde\sde -skx -- sample\test_util.exe 64bit vendor intel mmx mmx2 cmov sse sse2 sse3 ssse3 sse41 sse42 popcnt aesni rdtscp osxsave(xgetvb) pclmulqdq avx fma avx2 bmi1 bmi2 lzcnt enh_rep rdrand adx rdseed smap f16c movbe avx512f avx512dq avx512cd avx512bw avx512vl popcnt ok family=6, model=5, stepping=0, extFamily=0, extModel=5 display:family=6, model=55
それっぽく判定できてるようだ.
しかし命令が
VADDPD zmm1 {k1}{z}, zmm2, zmm3/m512/m64bcst{er}
こんな感じなのか. 一体どうやってC++で表現しよう.
AVX-512
勉強を始めようとマニュアルIntel Architecture instruction set extensions programming referenceをみたら1200ページ近くもあってびびる. まずはFoundationからかな. しかしAVXのときに一新したフォーマットVEX(vector extensions)じゃ対応しきれなくてEVEX(enhanced VEX)を導入するってどんなけその場凌ぎの拡張...
mcl
流石に簡単には速くできなくなってきた.2Mclk切るぐらいはすんなりいけるかなと思ったけど. 先に32bit対応(とARM)をするか.
プロファイル
mclは2.4Mclkでate-pairingとの差はあと2倍. 内訳はMiller loopが1.1Mclkでfinal expが1.3Mclk. ate-pairingはそれぞれ700Kclkと460Kclk.
final expの方がまだ無駄が多いな.
mclライブラリ続き
>ちなみにTEPLAは7.4Mclkだったのでまずはそれを越えないと.
5.5Mclkになって越えた. あと4～5倍.
もうちょっといじって3Mclk. あと2～3倍.
mclライブラリ
ようやく新しく作り直したナイーブなoptimal ate pairingが動いた. 一山越えた. ただがちがちに最適化したやつ(1.1Mclk)より24倍遅い(約27Mclk). 1.5倍程度の遅さ以内にはしたいとするとここから約20倍速くしないといけない.
ちなみにTEPLAは7.4Mclkだったのでまずはそれを越えないと.