読者です 読者をやめる 読者になる 読者になる

暗号計算機屋のブログ

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

nhira型モンゴメリ乗算アルゴリズム

とりあえず

この証明だけは、まだ、役に立つかも。というところ。

はじめに

RSA暗号を高速化するためのアルゴリズムとしてモンゴメリ乗算が有名だが、その改良型がC.C.Yang(IEEE 1998)によって論文発表された。それをもとにDesign Wave MagagineのRSA暗号コンテスト(2008年)で実装されたものが記事になっていた。それをヒントに、この参考文献の14.36のモンゴメリ乗算を改良したnhira型モンゴメリ乗算のアルゴリズムを説明する。自身が作ったアルゴリズムだがC.C.Yangと同じかもしれない。C.C.Yangの論文を取り寄せる前に自身で証明できてしまったため、C.C.Yangの論文を読む必要がなくなった。

ちなみにDesign Wave Magagineはモンゴメリリダクションを使ったモンゴメリ乗算にC.C.Yangを導入した感じだが、これは乗算とリダクションを同時に行うモンゴメリ乗算にC.C.Yangを導入した感じだと思っている。

nhira型モンゴメリ乗算

 {
INPUT : integers m  = (m_{n-1} … m_1 m_0)_b , x=(x_{n+1} x_n x_{n-1} … x_1 x_0)_b , \\
         y  = (y_{n+1} y_n y_{n-1} … y_1 y_0)_b \\
    with 0 ≦ x,y ≦ 2m-1 , R=b^{n+2} with gcd(m,b)=1 , \\
    and m' = -m^{-1} mod b \\
 \\
OUTPUT : xyR^{-1} mod m. \\
 \\
1. A←0.  (Notation : A = (a_n a_{n-1} … a_1 a_0)_b ) \\
 \\
2. For i from 0 to (n+1) do the following: \\
 2-1. u_i ← (a_0 + x_i y_0 ) m'  mod b. \\
 2-2. A ← (A + x_i y + u_i m) / b. \\
 \\
3. Return (A).  0≦A≦2m-1
}

【解説】 

xの範囲は2m-1までなので {x_{n+1}}は必ず0。 わかりやすくするためb=2として説明する。  Forループで0≦i≦n のとき0≦A<3m-1 とした場合、次のAの値の範囲は、2.2の式から (3m-2 + 2m-1 + m) / 2 = (6m-3) / 2 = 3m-1.5 以下なので0≦A <3m-1 の範囲になる。 Aの初期値は0なのでi=nのループ終了時のAは0≦A<3m-1 となる。Forループでi=n+1 のとき0≦A<3m-1なので、次のAの値の範囲は 2.2の式から(3m-2 + m ) / 2 = 2m-1 以下、すなわち0≦A≦2m-1

【解説2】(2017年2月26日 追加)

オリジナルの14.36に対してx,yがmを超える値となっていて正しい結果になる理由は、14.36は14.32のモンゴメリリダクションと同じ計算しかしていないということ。14.32で {a_i}を求めるタイミングまでにx×yを計算させて、最初からA=x×yであったのと同じ状態であればいい。14.32の定理でR= {b^{n}}をR= {b^{n+2}}とすればx,yが2m-1以下の範囲であれば正しい結果が得られる。

【参考文献】

Design Wave Magagineのコンテストの応募原稿にあった参考文献。(私は、読んでない) f:id:icf:20170225131209p:plain

(2017年2月27日 追加)

詳しく読んでないが、ここで説明するアルゴリズムを使って、やりたかったことが2008年のIEEEの論文になっていました。
↓↓↓これです↓↓↓

A New Modular Exponentiation Architecture for Efficient Design of RSA Cryptosystem - IEEE Xplore Document

HDD再生の研究

失敗は成功のもと

HDDを廃棄して都市鉱山の利用研究をするよりHDDを再生したほうが地球環境にいい。そんなことを考えている。そして再生したHDDによってお金の節約も狙う。個人でわかる範囲で調べているため、もっと詳しい人はいるだろうが、もっといい情報などあれば、コメントいただけるとうれしい。

利用環境はサーバーではなく普通の環境だ。BUFFALOのUSB HDD(320GB)の調子が悪くUSB 2.0で接続されているにもかかわらず転送速度が300KB/Secに落ちることがあるようになった。HDD-SCANを使って調べると32個のバッドセクタが見つかった。2007年のWestern Digital製なのだがバッドセクタは散在するため、パーティションでバッドセクタを避けると容量が激減するため、なんとかバッドセクタを修復する無料のソフトを探した。あるにはあったが実際にやってみると「これは修復できません」ということだった。そこでついにHDD Regeneratorを購入することにした。購入についての経験は次節で述べる。32個あったバッドセクタのすべてを修復したが、フォーマットをするとタイムアウトすることが多く結局、Western Digitalの320GBのHDDの再生には失敗した。 そこで先日、HDDレコーダから取り出した再生したHDDと交換した。

2016年9月10日 追加 HDD Regeneratorによるバッドセクタの修復(リカバリ)は32個所のすべてで成功している。過去にHDD Regeneratorの無料版で何度もバッドセクタの修復してHDDの再生に成功している。今回は、たまたまバッドセクタ以外の故障によるものと思われる。バッドセクタがあってもファイルシステムの機能を使ってHDDを利用可能にすることはできるが、HDDの運用を考えるとバッドセクタを修復してしまったほうが便利。 f:id:icf:20160909214208j:plain

HDD Regeneratorの購入についての経験

購入のページに移動すると自動的に円で購入できるが、左上にあるセレクトボックスでUSドルに変更する。相場にもよるのかもしれないが10800円が10200円で購入できたので約600円くらい安く買えた。Windows10対応みたいな表記があるが2011年版のようである。シリアル番号が送られてくるが、シリアル番号の名前とシリアル番号を入力しないとアクティベートできない。シリアル番号の名前はPayPalで購入したせいかもしれないが漢字が自動入力されて文字化けで困った。購入する前によく調べたほうがいいです。

世界一のRSA暗号LSI ICF3のアルゴリズム

アルゴリズムの解説ではなく僕が1998年にICF3の開発を開始したときに渡された資料の公開です。学生の頃は情報系の科目をメインに取っていましたが電気工学科の卒業ということもあり会社に入ってすぐの頃は電子回路シミュレーションをしていました。次にIBMLSIの配線データを日立のLSIの配線データに変換するエクセル作業、大型コンピュータの基板の配線長の計算、夜間勤務でクロック周波数のマージン測定、アメリカでFAXを受信するだけの仕事、暗号装置の線表管理と、雑用が続いた。その後、暗号LSI ICF1の改造、暗号LSI ICF2の改造をするのだが、論理設計の勉強など全くしないで改造作業をした。改造している時間は長かったが基本的なことを知らずに応用問題をやっていた感じで、技術が身についていくような感覚はなかった。このまま雑用係としていくのかと思えば1998年のある日、課長に呼ばれ「RSA暗号の論理設計をしてみないか?」と聞かれた「え、やっていいんですか?」と返答したのを覚えている。普通は論理の勉強のため7セグメントLEDとかのサンプル回路とか作り始めるのかもしれないが、そんな時間が与えられるはずもなく、いきなりRSA暗号の演算器の絵を描き始めた。CPUっぽいものだったので周囲からはポンチウムと笑われた。ポンチ絵+Pentiumの合成語。だが徐々に精度が上がっていき加算器4個を同時に動作させる演算器ができてきた。周囲には「4IPターボ」とか笑われた。IPとは、この部署固有の言葉でCPUの意味。本格的に4IPターボを作ろうとしたときに、日立製作所システム開発研究所から次のFAXが届いた。

f:id:icf:20160904202228p:plain

原書がインターネットでFreeで公開されています。FAXは14章の1ページです。

僕が数学科の卒業なら有限体を勉強していたはずで、もしかすると楽に読めたのかもしれないが、2日間、考えてFAXにある14.36のアルゴリズムからRSA暗号の計算方法を思いつくことに成功した。後から聞いた話、2日間で理解したのは評判になるくらい速かったらしい。 つい先日、原書がインターネットで公開されているという情報を教えてくれた人と話をして14.36を使う以外の方法でもモンゴメリを使ったRSA暗号の計算ができることを知った。ハードウェアで実装する場合と、ソフトウェアで実装する場合とで、若干、有利な方法が違うのだろうと思う。

廃棄するHDDレコーダからHDD 80GBを回収

日立製作所 2003年製のHDDレコーダを放置していたがHDDが足りなくなったので、HDDレコーダから回収することに。 中を開けてみると、でてきたHDDはSAMSUNG製(韓国)だった。

日立すら買わない日立のHDDを誰が買うのだと思いつつ、日本のHDD事業について調べてみた。縮小傾向だが東芝がまだやっているようだ。少し調べたくらいで、どうこう言えるほどではない。

しかし、一般的な話として、強く思うことはある。

経営者からしてみれば、海外の安いHDDを導入しHDDレコーダの価格を下げたほうが、いいのはわかる。しかし、その結果、将来、HDD事業を撤退することになった。

経営者に、将来日本が困るから、HDDを買い支えるよう、指示する仕組みは、あるのだろうか? なければ、同じことが、繰り返される。 日本人が愚かな民だと思いたくはない。

海外製のHDDに、リモートでHDD上のファイルを改ざんしてウィルスが発生する仕組みがあるかもしれない。 海外のHDDの工場にいって、1台1台製造に付き合っていたら、トータルで、割高になってしまうとか、なりそうで怖い。 なにか事業をしようとして、HDDの信頼性の問題に気が付き、割高くらいでは、すまなくなって後悔しまくることにならないだろうか?

f:id:icf:20160823230244j:plain

3、5、7で割った余りを計算する演算器

はじめに

1999年に世界一高速なRSA暗号のLSIを開発したが、その後、いろいろ研究していたときに見つけたのが7で割った余りを高速に演算する論理だった。 今では高校の数学で整数の性質で習う合同式を使えば簡単に証明できるが、それは2014年以降の話で、2000年頃にはまだそういうものはなかった。

3で割った余りが何の役に立つの?

LSIではノイズや宇宙線によってビットが反転し正しく動作しないことがある。パリティによってチェックすることもあるが乗算器のビット反転はパリティでは不可能だ。 乗算器を2重化すればいいが、トランジスタ数が増大する。そこでレシジュという方法が用いられていた。乗数と被乗数の3で割った余りを乗算して、乗算結果と比較するのだ。 3で割った余りどうしの乗算なので2bitの乗算器程度のトランジスタ数でいいためトランジスタ数の増大を防げる。 このときあった3で割った余りを計算する演算器を5、7に応用したのが「3、5、7で割った余りを計算する演算器」です。 3で割った余りを計算してくれる演算器に証明がついていれば良かったのだが、そういう便利なものはなく、自分で証明するしかなかった。ちなみに15で割った余りも可能です。 「すべてがFになる(計算機工学版)」がなぜ計算機工学なのかは、そういう理由からです。

3、5、7で割った余りが何の役に立つの?

3ではなく5や7で割った余りを使えば乗算器の誤りを検出する精度が上がる。が、それはおいておいて 素数判定をするにはミラー・ラビン素数判定法のような確率的アルゴリズムでも時間のかかる処理です。 そういった時間のかかる素数判定の前処理として「3、5、7で割った余りを計算する演算器」を使えば素数判定を高速化できる、、、と思っていた。 事実、GNUのGMPというライブラリでmpz_nextprime()という関数は3、5、7のような小さい素数で割り切れるかをテストしている。 しかし与えられた数よりも大きい素数を探すのに毎回、3、5、7で割れるかを判定しているわけではなく、前回の結果を記録しておいて、 増分を計算して判定するため、「3、5、7で割った余りを計算する演算器」の高速性は活かされないことが判明した。 VHDLなどの論理シミュレータでレシジュのシミュレーションを高速化できるなどの超マイナーな用途は、見つけたが、みなさんにも 考えてもらえるといいかなぁと。

7で割った余りを高速に演算する方法

f:id:icf:20160711083634p:plain f:id:icf:20160711083654p:plain f:id:icf:20160711083708p:plain

7で割った余りが演算できることの証明

計算機なので2進数です。3bitづつにまとめて8進数で考えます。 8進数は

 {
\begin{align}
A_n×8^n + A_{n-1}×8^{n-1} + … + A_2×8^{2} + A_1×8^{1} + A_0 
\end{align}
}

と表せます。

 {
A_n×8^n + A_{n-1}×8^{n-1} + … + A_2×8^{2} + A_1×8^{1} + A_0 (mod\;7)\\
= A_n + A_{n-1} + … + A_2 + A_1 + A_0\\
  \;\; + A_n×(8^n-1) + A_{n-1}×(8^{n-1}-1) + … + A_2×(8^{2}-1) + A_1×(8^{1} -1) (mod\;7) 
}

mが1以上の任意の整数に対して次が成り立ので

 {
\begin{align}
8^m-1 &= 7×8^{m-1} + 8^{m-1} - 1\\
 &= 7×8^{m-1}+7×8^{m-2}+8^{m-2} - 1\\
 &= 7×8^{m-1}+7×8^{m-2}+7×8^{m-3}+8^{m-3} - 1\\
 &= 7×8^{m-1}+7×8^{m-2}+7×8^{m-3}+7×8^{m-4}+…+8-1 \\
 &= 7×8^{m-1}+7×8^{m-2}+7×8^{m-3}+7×8^{m-4}+…+7 \\
 &= 7×(8^{m-1}+8^{m-2}+8^{m-3}+8^{m-4}+…+1)
\end{align}
}
 {
= A_n + A_{n-1} + … + A_2 + A_1 + A_0\\
  \;\; + A_n×(8^n-1) + A_{n-1}×(8^{n-1}-1) + … + A_2×(8^{2}-1) + A_1×(8^{1} -1) (mod\;7)\\
= A_n + A_{n-1} + … + A_2 + A_1 + A_0 (mod\;7)
}

あとはAnからA0を7で割った余りを求めながら加算していけば答えが計算されます。 これは7以外の3、15、31でも同じ方法で証明できます。

3,5で割った余りは、どうするの?

3だけなら7と同じ方法でいいのですが15で割った余りを計算してから3と5で割った余りを計算すると 一度に3,5が計算できます。

GMPよりも高速に計算できる?

GMPとはThe GNU MP Bignum Library の略で多倍長計算を非常に高速に行うライブラリです。 ハードウェア向きのアルゴリズムですがソフトウェアで実装するとどうなるか、やってみました。

mpz_fdiv_ui(p,3)

これを代替する関数をC言語で作ってみました。WORDSは8バイト単位で1024bitなら16です。 3固定のコードですがmpz_export()でデータを取り出さないといけない分、代替関数が不利です。

unsigned long mod3(mpz_t p,int WORDS) {
    static unsigned long long w[64];
    unsigned long W;
    size_t count = WORDS;       
    mpz_export (w, &count, -1, sizeof(unsigned long long), 0, 0, p);
    W = (unsigned long)(w[0]%15);
    for(int j=1 ; j < WORDS ; j++) W += (unsigned long)(w[j]%15);
    return W%3;
}

恐らく環境によって違ってくると思いますがCore i5 4690Tでは2048bitくらいまでは 代替関数が高速です。それより長いbit長ではGMPに負けます。 Core i5のAVX2命令を使ってもう少し高速化したものが、次です。 はじめてAVX2命令を使ったので十分に高速化できてないかもしれないですが C言語の代替関数よりも20%くらい高速化されました。

unsigned long mod3_avx(mpz_t p,int WORDS) {
    static unsigned long long w[64] __attribute__((aligned(32)));
    unsigned long long W;
    __m256i W0,W1,W2,W3,W4,W5,W6,W7;
    size_t count = WORDS;
    mpz_export (w, &count, -1, sizeof(unsigned long long), 0, 0, p);
    W = 0;        
    for(int j=0 ; j < WORDS/16 ; j++) {
        W0 = _mm256_load_si256((const __m256i *)(w+j*16+0));
        W1 = _mm256_load_si256((const __m256i *)(w+j*16+4));
        W2 = _mm256_slli_epi64(W0,32);
        W += (unsigned long)(w[j*16+8]%15);
        W4 = _mm256_slli_epi64(W1,32);
        W += (unsigned long)(w[j*16+9]%15);
        W3 = _mm256_srli_epi64(W0,32);
        W += (unsigned long)(w[j*16+10]%15);
        W5 = _mm256_srli_epi64(W1,32);
        W += (unsigned long)(w[j*16+11]%15);
        W6 = _mm256_srli_epi64(W2,32);
        W += (unsigned long)(w[j*16+12]%15);
        W7 = _mm256_srli_epi64(W4,32);
        W += (unsigned long)(w[j*16+13]%15);
        W0 = _mm256_add_epi64(W3,W6);
        W += (unsigned long)(w[j*16+14]%15);
        W1 = _mm256_add_epi64(W5,W7);
        W += (unsigned long)(w[j*16+15]%15);
        W2 = _mm256_add_epi64(W0,W1);
        W += _mm256_extract_epi64(W2,0);
        W += _mm256_extract_epi64(W2,1);
        W += _mm256_extract_epi64(W2,2);
        W += _mm256_extract_epi64(W2,3);
    }
    return W%3;
}

最後に

ここにあるような余りを求める計算など整数の性質が実社会で何の役に立つか? 僕の経験を話します。 IBMの暗号LSIを買ってきて日立の大型コンピュータの暗号装置として製品化することがあった。 国内向けの需要で暗号LSIから秘密鍵をバックアップできる必要があった。 しかしIBMが、その方法を教えてくれるはずもなかった。 なぜなら秘密鍵の盗難の方法を教えるのと同じだからだ。 別の暗号装置にリストアすることはできないが、同一の暗号装置ならバックアップできる方法はあった。 暗号装置にRSA公開鍵暗号を設定することで、それが可能になる。 説明は省略するがRSAの鍵は周期が大きくなるように作る。 そこに敢えて1とnを繰り返す周期の鍵を設定する。 暗号装置内部で生成される乱数によって1かnかになる。確率50%で1になる。 1であることがわかれば、それを利用してバックアップする方法が見つかった。 確率50%では、何度か繰り返せば、普通にバックアップが可能だ。

この話が整数の性質の学習の動機付けになればと思います。

すべてがFになる (計算機工学版)

【すべてがFになる】「すべてがFになる (計算機工学版)」イラスト/izuna [pixiv]