Movatterモバイル変換

Optimal Ateペアリングの実装詳細2013/7/3サイボウズ・ラボ光成滋生

目次 Optimal Ateペアリングの定義(さらっと) 今回はペアリングの話ではなく最適化全般のトピックが主 x64 CPUの概略 実行時間の計測 整数加算、減算の実装 整数乗算の実装 Haswell向けの改良 その他のトピック/ 282

BN曲線 𝐹𝑝上で定義される埋め込み次数12の楕円曲線 𝐸: 𝑦2 = 𝑥3 + 𝑏, 𝑏 ∈ 𝐹𝑝 𝑝 ≔ 𝑝 𝑧 = 36𝑧4+ 36𝑧3+ 24𝑧2+ 6𝑧 + 1𝑧が64bitなら𝑝は256bitぐらいの素数 𝑟 ≔ 𝑟 𝑧 = 36𝑧4 + 36𝑧3 + 18𝑧2 + 6𝑧 + 1 𝑡 ≔ 𝑡 𝑧 = 6𝑧2 + 1 #𝐸(𝐹𝑝) = 𝑟/ 283

記号 𝜋: 𝑥, 𝑦 → 𝑥 𝑝, 𝑦 𝑝 Frobenius写像 BN曲線に対してはtrace(𝜋 𝑝) = 𝑡 𝑓𝑠,𝑄: 𝐸上の有理関数(𝑠は整数𝑄は𝐸上の点) div 𝑓𝑠,𝑄 = 𝑠 𝑄 − 𝑠 𝑄 − 𝑠 − 1 𝒪 を満たすもの𝑠 𝑄 は 𝑄 の形式的な𝑠倍、 𝑠 𝑄は𝑄の 𝑠 倍された点を意味する 𝑙 𝑄1,𝑄2 𝑄1と𝑄2を通る直線/ 284

Optimal Ateペアリング 𝐺1 = 𝐸 𝑟 ∩ ker 𝜋 𝑝 − 1 = 𝐸 𝐹𝑝 𝑟 𝐺2 = 𝐸 𝑟 ∩ ker 𝜋 𝑝 − 𝑝 ⊆ 𝐸 𝐹𝑝12𝑟 𝐺3 = 𝜇 𝑟 ⊂ 𝐹𝑝12 ∗ 𝑒 ∶ 𝐺2 × 𝐺1 ∋ 𝑄, 𝑃 ⟼ 𝑚 𝑄, 𝑃𝑝12−1𝑟 ∈ 𝐺3 𝑚 𝑄, 𝑃 ∶= 𝑓6𝑧+2,𝑄 𝑃 ∙ 𝑔 𝑄 𝑃 𝑔 𝑄(𝑃) ≔ 𝑙 6𝑧+2 𝑄,𝜋 𝑝 𝑄 𝑃 ∙ 𝑙 6𝑧+2 𝑄+𝜋 𝑝 𝑄 ,−𝜋 𝑝2 𝑄 (𝑃)/ 285

ペアリングのアルゴリズム 1） 6𝑧 + 2 𝑄 と 𝑓6𝑧+2,𝑄 𝑃 を算出(Millerループ) 2） 𝑚 𝑄, 𝑃 = 𝑓6𝑧+2,𝑄 𝑃 ∙ 𝑔 𝑄 𝑃 を算出 3）𝑝12−1𝑟乗する(最終巾)/ 286

拡大体上の演算における戦略 𝐹 𝑝2上の乗算 x=a+bu, y = c+du, u^2 = -1 xy = (ac – bd) + ((a+b)(c+d) – ac – bd)u 従来 ac, bd, (a+b)(c+d)はFp:mulを使う Pairing2010における主要アイデア Fp:mul = mul256 + mod512mul256 : 256ビット整数乗算mul256 64bit乗算命令は速い(3clk, latency, 1clk throughput)mod512 : Montgomeryリダクション mul256の結果に対する加減算ac, bd, (a+b)(c+d)を512bit整数のまま加減算mod512の回数が3回から2回になる 512bit加減算は増える / 287

Aranhaらによる改良 𝐹 𝑝6などの拡大体にも容易に適用可能 拡大体の係数もより小さいものに 𝑏 = 2, z = −(262 + 255 + 1) 𝐹 𝑝2 = 𝐹𝑝 U / U2 − 𝛽 , 𝛽 = −1 ∈ 𝐹𝑝 𝐹 𝑝6 = 𝐹 𝑝2 V / V3− 𝜉 , 𝜉 = 1 + U ∈ 𝐹 𝑝2 𝐹 𝑝12 = 𝐹 𝑝6 W / W2 − V , 𝛾 = 𝑉 ∈ 𝐹 𝑝6 実装 最新の実装は上記を踏襲し，細部を改良 https://github.com/herumi/ate-pairing/0.35msec@Haswell(i7-4700MQ 3.4GHz)/ 288

x64 CPU概略 15個の汎用64bitレジスタ rax, rbx, rcx, rdx, rsi, rdi, rbp, r8, r9, ..., r15 フラグレジスタ 演算結果に応じて変わる1bitの情報群CF : 加算時に結果が64bitを超えた、減算でマイナスになったZF : 結果が0になったSF : 結果の最上位ビットが1だった 呼び出し規約 関数の引数に対応するレジスタ名WindowsとLinuxで異なる Windows : rcx, rdx, r8, r9 Linux : rdi, rsi, rdx, rcx 関数の中で壊してよいものと元に戻す必要のあるものLinux : r12, ..., r15, rbx, rbp, Win : 加えてrsi, rdi/ 289

算術演算 加減算 add x, y // x ← x + y; sub x, y // x ← x – y; carryつき加減算 adc x, y // x ← x + y + CF; 繰り上がりを加味 sbb x, y // x ← x – y – CF; 繰り下がりを加味 乗算 64bit x 64bit → 128bit mul x // [rdx:rax] ← x * rax (rax, rdxレジスタ固定) 除算 128bit / 64bit = 64bit あまり 64bit div x // [rdx:rax] / x ; 商 : rax, あまり : rdx/ 2810

条件比較 演算結果に応じてフラグが変わる フラグに応じて条件分岐する こういうコードはこんな感じ jg (jmp if greater), jge(jmp if greater or equal)などなど/ 2811if (x >= y) {Aの作業} else {Bの作業}cmp x, y // x-yの計算結果をCFに反映(CF = x >= y ? 0 : 1)jnc LABEL_A // jmp to LABEL_A if no carryBの作業jmp NEXTLABEL_A:Aの作業NEXT:

アセンブラの種類 gas, NASM, MASMなど 静的なアセンブラ マクロや条件式などの文法はそれぞれ独自構文 inline assembler おもにgcc(64bit Visual Studioでは非サポート) コンパイラが多少最適化してくれることも 記述が難しい LLVM 抽象度の高いアセンブラ/JIT可能 carryの扱いが難しく今回の用途では性能を出しにくい Xbyak(拙作) 抽象度は低い(gasやNASMと同じ)/JIT可能 C++の文法でアセンブラをかける/ 2812

実行時間の測り方 Vtune(Intel), CodeAnalyst(AMD)など CodeAnalystは無料Intel CPUでも使える perfコマンド(Linux only) perf listで測定したいパラメータを表示instructionsbranch-missessなどCPUによって様々なものがある perf stat –e L1-icache-load-misses 実行コマンド/ 2813

rdtsc CPUがもつカウンタ (2.8GHzなら1/2.8 nsec単位で)一つずつ増える Turboboostは切った方が周波数が固定になってよい駄目なら重たい処理を先に実行させてトップスピードにさせる マルチプロセス向けにrdtscpというのもある Xbyakではrdtscの薄いラッパークラスClockを提供 clk.begin(), clk.end()で測定したい関数をはさむ 最後にclk.getClock() / clk.getCount()で平均値を取得/ 2814Xbyak::util::Clock clk;for (int i = 0;i < N; i++) {clk.begin();some_function();clk.end();}printf("%.2fclk¥n", clk.getClock() / double(clk.getCount()));

256bit加算 記法 xi, yi, ziなどは64bitレジスタを表す [x3:x2:x1:x0]で256bit整数を表す(x0が最下位の64bit) 256bit整数z[]に256bit整数x[]を足すコードは次の通り 注意 z[], x[]が256bitフルに入ってると結果が257bitになる 今回はpを254bitに選んだため0 <= x, z < pならあふれない他にも様々な箇所で桁あふれがおきないため処理の簡略化が可能そのためセキュリティレベルが128bitではなく127bit/ 2815// [z3:z2:z1:z0] += [x3:x2:x1:x0]add z0, x0adc z1, x1 // carryつきadc z2, x2 // carryつきadc z3, x3 // carryつき

256bit加算を関数にする 呼び出し規約にしたがってレジスタを使う なかなか面倒 XbyakのStackFrameを使うとある程度抽象化、自動化可能LLVMはより汎用的にできる/ 2816//addNC(uint64_t z[4],const uint64_t x[4],const uint64_t y[4]);void gen_AddNC() {Xbyak::util::StackFrame sf(this, 3); //引数3個の関数const Xbyak::Reg64& z = sf.p[0]; // 一つ目の引数const Xbyak::Reg64& x = sf.p[1]; // 二つ目の引数const Xbyak::Reg64& y = sf.p[2]; // 三つ目の引数mov(rax, ptr [x]);add(rax, ptr [y]);mov(ptr [z], rax);for (int i = 1; i < 3; i++) {mov(rax, ptr [x + i * 8]);adc(rax, ptr [y + i * 8]);mov(ptr [z + i * 8], rax);} }

gen_addNCの結果 WindowsとLinuxのそれぞれに応じたコード生成 StackFrameはスタックを確保したり一時変数を使ったり、rcx, rdxを特別扱いする指定もできる自動的にレジスタの退避復元をおこなう/ 2817// Windows(引数はrcx,rdx,r8の順)mov rax,ptr [rdx]add rax,ptr [r8]mov ptr [rcx],raxmov rax,ptr [rdx+8]adc rax,ptr [r8+8]mov ptr [rcx+8],raxmov rax,ptr [rdx+10h]adc rax,ptr [r8+10h]mov ptr [rcx+10h],raxret// Linux(引数はrdi,rsi,rdxの順)mov rax,ptr [rsi]add rax,ptr [rdx]mov ptr [rdi],raxmov rax,ptr [rsi+0x8]adc rax,ptr [rdx+0x8]mov ptr [rdi+0x8],raxmov rax,ptr [rsi+0x10]adc rax,ptr [rdx+0x10]mov ptr [rdi+0x10],raxret

Fp::addの実装 addNCした結果zがz>=pならばpを引く if (z >= p) z -= p; アセンブラレベルでの比較の方法 z=[z3:z2:z1:z0]とx=[x3:x2:x1:x0]はどちらが大きいか 1. 頭から比較する 分岐がきわめて多くなる連続する分岐命令は好まれない 2. 引いてから考える 分岐は1回/ 2818cmp z3, x3ja z_gt_x // z3 > x3jb otherwise // z3 < x3cmp z2, x2 // here z3 == x3ja z_gt_x // z2 > x2jb otherwise // z2 < x2...z_gt_x:...otherwise:tmp_z = z // zの値を退避(mov x 4)subNC z, p // 引いてみる(z -= p)jnc .next // z >= 0ならnextへz = tmp_z // zの値を復元(mov x 4).next:

分岐しないFp::addの実装 CPUは分岐予測をする 当たると大体1clk 外れると大体20clk 一般のデータでは偏りがあるので結構精度よく当たる が、今回はランダムなので的中率は50%→平均10clk 分岐予測を排除する 条件移動命令cmovXX2clk latency1clk thrgouthput addの二つの実装分岐あり1.39Mclk, 分岐なし1.35MclkもちろんCPUによって異なる可能性あり(sandy, ivyで効果あり)単純ベンチだと分岐予測があたって分岐あり版が速くみえるかも/ 2819mov ti, zi x 4subNC z, pcmovc zi, ti ; 引きすぎてたら戻す

Fp::subの実装 subNCした結果が負ならpを足す addと違ってsubNCした時のCFを見ればよいので比較不要 分岐を使った実装 cmovを使った実装 0クリア cmov + メモリロード 加算結構命令数が多いので分岐に対してそれほどメリットがない cmovを使わない実装 命令数は同じだがcmovよりは速い@sandy/ 2820// z -= xの直後jnc .nextz[] += p[].next:t[] = 0cmovc t[] p[] //t[] = CF ? p[0] : 0z[] += t[]sbb t, t // t = CF ? -1 : 0and t[], p[] // t = CF ? p : 0z[] += t[]

256ビット加減算の命令順序 メモリから読んで演算する二つの方式 方式A(メモリまとめ読み) 方式B(メモリと演算を交互に) 実験によるとどちらが速いかCPUにより異なるOpteron, i7は方式Aが速い Westmereは方式Bが速かった out of orderだから関係ないと思ったが1%弱違った 実行時のCPU判別によりいずれかを選択上記方式はコード全般にわたって適用される/ 2821z0 ← x[0]z1 ← x[1]z2 ← x[2]z3 ← x[3]z0 ← z0 + y[0]z1 ← z1 + y[1] with carryz2 ← z2 + y[2] with carryz3 ← z3 + y[3] with carryz0 ← x[0]z0 ← z0 + y[0]z1 ← x[1]z1 ← z1 + y[1] with carryz2 ← x[2]z2 ← z2 + y[2] with carryz3 ← x[3]z3 ← z3 + y[3] with carry

256ビットx256ビット乗算(1/2) 256ビット整数を64ビット整数4個の組で表現する 64ビット→320ビット乗算4回と320ビット加算3回 筆算方式でmulするごとにaddを行う繰り上がり加算が3回余計に増える/ 2822𝑥3 𝑥2 𝑥1 𝑥0𝑦𝑥0 𝑦𝑥1 𝑦++（繰り上がり）𝑥2 𝑦++（繰り上がり）・・・1. 𝑥0 𝑦を計算2. 𝑥1 𝑦を計算3. 𝑥0 𝑦 𝐿 + 𝑥1 𝑦 𝐻t0 = 𝑥1 𝑦 𝐻 + 𝑐𝑎𝑟𝑟𝑦4. 𝑥2 𝑦を計算5. 𝑥2 𝑦 𝐿 + 𝑡0𝑡1 = 𝑥2 𝑦 𝐻 + 𝑐𝑎𝑟𝑟𝑦6. ...

256ビットx256ビット乗算(2/2) 乗算4回してから加算すると余計な加算が不要 ただし乗算結果を保持するワークエリアが必要 mulがCFを変更するためadcと同時に使えない 15個のレジスタを使い回して一時メモリを使わずに実装/ 28231. [𝑥𝑖 𝑦]（𝑖 = 0 … 3）を計算2. それらをまとめて加算↓加算は4回𝑥3 𝑥2 𝑥1 𝑥0𝑦𝑥0 𝑦𝑥1 𝑦𝑥2 𝑦𝑥3 𝑦加算が終わるまでどこかに保持する必要がある

256ビットx256ビット乗算 for Haswell HaswellではCFを変更しないmulxが導入された 加算(add, adc)しつつ乗算を繰り返しおこなえる 必要なレジスタ数が減る 退避、復元のためのmov命令が減る Montgomery reductionにも適用可能 ペアリング全体で13%の高速化 1.33Mclkから1.17Mclkへ(@Core i7 4700MQ 2.4GHz)/ 2824mov(a, ptr [py]); | ↓mul(x); | mul(x);mov(t0, a); | mov(t3, a);mov(t1, d); | mov(a, x);mov(a, ptr [py + 8]); | mov(x, d);mul(x); | mul(qword [py + 24]);mov(t, a); | add(t1, t);mov(t2, d); | adc(t2, t3);mov(a, ptr [py + 16]);| adc(x, a);↓ | adc(d, 0);mov(d, x);mulx(t1, t0, ptr [py]);mulx(t2, a, ptr [py + 8]);add(t1, a);mulx(x, a, ptr [py + 16]);adc(t2, a);mulx(d, a, ptr [py + 24]);adc(x, a);adc(d, 0);

記述の簡便さのための手法 各種2項演算はsrc x 2 + dstのglobal関数を作る Fp::add(z, x, y); // z = x + yなど &z == &x == &yなどのときでも正しく動くように注意 演算子オーバーロード Fp operator+(const Fp&, const Fp&)などをFp::addを使って定義するz = x + y;などとかける。 Fp2, Fp6, Fp12などの拡大体でも同様に作る コピペばかりになって間違いやすい/ 2825

CRTPによる半自動的生成手法 add, subなどを使ってoperator+, operator-を定義するtemplateクラス Fp, Fp2などはadd, subさえつくればaddsubmulを継承することでoperator+が使えるようになる virtual継承ではないので呼び出し時のコストは(通常)ない/ 2826template<class T, class E = Empty<T> >struct addsubmul : E {template<class N>T& operator+=(const N& rhs) {T::add(static_cast<T&>(*this), static_cast<T&>(*this), rhs);return static_cast<T&>(*this); }...strut Fp : addsubmul<Fp>{static void add(Fp&, const Fp&, const Fp&);};

記法の簡便さと演算性能 z = x + y;とFp::add(z, x, y); 一般的に前者の方が書きやすく可読性も高い しかし隠れた一時変数の生成とコピーに注意する x + yの結果をtmpに保存してz = tmpを実行 方針 最初は数式を書きやすい前者で始める 動くことがわかったら一時変数や移動を減らすように後者に置き換える式Templateによる一時変数除去テクニックはあるが正直使いにくい、挙動を把握しにくいため勧めない/ 2827// Fp::add(z, x, y);lea r8,[y]lea rdx,[x]lea rcx,[z]call [mie::Fp::add]// z = x + y;lea r8,[rbp+7]lea rdx,[rbp-19h]lea rcx,[rbp-39h]call [mie::Fp::add]movaps xmm0,[rbp-39h] //データ移動movaps [rbp+37h],xmm0movaps xmm1,[rbp-29h]movaps [rbp+47h],xmm1

 Fp6などの演算は基礎体のmulやaddを呼び出す mulはレジスタをフルに使うため関数の中でレジスタの退避と復元をおこなっている 連続してmulを呼び出すならその間の復元と退避は除去可能 退避復元をしない専用関数を用意する呼び出し規約からの逸脱 コンパイラの関知しないところのため手作業が必要LLVMがこの分野で使えるならoptに任せることも可能になるか メリット 速度向上 デメリット デバッグが難しいかもしれないFp2_mul:call Fp_mulcall Fp_mulretFp_mul:レジスタの退避本体レジスタの復元レジスタの退避・復元の省略の一般論/ 2828Fp2_mul:call in_Fp_mulcall in_Fp_mulret// Cからは呼べないin_Fp_mul:本体

Movatterモバイル変換

Change Language

optimal Ate pairing

Embed presentation

Recommended

More Related Content

What's hot

Similar to optimal Ate pairing

More from MITSUNARI Shigeo

optimal Ate pairing