繰り返し自乗法
10月 8th, 2012 | 
Author: [latexpage]
繰り返し自乗法とは、ある自然数$a$の$m$を法としたべき乗、つまり、
\begin{eqnarray}
a^k \bmod m
\end{eqnarray}
を高速に求める計算法です。
この計算は、べき乗数$k$が小さい場合は何ともありませんが、たとえば、
\begin{eqnarray}
13^{400} \bmod 31
\end{eqnarray}
のような場合は、単純にべき乗していくと値が指数爆発して簡単には求まりませんし、毎回積のmodをとっていくにしても計算量が多くなってしまいますね。この場合は法$m$が素数ですから、フェルマーの小定理によって計算が楽になりますし、そうでない場合も法$m$のオイラー関数$\phi (m)$が計算できれば、オイラーの定理が使えます。しかしながら、この考え方は、法$m$の素因数分解が出来ることが前提ですし、何より一般的ではありません。
指数のべき乗の余りなんて、どこで使うんだと思う方もいると思います。具体的な例として、RSA暗号が挙げられます。(参考：RSA暗号の数理 http://nsplat.wordpress.com/2012/02/12/)
暗号化と復号化は、まさにこの巨大数のべき乗です。
そこで本記事では、一般の法$m$に対して用いることが出来る繰り返し自乗法という方法を考えます。
繰り返し自乗法の手順は3ステップに分けられます。
1. 指数$k$を二進展開する
2. 指数の底$a$の自乗を計算し、これを繰り返し自乗していく
3. 1. で求めた二進数の{0, 1}に対応して、1のときだけ、2.の値を掛け合わせる
では、式(1)に関して計算してみましょう。
1. 指数$k$を二進展開する
$k$を二進展開します。
\begin{eqnarray}
k=u_0\cdot 2^0+u_1\cdot 2^1+u_2\cdot 2^2+\cdots
\end{eqnarray}
となります。
係数$u_i$は二進数に展開したときの第$i$位の値です。
これを指数の肩に載せてみましょう。
\begin{eqnarray}
a^k=a^{u_0\cdot 2^0+u_1\cdot 2^1+u_2\cdot 2^2+\cdots}
\end{eqnarray}
2. 指数の底$a$の自乗を計算し、これを繰り返し自乗していく
次に$a$を$A_0$とし、$A_0$の自乗$a^2$を$m$で割った余りを計算し、これを$A_1$とします。つまり
\begin{eqnarray}
a\equiv A_0 \bmod m
\end{eqnarray}
\begin{eqnarray}
\left(A_0\right)^2\equiv A_1 \bmod m
\end{eqnarray}
この得られた$A_1$をさらに自乗し、$m$で割った余りを$A_2$とします。これを繰り返し、$A_2, A_3, \cdots $としていきます。
\begin{eqnarray}
\left(A_1\right)^2\equiv A_2 \bmod m
\end{eqnarray}
\begin{eqnarray}
\left(A_2\right)^2\equiv A_3 \bmod m
\end{eqnarray}
この$A_i$は結局、$a^{2＾i}$を$m$で割った余りに相当することがわかると思います。
すなわち、
\begin{eqnarray}
a^{2^i}\equiv A_i \bmod m
\end{eqnarray}
3. 1. で求めた二進数の{0, 1}に対応して、1のときだけ、2.の値を掛け合わせる
さて、式(4)を指数の性質から積に分解すると、
\begin{eqnarray}
a^k=a^{u_0\cdot 2^0}\cdot a^{u_1\cdot 2^1}\cdot a^{u_2\cdot 2^2}\cdot \cdots}
\end{eqnarray}
と表せます。
これを法$m$で考え、式(8)を適用すると、
\begin{eqnarray}
a^k&=&a^{u_0\cdot 2^0}\cdot a^{u_1\cdot 2^1}\cdot a^{u_2\cdot 2^2}\cdot \cdots\\
&\equiv& {A_0}^{u_0}\cdot {A_1}^{u_1}\cdot {A_2}^{u_2}\cdot \cdots\bmod m
\end{eqnarray}
したがって、二進展開した$k$の二進数における各位$u_i$が$1$に相当する$A_i$の積をとれば、求める指数が得られることがわかります。
これが高速であることを示すために、式(2)の値を具体的に計算してみましょう。
まず、400を二進展開します。
\begin{eqnarray}
400 = 0\times 2^0+0\times 2^1+0\times 2^2+0\times 2^3+1\times 2^4\\+0\times 2^5+0\times 2^6+1\times 2^7+1\times 2^8
\end{eqnarray}
次に、$13$の自乗を繰り返し計算します。
\begin{eqnarray}
\left(13^{2^0}\equiv \right)13 \equiv 13 \bmod 31\\
\left(13^{2^1}\equiv\right)13^2 \equiv 14 \bmod 31\\
\left(13^{2^2}\equiv\right)14^2 \equiv 10 \bmod 31\\
\left(13^{2^3}\equiv\right)10^2 \equiv 7 \bmod 31\\
\left(13^{2^4}\equiv\right)7^2 \equiv 18 \bmod 31\\
\left(13^{2^5}\equiv\right)18^2 \equiv 14 \bmod 31\\
\left(13^{2^6}\equiv\right)14^2 \equiv 10 \bmod 31\\
\left(13^{2^7}\equiv\right)10^2 \equiv 7 \bmod 31\\
\left(13^{2^8}\equiv\right)7^2 \equiv 18 \bmod 31\\
\end{eqnarray}
最後に、式(11)の指数に合わせて積をとると、
\begin{eqnarray}
13^{400} &=& 13^{0\times 2^0+0\times 2^1+0\times 2^2+0\times 2^3+1\times 2^4+0\times 2^5+0\times 2^6+1\times 2^7+1\times 2^8}\\
&\equiv & 13^0\cdot 14^0\cdot 10^0\cdot 7^0\cdot 18^1 \cdot 14^0 \cdot 10^0 \cdot 7^1 \cdot 18^1 \bmod 31\\
&\equiv & 18^1\cdot 7^1 \cdot 18^1 \bmod 31\\
&\equiv & 5 \bmod 31
\end{eqnarray}
となって計算できました。
どれぐらい簡単になったか、積の回数で比較しましょう。
普通に13を400回かけるとすると、積の回数は399回ですね。
一方、繰り返し自乗法では、13の自乗を繰り返し計算するフェーズ2.で8回、それらの積を計算するフェーズで2回の積、合計10回で済みます。
ずいぶんと計算が簡略化されましたね！
一般的に考えると、繰り返し自乗法では、指数$k$の二進展開の桁数を$d$すなわち、$\log_2 k = d$としたとき、その積の回数は$\mathcal{O}(d)$となることが知られています。べき乗をただ計算する計算法のオーダーが$\mathcal{O}(2^d)$となることを考えると、非常に効率的なアルゴリズムであることがわかるでしょう。
アルゴリズムの実装は多少難しいように思えるかもしれませんが、実は結構簡単です。
以下はJavaでの実装ですが、次のように、非常に簡潔に書くことが出来ます。
public static BigInteger power(BigInteger x, BigInteger k, BigInteger m){
BigInteger v = BigInteger.ONE;
if(k.equals(BigInteger.ZERO))
return v;
BigInteger p = x.remainder(m);
BigInteger two = BigInteger.ONE.add(BigInteger.ONE);
while(k.compareTo(BigInteger.ZERO) > 0){
if(k.remainder(two).equals(BigInteger.ONE))
v = (v.multiply(p)).remainder(m);
p = (p.multiply(p)).remainder(m);
k = k.shiftRight(1);
}
return v;
}
参考文献
[1] ジョセフ・H・シルヴァーマン著、鈴木治郎訳、はじめての数論、ピアソン・エデュケーション
Posted in 数学