楕円曲線上の離散対数問題 (ECDLP) への攻撃手法と SageMath による実装のまとめです。 前回の離散対数問題の攻撃手法の記事では「有限体上の離散対数問題 (DLP)」への攻撃手法を扱いましたが、ここでは楕円曲線上の離散対数問題について説明します。
楕円曲線上の離散対数問題 (ECDLP)
有限体 $F_p$ 上の楕円曲線 $E$ について、$E$ 上のある点 $G$ と整数 $x$ から、点 $Y = xG$ (スカラー倍算) を満たす整数 $x$ を探す問題のことを、楕円曲線上の離散対数問題といいます。 英語では ECDLP (Elliptic Curve Discrete Logarithm Problem) といいます。
ここでは次の3つの楕円曲線上の離散対数問題への攻撃手法について説明します。
- Baby-step Giant-step法
- Pollard's rho法 (ポラード・ロー)
- Pohlig–Hellman法 (ポーリッヒ・ヘルマン)
Baby-step Giant-step法
楕円曲線 $E$ 上の離散対数問題 $Y = xG$ の $x$ を求める方法の1つである Baby-step Giant-step法 (BSGS) は数え上げ法よりも少ない回数の演算でできますが、多くの格納領域を必要とします。
まず、楕円曲線 $E$ の位数を \(\#E\)、\(m = \lceil \sqrt{\#E - 1}\rceil\) とおき、$x$ を $m$ で割ったときの商 $q$ と余り $r$ の式を作ります。
\[x = qm + r \;\;\;(0 \le r < m)\]この $q, r$ をBSGSで求めることで、離散対数問題の $Y = xG$ を満たす最小の $x$ が求まります。 まず、以下のように式変形をします。
\[\begin{aligned} Y &= xG \\ Y &= (qm+r)G \\ Y - rG &= qmG \end{aligned}\]次にBaby-stepの処理として、左辺の計算を事前にし、集合$B$を求めます。なお $r$ は余りなので $m$ 未満の整数となります。
\[B = \{ (Y - rG, r) \;\vert\; 0 \le r < m \}\]もし、集合 $B$ の中に $(1,r)$ があれば、$Y - rG = 1$ より $Y = rG$ で $x = r$ なので離散対数問題が解かれます。 それ以外のときは、Giant-stepの処理として、右辺の計算をします。$q = 1,2,…,m-1$ に対して $qmG$ が集合 $B$ の中に含まれるかを調べ、集合に含まれるときは $q$ と $B$ に保存した $r$ から $x = qm + r$ で離散対数問題が解かれます。
実装は以下の通りです（SageMath 9.0, Python 3.8）
def babystep_giantstep(G, Y, E):
m = int((E.order()-1)**0.5 + 0.5)
# Baby step
table = {}
YrG = Y # Y-r*G
for r in range(m):
table[YrG] = r
YrG -= G
# Giant step
mG = m * G # m*G
qmG = mG # qm*G
for q in range(1, m):
if qmG in table: # 左辺と右辺が一致するとき
return q * m + table[qmG]
qmG += mG
return None
p = 240556067
F = GF(p)
E = EllipticCurve(F, [0, 486662, 0, 1, 0]) # Curve25519
G = E(103666880, 133544401)
Y = E(220898463, 208070124)
x = babystep_giantstep(G, Y, E)
print(x)
print(Y == x*G)
なお、SageMath には discrete_log という離散対数問題を解く関数があるので、実装したくない人
はこれを使うのが一番手っ取り早いです。
内部的には Baby-step Giant-step法とPohlig–Hellman法を使っているようです。
x = G.discrete_log(Y)
print(x)
Pollard's rho法
楕円曲線 $E$ 上の離散対数問題 $Y = xG$ の $x$ を求める方法の1つである Pollard's rho法（ポラードのローアルゴリズム）はBSGSと同じ計算量ですが、使用するメモリ空間は定数量のみです。
まず、関数 $f$ を次のように定義します。
\[f(X) = \begin{cases} Y+X & (X \in P_1 のとき) \\ X+X & (X \in P_2 のとき) \\ G+X & (X \in P_3 のとき) \end{cases}\]次に、数列 $X_i$ について考えます。乱数 $a_0$ を選び、$X_0 = a_0 G$ とするときの数列 $X_i$ を漸化式 $X_{i+1} = f(X_i)$ で計算すると、この数列の項は $X_i = a_0 G + b_0 Y$ となります。 ここで、$b_0 = 0$、$p$ を楕円曲線 $E$ の位数とすると、$a_{i+1}, b_{i+1}$ は関数 $f$ の定義から次のようになります。
\[a_{i+1} = \begin{cases} a_i & (X_i \in P_1 のとき) \\ 2a \pmod{p} & (X_i \in P_2 のとき) \\ a_i + 1 \pmod{p} & (X_i \in P_3 のとき) \end{cases}\] \[b_{i+1} = \begin{cases} b_i + 1 \pmod{p} & (X_i \in P_1 のとき) \\ 2b \pmod{p} & (X_i \in P_2 のとき) \\ b_i & (X_i \in P_3 のとき) \end{cases}\]ポラード・ローのアルゴリズムでは、$Y = xG$ を求めるために、まず $a_iG + b_iY = A_i G + B_i Y$ を満たす $(a_i, b_i), (A_i, B_i)$ を求めます。$G$ は巡回群の生成元なので、$X_{i+k} = X_i$ となる $i, k$ が存在します。つまり、
\[\begin{aligned} a_i G + b_i Y &\equiv a_{i+k}G + b_{i+k}Y \\ (a_i - a_{i+k})G &\equiv (b_{i+k} - b_i)Y \\ (a_i - a_{i+k})G &\equiv (b_{i+k} - b_i)xG \\ (a_i - a_{i+k}) &\equiv x(b_{i+k} - b_i) \\ \end{aligned}\]が成立します。よって、次の式から離散対数 $x$ を求めることができます。
\[x \equiv (a_i - a_{i+k}) (b_{i+k} - b_i)^{-1} \pmod{p}\]${a_i}G + {b_i}Y = {A_i}G + {B_i}Y = {a_{i+k}}G + {b_{i+k}}Y$ を満たす $(a_i, b_i), (a_{i+k}, b_{i+k})$ の探し方は、巡回群の数列 $x_i$ は周期的であるので、フロイドの循環検出法（ウサギとカメのアルゴリズム）を使って求めます。
有限体上の離散対数問題 (FFDLP) のときとは違う点として、楕円曲線の位数が素数ではなく合成数のときがあります。 位数が合成数の時は、まず位数を素因数分解して ($p = q_1 q_2 \cdots q_n$)、それぞれの素因数を位数としてPollard's rho法で離散対数を解き、最後に中国人剰余定理(CRT)で各式を全て満たす一つの解を求めます。 ただし、楕円曲線の位数が合成数のときは、次のPohlig–Hellman法が有効なので、ここでは素数位数の場合の実装のみ示します。 通常、暗号技術を安全に運用する時は素数位数となるように楕円曲線のパラメータを設定します。
実装は以下の通りです（SageMath 9.0, Python 3.8）
# Pollard's rho法
def pollards_rho(G, Y, E):
q = E.order()
print('order =', q)
e = E(0,1,0) # 単位元
def new_xab(x, a, b, g, y, q):
try:
subset = mod(x.xy()[0], 3)
except ZeroDivisionError:
subset = 2
if subset == 0:
return (x+x, (a*2) % q, (b*2) % q)
if subset == 1:
return (x+g, (a+1) % q, b )
if subset == 2:
return (x+y, a , (b+1) % q)
x, a, b = e, 0, 0
X, A, B = x, a, b
for i in range(1, p):
x, a, b = new_xab(x, a, b, G, Y, q)
X, A, B = new_xab(X, A, B, G, Y, q)
X, A, B = new_xab(X, A, B, G, Y, q)
if x == X:
break
res = ((a - A) * inverse_mod(B - b, q)) % q
if G * res == Y:
return res
return None
p = 47212873
E = EllipticCurve(GF(p), [3, 30])
G = E(47204974, 33464750)
Y = E(41264829, 15780002)
x = pollards_rho(G, Y, E)
print(x)
print(Y == x*G)
Pohlig–Hellman法
Pohlig–Hellman法（ポーリッヒ・ヘルマンのアルゴリズム）は楕円曲線 $E$ の位数 \(\#E\) が $p_1, p_2, …$ に因数分解できるとき、$E$ 上の離散対数 $Y = xG$ の計算問題を、素数位数 $p_i$ での離散対数問題に帰着させることで問題を小さくして解き、得られた複数の結果から、最後に中国人剰余定理で一つの答えを求める方法です。
まず、楕円曲線の位数 \(\#E\) が素数 $q_i$ に素因数分解できるとします。
\[\#E = q_1 q_2 \cdots q_k\]離散対数 $x$ は任意の自然数なので、商 $b_i$ と余り $a_i$ を使って $x = a_i + b_i q_i$ と書くことができます。 次に $Y = xG$ の両辺を \(\frac{\#E}{q_i}\) 倍すると、以下のように式変形できます。 なお、途中で生成元 $G$ を位数倍すると無限遠点 $O$ になることを使っています。
\[\begin{aligned} {\frac{\#E}{q_i}} Y &= \frac{\#E}{q_i} xG \\ &= \frac{\#E}{q_i} (a_i + b_i q_i) G \\ &= {\frac{\#E}{q_i}}a_i G + {b_i (\#E)}G \\ &= {\frac{\#E}{q_i}}a_i G \\ \end{aligned}\]見やすくするために、\(Y_i = {\frac{\#E}{q_i}}Y,\; G_i = {\frac{\#E}{q_i}}G\) とおくと、
\[Y_i = {a_i}G_i \pmod{p}\]となり、最初に示した離散対数問題の形になりますが、余り $a_i$ は割る数 $q_i$ よりも小さいという事実から、より高速に離散対数 $a_i$ を求めることができます。
そして \({q_i}G_i = {\frac{\#E}{q_i}}{q_i}G = 1\) より $G_i$ の位数は $q_i$ です。ここから、全ての素因数での離散対数 $a_i$ をまとめると次のようになります。
\[\begin{aligned} x &\equiv a_0 \pmod{q_0} \\ x &\equiv a_1 \pmod{q_1} \\ &\;\;\vdots \\ x &\equiv a_k \pmod{q_k} \\ \end{aligned}\]$q_0, …, q_k$ は互いに素なので、中国人剰余定理(CRT)を使って、これらを満たす離散対数 $x$ を求めることができます。
実装は以下の通りです（SageMath 9.0, Python 3.8）
def babystep_giantstep(G, Y, E):
m = int((E.order()-1)**0.5 + 0.5)
# Baby step
table = {}
YrG = Y # Y-r*G
for r in range(m):
table[YrG] = r
YrG -= G
# Giant step
mG = m * G # m*G
qmG = mG # qm*G
for q in range(1, m):
if qmG in table: # 左辺と右辺が一致するとき
return q * m + table[qmG]
qmG += mG
return None
# Pohlig–Hellman法
def pohlig_hellman_ECDLP(G, Y, E):
crt_moduli = []
crt_remain = []
for q, _ in factor(G.order()):
print('q =', q)
G0 = G * (E.order() // q)
Y0 = Y * (E.order() // q)
res = babystep_giantstep(G0, Y0, E)
if (res is None) or (res <= 1):
continue
crt_moduli.append(q)
crt_remain.append(res)
ans = crt(crt_remain, crt_moduli)
if ans is None:
return None
if G * ans == Y:
return ans
return None
p = 240556067
F = GF(p)
E = EllipticCurve(F, [0, 486662, 0, 1, 0]) # Curve25519
g = E(103666880, 133544401) # E.gen(0)
y = E(220898463, 208070124)
x = pohlig_hellman_ECDLP(g, y, E)
print(x)
print(y == x*g)
以上です。