离散对数

在整数中,离散对数(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的对数;ordbase的阶,可以缺省,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