暗号計算機屋のブログ

なにか思いついたことを不定期に更新。

わかりやすいICF3-Fのモンゴメリ乗算器の説明

はじめに

米アマゾンとイーサリアム財団らがFPGAコンテストをするそうです。簡単に言うとA×A mod N (A,N: 1024bit)をできる限り速く演算するというコンテスト。コンテストのサイトで紹介されているサバンチ大学の論文「Low-Latency Modular Multiplication Algorithm - Erdinc Ozturk」にはXilinxのFPGAで剰余乗算器、A×B mod N(A,B,N: 512bit)を実装した結果がある。巨大な剰余乗算器で、ICF3-Fの剰余乗算器(=モンゴメリ乗算器)よりも、かなり高速な演算器です。コンテストで1024bitのものが実装できれば、ですが。コンテストで1024bitの性能の結果が出たようです。詳しくはコンテストのサイトほうへ。

ご注意！　ICF3-Fはオープンソースではありません。

目的

サバンチ大学の論文には128bit、256bit、512bitの結果があるので、そこから1024bitのDSPの数を予想すると4225個。一方、ICF3-Fは44個×2です。 サバンチ大学の論文の方法で1024bitのものが実装されてもICF3-Fは、RSAのSSLアクセラレータの用途など、実際に利用される用途では負けていないということを言いたいのですが、 あまり細かいことを書きたくなかったので、わかりやすいICF3-Fのモンゴメリ乗算器の説明になりました。ICF3-Fとの違いで大きなものはNの入れ替え。 ICF3-Fは演算毎にNを入れ替えられるが、Nの入れ替えに時間がかかると思われます。つまりサバンチ大学の論文では複数のSSL証明書を運用したい場合、運用が難しい。周波数も考慮しないといけないのですがICF3-Fの約100倍のDSPを使うので100倍以上の性能が出ないとコストパフォーマンスが悪い。100倍以上としたのは、ICF3-FはリダクションもDSPで演算していますが、サバンチ大学の論文では巨大なLUTとBRAMを必要としている点があります。サバンチ大学の論文では小型のSSLアクセラレータを作れないであろうとか。

ICF3-Fのモンゴメリ乗算器(1024bit版 DSP 44個)

f:id:icf:20190812212327p:plain
ICF3-Fのモンゴメリ乗算器
1ループ=3サイクルで17bit分のモンゴメリ乗算。1024×1024bitをするには61ループ(183サイクル)。DSP内のレジスタに入力データを設定、ループ終了後、DSP内のレジスタに答えが格納されている。 uのブロッドキャスト以外は、隣接するDSPしかデータ転送がなくデータのローカリティが高いアーキテクチャ。k番目のDSPは、k-1番目のDSPとk+1番目のDSPと通信するのみで配線実装しやすい。 3サイクル中、乗算は2回行われる。乗算のうち、ひとつはuの値がないと演算できないためDSPが2倍あっても、あまり高速化できない。1つのDSPでuの値が不要な乗算を先に行い、その後、uの値を使う乗算を行う。DSPは乗数を17bitづつシフトする用途にも使われている。サバンチ大学の論文では乗算にRNSを使っているためXilinxのDSPの24×17の乗算器を16×16として使っているがICF3-Fでは乗算にRNSは使っていないので24×17として使えている。つまり(24×17)÷(16×16) = 1.59倍、効率良く乗算器を使っている。

このようにループさせるだけでモンゴメリ乗算が正しく演算できることの証明はモンゴメリ乗算の累積加算における分割加算の証明にあります。

大きなモンゴメリ乗算器が実装できる仕組み

f:id:icf:20190812212444p:plain
大きなモンゴメリ乗算器が実装できる仕組み
乗算器が大きくなっていくとuのディレイが間に合わなくなってくる。そこでuのブロッドキャストに中継FFを入れる。 kからk+1のデータはFFで中継。k+1からkのデータを1サイクル早く作る。ここがポイント。このときディレイが心配になるが、それでもクリティカルパスはu生成になると思うので、ここはクリティカルパスにならない。このため、途方もなく大きなモンゴメリ乗算器を実装できる。途方もなく大きな乗算器になってもクリティカルパスが、長くなることはない。つまりDSP当たりの演算効率は、ほとんど下がらない。

ICF3-Fの性能を2倍にする

RSA暗号のSSLアクセラレータにおいて、RSA暗号の鍵長が大きくなるとサバンチ大学の論文では実装できなくなるし、ICF3-Fも実装できても性能が遅すぎる問題がでてくる。そこで2個のICF3-Fを連結して中国人剰余定理を使って、性能をほぼ２倍にする方法が、現実可能だということです。

僕の分割加算と準同型暗号

目的

僕の考えた分割加算が準同型暗号を使ったシステムで役に立ちそう。

分割加算とは

大きな数のモンゴメリ乗算を高速に演算するアルゴリズム。僕が考えたもの。 詳しくは「モンゴメリ乗算の累積加算における分割加算の証明」。

準同型暗号とは

暗号化されたまま加算や、乗算が行える暗号。RSA暗号など、いくつかの公開鍵暗号がある。僕はアルゴリズムや安全性については詳しくない。 電子投票や電子マネーでの応用が期待されている様子。具体的に社会実装するという動きもあるようです。

分割加算は使えるか？

RSA暗号はもちろん。Paillier暗号でも、使えると思われます。 分割加算では基数2n(nは、17～64程度)なので 使えるかどうかは暗号アルゴリズムにあるmod xのxがgcd(2,x)=1であること。すなわちxが奇数かどうかになります。 準同型暗号を使ったシステムで、どのくらいモンゴメリ乗算を使うかによって、実際に分割加算を実装したFPGAなどのハードデバイスが役に立つかが決まります。

分割加算はFPGAに向いている

分割加算を使った暗号プロセッサICF3-FはASICのみならずFPGAでも効率的に実装できる。 ICF3-Fの元となっているICF3は、オープンソースハードウェアとして公式サイトがある。 ただし、まだライセンスの整備が十分でない。 構想の話ですが、ICF3-FとICF3の命令セットを、なるべく共通化して、オープンソースのICF3で広く普及させ、 真に性能が必要なシステムでは、ICF3-Fが売れる。その利益によってオープンソースのICF3が整備されるというような、、、。

仮想マシンの加速支援機構つきの新型8bit CPU

2019年12月29日 修正 XilinxのFPGA(XC7A35TICSG324-1L)を搭載したArtyに実装した結果の更新

概要

新型の8bit CPU(ICF3-Z)を設計してFPGAに実装。仮想マシンの加速支援機構を使った、簡単なスタックマシンで、サンプルプログラムがXilinx FPGA Artix-7(-1)の実機上で動作した。たった約400LUTの小さいCPUで。Intel 8051マイコンのFREEな実装にlight52があります。軽量なマイコンですがCPU部だけで約800LUTあるようです。light52は1命令、6サイクルですが、このICF3-Zは1命令、1サイクルが基本で、時に1サイクルに複数命令の実行ができます。周波数比較もICF3-Zが、おおよそ２倍のようで、圧倒的に勝ったのかも。他にあまり比較できるものがなかったという話もありますが。

仮想マシンの加速支援機構とは

新型8bit CPU(ICF3-Z)なのですが命令コードが32bitで、一般的な8bit CPUと比べてメモリ効率が悪い。そこで16bitの圧縮命令の機能を追加した。この機能が仮想マシンの加速支援機構なのです。 圧縮命令のために作られたもので、加速支援機構として、最終的に成立するのかは、わからないですが、需要があって無理がなければ、加速支援機構として考えていく方向です。

加速支援機構が何の役に立つのか？

超低スペックなCPUでは、ソフトウエアで仮想マシンを実装しても性能的な問題で用途が狭められることがあります。この加速支援機構があれば実用の範囲が広がります。

2019年11月23日 追記 仮想マシンの命令列を大容量の安価なプログラムメモリに置けるので、容量の少ない高価なデータメモリを圧迫しないということも利点です。

なぜ低スペックなCPUで仮想マシンなのか?

仮想マシンをJavaのようなスタックマシンにしたり、レジスタ１個のマシンにしたりできます。 CPUのアーキテクチャの違いを吸収して、いろいろな言語のコンパイラの開発を容易にします。 非C言語が利用できるようになることで、より多くの人がIoTデバイスなどを開発できるようになります。

試作プログラム

プログラムの上部で簡単なスタックマシンの命令を実装して、仮想マシンの命令列をFPGAの実機で動作させたものです。独自のアセンブラなので多少、わかりにくかもしれません。自作したアセンブラでバイナリを生成して、FPGAのプログラムメモリに、置きます。 即値の5をPUSH、メモリの0番の値をPUSH、メモリの1番の値をPUSH、乗算、加算、LEDに結果を出力。 メモリ0番には、最初のところで2を書き込んでいます。メモリ1番には3。つまり3*2+5 = 11がLEDに出力されます。(確認しました)

########################################################
# vm (Virtual Machie)
#
# C Register : SP(Stack Pointer)
# D Register : PRINT (LED x 8)
# MEM[0-255] : Memory
#
# PUSHI: SP=SP-1 , MEM[SP]=nn                     (5cyc)
# PUSHR: SP=SP-1 , MEM[SP]=MEM[nn]                (5cyc)
# POP  : MEM[nn]=MEM[SP] , SP=SP+1                (4cyc)
# ADD  : MEM[SP+1]=MEM[SP]+MEM[SP+1] , SP=SP+1    (7cyc)
# MUL  : MEM[SP+1]=MEM[SP]*MEM[SP+1] , SP=SP+1    (16cyc)
# PRINT: PRINT MEM[SP]      
########################################################
# Initial MEM[0]=0x02 , MEM[1]=0x03
########################################################
    ADDL=N;N=0x02;A=ANS
    STORE;Z=0
    ADDL=N;N=0x03;A=ANS
    B(START)
    STORE;Z=1;C=ANS;D=ANS
########################################################
%NOP:
    JRETURN
%PUSHI:
    ADDR=C;ADDL=N;N=0xFF;C=ANS
    RETURN;COPR;ADDL=N;A=ANS
    ADDR=C;R(ADDR);STORE
%PUSHR:
    ADDR=C;ADDL=N;N=0xFF;C=ANS
    RETURN;R(N);COPR
    MOV;ADDR=C;R(ADDR);STORE
%POP:
    RETURN;R(ADDR);ADDR=C;ONE;C=ANS
    MOV;COPR;R(N);STORE
%ADD:
    R(ADDR);ADDR=C;ONE;C=ANS
    B=REG;R(ADDR);ADDR=C
    B=REG;ADDR=B;A=ANS
    RETURN;ADDL=A;ADDR=B;A=ANS
    R(ADDR);ADDR=C;STORE
%MUL:
    R(ADDR);ADDR=C;ONE;C=ANS;D=ANS
    C=REG;R(ADDR);ADDR=C
    C=REG;ADDR=C;A=ANS
    B=ANS;I=X;X=7
    MUL;ADDL=A;ADDR=B;B=RSH;C=RSH;WAIT
    RETURN;ADDR=C;A=ANS
    R(ADDR);ADDR=D;STORE;C=ANS
%PRINT:
    RETURN;R(ADDR);ADDR=C
    D=REG
START:
    _%PUSHI,0x05,%PUSHR,0x00
    _%PUSHR,0x01,%MUL,0x00
    _%ADD,0x00,%PRINT,0x00
    EXIT(0);@D
END:
    B(END)
    NOP

ちなみに

割込みの実装もついていて約400LUTです。2つの割込みを受けられます。

2019年12月29日更新

XilinxのFPGA(XC7A35TICSG324-1L)に搭載したArtyに実装して面積を確認しています。 論理合成のオプションで結果が違うので、面積重視、性能重視の結果です。

合成方針 LUT数 FF数 BRAM数 LUT-RAM 周波数
面積重視 390 197 1.5 10 150MHz
面積重視 406 197 1.5 10 160MHz
性能重視 481 197 1.5 10 175MHz

BRAM 1.5の内訳はプログラムメモリ 4KBとデータメモリ 2KB

データメモリの先頭32バイトがレジスタ、先頭256バイトがスクラッチパッドメモリ。この合成結果はデータのアドレス線を8bitにしたものです。11bitにすればBRAMを消費することなく2KBのデータメモリをすべて使えます。データのアドレス線を16bitまで設定可能ですがBRAMを消費します。圧縮命令(仮想マシン)を利用する場合はプログラムメモリの先頭の領域を消費します。圧縮命令の命令数によって消費する領域が変化します。

割込みを使ったデバッグ機能の検証をしました。 プログラムのコード中に、ブレークポイントを設定して、そこに、たどり着いた回数を計算して、設定した値(パスカウント)のとき、ブレークさせるみたいな機能が実現できます。 具体的には毎サイクル割込みを入れて、アドレスを確認して、設定されたアドレスなら、指定された回数に到達したかをチェック。 指定された回数の場合は、LEDを点灯させるというプログラムを作成して、FPGAの実機でテストしました。ちゃんと動いているようです。 ただプログラムコード中の一部の命令ではパスカウントを正しく計算できません。遅延分岐で遅延スロットをキャンセルする命令などです。 分岐命令の直後の命令にブレークポイントを設定した場合は、パスカウントは正しくありませんという仕様にしようかと考えています。 割込み処理中でキャンセル信号をチェックしてもいいのですが、 デバッグ機能のためにキャンセル信号をチェックする論理を追加するのは面積最小を追求するプロセッサでは、あまり良くないと考えてのことです。 (面倒とかではなくて)

約400 LUTの面積の小さいCPUで、追加のデバッグモジュールなしにパスカウントを 使ったブレークポイントを設定できるのは、良くできるほうなんだろうと思うが、 他のCPUの面積とかわからないから、なんとも。

講演会に参加しました「IEEE Computer Society 2018会長を経験して」

2018年12月8日(土) 早稲田大学 グリーン・コンピューティング・システム研究開発センターで昨年からIEEEの会長をされている笠原先生の講演会「IEEE Computer Society 2018会長を経験して」に参加しました。 正しくは成田/笠原/木村研究室のOB会で身内だけで開催されました。僕は笠原研のOBで第一期生なのです。定義にもよりますが僕より上の学年は成田研から笠原研に異動しているので僕の学年が第一期かと。 IEEEは、とても有名な学会と思っていましたが、北米から初めて日本に、そして次期会長は、イタリアの方になるようです。ますます有名になっていくように思いました。 会長の仕事には、選挙結果の不平の仲裁みたいなものがあるとか、いろいろあったように思いますが、早稲田の今後の展開についてなどの話があって、あれと思ったら笠原先生からいただいた名刺に「副総長」とあり、なるほどと。早稲田の国際的な、さらなる知名度向上に向けての力強い話が良かったように思いました。

久々に研究室時代の人と話すのですが、僕の同期はOB会に来ないことが多く、１年下も来てませんでした。 ２年下が１人来ていたので話せて、インターネットのシステムに詳しいようでした。 研究室時代の研究上の先輩にRFCとか論文とか書いている人がいてDNSに詳しくSSLアクセラレータの話をしても「ECDSAだよね」とか。^^; 楽しい１日でした。

モンゴメリ乗算の累積加算における分割加算の証明

はじめに

今年の4月にFPGAをはじめたときにFPGAにDSP(乗算器)が多数あることを知った。 これでモンゴメリ乗算が高速化できるわけだがアルゴリズムをそのまま実装するとビット長の大きな累積加算が、ビット長とともに周波数が下がり性能がでない。 そこで分割して加算しても、正しい加算結果になるような方法を考えた。これを使ったSSLアクセラレータ ICF3-F(商用版)の実装を急いでいるところだが、証明もしておかないと、いけないかと思って、急いで公開することにしました。産業スパイによる海賊版を抑止する目的です。ちなみに他にも応用ができる方法ではないかと思います。 この証明にコメントしたい方が、いらしゃるのかわかりませんが、もしあれば公開日(2018年9月26日)の２日後、28日の午後12時からでお願いします。とりあえず、その記録は残すようにします。

問題

非負の整数の変数Aを2進数表現すると、最下位ビットからrビットづつ区切ってs個のブロックに分割して表現できる。

\displaystyle A = \sum_{k = 0}^{s-1} a_k \cdot 2^{r k}

a_k がrビットの場合はAはr×sビットの数を表現可能だが、 a_kが(r+1)ビットあってもr×sビットの数を表現できる。 最上位ビットを e_k、それ以外を f_kとすると a_k = e_k \cdot 2^r + f_kになるが

　\displaystyle A = \sum_{k = 0}^{s-1} e_k \cdot 2^{r (k+1)} + \sum_{k = 0}^{s-1} f_k \cdot 2^{r k}

として全加算してやればAは冗長性がなくなりr×s+2個のビットの数になる。 このような冗長性をもったAを考える。

A = (A + t \cdot B + u \cdot C ) \div 2^{n} --- (1)

r \geqq n + 2 --- (2)

t,uはnビットの数。B、Cはr×s-2ビットの数。 \displaystyle B = \sum_{k = 0}^{s-1} b_k \cdot 2^{r k}、 \displaystyle C = \sum_{k = 0}^{s-1} c_k \cdot 2^{r k}と表す。

式(1)の右辺を計算した結果、任意のA、B、C、t、uにおいて、 同じ冗長性、すなわち a_kが(r+1)ビットある左辺Aに正しく格納できることを証明せよ。

証明

式(1)右辺のカッコ内は

A + t \cdot B + u \cdot C

\displaystyle = \sum_{k = 0}^{s-1} a_k \cdot 2^{r k} + t \times \sum_{k = 0}^{s-1} b_k \cdot 2^{r k} + u \times \sum_{k = 0}^{s-1} c_k \cdot 2^{r k}

\displaystyle = \sum_{k = 0}^{s-1} a_k \cdot 2^{r k} + \sum_{k = 0}^{s-1} t \cdot b_k \cdot 2^{r k} + \sum_{k = 0}^{s-1} u \cdot c_k \cdot 2^{r k}

\displaystyle = \sum_{k = 0}^{s-1} ( a_k \cdot 2^{r k} + t \cdot b_k \cdot 2^{r k} + u \cdot c_k \cdot 2^{r k} )

\displaystyle = \sum_{k = 0}^{s-1} ( a_k + t \cdot b_k + u \cdot c_k ) \cdot 2^{r k}

ここで p_k = a_k + t \cdot b_k + u \cdot c_kとする。t,uがnビットであるから p_kは最大でもr+n+2ビット。 p_kの最上位ビットから2bitを px_k、最上位ビットから2ビット、最下位ビットからnビットを除いたrビットを py_k、 最下位ビットからnビットを pz_kとする。

p_k = px_k \cdot 2^{r+n} + py_k \cdot 2^{n} + pz_k (0 \leqq k \leqq s-1) --- (3)

(1)の右辺は

\displaystyle \sum_{k = 0}^{s-1} p_k \cdot 2^{r k} \div 2^{n}

(3)より

\displaystyle = \sum_{k = 0}^{s-1} ( px_k \cdot 2^{r+n} + py_k \cdot 2^{n} + pz_k ) \cdot 2^{r k} \div 2^{n}

k=0の項をシグマの外に出して整数にすると

\displaystyle = px_0 \cdot 2^{r} + py_0 + \sum_{k = 1}^{s-1} ( px_k \cdot 2^{r+n} + py_k \cdot 2^{n} + pz_k ) \cdot 2^{r k} \div 2^{n}

\displaystyle = px_0 \cdot 2^{r} + py_0 + \sum_{k = 1}^{s-1} ( px_k \cdot 2^{r} + py_k ) \cdot 2^{r k} + \sum_{k = 1}^{s-1} pz_k \cdot 2^{r k - n}

\displaystyle = px_0 \cdot 2^{r} + py_0 + \sum_{k = 1}^{s-1} ( px_k \cdot 2^{r} + py_k ) \cdot 2^{r k} + \sum_{k = 1}^{s-1} pz_k \cdot 2^{r(k-1) + r-n}

\displaystyle = py_0 + pz_1 \cdot 2^{r-n} + \sum_{k = 1}^{s-2} ( px_{k-1} + py_k + pz_{k+1} \cdot 2^{r-n} ) \cdot 2^{r k} + ( px_{s-2} + py_{s-1} ) \cdot 2^{r (s-1)} + px_{s-1} \cdot 2^{r s}

ただしB、Cがr+n-2ビットなので必ず px_{s-1} =0。

\displaystyle = py_0 + pz_1 \cdot 2^{r-n} + \sum_{k = 1}^{s-2} ( px_{k-1} + py_k + pz_{k+1} \cdot 2^{r-n} ) \cdot 2^{r k} + ( px_{s-2} + py_{s-1} ) \cdot 2^{r (s-1)}

ここで式(2)の条件からr-nは2以上。 w_k = px_{k-1} + pz_{k+1} \cdot 2^{r-n} (1 \leqq k \leqq s-2)とすると w_kはrビットの記憶素子を持つ1つの変数としてまとめられる。

\displaystyle = py_0 + pz_1 \cdot 2^{r-n} + \sum_{k = 1}^{s-2} ( py_k + w_k ) \cdot 2^{r k} + ( px_{s-2} + py_{s-1} ) \cdot 2^{r (s-1)}

これは(1)左辺のAとして次にようになる。

\displaystyle A = \sum_{k = 0}^{s-1} a_k \cdot 2^{r k}

a_k = \left\{ \begin{array}{} py_0 + pz_1 \cdot 2^{r-n} ( k = 0) \\ py_k + w_k ( 1 \leqq k \leqq s-2 ) \\ ( px_{s-2} + py_{s-1} ) ( k = s - 1) \end{array} \right.

py_k、 w_kともにrビットの数なので加算してもr+1ビットの数である。k=0、k=s-1も同様にr+1ビット。 式(1)の右辺を計算した結果、同じ冗長性、すなわち a_kが(r+1)ビットある左辺Aに正しく格納できる。

補足

実際の実装では最後に、全加算を行って冗長性をなくしてr×sビットのレジスタに格納します。 このとき桁あふれが発生する可能性がありますが、一般のCPUのレジスタと同じで桁あふれはアルゴリズム側の責任。 モンゴメリ乗算では最終的な計算結果はCの2倍以下になるので全加算後、r×sビットのレジスタに必ず格納できます。 計算過程においても、上記証明によって、桁あふれすることなく計算できていることがわかります。 この問題では積和の数が２個の場合ですがn個にした証明など一般化できるように思います。いづれまた。