新宿は3年半振りぐらい. ビルがでかい.
最近将来プログラム書いてて食ってけるのか心配になることも多いので, サイボウズの青野さんの講演「日本のITエンジニアは幸せになれるのか？」と首藤さんの対談に興味があって行ってみた.
やっぱり得意な分野を突き進めるのは大事. ただその技術を使う側の人のことを考えるのも同じぐらい大事とのこと.
基礎技術ならその重要性をきちんを説明することかな. まあそう書いてしまうと当たり前みたいだけど.
午後の発表ではみんな型推論好きだなーと思いつつ聞いていた.
フリーディスカッションではこの事業をどう存続していくかという話が出た. 一人一人が周りに地道に知らせていくという意見もあったが, やはりお上(お金を出してくれるところ)に成果をきちんとアピールすること(個別にではなく母体として --- それはたとえばとりあえずは各自の成果やblogへの"公式"ポータルを作るでもいいだろう ---)が必要じゃなかろうか. というかそれはまさに青野さんが仰ってたことに当てはまるのではなかろうか. と, 部外者が勝手に思ってみた.
交流会では名刺を持ってる中学生がいた. 私が初めて持ったのは彼の2倍の年だ. うひゃー.
○ src/zm1.cppのmakeMul4x1()
256bit整数pyと64bit整数xを掛けて320bit整数にする. 筆算方式でやる. ただ掛ける度に加算してるとトータルの加算回数が6回になってしまう. py[i] * xを全て求めてから足すと4回ですむ. その代わり今度は途中のレジスタがたくさん必要になるのでうまくレジスタを入れ換えてやりくりする. 後者の場合はmul + adcのレイテンシも多少隠蔽できるかもしれない.
const Reg64& a = rax; const Reg64& d = rdx; mov(a, ptr [py]); mul(x); // [d:a]← a * x mov(t0, a); mov(t1, d); mov(a, ptr [py + 8]); mul(x); mov(t, a); mov(t2, d); mov(a, ptr [py + 8 * 2]); mul(x); xchg(a, x); // swap(a, x) mov(t3, d); mul(qword [py + 8 * 3]); add(t1, t); adc(t2, x); adc(t3, a); adc(d, 0);
○ src/zm1.cppのvoid makeMontgomery1()
Montgomery乗算の1ライン分を行う. 擬似コードは次の通り: ppはpの逆数で初期化時に一度だけ計算しておく.
c = [c4:c3:c2:c1:c0] c += x[3..0] * y q = uint64_t(c0 * pp) c = (c + q * p) >> 64
二段目以降のみ前回の分のを足す必要があるので切り替える.
if (isFirst) {
makeMul4x1(c3, c2, c1, c0, px, y, c4);
mov(c4, d);
// [c4:c3:c2:c1:c0] = px[3..0] * y
} else {
makeMul4x1(t3, t2, t1, t0, px, y, t4);
// [d:t3:t2:t1:t0] = px[3..0] * y
add_rr(c3, c2, c1, c0, t3, t2, t1, t0);
adc(c4, d);
}
その続き
mov(rax, pp_); mul(c0); // q = a mov(y, a); makeMul4x1(t3, t2, t1, t0, p, y, t4); add(c0, t0); // A mov(c0, 0); // *** adc(c1, t1); // B ...
***は何故xor(c0, c0)ではなくmov(c0, 0)なのか. x64では前者の命令長は3byte, 後者は7byteと大きいのに. 答えはxorはcarryをクリアするけどmovはしないから. Aで作ったcarryをBで使いたいのでこのタイミングでクリアするためにxorは使えない. c0の0クリアをadd()より前でできるならxorでよいが無かったのでこうした.
ちなみにnasmでmov rbx, 0とするとなんと10byteにエンコードされる. mov rbx, dword 0とすると7byteに減る. Xbyakではmov(rbx, 0);で後者のエンコードにするようにした.
○ src/zm1.cppのvoid makeNegEntry()
256bit整数xに対して-xを得る. xが0のときは0, それ以外はp - xを求める. この場合, x == 0である確率は極めて小さいと考えられるので分岐予測がヒットするはず. だから条件分岐を使うのがよい.
inLocalLabel();
load_rm(z3, z2, z1, z0, x);
mov(rax, z0);
or(rax, z1);
or(rax, z2);
or(rax, z3);
jz(".skip");
mov(rax, (uint64_t)&s_pTbl[1]);
load_rm(z3, z2, z1, z0, rax);
sub_rm(z3, z2, z1, z0, x);
L(".skip");
store_mr(z, z3, z2, z1, z0);
ret();
outLocalLabel();
inLocalLabe(), outLocalLabel()はXbyakの独自命令関数. この間にあるピリオドで始まるラベルはこの中でのみ有効となる. makeNegEntry()を複数回呼び出した場合, 複数の場所にある".skip"が混在しないようにするための仕様.
256bit整数xが0かどうかはx[3], x[2], x[1], x[0]の全てorをとってから判断することで分岐を減らす. 常套手段やね. 64bitのアドレッシングでは64bit即値は扱えないので一度raxにポインタを代入して使っている.
○ src/zm1.cppのvoid mul4x4()
256bit整数px[]と256bit整数py[]の積, 512bit整数pz[]を求める. 流石に16個のレジスタでもやりくりは苦しい. それでも同じメモリを複数回readすることはあってもtemporaryメモリは使わないようにした. そのおかげもあってかなり速い. mul命令は16回.
まず128bit([a:b])×128bit([c:d])→256bitではKaratsuba法によるmul3回よりも筆算方式mul4回が速い. なぜならad + bc = (a + b)(c + d) - ac - bdとしてmulを減らす場合, a + b, c + dが途中に必要. これらの値は65(= 64 + 1)bitなので(a + b)(c + d)は65bit×65bit→130bitの計算. その計算に"mul"は1回でよいが, carryをともなる計算が何回か必要になりそのコストはmul1回分よりも大きい.
次の可能性として128bit×128bit→256bitを筆算方式で行い, それを用いてKaratsuba法で256bit×256bit→512bitを求める方法がある. が, やっぱりそれでもコストは大きかった. 結局現時点でのCPUでは16回のmulを使う筆算方式が最速だった.
○ src/zm1.cppのvoid mul5(), dblMul5()
256bit×5, 512bit×5を求める. 前者の場合, (x << 2) + xとして求めたが後者ではmul命令を使った. シフト + 加算で5倍を実装するとtemporary変数が必要になり512bitではそれらをレジスタにおく場所が無い. また速度的にもあまり変わらないようであった. むしろmulの方がいいかもしれない. いずれにせよこれらの関数は全体の中でたいした時間を占めてないのであまり気にしない.
○ src/zm1.cppのvoid modp()
小さい値x(x < 14p)に関するx mod pを求める. 割とよく出てくるのでなんとか速くできないかと考えていたら, 今回選択したp特有の性質だけれども, xが何倍のpであるかをx[4] >> 61で概算できることに気がついた. でその概算値を予めテーブルtblに入れておき引き算する. 概算値の誤差は最大1なので大きめに引いておいて, もし負だったら引き過ぎということでpを足す. その部分のコード:
// x -= tbl[x>>61]; // if (x < 0) x += p; sub_rm(x3, x2, x1, x0, a + x4); // [x3:x2:x1:x0] -= tbl[x>>61] sbb(x4, x4); // x4 = x < 0 ? -1 : 0 load_rm(t3, t2, t1, t0, a + 32); // [t3:t2:t1:t0] ← p and(t0, x4); // x > 0ならt0 ← 0 and(t1, x4); // x > 0ならt0 ← 0 and(t2, x4); // x > 0ならt0 ← 0 and(t3, x4); // つまり[t3:t2:t1:t0] = x < 0 ? p : 0 add_rr(x3, x2, x1, x0, t3, t2, t1, t0); // x[] += t[]
ややトリッキー. sub_rmの時点でcarryが立っていたら, 引きすぎたということでsbb(x4, x4)により0xffffffff, そうでなければ0というマスクを生成する. load_rm()で値を読んで, andをすればcarryが立っているなら値が入って, なければ0になる.
もちろんt0, ..., t3を事前に0初期化しておきcmovcを使ってcarryがあるときのみメモリから値を読んでもいい. マイクロベンチマークだとandの方が気持ち速かった.
○ FizedZpCode()::init()
コード生成初期化部分. 関数ごとにクラスを作るとそれぞれごとにページサイズ分(4096byte)確保されて無駄が多い. だから予め12KBぐらいと多めに取っておいて, 現在のコード位置を取得して関数を生成し, align(16)を呼び出して次の関数を生成するというサイクルを繰り返す. そうすれば関数を詰めて生成できる.
mie::ZmZ_1::add = (void (*)(ZmZ_1& out, const ZmZ_1& x, const ZmZ_1& y))getCurr(); makeAddEntry(); align(16); mie::ZmZ_1::sub = (void (*)(ZmZ_1& out, const ZmZ_1& x, const ZmZ_1& y))getCurr(); makeSubEntry(); ...
○ src/zm1.cppのmakeAddEntry()
void makeAddEntry()
{
#ifdef XBYAK64_WIN
mov(ptr [rsp + 8 * 1], rbx);
mov(ptr [rsp + 8 * 2], rsi);
mov(ptr [rsp + 8 * 3], rdi);
const Reg64& z = rcx;
const Reg64& x = rdx;
const Reg64& y = r8;
const Reg64& z2 = r9;
const Reg64& z1 = r10;
const Reg64& z0 = r11;
const Reg64& t0 = rbx;
const Reg64& t1 = rsi;
const Reg64& t2 = rdi;
#else
sub(rsp, 8 * 1);
mov(ptr [rsp], r12);
const Reg64& z = rdi;
const Reg64& x = rsi;
const Reg64& y = rdx;
const Reg64& z2 = rcx;
const Reg64& z1 = r8;
const Reg64& z0 = r9;
const Reg64& t0 = r10;
const Reg64& t1 = r11;
const Reg64& t2 = r12;
#endif
...
これに限らず関数生成コード全般に, win64とlinux64を吸収するための#ifdef XBYAK64_WIN ... #else ... #endifがある. どうみてもかっこ悪い. jitasmみたいに抽象化しようと何度も仕様を検討した. ただjitasmはebpを使ったスタックフレームを作ってるけど私はそういうことはしたくないし, レジスタの使い回しのところでmul命令を使うときと使わないときでrdxを退避するのかどうかraxはつぶしていいのかとか, の指定は必須になり使う側もややこしい(覚えることが増える). 結局諦めてベタベタ書くことにした.
上のコードについてコメント. Win64でrspを動かさずにいきなりrbx, rsi, rdiを保存してるが, これはWindowsでのスタックに書いたように呼び出し元で確保された領域なので大丈夫. Linuxもそういう領域(red zone)はあるみたいだけど, -mno-red-zoneなんてオプションもあるし(挙動未調査)安全側に倒した.
今回はSSEを使わない実装なのでスタック退避の代わりにxmm0とかに保存してるときもある. まあ気持ちの問題.
#そのあたり, いろいろ試したかったというのもラップしなかった理由の一つ.
const Reg64&参照でレジスタの別名をつける. これによりこれ以降z, x, yとみやすいコードで記述できる.
void load_rm(const Reg64& z3, const Reg64& z2, const Reg64& z1, const Reg64& z0,
const Reg32e& m)
{
mov(z0, ptr [m + 8 * 0]);
mov(z1, ptr [m + 8 * 1]);
mov(z2, ptr [m + 8 * 2]);
mov(z3, ptr [m + 8 * 3]);
}
mで指定されたアドレスから256bitデータを[z3:z2:z1:z0]のレジスタに読み込む. const Reg32e&はアドレス生成を保持する特別なクラス. このクラスで受けると呼び出すときにレジスタ指定だけでなくて
load_rm(x, z2, z1, z0, x); load_rm(t3, t2, t1, t0, px + 8 * 4); load_rm(z3, z2, z1, z0, rbp + rax * 4 + 16);
みたいに書けて便利. add_rm(), sub_rr()なども同様に用意しておく.
さて, makeAddEntry()で作る関数の話. 256bit整数xと256bit整数yを足して256bitの固定素数pで割る. 今回p<N/2, N=2256なのでx + yをした時点で256bitを超えることは無い. これにより最後のcarryを保存しておく必要がないので少し簡単になる.
z = x + y; if (z >= p) z -= p;
if()をどう処理するかで悩む. zがpより大きいかはアセンブリレベルではz[3]のほうから順次比較するか, 引ききってみて最後にcarryがたったかのどちらかで判断することになる. 前者の場合, 実際にzがpより大きいと判断した時点で引き算をすることになる.
試した感じでは, cmovが速いAMD系は後者, jmpが遅いCore2Duoも後者, i7は微妙に前者なときもある?って感じ.
将来Intel系ももっとcmovが速くなってくだろうことを期待して後者のロジックで実装した(まあjmpは汚いというイメージもある).
もちろんJITなのでcpu判別して切り替えてもいい. 今回は有為な差で無かったのでそうしなかった.
擬似コードはこう. 演算は全て256bit整数:
x += y // x + yを計算してxに代入
tmp ← x // xの値をtmpに保存
x -= p // x - pを計算してxに代入
if (carry) { // carryがたつ = x < pだったのでxをpを引く前の値に戻す
x ← tmp
}
最後のifのところでcmovを使う.
sub_rm(x, z2, z1, z0, rax); // [x:z2:z1:z0] -= ptr rax[3:2:1:0] cmovc(z0, t0); // carryがあれば[x:z2:z1:z0] ← [y:t2:t1:t0] cmovc(z1, t1); cmovc(z2, t2); cmovc(x, y);
; nasm call $ + 5 E800000000 ; Xbyak call (getCode() + 5 ); E800000000 ; yasm call $ + 5 E8FBFFFFFF$やgetCode()は現在値のアドレスを取得する命令だけど, それをcallの命令を構築する前に確定するのか, 後に確定するのかという話. nasmとyasmは殆ど互換性があるだけに$を含んでるとはまりそうだなあ.
struct uint256_t {
uint64_t data_[4];
static void (*add)(uint256_t&, const uint256_t&, const uint256_t&);
};
としてつくった場合に(昨日はaddをXbyakで生成したけど本質じゃないので略), やっぱりa += bとかしたいので次のようなテクニックを使う.
template<class T>
struct Addable {
T& operator+=(const T& rhs) {
T::add(static_cast<T&>(*this), static_cast<T&>(*this), rhs);
return static_cast<T&>(*this);
}
friend T operator+(const T& a, const T& b) {
T c; T::add(c, a, b); return c;
}
};
こういうtemplateを用意しておけば, CRTPを使って
struct uint256_t : public Addable<uint256_t> {...};
とすることでadd()を持つクラスから余計なコストを払わずに自動的にa += bを持つクラスを生成できる(これもメタプログラミングと言えるかな).
struct uint256_t : public Addable<uint256_t>, Subable<uint256_t> {...};
とすると何故か遅くなったり, 場合によっては動かなくなったりすることがある. それは何故か?
というのがクイズでした.
// 32bitVCでの結果
struct uint256_t {};の場合sizeoは32
struct uint256_t : Addable<A> {}; の場合sizeofは32
struct uint256_t : Addable<A>, Subable<A> {}; の場合sizeofは40
これはC++の規格上増えても増えなくてもよいので仕方がないようだ. でも余計に増えた8byteのせいでメモリレイアウトがくるって&uint256_t == uint256_t.data_でなくなってしまうかもしれない.
単一継承なら普通問題ないのでこれは次のようなEmptyクラスを用意して回避するのが常套手段.
template<class T> struct Empty {};
template<class T, class E = Empty<T>>struct Addable : E {上と同じ};
template<class T, class E = Empty<T>>struct Subable : E {...};
struct uint256_t : public Addable<uint256_t,Subable<uint256_t>> {...};
こうするとsizeofは32のまま.
まあうっかり二つで継承してしまっても動くことも多いので気づかないかもしれないのだけど,
asmな関数と連係するときはメモリレイアウトが狂うと終わりなので大事なんだということが言いたかった.
詳しくはboost::operatorsの注意書き参照.