離散対数問題に対する攻撃手法 (Python & SageMath)

離散対数問題 (DLP) への攻撃手法と Python & SageMath による実装のまとめです。暗号技術として、Diffie-Hellman鍵共有などの安全性は「離散対数問題」に依存しています。今年のセキュリティキャンプ2020の暗号解読ゼミでは、離散対数問題をテーマにしている方がいたので、話についていくために独学で勉強したときのまとめです。なお、楕円曲線上の離散対数問題について勉強したときのまとめは次の記事で説明しています。

離散対数問題 (DLP)

有限体 $F_p$ について、素数 $p$ と生成元 $g$ 、 $y \in \{1, ..., p-1\}$ が与えられたとき、 $y = g^x$ を満たす $x$ を探すことを離散対数問題といいます。英語では DLP (Discrete Logarithm Problem) と言ったり、有限体上を強調するために FFDLP (Finite Field DLP) と書いたりします。

ここでは次の4つの離散対数問題への攻撃手法について説明します。

Baby-step Giant-step法
Pollard's rho法 (ポラード・ロー)
Pohlig–Hellman法 (ポーリッヒ・ヘルマン)
指数計算法

Baby-step Giant-step法

離散対数問題 $y = g^x \mod{p}$ の $x$ を求める方法の1つである Baby-step Giant-step法 (ベイビーステップ・ジャイアントステップ; BSGS) は数え上げ法よりも少ない回数の演算でできますが、多くの格納領域を必要とします。

まず、 $m = \lceil \sqrt{p - 1}\rceil$ とおき、 $x$ を $m$ で割ったときの商 $q$ と余り $r$ の式を作ります。

x = qm + r \;\;\;(0 \le r < m)

この $q, r$ をBSGSで求めることで、離散対数問題の $y = g^x$ を満たす最小の $x$ を求まります。まず、次のように式変形をします。

\begin{aligned} y &= g^x \\ y &= g^{qm+r} \\ y g^{-r} &= (g^m)^q \end{aligned}

次に、Baby-stepの処理として、左辺の計算を事前にし、集合 $B$ を求めます。なお $r$ は余りなので $m$ 未満の整数となります。

B = \{ (y g^{-r}, r) \;\vert\; 0 \le r < m \}

もし、集合 $B$ の中に $(1,r)$ があれば、 $yg^{-r} = 1$ より $y = g^r$ で $x = r$ なので離散対数問題が解かれます。それ以外のときは、Giant-stepの処理として、右辺の計算をします。 $q = 1,2,…,m-1$ に対して $g^{qm}$ が集合 $B$ の中に含まれるかを調べ、含まれるときは $q$ と $B$ に保存した $r$ から $x = qm + r$ で離散対数問題が解かれます。

実装は以下の通りです（Python 3.8 以降）

def babystep_giantstep(g, y, p):
    m = int((p-1)**0.5 + 0.5)
    # Baby step
    table = {}
    gr = 1  # g^r
    for r in range(m):
        table[gr] = r
        gr = (gr * g) % p
    # Giant step
    gm = pow(g, -m, p)  # gm = g^{-m}
    ygqm = y            # ygqm = y * g^{-qm}
    for q in range(m):
        if ygqm in table: # 右辺と左辺が一致するとき
            return q * m + table[ygqm]
        ygqm = (ygqm * gm) % p
    return None

g = 7
y = 765686981
p = 35808104999
x = babystep_giantstep(g, y, p)
print(x)
print(pow(g, x, p) == y)

Pollard's rho法

離散対数問題 $y = g^x \mod{p}$ の $x$ を求める方法の1つである Pollard's rho法（ポラード・ローのアルゴリズム）はBSGSと同じ計算量ですが、使用するメモリ空間は定数量のみです。なので、BSGSよりも少ないメモリ空間で離散対数を計算することができます。

まず、関数 $f$ を次のように定義します（これはランダムに移動するための関数です）。

f(x) = \begin{cases} yx & (x \in G_1 のとき) \\ x^2 & (x \in G_2 のとき) \\ gx & (x \in G_3 のとき) \end{cases}

次に、数列 $x_i$ について考えます。乱数 $a_0$ を選び、 $x_0 = g^{a_0}$ とします。このときの数列 $x_i$ を漸化式 $x_{i+1} = f(x_i)$ で計算すると、この数列の項は $x_i = g^{a_0} y^{b_0}$ となります。ここで、 $b_0 = 0$ で、 $a_{i+1}, b_{i+1}$ は関数 $f$ の定義から次のようになります。

a_{i+1} = \begin{cases} a_i & (x_i \in G_1 のとき) \\ 2a \pmod{p} & (x_i \in G_2 のとき) \\ a_i + 1 \pmod{p} & (x_i \in G_3 のとき) \end{cases}

b_{i+1} = \begin{cases} b_i + 1 \pmod{p} & (x_i \in G_1 のとき) \\ 2b \pmod{p} & (x_i \in G_2 のとき) \\ b_i & (x_i \in G_3 のとき) \end{cases}

ポラード・ローのアルゴリズムでは、 $y = g^x$ を求めるために、まず $g^{a_i}y^{b_i} = g^{A_i}y^{B_i}$ を満たす $(a_i, b_i), (A_i, B_i)$ を求めます。 $G$ は $g$ が生成する巡回群なので、 $x_{i+k} = x_i$ となる $i, k$ が存在します。つまり、

\begin{aligned} g^{a_i}y^{b_i} &\equiv g^{a_{i+k}}y^{b_{i+k}} \pmod{p} \\ g^{a_i - a_{i+k}} &\equiv y^{b_{i+k} - b_i} \pmod{p} \\ g^{a_i - a_{i+k}} &\equiv g^{x(b_{i+k} - b_i)} \pmod{p} \\ (a_i - a_{i+k}) &\equiv x(b_{i+k} - b_i) \pmod{p-1} \\ \end{aligned}

が成立します。よって、次の式から離散対数 $x$ を求めることができます。

x \equiv (a_i - a_{i+k}) (b_{i+k} - b_i)^{-1} \pmod{p-1}

$g^{a_i}y^{b_i} = g^{A_i}y^{B_i} = g^{a_{i+k}}y^{b_{i+k}}$ を満たす $(a_i, b_i), (a_{i+k}, b_{i+k})$ の探し方は、巡回群の数列 $x_i$ は周期的であるので、フロイドの循環検出法（ウサギとカメのアルゴリズム）を使って求めます。

実装は以下の通りです（Python 3.8 以降）

def pollard_rho(g, y, p):
    q = (p-1) // 2
    # ランダムに移動するための関数
    def new_xab(x, a, b,  g, y, p, q):
        subset = x % 3
        if subset == 0:
            return ((x*x) % p, (a*2) % q, (b*2) % q)
        if subset == 1:
            return ((x*g) % p, (a+1) % q, b        )
        if subset == 2:
            return ((x*y) % p, a        , (b+1) % q)
    # フロイドの循環検出法
    x, a, b = 1, 0, 0
    X, A, B = x, a, b
    for i in range(1, p):
        x, a, b = new_xab(x, a, b,  g, y, p, q)
        X, A, B = new_xab(X, A, B,  g, y, p, q)
        X, A, B = new_xab(X, A, B,  g, y, p, q)
        if x == X:
            break
    res = ((a - A) * pow(B - b, -1, q)) % q
    if pow(g, res, p) == y:
        return res
    if pow(g, res + q, p) == y:
        return res + q
    return None

g = 7
y = 765686981
p = 35808104999
x = pollard_rho(g, y, p)
print(x)
print(pow(g, x, p) == y)

Pohlig–Hellman法

Pohlig–Hellman法（ポーリッヒ・ヘルマンのアルゴリズム）は巡回群 $G$ の位数 $\lvert G\rvert$ が $p_1, p_2, …$ に因数分解できるとき、 $G$ での離散対数 $y = g^x \pmod{p}$ の計算問題を、素数位数 $p_i$ での離散対数問題に帰着させることで問題を小さくして解き、得られた複数の結果から、最後に中国人剰余定理で答えを求める方法です。

まず、 $p-1$ が素数 $q_i$ に素因数分解できるとします。

p-1 = q_1 q_2 \cdots q_k

離散対数 $x$ は任意の自然数なので、商 $b_i$ と余り $a_i$ を使って $x = a_i + b_i q_i$ と書くことができます。次に $y = g^x \pmod{p}$ の両辺を $\frac{p-1}{q_i}$ 乗すると、以下のように式変形できます。なお、途中でオイラーの定理 $g^{p-1} \equiv g^{\varphi(p)} \equiv 1 \pmod{p}$ を使っています。

\begin{aligned} y^{\frac{p-1}{q_i}} &\equiv g^{x \,\frac{p-1}{q_i}} \pmod{p} \\ &\equiv (g^{a_i + b_i q_i})^{\frac{p-1}{q_i}} \pmod{p} \\ &\equiv g^{\frac{a_i (p-1)}{q_i}} g^{b_i (p - 1)} \pmod{p} \\ &\equiv g^{\frac{a_i (p-1)}{q_i}} \pmod{p} \\ \end{aligned}

見やすくするために、 $y_i = y^{\frac{p-1}{q_i}},\; g_i = g^{\frac{p-1}{q_i}}$ とおくと、

y_i \equiv g_i^{a_i} \pmod{p}

となり、最初に示した離散対数問題の形になりますが、余り $a_i$ は割る数 $q_i$ よりも小さいという事実から、より高速に離散対数 $a_i$ を求めることができます。

そして $(g_i)^{q_i} = (g^{\frac{p-1}{q_i}})^{q_i} = 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 以降）

# Baby-step Giant-step法
def babystep_giantstep(g, y, p, q=None):
    if q is None:
        q = p - 1
    m = int(q**0.5 + 0.5)
    # Baby step
    table = {}
    gr = 1  # g^r
    for r in range(m):
        table[gr] = r
        gr = (gr * g) % p
    # Giant step
    try:
        gm = pow(g, -m, p)  # gm = g^{-m}
    except:
        return None
    ygqm = y                # ygqm = y * g^{-qm}
    for q in range(m):
        if ygqm in table:   # 左辺と右辺が一致するとき
            return q * m + table[ygqm]
        ygqm = (ygqm * gm) % p
    return None

# Pohlig–Hellman法
def pohlig_hellman_DLP(g, y, p):
    crt_moduli = []
    crt_remain = []
    for q, _ in factor(p-1):
        x = babystep_giantstep(pow(g,(p-1)//q,p), pow(y,(p-1)//q,p), p, q)
        if (x is None) or (x <= 1):
            continue
        crt_moduli.append(q)
        crt_remain.append(x)
    x = crt(crt_remain, crt_moduli)
    return x

g = 2
y = 1094511311619717224471473901707
p = 2 * 32803 * 196159 * 1981991353 * 47814426923 + 1  # 素数p
x = pohlig_hellman_DLP(g, y, p)
print(x)
print(pow(g, x, p) == y)

指数計算法

離散対数問題 $y = g^x \mod{p}$ の $x$ を求める方法の1つである指数計算法（Index Calculus Algorithm）は、BSGSやロー法 ( $\rho$ 法) よりも効率の良いアルゴリズムです。

合同式 $y \equiv g^x \pmod{p}$ を解くために、一つの上界 $B$ を決めます。任意の整数 $\beta$ の全ての素因数 $q \in \mathbb{P}$ が $B$ 以下のとき、この集合（因子基底）は $F(B) = \{ q \in \mathbb{P} \;\vert\; q \le B \}$ と書き、このとき整数 $\beta$ は $B$ スムーズ（ $B$ -smooth）といいます。

まず、因子基底の全ての元に対して離散対数を計算します。乱数 $k \in \{1,...,p-1\}$ を選び、 $g^k \;\mathrm{mod}\;p$ が $B$ スムーズな数であるかを調べます。 $B$ スムーズな数であれば、その素因数分解を計算します ( $g^k \;\mathrm{mod}\; p = p_1^{e_1} p_2^{e_2} \cdots p_n^{e_n}$ )。これが十分な数 $C$ 個の合同式が集まるまで繰り返すと、次の式が得られます。

\begin{aligned} g^{k_1} &\equiv p_1^{e_{11}} p_2^{e_{21}} \cdots p_n^{e_{n1}} \pmod{p} \\ g^{k_2} &\equiv p_1^{e_{12}} p_2^{e_{22}} \cdots p_n^{e_{n2}} \pmod{p} \\ &\;\;\vdots \\ g^{k_C} &\equiv p_1^{e_{1C}} p_2^{e_{2C}} \cdots p_n^{e_{nC}} \pmod{p} \\ \end{aligned}

これらの合同式は、以下のように式変形することができます ( $p_i^{e_{ij}} = (g^{\log_g p_i})^{e_{ij}}$ と、巡回群の位数は $p-1$ なので)。

\begin{aligned} k_1 &\equiv e_{11} \log_g p_1 + e_{21} \log_g p_2 + \cdots + e_{n1} \log_g p_n \pmod{p-1} \\ k_2 &\equiv e_{12} \log_g p_1 + e_{22} \log_g p_2 + \cdots + e_{n2} \log_g p_n \pmod{p-1} \\ &\;\;\vdots \\ k_C &\equiv e_{13} \log_g p_1 + e_{2C} \log_g p_2 + \cdots + e_{nC} \log_g p_n \pmod{p-1} \\ \end{aligned}

次の目標は離散対数 $\log_g p_i$ を求めることです。この値はすべての合同式の各項で共通なので、行列の掛け算の形にすることができます。

\begin{pmatrix} k_1 \\ k_2 \\ \vdots \\ k_C \end{pmatrix} = \begin{pmatrix} e_{11} & e_{21} & \ldots & e_{n1} \\ e_{12} & e_{22} & \ldots & e_{n2} \\ \vdots & \vdots & \ddots & \vdots \\ e_{1C} & e_{2C} & \ldots & e_{nC} \end{pmatrix} \begin{pmatrix} \log_g p_1 \\ \log_g p_2 \\ \vdots \\ \log_g p_n \end{pmatrix} \pmod{p-1}

行列の形にしたことで、ガウスの消去法により、 $(\log_g p_i)$ を求めることができます。ここまでの事前計算が終わったら、最後に、求める離散対数 $x = \log_g y$ を計算します。まず、乱数 $\kappa \in \{1,...,p-1\}$ を選び、以下の式を計算します。

y g^\kappa \;\mathrm{mod}\;p

この値が $B$ スムーズになるまで乱数を選びなおします。もし $B$ スムーズであれば因数分解します。

y g^\kappa \equiv p_1^{c_1} p_2^{c_2} \cdots p_n^{c_n} \pmod{p}

この合同式は以下のように式変形できます。

\log_g y + \kappa \equiv c_1 \log_g p_1 + c_2 \log_g p_2 + \cdots + c_n \log_g p_n \pmod{p-1}

ここで $\log_g y$ 以外は全てわかっているので、離散対数 $\log_g y = x$ は以下で求めることができます。

x = \left( \sum_{i} c_i \log_g p_i - \kappa\right) \;\mathrm{mod}\;(p-1)

実装は以下の通りです（SageMath 9.0, Python 3.8）

# nがBスムーズな数であるかを調べる
def is_Bsmooth(b, n):
    factors = list(factor(int(n)))
    if len(factors) != 0 and factors[-1][0] <= b: 
        return True, dict(factors)
    else:
        return False, dict(factors)

# 連立合同式を求める
def find_congruences(B, g, p, congruences=[]):
    unique = lambda l: list(set(l))
    bases = []
    max_equations = prime_pi(B)
    while True:
        k = randint(2, p-1)
        ok, factors = is_Bsmooth(B, pow(g,k,p))
        if ok:
            # TODO: 線形独立のときだけ追加する
            congruences.append((factors, k))
            if len(congruences) >= max_equations:
                break
    bases = unique([base for c in [c[0].keys() for c in congruences] for base in c])
    return bases, congruences

# 連立合同式を行列に変換する
def to_matrices(R, bases, congruences):
    M = [[c[0][base] if base in c[0] else 0 \
            for base in bases] for c in congruences]
    b = [c[1] for c in congruences]
    return Matrix(R, M), vector(R, b)

# 指数計算法
def index_calculus(g, y, p, B=None):
    R = IntegerModRing(p-1)
    if B is None:
        B = ceil(exp(0.5*sqrt(2*log(p)*log(log(p)))))
    bases = []
    congruences = []
    # 合同式を満たす行列ができるまで繰り返す。
    for i in range(100):
        bases, congruences = find_congruences(B, g, p, congruences)
        M, b = to_matrices(R, bases, congruences)
        # Mx = b を満たす行列xを求める
        try:
            exponents = M.solve_right(b)
            break
        except ValueError:
            # matrix equation has no solutions
            continue
    else:
        return None
    # ag^y mod p がBスムーズである指数kを決定する
    while True:
        k = randint(2, p-1)
        ok, factors = is_Bsmooth(B, (y * pow(g,k,p)) % p)
        if ok and set(factors.keys()).issubset(bases):
            print('found k = {}'.format(k))
            break
    print('bases:', bases)
    print('q:', factors.keys())
    dlogs = {b: exp for (b,exp) in zip(bases, exponents)}
    x = (sum(dlogs[q] * e for q, e in factors.items()) - k) % (p-1)
    if pow(g, x, p) == y:
        return x
    return None

g = 2
y = 330456869054588
p = 924614919573299
x = index_calculus(g, y, p)
print(x)
print(pow(g, x, p) == y)

数体ふるい法

指数計算法を改良したものに「数体ふるい法」があります。いろいろ論文を見ながら特殊数体ふるい法の実装を試みたのですが、特定の値のときしか離散対数問題を解くことができないプログラムができてしまいました。一応、数体ふるい法のPython実装は検索してもあまりないので、誰かの役に立てばいいなと思い、残骸をGistに残しておきます。

　[WIP] Special Number Field Sieve (SageMath) – Gist

数体ふるい法で参考にした論文：

Worked Example for the Special Number Field Sieve … 計算途中の値があるので実装時のテストに使えます
An Introduction to the General Number Field Sieve … 5.5 Sieving や 5.6 Forming the Matrix の行列の作り方などが参考になります
離散対数問題の困難性に関する計算量についての調査・研究報告書 (CRYPTOREC) … 2.2 有限体上の離散対数問題へのNFS とかが参考になる (数少ない日本語での解説)
smallnfs/smallnfs.sage at master · lgremy/smallnfs … 数体ふるい法のSageMathによる実装らしきもの。正しく動作するか未確認

以上です。

🎄この記事はセキュリティキャンプ Advent Calendar 2020 - Adventar の9日目です。明日は TumoiYorozu さんで「キャンプの講義の記事」です。お楽しみに🎄

参考文献

J.A.ブーフマン著、林芳樹訳『暗号理論入門原書第3版』丸善出版 2012
Douglas R. Stinson 著、櫻井幸一監訳『暗号理論の基礎』共立出版 1994
代数曲線暗号とその安全性 - 第 15 回整数論サマースクール報告集, pp.223-238

晴耕雨読

working in the fields on fine days and reading books on rainy days