离散对数

crypto 评论

​

在整数中，离散对数（Discrete Logarithm, DL）是一种基于同余运算和原根的一种对数运算。

在任何群 $G$ 中可为所有整数 $k$ 定义一个幂数为 $b^k$，而离散对数 $\log_ba$ 是指使得 $b^k = a$ 的整数 $k$。 离散对数在一些特殊情况下可以快速计算。然而，通常没有具非常效率的方法来计算它们。公钥密码学中几个重要算法的基础，是假设寻找离散对数的问题解，在仔细选择过的群中，并不存在有效率的求解算法。

​

基本定义

生成元

在一个群 $G$ 中，如果 $g$ 是 $G$ 的生成元，即所有 $G$ 中的所有元素都可以被表示成 $y=g^x$，此时的 $x$ 称为 $y$ 在 $G$ 中的对数。

阶

设 $m \ge 1,\gcd(a,m)=1$，使得 $a^d \equiv 1 \pmod m$ 成立的最小正整数 $d$ 称为 $a$ 对模 $m$ 的指数或者阶。一般将其计为 $\delta_m(a)$。

满足 $a^d \equiv 1 \pmod m$ 的 $d$ 一定满足 $d \mid \varphi(m)$。

原根

当 $\delta_m(a)=\varphi(m)$ 时，称 $a$ 是模 $m$ 的原根，简称 $m$ 的原根。

只有 $m=2,4,p^\alpha,2p^\alpha$（$p$ 为奇素数，$\alpha$ 为正整数）时，模 $m$ 的剩余系存在原根。（充要条件）

​

常规Sage函数

参数说明：求解以base为底，a的对数；ord为base的阶，可以缺省，operation可以是+与*，默认为*；bounds是一个区间(ld,ud)，需要保证所计算的对数在此区间内。

  • discrete_log(a,base,ord,operation)

    通用的求离散对数的方法。

  • discrete_log_rho(a,base,ord,operation)

    求离散对数的Pollard-Rho算法。

  • discrete_log_lambda(a,base,bounds,operation)

    求离散对数的Pollard-kangaroo算法（也称为lambda算法）。

  • bsgs(base,a,bounds,operation)

    小步大步法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#生成64位的素数p作为模数，int为32位，超过int要在数字后加L
p=random_prime(2L**64)
 
#定义有限域GF(p)
G=GF(p)
 
#找一个模p的原根
gp ('znprimroot('+str(p)+')')
#输出Mod(rt,p)，则x是模p的原根
g=G(rt)
 
#生成私钥
x=G(ZZ.random_element(p-1))
 
#公钥y=g^x mod p，由于已经定义在GF(p)上，因此g**x就是g^x mod p
y=g**x
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#计算DLP的通用方法
discrete_log(y,g)==x
#n为合数（Pohlig-Hellman）
x = discrete_log(mod(b,n),mod(a,n)) 
#n为质数或质数幂（线性筛Index Calculus）
R = Integers(99)
a = R(4)
b = a^9
b.log(a)
#或
x = int(pari(f"znlog({int(b)},Mod({int(a)},{int(n)}))"))
x = gp.znlog(b, gp.Mod(a, n))
 
#计算离散对数的lambda方法
discrete_log_lambda(y,g,(floor(ZZ(x)/2),2*ZZ(x)))==x
 
#小步大步法计算DLP
bsgs(g,y,(floor(ZZ(x)/2),2*ZZ(x)))==x

#GF(p)下计算DLP
F = GF(p)
ax = F(b + (a - 1) * A) / F(a * s - s + b)
x = discrete_log(F(ax), F(a))

自改DLP：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
multiplication_names = ( 'multiplication', 'times', 'product', '*')
addition_names       = ( 'addition', 'plus', 'sum', '+')
def old_discrete_log(a, base, ord=None, operation='*', identity=None, inverse=None, op=None):
     b = base

     from operator import inv, mul, neg, add
     Z = Integers()

     if operation in multiplication_names:
         identity = b.parent()(1)
         inverse  = inv
         op = mul
         if ord==None:
             ord = b.multiplicative_order()
     elif operation in addition_names:
         identity = b.parent()(0)
         inverse  = neg
         op = add
         if ord==None:
             ord = b.order()
     else:
         if ord==None or identity==None or inverse==None or op==None:
             print(ord, identity, inverse, op)

     if ord < 100:
         c = identity
         for i in range(ord):
             if c == a:        # is b^i
                 return Z(i)
             c = op(c,b)

     m = ord.isqrt()+1  # we need sqrt(ord) rounded up
     table = dict()     # will hold pairs (b^j,j) for j in range(m)
     g = identity       # will run through b**j    
     for j in range(m):
         if a==g:
             return Z(j)           
         table[g] = j
         g = op(g,b)

     g = inverse(g)     # this is now b**(-m)
     h = op(a,g)        # will run through a*g**i = a*b**(-i*m)
     for i in range(1,m):
         j = table.get(h)
         if not j==None:  # then a*b**(-i*m) == b**j
             return Z(i*m + j)
         if i < m-1:
             h = op(h,g)

def d_log(q): 
    print("dlogging: ", q)
    dlogs = []
    for f in factors:
        t = order//f
        qt = G_mul(q,t)
        gent = G_mul(gen, t)
        dlog = old_discrete_log(qt, gent, ord=f, operation='NONE', op=G_add, identity=zero, inverse=inverse)
        dlogs.append(dlog)
        if None in dlogs:
            raise ValueError("oh no")
    l = CRT_list(dlogs, factors)
    return l

gen = (,)
c = (,)
k = d_log(c)
print(k)

​

大步小步算法（BSGS）

大步小步算法常用于求解用于解决解高次同余方程的问题，问题形式如：有同余方程$a^x \equiv b \pmod p$，$p$ 为质数，求最小非负整数解 $x$ 使得原方程成立。这类问题也称为离散对数问题。该算法的复杂度可以达到$O(\sqrt{p}\log{n})$甚至更低。

1
2
#Sage
bsgs(base,a,bounds,operation)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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 = 
y = 
p = 
x = babystep_giantstep(g, y, p)
print(x)
print(pow(g, x, p) == y)

​

Pollard Rho算法

Pollard Rho算法是一种随机性的概率型算法，基于PR算法我们可以加速大整数的因子分解，这也是1975年Pollard提出该算法时的应用方面，这可以用来攻击RSA的公钥密码体制，之后在1978年Pollard又利用循环群 $ρ$ 形的特点提出了用于解决离散对数问题的PR算法，后来也扩展到了椭圆曲线上，这样PR算法便成功威胁到了整个公钥密码体系。

解决离散对数问题的Pollard Rho算法，可以帮助我们在有限的循环群上解决离散对数问题。简单来看我们取一个大素数 $P$，设 $G$ 为一个乘法的循环群，其生成元为 $g$，则该群的阶就是 $N=P-1$。此时在 $G$ 上的离散对数问题就是对于 $G$ 中的元素 $q$，找到 $x$ 使得 $g^x=q$，而Pollard Rho算法主要就是利用了循环群的生成序列呈现一个 $ρ$ 字形状。

使用Pollard Rho算法后可以将求解的时间复杂度降到 $O(\sqrt{N})$，至于空间复杂度则可以忽略不计，因为该算法是个概率型算法，每次只随机选取一对进行比较，所以不会占用存储，这与大步小步算法相反，其实PR算法也可以看作是大步小步算法的随机化版本，继承令时间复杂度的有点，同时还减少了空间的消耗，因为大步小步算法需要占用大量的存储空间。

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

参考：https://xz.aliyun.com/t/2780

1
2
#Sage
discrete_log_rho(a,base,ord,operation)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
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 = 
y = 
p = 
x = pollard_rho(g, y, p)
print(x)
print(pow(g, x, p) == y)

​

Pohlig-Hellman算法

Pohig-Hellman算法是一种求解光滑阶循环群上的离散对数的方法。

Pohlig-Hellman算法的复杂度在一般情况下比BSGS算法高，但是在特殊情况下（循环群的阶是光滑数，即可以因子分解成较小的数的乘积），使用Pohlig-Hellman能取得好的效果。而且有些时候，尽管BSGS能够将复杂度降至$\sqrt{p}$，但是这个数依然很大，所以不能用。这时可以考虑Pohlig-hellman方法能不能起作用。

  • 算法思想

    考虑DLP问题：$a^x \equiv b \pmod p$，因为 $p$ 是大素数，模 $p$ 的循环群的阶是 $p−1$。假设模 $p$ 的最小的本原元是 $g$（本原元是可以求的），那么有

    $\left\{\begin{array}{c} a \equiv g^{a’} \pmod p \\ b \equiv g^{b’} \pmod p \end{array}\right.$

    进一步有 $a^x \equiv b \pmod p \Longleftrightarrow g^{a’x} \equiv g^{b’} \pmod p$，

    有 $a’x \equiv b’ \pmod {p-1}$，

    如果求出了满足上式的 $a’$ 和 $b’$ ，通过扩展gcd方法可以求一次同余方程的解得到 $x$ 。

    问题归结成如何求 $a’$ 和 $b’$ ，即原本的一个离散对数问题，现在变成两个离散对数问题：

    $\left\{\begin{array}{c} a \equiv g^{a’} \pmod p \\ b \equiv g^{b’} \pmod p \end{array}\right.$

    以求 $a’$ 为例，解DLP问题：$g^x \equiv a \pmod p$：

    1. 将 $p-1$ 进行标准的素因子分解，即 $p-1=\prod\limits_{i=1}^{m} p_i^{k_i}=p_1^{k_1}p_2^{k_2}p_3^{k_3} \cdots p_m^{k_m}$；

    2. 对每个素因子 $p_i$，将 $x$ 表示成 $p_i$ 进制，有

      $x=a_0+a_1p_1+a_2p_2^2+a_3p_3^3+\dots+a_{k_i-1}p_i^{k_i-1} \pmod {p_i^{k_i}}$，

      这样的 $p_i$ 进制表示，系数 $a_i$ 自然是小于 $p_i$；

    3. 令 $r=1$，有 $(g^x)^{\frac{p-1}{p_i^r}} \equiv a^{\frac{p-1}{p_i^r}} \pmod p$，展开 $x$ 有

      $(g^{a_0+a_1p_1+a_2p_2^2+a_3p_3^3+\dots+a_{k_i-1}p_i^{k_i-1}})^{\frac{p-1}{p_i^r}} \equiv a^{\frac{p-1}{p_i^r}} \pmod p$

      $g^{a_0{\frac{p-1}{p_i}}} \cdot g^{a_1(p-1)} \cdot g^{a_2(p-1)p_i} \cdots g^{a_{k_i-1}(p-1)p_i^{k_i-2}} \equiv a^{\frac{p-1}{p_i}} \pmod p$

      注意从第二项开始，每一项指数都包含 $p−1$（由费马小定理知 $g^{p-1} \equiv 1 \pmod p$），所以式子变成

      $g^{a_0{\frac{p-1}{p_i}}} \equiv a^{\frac{p-1}{p_i}} \pmod p$

      这个式子中只有 $a_0$ 是未知的，因为 $a_0 \in [0,p_i−1]$，所以可以穷举得到 $a_0$ 的值。

    4. 再令 $r=2,3,4,\cdots,k_i$，重复步骤3，依次穷举求出 $a_1,a_2,\cdots,a_{k_i-1}$，整个的时间复杂度是 $\mathrm{O}(p_ik_i)$。可以得到

      $x=a_0+a_1p_1+a_2p_2^2+a_3p_3^3+\dots+a_{k_i-1}p_i^{k_i-1} \pmod {p_i^{k_i}}$

    5. 重复上述过程，得到 $m$ 个关于 $x$ 的式子，利用中国剩余定理（CRT），可以计算出 $x$ 的值。

    利用这个方法求出 $a’$ 和 $b’$ 后，就可以得到原DLP问题的解。

    注：

    Pohlig-Hellman算法最重要的点是利用了原根的性质，它只能解决 $a \equiv g^x \pmod p$（$g$ 是原根）的问题，对于 $g$ 不是原根的情况，需要利用原根将原方程转化成两个关于原根的离散对数问题。

    注意Pohlig-hellman只能用于求解光滑阶群，也就是 $p-1$ 可以分解成小的素因子乘积。否则，穷举 $a_i$ 的时间复杂度依旧很高。另外可以考虑在穷举 $a_i$ 时利用小步大步算法，进一步优划算法复杂度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
# 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 = 
y = 
p = 
x = pohlig_hellman_DLP(g, y, p)
print(x)
print(pow(g, x, p) == y)

​

指数计算法 / Index Calculus Algorithm

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
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:
            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)
        try:
            exponents = M.solve_right(b)
            break
        except ValueError:
            # matrix equation has no solutions
            continue
    else:
        return None
    # ag^y mod p
    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 = 
y = 
p = 
x = index_calculus(g, y, p)
print(x)
print(pow(g, x, p) == y)

​

多项式DLP

参考：InCTF 2020 - DLPoly

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
#Sage
def brute_dlp(gi, ci, n, lim):
	bi = gi
	for i in range(1, lim+1):
		if bi == ci:
			return i
		bi = (bi * gi) % n
	print("[-] NOT in the range")
	print("[-] Something's Wrong, you gotta check the range", lim)

def pohlig_hellman(g, c, s, n, factors):
	res = []
	modulus = []
	for q, e in factors:
		assert pow(g, s//(q**e), n) != 1
		gi = pow(g, s//(q**e), n)
		ci = pow(c, s//(q**e), n)
		dlogi = brute_dlp(gi, ci, n, q**e)
		print("[+] dlog modulo {0} == {1}".format(q**e, dlogi))
		res.append(dlogi)
		modulus.append(q**e)
	print("\n[*] res = ", res)
	print("[*] modulus = ", modulus)
	dlog = CRT(res, modulus)
	print("\n[+] dlog modulo {0} == {1}".format(prod(modulus), dlog))
	return dlog

p = 
P = PolynomialRing(Zmod(p), name='x')
x = P.gen()
n = 
g = 
nfactors = n.factor()
s = 1
for i in nfactors:
	s *= p**(i[0].degree()) - 1

print(s)
print(factor(s))
qe = 
dlog = pohlig_hellman(g, c, s, n, qe)

flag = bytes.fromhex(hex(dlog)[2:])
print("\n[*] flag = ", flag.decode())

​

矩阵DLP

给定 $n$ 阶矩阵 $G$，求满足 $G^x=H$ 的 $x$ 值。

通常矩阵的乘法是复杂的，而矩阵的幂就更是这样了。注意到 $A^n=P(P^{-1}AP^n)P^{-1}$，

如果能找到合适的矩阵 $P$ 使得 $P^{−1}AP$ 的幂是简单的，就找到了计算矩阵的幂乃至多项式的简便方法。例如计算对角矩阵的幂只需对分量求幂，这就非常简单。

但是并非所有的 $n$ 阶矩阵都可以对角化，经过复杂的论证，得知 $n$ 阶矩阵一定有 Jordan 标准型，也就是对于 $n$ 阶矩阵 $A$，构造可逆矩阵 $P$ 使得 $P^{-1}AP=\text{diag}\{J_1,J_2,\cdots,J_s\}$，其中 $J_1,J_2,\cdots,J_s$ 是 Jordan 块。所谓特征值为 $\lambda$ 的 $n$ 阶 Jordan 块，就是 $n$ 阶矩阵

$J_n(\lambda)=\left(\begin{matrix} \lambda & 1 & 0 & \cdots & 0 & 0 \newline 0 & \lambda & 1 & \cdots & 0 & 0 \newline 0 & 0 & \lambda & \cdots & 0 & 0 \newline \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \newline 0 & 0 & 0 & \cdots & \lambda & 1 \newline 0 & 0 & 0 & \cdots & 0 & \lambda \newline \end{matrix}\right)$

特别地，对于矩阵 $A$ 可对角化的情形，相应的 Jordan 块 $J_1,J_2,\cdots,J_s$ 均为一阶的。

计算 $(J_n(\lambda))^k$. 注意到 $J_n(\lambda)=\lambda E+J_n(0)$, 应用二项式定理，得到

$(J_n(\lambda))^k=\left(\begin{matrix} \lambda^k & k\lambda^{k-1} & \cdots & \mathrm{C}_k^{n-1}\lambda^{k-n+1} \newline 0 & \lambda^k & \cdots & \mathrm{C}_k^{n-2}\lambda^{k-n+2} \newline \vdots & \vdots & \ddots & \vdots \newline 0 & 0 & \cdots & \lambda^k \newline \end{matrix}\right)$​

推导：

对于 $v=A^nu$，求出 $A=PJP^{-1}$，有 $A^n=(PJP^{-1})^n=PJ^nP^{-1}$。

所以 $v=A^nu$ 可以化为 $v’=J^nu’$，其中 $v’=P^{-1}v$ 和 $u’=P^{-1}u$​。

计算 $J^n$ 可借助 WolframAlpha 计算通项公式。

1
2
3
4
5
6
7
8
9
10
n = 
p = 
A = []
B = []

G = matrix(GF(p), n, n, A)
H = matrix(GF(p), n, n, B)

G_Jor, P = G.jordan_form(transformation=True)
H_Jor = ~P * H * P

其他：

Pohlig–Hellman算法（ $p$ 光滑）

discrete_log(H, G,algorithm='lambda')

参考：

Sharif CTF 8 - ElGamat

AIS3 EOF CTF Quals 2021 - notRSA