RSAのフェルマー攻撃
訳あって、春頃に話題になっていたRSAのフェルマー攻撃1に関して調べていました。フェルマー法2については各所で解説があるのでここでは詳しくは触れませんが、素数 p, q の積がN であるときに、Nから p, q の値を求める計算手順（アルゴリズム）です。RSAの場合、Nは公開鍵に、p, qは秘密鍵に含まれる値であり、Nからp, q が計算で求められるということは公開鍵から秘密鍵がバレることに等しいわけで、そんな方法があるのかと話題になったわけです。
脆弱なｐ、ｑの値について
ただ、フェルマー法による計算が計算機で現実的な時間内に終わるためには条件があり、「p, q が近い値である」必要があります。p, q の値が近いということは、その積であるNの平方根の値とp, qの距離も近いということが言えるので、計算機によるフェルマー法でわりと簡単に見つけることができます。
では、どれくらい「近い」とダメなのか。上記の脆弱性を発見された方のブログでは、pとqが最大2517まで異なる数を確実に因数分解したとしており、下位64バイトしか差がない素数は脆弱、としています。これはちょっと言ってることが分かりにくいのですが、鍵長2048ビットのRSA鍵の場合、p, q の値はそれぞれ1024ビット(128バイト)になります。その半分である上位512ビット(64バイト)が一致する場合は脆弱だと言っています。「2517 の差までは因数分解した」とのことなので正確には507ビット以上が一致しても脆弱ということです。
攻撃のラウンド数
この説明での仮定は「100ラウンドのフェルマーアルゴリズム」ということなので、ラウンド数が増えればもっと差が大きいp, q も計算できます。ただし、2514の差までは、大抵最初のラウンドで因数分解できるということで、100ラウンドは過剰とのこと。このあたりは直観的には分かりにくく実際にやってみないとどういうことなのかピンとこない人も多いかもしれません。
そこで、もう少し分かり易くするために実験をしてみました。フェルマー法を試した実験はQiitaの記事3でpythonコードが公開されていたのでこちらを拝借しつつも、このコードはp, q の差を2バイト程度しか作らないので、差を可変で下位Xビット分にして、p, qを発見するまでの時間を計測するコードに改造し、平均的な処理時間を計測してみました。（記事最後に掲載）
環境はAMD Ryzen7 2700X 3.7GHzです。
以下は鍵長2048ビット、p, qが1024ビットの場合の結果です。
| X bits | round | elapsed(s) |
|---|---|---|
| 128 | 1 | 0.000256 |
| 256 | 1 | 0.000252 |
| 384 | 1 | 0.000395 |
| 512 | 1 | 0.000737 |
| 513 | 1 | 0.000926 |
| 514 | 2 | 0.001262 |
| 516 | 41 | 0.021497 |
| 518 | 1074 | 0.534091 |
| 519 | 2303 | 1.162753 |
| 520 | 18304 | 9.241641 |
| 521 | 37168 | 18.735747 |
| 522 | 183703 | 94.061339 |
513ビットまではラウンド数は1回であり処理時間も1秒もかかりません。ところが514ビットを越えたあたりからラウンド数が増え始めます。516ビット以降はラウンド数も処理時間も顕著な増え方をしています。ラウンド数も30回程度の試行の中から平均で取っているため、513ビットにはわずかにラウンド数2回が含まれていますし、514ビットにはラウンド数１回もかなり含まれています。「2514までは大抵最初のラウンドで因数分解できる」というのは本当のようです。
では鍵長4096ビット、p, q が2048ビットの場合はどうなるのでしょうか？こちらは素数発見に時間がかかるため範囲と試行回数を減らして試しました。ですので、あまりあてになりませんが傾向はなんとなく分かります。
| X bits | round | elapsed(s) |
|---|---|---|
| 512 | 1 | 0.000716 |
| 768 | 1 | 0.001538 |
| 1024 | 1 | 0.002580 |
| 1025 | 1 | 0.003465 |
| 1026 | 4 | 0.010189 |
| 1027 | 3 | 0.009397 |
| 1028 | 7 | 0.018712 |
今度は1025ビットまではラウンド数1回で発見できています。つまり、鍵長の1/4＋1ビットまではほぼ1ラウンドで発見できるが、そこから先はラウンド数が増加し始めるという傾向が分かります。
「差が半分以下」の意味するところ
さて、半分を境にたかが数ビットの違いで発見までのラウンド数が全く違ってくるというのはあまり直観的ではありませんし、鍵長が増えるとラウンド数の増加する差分も変わるというのはなんとなく不思議です。もう少し分かり易くするために小さい値で考えてみます。
まず、全体の半分の上位ビットが一致する8ビットのp, q を考えてみます。まずpですが、8ビットで取り得る最大の素数として 251 = 1111 1011 を採用します。ここでp=251との距離を最大に取る素数をqとするなら 241 = 1111 0001 が素数として採用できます。あれ？もっと差を拡げられないかな？と思っても、これより小さい素数は上位4ビットが変化してしまうので条件を満たしません。しかも241の次の素数は251となりpと同じです。つまりqの選択肢は241の１つしかありません。
$$ p=251、q=241のケースで、Nと\sqrt{N}を計算すると\\
N=251\times 241=60491, \sqrt{N}=\sqrt{60491}=246\\
フェルマー法でN=(a-b)(a+b), p=a+b, q=a-b と書けるので\\
a=\sqrt{N}=246だとすると、b=5 が1回目で計算できる $$
では、上位２ビットが一致するp, q の場合はどうでしょうか？同じp = 251 = 1111 1011 として上位2ビットが一致する素数は 193 = 1100 0001、197 = 1100 0101、199 = 1100 0111、211 = 1101 0011、223 = 1101 1111、227 = 1110 0011、229 = 1110 0101、233 = 1110 1001、239 = 1110 1111、241 = 1111 0001 と選択肢がかなり広がりました。
これは当たり前の結果で、一致する上位ビットが少ないことは採用できる数値の幅が広くなるので、選択可能な素数が増えることは理に適っています。ただ、上位4ビット(半分)の一致と2ビットの一致の差（数値の幅としては4倍程度）が随分大きいようにも感じます。どうも素数の分布には「偏り」があるようです。では素数は実際どういう分布になっているのでしょうか。
素数定理
ある任意の桁を上限とする自然数の中にいくつ素数が含まれているのか、その割合や分布についてはかなり昔から研究されてきました。素数の分布は不規則かつ複雑で未知であると言われていて、確定的な予測はできないと言われていますが、ガウスの予想した近似式である「素数定理」は証明されている定理です。
素数定理をざっくりと表現すると、「素数は大きい数になるほどまばらに分布している」ということになります。
x以下の素数の個数を \(\pi(x)\) とした場合に \(\pi(x)\sim \frac{x}{\log x}\) が成立する、というのが素数定理ですが、「\(\sim\)」というのは両辺の比の極限が１に収束することを意味しています。つまり、xが大きくなるほどほぼ等しくなるよということです。\(\log x\) は自然対数であり、素数の数の割合は対数関数の逆数に近似した分布に近くなるということです。自然数x以下に存在する素数の数の割合\(p(x)\)は \(\frac{1}{\log x}\) に近似する、とも言えます。
素数定理によれば1024ビットの自然数と2048ビットの自然数の中に分布する素数の密度を比べると1024ビット（\(p(x) ≒ 1 / 709\)）より2048ビット（\(p(x) ≒ 1 / 1419\)）の方が疎である、ということも分かります。pythonによる実験で鍵長を増やしても傾向が変わらなかったのは、鍵長が増えた分だけ素数の密度も低くなっているからということが言えそうです。
桁が大きい数値ほど素数の密度が疎になるということは、ある特定の自然数以下における素数の分布もグラデーションになっているであろうことが分かります。つまり素数定理の密度は全体に万遍なく適用できる密度ではなく、小さい桁の方に偏っていることが分かります。その傾向は小さい桁の数値で見ても明らかです。
0000 0001 … 0000 1111 の中に存在する素数：６個
1111 0001 … 1111 1111 の中に存在する素数：２個
となると最上位ビットが1で、かつその上位半分のビットが一定の数値の範囲に出現する素数は素数定理の密度よりもかなり少なくなくなるであろうことは容易に予想できます。素数定理においてはそのような範囲のxにおける素数の数\(\pi(x)\)には差がほぼ生じません。定理における素数の密度は対数関数の逆数に沿って変化するので、1024ビットの下位半分のビット、2512分の数値が変化しても、密度は誤差範囲の差しか生じないのです。十分に巨大な自然数において上位半分のビットが一致する範囲における素数の数はとても少ない（稀）であろうということ、自然数が大きければ大きいほどその傾向も強くなることが予想できます。
もちろん素数定理は近似式なのでその範囲に素数が存在しないわけではないですが、実際希薄です。そのため、この範囲に存在する素数をp, q に選んでも大抵はフェルマー法の1ラウンド目であっという間に発見されてしまうということになります。
素数定理についての詳細はネット上に沢山解説があるので興味のある人は是非調べてみてください。
実験コード
import random as rnd
from re import X
import sympy
import time
key_length = 1024
# key_length = 2048
distance = 100000000
## 素数生成
def gen_prime(length):
p = 0
q = 0
# 乱数xを生成する
x = rnd.randrange(2, pow(2, key_length))
# lengthビットの反転ビット
LB = 2 ** length -1
# xの次の素数を探索し、それをpとする
for i in range(distance):
if sympy.isprime(x+i):
if p == 0:
p = x + i
x = p
# xの下位lengthビットを反転
x ^= LB
break
# pの下位lengthビットを反転した値から素数を探索し、それをqとする
for i in range(distance):
if sympy.isprime(x+i):
if q == 0:
q = x + i
break
N = p * q
if p > 0 and q > 0:
print('p:', p)
print('q:', q)
print('N:', N)
else:
print('Not found')
quit()
return p, q, N
## 平方根
def isqrt(n):
x = n
y = (x + 1) // 2
while y < x:
x = y
y = (x + n // x) // 2
return x
## フェルマー法
def fermat_method(p, q, N):
# 素数が近ければ、aはNの平方根に近い
a = isqrt(N)
b = 0
time_start = time.perf_counter()
for i in range(distance):
# このため、Nの平方根から初めて、一周するごとに1ずつ推測を増やしていけば、aの値を推測することができる
a = a + 1
# それぞれの推測にあ値して、b^2=a^2-Nを計算することができる
b2 = a*a - N
t = isqrt(b2)
# その結果が平方数であれば、aを正しく推測したことが分かる
if t*t == b2:
b = t
round = i + 1
print('found it (round=%d)' % round)
break
time_end = time.perf_counter()
elapsed = time_end - time_start
if b > 0:
print('a:', a)
print('b:', a)
# ここから、p=a-b と q=a+b を計算することができる
pp = a - b
qq = a + b
if p == pp:
print('pp:', pp)
if q == qq:
print('qq:', qq)
NN = pp * qq
print('NN:', NN)
print('Elapsed time:', elapsed)
else:
print('Not found')
return elapsed, round
## main
def main():
# 結果
elapsed_avg = {}
round_avg = {}
# 試行回数
trials = 30
# p, q の差を作る下位ビット数
dif_length = [128, 256, 384, 512, 514, 516, 518, 519, 520, 521, 522]
#dif_length = [512, 768, 1024, 1025, 1026, 1027, 1028]
for bits in dif_length:
trial_elapsed = []
trial_round = []
print()
print('===========================')
print('****** Different LSB : ', bits)
for tr_no in range(0, trials):
print('--------------------------')
print('trial no : ', tr_no+1)
p, q, n = gen_prime(bits)
elapsed , round = fermat_method(p, q, n)
trial_elapsed.append(elapsed)
trial_round.append(round)
elapsed_avg[str(bits)] = sum(trial_elapsed) / len(trial_elapsed)
round_avg[str(bits)] = sum(trial_round) / len(trial_round)
print()
print('===========================')
print('****** result')
print('bits : round : elapsed')
print('-------------------------')
for bits in dif_length:
print(str(bits)+' : %8d : %3.6f' % (round_avg[str(bits)],elapsed_avg[str(bits)]))
print('-------------------------')
print(round_avg)
print(elapsed_avg)
if __name__ == "__main__":
main()
コメント