月別アーカイブ: 2016年10月

モンゴメリ乗算

#このページがGoogle検索の上の方に引っかかってるみたいなので,慌てて内容を充実させた.何かおかしな点やわかりづらい点があればコメントを頂ければと思う

#2017/9/30 わかりやすいように編集

\(N > 0\)を法とした合同式を計算する場合,\(N\)による割り算は避けては通れない問題となる.四則演算の中で,割り算が最も遅いことを考えれば,この剰余を求める操作を高速に行うのが簡単ではないことは明らかだ.
しかし,世の中には不思議なことにこれを比較的高速に演算する方法がある.その一つとして今回はモンゴメリ乗算を紹介しよう.

モンゴメリ乗算の形式的説明

まずモンゴメリ乗算を説明するための前準備として,次のモンゴメリ写像\(f\)とモンゴメリリダクション\({\rm MR}\)を導入する.

$$f(x)=xR \mod N \\
{\rm MR}(x) = xR^{-1} \mod N$$

このとき,\(N\)は\(R\)と互いに素であり,かつ\(R > N\)が成り立つとする.
\(R^{-1}\)は,\(R\times R^{-1} = 1 \mod N\)となるような整数である.これは,例えば5が法の場合は,\(R=2\)であれば,\(R^{-1}=3\)となる.また,上式の関係から,モンゴメリ写像とモンゴメリリダクションは逆写像の関係にあるといえる.

次に\(N\)より小さな整数,\(a\),\(b\)を用意し,その積の剰余\(ab \mod N\)を考えたい.ここではモンゴメリリダクションを使って,うまいこと表現できないかを考える.そのためにまず

$$
\begin{align}
f(ab) &= abR \mod N \\
&= aR \times bR \times R^{-1} \mod N \\
&= f(a)f(b)R^{-1} \mod N
\end{align}
$$

であることに注意する.今,明らかに\(f(a) = {\rm MR}(aR^2)\)が成り立つので,

$$
f(ab) = {\rm MR}({\rm MR}(aR^2){\rm MR}(bR^2))
$$

両辺に対して,モンゴメリリダクションを取れば

$$
ab \mod R = {\rm MR}({\rm MR}({\rm MR}(aR^2){\rm MR}(bR^2)))
$$

さて,ここまでは単純に式の形を変えたにすぎない.当初の目的は如何に\(\mod N\)を高速に計算するかだった.それをモンゴメリリダクションを用いた形に直したところで,その中に剰余演算があるのでは意味がないように思われる.ところがこのモンゴメリリダクションを高速に計算するアルゴリズムがある.すなわち

$${\rm MR} (T) = \begin{cases} g(T) & (g(T) < N) \\ g(T)-N & (otherwise) \end{cases}$$

$$g(T)=\frac{T+(TN’ \mod R)N}{R}$$

ただし$$NN’\equiv -1 \mod R$$が成り立つ.

モンゴメリリダクションの計算方法を見ると,今までネックだった\(\mod N\)がいつの間にかなくなり,代わりに\(\mod R\)だけになっている.他にも\(R\)による除算が増えているが,結局どちらも\(R\)を恣意的に選べることから,\(R\)を2のべき乗にすれば,アンド演算とシフト演算のみでそれぞれ計算できる.残りは除算以外の演算なので,それほど重くない.よって,モンゴメリリダクションで\(ab\)を高速に計算する方法がわかったことになる.

さて,ここまではモンゴメリ乗算についての形式的な説明である.モンゴメリ写像から始まり,モンゴメリリダクションについて紹介した.また,その高速な計算方法を天下り的に与えたので,上の話は多くの人にとって納得できるものではないだろう.そこでここからは,上記の内容について1つずつ整理して説明していく.

モンゴメリ乗算の意味

まず,そもそもモンゴメリ乗算とは,ある数の積をモンゴメリ写像したものを指す.例えば\(N\)未満の数,\(a,b\)を用意して,
$$
\begin{align}
f(ab) &= abR \mod N \\
&= (aR)(bR)R^{-1} \mod N \\
&= f(a)f(b)R^{-1} \mod N
\end{align}
$$
はモンゴメリ乗算である.この定義から明らかだが,別にモンゴメリ写像された2つの変数の積は,モンゴメリ乗算ではない.その後, \(R^{-1}\)を掛けてから\(\mod N\)を取って,はじめてモンゴメリ乗算となる.そしてその「\(R^{-1}\)を掛けてから\(\mod N\)」を取る演算がちょうどモンゴメリリダクション\(\rm MR\)の定義になっている.

また別の言い方をすれば,モンゴメリ表現(モンゴメリ写像されたもののこと)された2つの変数の積をモンゴメリ表現(モンゴメリ乗算の結果)にするには,積に対してモンゴメリリダクションを取ればいい.よって,例えば\(abc \mod N\)を求めたいとすれば,まず\(a,b,c\)それぞれのモンゴメリ表現\(A,B,C\)を用意し,次に
$$
D = {\rm MR}(AB)
$$
を求める.\(D\)はモンゴメリ表現なので,\(D\)と\(C\)のモンゴメリ乗算を考えることが出来て,
$$
E = {\rm MR}(DC)
$$
となる.モンゴメリ写像の逆変換はモンゴメリリダクションなので,結局\(abc \mod N\)は
$$
{\rm MR}(E)={\rm MR}(f(abc))=abc
$$
で与えられる.

このように,与えられた変数をすべてモンゴメリ表現にしてから,掛けてはモンゴメリリダクションを行っていけば,最終的にすべての変数の積のモンゴメリ表現を得ることが出来る.あとはそれをモンゴメリリダクションすることで知りたい答えが得られるというのがモンゴメリ乗算の全容である.

そうすると,モンゴメリ乗算は,モンゴメリリダクションが肝であることがわかる.もしこのモンゴメリリダクションが高速に行えればモンゴメリ表現を素早く求められることになるので,結果的に剰余演算が高速できることになる.そして,そのモンゴメリリダクションを高速に計算する方法が先程示した計算法だったわけだ.ではそこにある仕組みは一体何なのだろうか.それを探るためにもう一度モンゴメリリダクションの計算式を以下に示す.

$${\rm MR} (T) = \begin{cases} g(T) & (g(T) < N) \\ g(T)-N & (otherwise) \end{cases}$$

$$g(T)=\frac{T+(TN’ \mod R)N}{R}$$

ただし$$NN’\equiv -1 \mod R$$が成り立つ.

まずは \(g\)の式に注目しよう.こいつの分子は与えられた引数\(T\)に\(N\)の定数倍を足し合わせている.したがって,この分子の式自体は\(\mod N\)において,値を保存することは確かだろう.次に,面白い事実として\((TN’ \mod R)N\)は\(T\)の\(R\)以下のビットの桁をすべて繰り上げて0にするような数となっている.これは例えば\(R\)が128であれば,\(T+(TN’ \mod R)N\)は下位7ビットまでがすべて0となるということだ.試しにこのことを確認するためのC++のコードを次に示す(一部wikipediaのパクリ).

上記のコードのT,R,Nの値を色々と変更してみてほしい.どのパターンでも述べたとおりになっているだろう.要するに,これは\(T\)の下位桁のビットがすべて0になるような数を足し合わせることで,下位桁の値を上位桁に寄せ集めている事になっている.この寄せ集める操作は,あくまで\(\mod N\)においては保存するので,結果的にビットの広がりを小さくする効果のみを持つ.最後に寄せ集めた数を\(R\)で割って,右にシフトさせることで,\(N\)までの範囲に収めようとしているのである.ただし,\(T\)の値によっては,\(g(T)\)が\(N\)を少しだけはみ出すことがあり得るので,それを最後の条件式で補正している.

以上がモンゴメリリダクションの仕組みである.

内容についてどう思いますか?
  • いいね (15)
  • だめ (0)

sptkのwav2rawのエラー

wav2rawを使おうとしたら

uid: Undefined variable.

というエラーが出た.

よくわからないがどうやらuidという変数がないことを怒られていることがわかるので

export uid=$UID

をしたら動作した.

備忘録

内容についてどう思いますか?
  • いいね (0)
  • だめ (0)