べき乗とは
$$ z=x^y $$
義務教育で習うと思いますが、xのy乗を計算することです。
公開鍵暗号暗号では重要な処理になっています。
（正確に言うと、公開鍵暗号では 有限体の演算 であるため、$z=x^y\bmod{p}$ のように$\bmod{p}$ が必要です。）
例
$ x=2, y=5$ という値ならば、
$$ x^y = 2^5$$
$$ 2\times 2=4, \ 4 \times 2=8, \ 8 \times 2=16, \ 16\times 2=32 $$
と計算すれば良いですが、数が大きくなると太刀打ちできません。
$ y=11957708941720303968251$ の場合、一人の人生をかけるほどの時間が掛かりそうです。
べき乗算アルゴリズムと公開鍵暗号
公開鍵暗号のように256bitや2048bitなどの大きな数の場合、宇宙が誕生してから現在までの時間をかけても計算することができません。
そこで、高速にべき乗を計算するアルゴリズムが提案されており、現代ではほとんどすべてのコンピュータで実行されています。
DH鍵共有では$g^a$と$g^b$というような2つのそれぞれの秘密鍵を利用して鍵を共有します。
最近良く利用される楕円曲線暗号では、楕円曲線上の点のスカラー倍$Q=kP$という処理にべき乗算と同じアルゴリズムを利用します。
つまり、現代の情報社会の安全性を支える公開鍵暗号はべき乗算アルゴリズムが重要なのです！！！😃
高速べき乗アルゴリズム
上記の一回ずつ計算する方法を直感的な計算方法とすると、
$y$が256bitならば$2^{256} \approx 1.16\times 10^{77}$です。
宇宙に存在する素粒子の数が約$10^{80}$個らしいので、いかに巨大な数値であるかわかるでしょう。
これでは非常に遅いので効率的な方法が存在します。
バイナリ法
べき乗の値$y$を2進数で表記し、1ビットごとにチェックします。
初期値として$z=1$を設定します。
チェックしたビットの値で計算を変更するため、バイナリ法と呼ばれています。
0ならば$x=x^2$を計算1ならば$z=zx$と$x=x^2$
というようにべき乗の値のビットを手がかりにべき乗を計算します。
例
$y=11957708941720303968251$は2進数にすると
10100010000011101010001100000111111101011001011110111000000001111111111011
になります。
計算量
これは74bitで、74回の2乗算と40回の乗算($74S+40M$)で計算可能です。
もう少し発展させて、256bitのべき剰余算を行う場合、256回の2乗算と平均して128回の乗算で計算できます。
C/C++プログラム(GMP利用)
#include "gmpxx.h"
void exp_binary(mpz_class &z, const mpz_class& x, const mpz_class& y, const mpz_class& mod)
{
z = 1;
size_t bit_length = mpz_sizeinbase(y.get_mpz_t(), 2)
for (int i = bit_length - 1; i >= 0; --i) {
z*=z;
z%=mod;
if (mpz_tstbit(y.get_mpz_t(), i) == 1)
z *= x;
z%=mod
}
}
Pythonプログラム
def binary_exp(x, y, mod):
z=1
bin_exp = bin(y)[2:][::-1]
bit_length = len(bin_exp)
for i in range(bit_length - 1, -1, -1):
z = (z * z)%mod
if bin_exp[i] == '1':
z = (z * x)%mod
return z
Sliding window法
バイナリ法を改良したものとして、高速な暗号ライブラリで用いられています。
バイナリ法ではすべてのビットを一つずつチェックところ、この手法では複数ビットをひとまとめに見て、計算します。
ひとまとめにするビット幅を窓幅$k$とします。
窓幅に応じて事前計算を行い、その値を計算で利用します。
これは事前計算で計算コストが増えますが、2048bitなどの巨大な数値でバイナリ法と比較して効果を発揮します。
事前計算は次のような$x$の奇数乗を$k$の値にしたがって求めます。 $$x^3,\ x^5, \ … \ x^{2^k-1}$$
例えば$x=2, \ k=3$だとテーブルは $$[ 2, 8, 32, 128 ]$$ となりますが、 下記のプログラムでは次のようにしています。 $$[ 0, 2, 0, 8, 0, 32, 0, 128 ]$$
これはテーブル参照の手間を削減するためであり、実用的な暗号ライブラリでも利用されています。
メモリがもったいないという方はテーブル参照の処理を書いて配列数を減らすと良いでしょう()。
$y$の値のビットをチェックし、0ならば2乗算を実行します。
1ならば、窓幅$k$をチェックします。
$k=3$のとき、ウィンドウとして観測できるものとしては、
111101111
という奇数が出てきます。
観測した窓幅の数回の2乗算を実行後、ウィンドウの値からテーブルを参照して乗算を実行します。
例えば、11が観測されたら、table[3]を読み出し、$z=z\times 8$を計算します。
処理した窓幅分のビットを進めて最後まで処理を繰り返します。
例
先程の$y=11957708941720303968251$を使います。
窓幅$k=3$とすると、次のようになります1。
$$ \underset{5}{\underline{101}}000\underset{1}{\underline{1}} 00000 \underset{7}{\underline{111}}0\underset{5}{\underline{101}}000\underset{3}{\underline{11}}00000 \underset{7}{\underline{111}} \ \underset{7}{\underline{111}} \ \underset{5}{\underline{101}}0\underset{3}{\underline{11}}00\underset{5}{\underline{101}} \ \underset{7}{\underline{111}}0 \underset{7}{\underline{111}}00000000\underset{7}{\underline{111}} \ \underset{7}{\underline{111}} \ \underset{7}{\underline{111}} \ \underset{5}{\underline{101}} \ \underset{1}{\underline{1}} $$
左側の$101$から処理を行います。
ビット幅は3なので3回の2乗算を実行し、table[5]の値を乗算します。
次のビットは0ですので、2乗算のみ実行します。
というような処理を最後まで行います。
計算量
テーブル生成の処理を含めて75回の2乗算と20回の乗算($75S+20M$)で計算可能です。
コンピュータでは乗算を筆算のようなアルゴリズムで計算します。 そのため、2乗算は乗算よりも乗算の回数を減らすことができます。
例えば、$12 \times 15$という計算は $2\times 5$、$1 \times 5$、$1 \times 2, 1 \times 1$ を計算し、平行四辺形を加算することで計算できます。 このとき、掛け算は4回必要です。
2乗算の例$15\times 15$では、 $5 \times 5, 1 \times 5$、$5 \times 1, 1 \times 1$ を計算します。 ここで、$1 \times 5$と$5 \times 1$ は結果が同じなので、乗算の回数を減らすことができます。
コストを比較する際、$1S=0.8M$ とすることがあり、 その場合の計算コストは$75S+20M=80.0$ となります。
バイナリ法の場合は$74S+40M=97.6$ですので、 Sliding window法を利用することで計算量を減らすことができます。
Cプログラム
#include <iostream>
#include <gmpxx.h>
// Window幅
const size_t k = 3;
int max(const int a, const int b) {
if (a >= b) {
return a;
} else {
return b;
}
}
void exp_sliding_window(mpz_class& z, const mpz_class& x, const mpz_class& y, const mpz_class& mod) {
// テーブルの生成
mpz_class table[1 << k];
mpz_class aa;
for (int i = 0; i < (1 << k); ++i) {
table[i] = 0;
}
table[1] = x;
aa = (x * x) % mod;
for (int i = 1; i < (1 << (k - 1)); ++i) {
table[2 * i + 1] = (table[2 * i - 1] * aa) % mod;
}
// exponentの2進数表現を取得
std::string bin_exp = y.get_str(2);
std::reverse(bin_exp.begin(), bin_exp.end());
size_t bit_length = mpz_sizeinbase(y.get_mpz_t(), 2);
z = 1;
int i = bit_length - 1;
while (i >= 0) {
if (bin_exp[i] == '0') {
z = (z * z) % mod;
i -= 1;
} else {
int s = max(i - k + 1, 0);
while (bin_exp[s] == '0') {
s += 1;
}
for (int j = 0; j < i - s + 1; ++j) {
z = (z * z) % mod;
}
std::string window_bits = bin_exp.substr(s, i - s + 1);
std::reverse(window_bits.begin(), window_bits.end());
mpz_class window = mpz_class(window_bits, 2);
z = (z * table[window.get_ui()]) % mod;
i = s - 1;
}
}
}
Pythonプログラムの例
def exp_sliding_window(x, y, mod):
# Window幅
k = 3
#テーブルの生成
table = [0] * (1 << k)
table[1]=x
aa=x*x % mod
for i in range(1, 1<<(k-1)):
table[2*i+1] = (table[2*i-1] * aa)%mod
print(table)
bin_exp = bin(y)[2:][::-1]
length = len(bin_exp)
print(bin_exp)
result = 1
i = length-1
while i >= 0:
if bin_exp[i] == '0':
result = (result * result)%mod
i -= 1
else:
s = max(i-k+1, 0) # s is the start of the window
while bin_exp[s] == '0':
s += 1
for _ in range(i-s+1):
result = (result * result)%mod
window_bits = bin_exp[s:i+1][::-1]
window = int(window_bits, 2)
result = (result * table[window])%mod
i = s - 1
return result
window法
Sliding window法よりも、テーブルの生成や乗算の回数が増えますが、奇数ビット列以外のウィンドウも作っておく方法があります。
テーブル参照が面倒な場合には良いかもしれませんが、ウィンドウサイズで割り切れるようにするために左に'0’を追加する必要があります。
def window(base, exponent, mod):
k = 3
table = [0] * (1 << k)
table[0]=1
table[1]=base
# 2^k-1までループ
for i in range(1, (1<<k)-1):
table[i+1] = (table[i] * base)%mod
mul+=1
bin_exp = bin(exponent)[2:][::-1]
length = len(bin_exp)
result = 1
i = length-1
# kで割り切れるように'0'を追加
while (i+1) % k != 0:
bin_exp += '0'
i += 1
# kビットずつ計算
while i >= 0:
s = i-k+1
for _ in range(k):
result = (result * result)%mod
window_bits = bin_exp[s:i+1][::-1]
window = int(window_bits, 2)
result = (result * table[window])%mod
i = s - 1
return result
※広告
素数の本ですが多倍長整数の計算アルゴリズムなども記載されており、辞書代わりによく読む本です。
数値を入力するとLaTeXで出力してくれるPythonコード( https://tweet.mimoex.net/m/58 ) ↩︎