离散对数问题(Discrete logarithm Problem,DLP)：

已知$g,p,y$，对于方程$y \equiv g^x \pmod p$，求解x。

其中：p为素数，g为p的生成元

注：选择p的生成元作为g以最大程度上降低被爆破的可能性，更加安全，否则会提前进入循环

基本概念

生成元：

设 G 是一个群（满足封闭性、结合律、存在单位元、每个元素有逆元的集合），若存在子集 $S \subseteq G$，使得 G 中每一个元素都可以表示为 S 中元素及其逆元的有限次乘积，则称 S 是群 G 的生成集，S 中的每个元素称为 G 的生成元。

元素的阶：

设 G 是一个群，$g \in G$ 是群中一个元素，元素 g 的阶定义为：

满足$g^k = e$（e 是群的单位元，$g^k$ 表示 g 自身运算 k 次，如加法群中为 kg）的最小正整数 k，记为 $\text{ord}(g)$

若不存在这样的正整数 k（即 $g^k \neq e$ 对所有 $k > 0$ 成立），则称 g 的阶为无穷大。

原根：

模 m 的原根：设 m 是正整数，a 是整数，若满足：

1. a 与 m 互素(即 $\gcd(a, m) = 1$)；
2. a模m的阶(在乘法群($\mathbb{Z}_m^*, \times$)中)等于该群的阶$\phi(m)$(欧拉函数，即小于m且与m互素的正整数个数)，则称a是模 m 的原根。

求解方式

暴力破解

枚举x进行求解，通常在p很小的时候如此进行

小步大步法(Baby-step giant-step,BSGS)算法

使用的是中间相遇攻击的思想，适用于x的范围较小(修改m为$int(\sqrt {x'})$，x'为x的最大范围)或者p较小的情况。

令$x = im + j$，其中$m = \sqrt{n}$，那么整数$i$和$j$都在$0 \sim m$的范围内

有

$$y = g^x = g^{im+j}$$

$$y(g^{-m})^i = g^j$$

枚举所有的j并进行计算，将结果储存到一个集合S中，再枚举i，计算$y(g^{-m})^i$，一旦发现结果在S中，则说明得到一组i和j满足需求。

import math

"""快速幂运算"""
def spow(a, b, mod):
    result = 1
    a = a % mod
    while b > 0:
        if b % 2 == 1:
            result = (result * a) % mod
        a = (a * a) % mod
        b = b // 2
    return result

"""BSGS算法"""
def BSGS(g, y, p, order=None):
    if y == 1:
        return 0

    if order==None:
        order = p - 1

    m = int(math.isqrt(order)) + 1
    table = {}

    current = 1
    for r in range(m):
        if current not in table:
            table[current] = r
        current = (current * g) % p

    gm_inv = spow(g, m * (p-2), p)
    current = y % p

    for q in range(m):
        if current in table:
            return q * m + table[current]
        current = (current * gm_inv) % p

    return None

g = 
y = 
p = 

x = BSGS(g,y,p)
assert pow(g,x,p) == y
print(x)

Pollard Rho算法(PR)

适用于生成元的阶的素因子都是大数的情形。

关键是构造$g^{a_1}h^{b_1} = g^{a_2}h^{b_2} \pmod p$，进而推导出x：

序列构造函数定义：

$$f(x,a,b) = \begin{cases} (x^2 \mod p,2a \mod (p-1),2b \mod (p-1)) \quad x \mod 3 = 0 \\ (x*g \mod p,a+1 \mod (p-1),2b \mod (p-1)) \quad x \mod 3 = 1\\ (x*h \mod p,a \mod (p-1),b+1 \mod (p-1)) \quad x \mod 3 = 2\\ \end{cases}$$

这样可以保证$x_k = g^{a}h^{b}$一直成立

碰撞过程：一个一次一步(龟速迭代)，另外一个一次两步(兔速迭代)

当两者相等时：$g^ah^b \equiv g^ch^d \mod p$，由于g是生成元，其幂次在模n意义下唯一，因此指数必须满足：

$a + bx \equiv c + dx \mod (p-1)$ （$h = g^x$代入)

随后求解关于x的线性同余方程组

"""快速幂运算"""
def spow(a, b, mod):
    result = 1
    a = a % mod
    while b > 0:
        if b % 2 == 1:
            result = (result * a) % mod
        a = (a * a) % mod
        b = b // 2
    return result

"""扩展欧几里得算法"""
def exgcd(a, b):
    if b == 0:
        return (a, 1, 0)
    else:
        g, x, y = exgcd(b, a % b)
        return (g, y, x - (a // b) * y)

"""Pollard Rho算法"""
def pollards_rho(g, h, p):
    if h == 1:
        return 0
    if g == h:
        return 1

    n = p - 1

    def f(x, a, b):
        if x % 3 == 0:
            return ((x * x) % p, (2 * a) % n, (2 * b) % n)
        elif x % 3 == 1:
            return ((x * g) % p, (a + 1) % n, b)
        else:
            return ((x * h) % p, a, (b + 1) % n)

    x, a, b = 1, 0, 0
    y, c, d = x, a, b

    while True:
        x, a, b = f(x, a, b)
        y, c, d = f(y, c, d)
        y, c, d = f(y, c, d)

        if x == y:
            da = (d - b) % n
            bc = (a - c) % n

            gcd_val, inv_da, _ = exgcd(da, n)

            if bc % gcd_val != 0:
                x, a, b = 1, 0, 0
                y, c, d = x, a, b
                continue

            inv = (inv_da * (bc // gcd_val)) % (n // gcd_val)

            for k in range(gcd_val):
                candidate = (inv + k * (n // gcd_val)) % n
                if spow(g, candidate, p) == h % p:
                    return candidate

            return None

g =
y = 
p =
x = pollards_rho(g,y,p)
assert pow(g,x,p) == y
print(x)

Pohlig-Hellman算法(PH)

适用于求解光滑阶循环群上的离散对数，即：

$$p-1 = \prod_{i=1}^n p_i^{k_i} = p_1^{k_1}*p_2^{k_2} \cdots$$

根据拉格朗日定理，G存在阶为$p_i^{k_i}$的子群，对于每个$p_i^{k_i}$，构造子群元素：$g_i = g^{n/{p_i^{k_i}}}$和$h_i = h^{n/{p_i^{k_i}}}$，此时有$g_i^{x_i} = h_i \mod p$，$x_i = x \mod p_i^{k_i}$：

$$h_i = h^{n/{p_i^{k_i}}} = (g^x)^{n/{p_i^{k_i}}} = (g^{n/{p_i^{k_i}}})^x = g_i^x$$

$$g_i^{x_i} = g_i^x \mod p$$

由于$g_i$的阶为$p_i^{k_i}$，所以指数必满足$x_i = x \mod p_i^{k_i}$

(子群的幂次满足模阶律：$g_i^a \equiv g_i^b \mod p \iff a ≡ b \mod p_i^{k_i}$)

因为$p_i^{k_i}$较小，可以通过BSGS求出$x_i$，最后多组$(x_i,p_i^{k_i})$可以通过CRT求出x。

from sympy import factorint 
from sympy.ntheory.modular import crt
import math

"""快速幂运算"""
def spow(a, b, mod):
    result = 1
    a = a % mod
    while b > 0:
        if b % 2 == 1:
            result = (result * a) % mod
        a = (a * a) % mod
        b = b // 2
    return result

"""扩展欧几里得算法"""
def exgcd(a, b):
    if b == 0:
        return (a, 1, 0)
    else:
        g, x, y = exgcd(b, a % b)
        return (g, y, x - (a // b) * y)

"""BSGS算法"""
def BSGS(g, y, p, order=None):
    if y == 1:
        return 0

    if order==None:
        order = p - 1

    m = int(math.isqrt(order)) + 1
    table = {}

    current = 1
    for r in range(m):
        if current not in table:
            table[current] = r
        current = (current * g) % p

    gm_inv = spow(g, m * (p-2), p)
    current = y % p

    for q in range(m):
        if current in table:
            return q * m + table[current]
        current = (current * gm_inv) % p

    return None

"""Pohlig-Hellman算法"""
def pohlig_hellman(g, h, p):
    if math.gcd(g, p) != 1:
        raise ValueError("g必须与p互质")

    n = p - 1
    factors = factorint(n)

    remainders = []
    moduli = []

    for prime, exp in factors.items():
        pe = pow(prime, exp)
        g_i = spow(g, n // pe, p)
        h_i = spow(h, n // pe, p)

        x_i = BSGS(g_i, h_i, p, pe)
        if x_i is None:
            return None

        remainders.append(x_i)
        moduli.append(pe)

    x, _ = crt(moduli, remainders)
    return x

g = 
p = 
h = 

x = pohlig_hellman(g, h, p)
assert pow(g,x,p) == h
print(x)

参考

• https://lazzzaro.github.io/2020/05/07/crypto-%E7%A6%BB%E6%95%A3%E5%AF%B9%E6%95%B0/
• https://ctf-wiki.org/crypto/asymmetric/discrete-log/discrete-log/