勉強したメモ。
大きい素数でちゃんと動くRSAの実装の記事があんまりなかったので。
目的
256bitRSA暗号を最小構成で実装する。
RSAのやりかた
詳細はwikipediaでも見て頂いて、流れだけ書きます。
- bit数を決める。(今回は256bit)
- 大きな素数$p,q$を作る。(128bit。256の半分)
- $n=pq$
- $(p-1)(q-1)$より小さく、互いに素な数$e$を決める。
- $ed\equiv 1 \pmod n$である$d$を求める。
これで準備ができました。暗号化は、
- 平文$m$を作る($m<n$)
- 暗号文$c$は、$m^e \pmod n$
復号は、$m=c^d \pmod n$でOKです。
そんなに難しい計算はないですが、ポイントは
- べき乗計算が重い
- 大きな素数$p,q$をどうやって作る?
- $e$をどうやって決める?
- $d$は?
べき乗計算が重い
RSAでは数百桁の数でべき乗するので、素直に指数回掛け算などやってられません。
高速指数演算というのを使うと速くなります。
さらに、求めたいのは冪の剰余であって冪そのものではないので、この点でも最適化できます。
$a$の$b$乗を$n$で割ったあまりを計算するとしましょう。まず$b$を2進数にします。2進数の各桁を$b_i$として、
$$
a^b = a^{\sum 2^{i}b_i}=\prod a^{2^ib_i}
$$
となります。
ここで、
$$
a^{2^{i+1}}=a^{2^i}*a^{2^i}
$$
なので、掛け算一回で次の桁を計算できます。
さらに、計算したいのはべき乗の値を$n$で割ったあまりなので、掛け算の各ステップで随時あまりを取ってあげれば、$n$より小さい値を扱うだけで済みます。
例として、$a=3,b=11,n=21$とすると、まず$b=(1011)_2$つまり、
$$
b=2^3+2^1+2^0
$$
として、
$$
a^b=3^{2^3+2^1+2^0}=3^{2^3}*3^{2^1}*3^{2^0}
$$
$$
3^{2^0}=3
$$
$$
3^{2^1}=3^{2^0}*3^{2^0}=9
$$
$$
3^{2^2}=3^{2^1}*3^{2^1}=81
$$
$$
3^{2^3}=3^{2^2}*3^{2^2}=6561
$$
それぞれ$21$で割ったあまりにして、計算中も21を超えたらあまりを取ります。
$$
396561\equiv (3*9)9\equiv 69\equiv 12
$$
def modular_exp(a, b, n):
res = 1
while b != 0:
if b & 1 != 0:
res = (res * a) % n
a = (a * a) % n
b = b >> 1
return res
大きな素数をどうやって作るか
どうやら直接素数を生成するのは結構難しいようです。
実用的には、運良く素数に当たるまで乱数を生成するというのが簡単らしい。
乱数生成
bit数が指定されるので、その数だけ0,1をランダムに並べればOK。
def gen_rand(bit_length):
bits = [random.randint(0,1) for _ in range(bit_length - 2)]
ret = 1
for b in bits:
ret = ret * 2 + int(b)
return ret * 2 + 1
最大ビットと最小ビットはなからず1になるようにしました。
最大ビットは桁数をが減らないように、最小ビットは奇数にするためです。
なぜ奇数かと言えば、今は素数がほしくて、偶数は素数ではないからです。
素数判定
巨大な乱数$p$が素数かどうか判定したい。
$p$はとても大きいので、順に割っていくような方法ではおそすぎる。
**ミラー–ラビン素数判定法**というのがよく使われるらしい。アルゴリズムはwikipediaにも乗ってるのでその通りに実装。
簡単に言えば、素数の必要条件を何回もチェックして、全部通ったらたぶん素数、一つでも通らなかったら合成数だと判定してます。
これは確率的素数判定のアルゴリズムなので、合成数を素数と判定してしまうことがあります。
ただ、その誤検知は1回のチェックで$\frac{1}{4}$程度だそうなので、何十回かやれば実用上ほぼ0になるので問題ないということらしいです。
def mr_primary_test(n, k=100):
if n == 1:
return False
if n == 2:
return True
if n % 2 == 0:
return False
d = n - 1
s = 0
while d % 2 != 0:
d /= 2
s += 1
r = [random.randint(1, n - 1) for _ in range(k)]
for a in r:
if modular_exp(a, d, n) != 1:
pl = [(2 ** rr) * d for rr in range(s)]
flg = True
for p in pl:
if modular_exp(a, p, n) == 1:
flg = False
break
if flg:
return False
return True
ということで、素数生成はこんな感じです。
def gen_prime(bit):
while True:
ret = gen_rand(bit)
if mr_primary_test(ret):
break
return ret
eをどうやって決める?
$(p-1)(q-1)$と互いに素ならなんでもいいので、素数にするのが手早いです。
p = gen_prime(128)
q = gen_prime(128)
e = gen_prime(128)
とかで問題ないのではないかと思います。
wikipediaには$65537$がよく使われると書いてました。決め打ち定数でもいいみたいですね。
dは?
拡張ユークリッドの互除法というやつが使えます。
もともと自然数$a,b$に対して、
$$
ax+by=gcd(a,b)
$$
となるような$x,y$を求めるアルゴリズムです。(gcdは最大公約数)
これを、$e,(p-1)(q-1)$に対して使います。最大公約数は$e$の定義から1なので、
$$
ex+(p-1)(q-1)y=1
$$
両辺を$(p-1)(q-1)$で割ると、
$$
ex\equiv 1
$$
となるので、この$x$が$d$です。
拡張ユークリッドの互除法はここにあったのでコピペしました。$x_0$が$d$です。
def xgcd(b, n):
x0, x1, y0, y1 = 1, 0, 0, 1
while n != 0:
q, b, n = b // n, n, b % n
x0, x1 = x1, x0 - q * x1
y0, y1 = y1, y0 - q * y1
return b, x0, y0
def gen_d(e, l):
_, x, _ = xgcd(e, l)
return x % l
実行
あとは平文を用意すれば動きます。
bit_length = 128
p = gen_prime(bit_length)
q = gen_prime(bit_length)
e = gen_prime(bit_length)
d = gen_d(e, (p - 1) * (q - 1))
n = p * q
m = 123456789
c = modular_exp(m, e, n) # 暗号文
m_ = modular_exp(c, d, n) # 123456789
1024bitでも一応動きました。(20秒くらいかかった)
実際は、同じ平文を何度も暗号化するのは危険なので、末尾に乱数を追加したり(復号時に削除)、mが小さすぎたらパディングしたりするそうです。