目录

Implementation of discrete logarithm algorithm

Notes for $An\ Introduction\ to\ Mathematical\ Cryptography$

在整数中,离散对数(英语:Discrete logarithm)是一种基于同余运算和原根的一种对数运算。而在实数中对数的定义 $\log_ba$ 是指对于给定的 $a$ 和 b,有一个数 $x$,使得 $b^{x} = a$。相同地在任何群 $G$ 中可为所有整数 $k$ 定义一个幂数为 $b^{k}$,而离散对数 $\log_ba$ 是指使得 $b^{k} = a$ 的整数 $k$。 离散对数在一些特殊情况下可以快速计算。然而,通常没有具非常效率的方法来计算它们。公钥密码学中几个重要算法的基础,是假设寻找离散对数的问题解,在仔细选择过的群中,并不存在有效率的求解算法。

Shanks’s Babystep-Giantstep Algorithm

对于同余式 $g^x ≡ h\ (mod\ p)$ ,$g$ 为质数 $p$ 的一个原根,这种情况便可以使用 Shanks 算法。令 $N = p-1$,相比于暴力枚举法 $O(N)$ 的时间复杂度,Shanks 算法可以将时间复杂度降到 $O(\sqrt{N})$ .

  1. 计算 n = $\lfloor {\sqrt{N}} \rfloor + 1$.
  2. 构造出两个列表:
    • List1 = $[1,g,g^2,g^2,\dots ,g^n]$
    • List2 = $[h,h\cdot g^{-n},h \cdot g^{-2n},h \cdot g^{-3n},\dots ,h \cdot g^{-n^2}]$
  3. List1List2 中找到两个相同的元素,即 List1[i] == List2[j] ,那么就会有 $g^i = hg^{-jn}$.
  4. 那么 $x = i + jn$ 便是 $g^x = h$ 的一个解.
 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
#调用函数 sDLP(g,h,p) 返回 g^x≡h (mod p) 的一个解
#Shanks's Babystep-Giantstep Algorithm
from gmpy2 import invert,iroot
from Crypto.Util.number import getPrime

class node:
    def _init_(self):
        self.vue=0
        self.num=0
def cmp(a):
      return a.vue
def init_list(first,g,n,p):
      List=[]
      temp=node()
      temp.vue,temp.num=first,0
      List.append(temp)
      for i in range(1,n+1):
            temp=node()
            temp.num = i
            temp.vue = List[i-1].vue * g % p
            List.append(temp)
      List.sort(key=cmp)
      return List
def sDLP(a,b,p):
    ans=p
    n=iroot(p,2)[0]+1
    L1=init_list(1,a,n,p)
    aa=pow(invert(a,p),n,p)
    L2=init_list(b,aa,n,p)
    i = 0
    j = 0
    while True :
        if (i>=n or j>=n): break
        while (L1[i].vue < L2[j].vue and i<n): i += 1
        while (L1[i].vue > L2[j].vue and j<n): j += 1
        if L1[i].vue == L2[j].vue :
            x=L1[i].num+L2[j].num*n
            return int(x)
p = 552022109
g = 520158203
h = 525148510
print(sDLP(g,h,p))

Pohlig-Hellman Algorithm

Chinese remainder theorem

其实这里应该介绍 Pohlig–Hellman 离散对数算法了,但是其主要过程利用了中国剩余定理(CRT),相关计算和证明过程在百度百科就很容易找到。

这里引用 An Introduction to Mathematical Cryptography 书中的一段话:

In addition to being a theorem and an algorithm, we would suggest to the reader that the Chinese remainder theorem is also a state of mind.

CRT更是一种解决问题的思想,把一个大的问题分解为若干个小的问题分别求解,CRT使得若干个子问题的解可以再组合成原问题的解(同余方程组)。这种算法应用到解决离散对数问题显得非常巧妙。这也是 Pohlig-Hellman 离散对数算法的思想。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
# m 和 a 为两个列表,表示同余方程组 x mod m = a (m1,a1;m2,a2;...)
from functools import reduce
from gmpy2 import invert

def CRT(m,a):
    Num=len(m)
    M=reduce(lambda x,y: x*y, m)
    Mi=[M//i for i in m]
    t=[invert(Mi[i], m[i]) for i in range(Num)]
    x=0
    for i in range(Num):
        x+=a[i]*t[i]*Mi[i]
    return x%M

Pollard’s rho Algorithm

同样是为了实现 Pohlig-Hellman 离散对数算法,我们还需要了解一种分解质因数算法——Pollard’s rho 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
#通过调用 Factor(n) 进行质因数分解,返回值是因数列表。
from Crypto.Util.number import isPrime
from math import gcd

def f(x):
    return x**2 + 1

def pollard_rho(N):
    xn = 2
    x2n = 2
    d = 1
    while d == 1:
        xn = f(xn) % N
        x2n = f(f(x2n)) % N
        abs_val = abs(xn - x2n)
        d = gcd(abs_val, N)
    return d

def Factor(n):
    ans=[]
    while True:
        temp=pollard_rho(n)
        ans.append(temp)
        n=n//temp
        if n==1:return ans
        if isPrime(n):
            ans.append(n)
            return ans
'''
n=12345678754345678765456789876587654567899876
print(Factor(n))
output:[4, 3109, 3553454208763, 279372423577347576184497407]
'''

Pohlig-Hellman Algorithm

若 p-1 是一个 B-smooth number ,Pohlig-Hellman 离散对数算法只有在 B 较小的时候表现良好,所以对于解决离散对数问题(DLP)来说,其仍是微不足道的。同时也说明了如果在加密过程中选择了 p-1 的最大质因子较小的质数 p,那么此时的 DLP 可能容易受到Pohlig-Hellman 算法的攻击。

In number theory, a n-smooth (or n-friable) number is an integer whose prime factors are all less or equal to n.–Wikipedia

对于离散对数问题 $g^x ≡ h\ (mod\ p)$ , $N = p-1$ 为 $G$ 的秩,根据 唯一分解定理 : $N = q_1^{e_1}\cdot q_2^{e_2}\cdots q_t^{e_t}$

  1. 对 $N$ 进行分解,得到因数列表 $List\ q^e$ . 即 $qe = [q_1^{e_1},q_2^{e_2},\cdots ,q_t^{e_t}]$
  2. 用已知的每一项 $q_i^{e_i}$ 计算 $List\ g^{(p-1)/q^e}$ 和 $List h^{(p-1)/q^e}$ .
  3. 对于每一组 $g^{(p-1)/q^e}$ 和 $h^{(p-1)/q^e}$ ,计算出 $($g^{(p-1)/q^e})^x = h^{(p-1)/q^e}$ 中的 $x$ ,并存入 $List\ x$ .
  4. 最后一步利用CRT,对于 List 中对应的各 $i$ 项,$x =\ List\ x[i] (mod\ List\ q^e[i])$ 构成同余方程组 .

有了前面的 CRT 和 Pollard’s rho 质因数分解算法的铺垫,实现 Pohlig-Hellman 离散对数算法就会容易很多。

 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
from Crypto.Util.number import long_to_bytes
from functools import reduce
from gmpy2 import gcd,invert
from ShanksDLP import sDLP
from PollardRhoFactor import Factor
import time
#g^x = h (mod p)

def CRT(m,a):
      Num=len(m)
      M=reduce(lambda x,y: x*y, m)
      Mi=[M//i for i in m]
      t=[invert(Mi[i], m[i]) for i in range(Num)]
      x=0
      for i in range(Num):
            x+=a[i]*t[i]*Mi[i]
      return x%M
def BruteForceDLP(A,B,P):
      for i in range(P):
            if pow(A,i,P)==B:
                  return int(i)
def PohligHellman(g,h,p):
      qe=Factor(p-1)
      assert reduce(lambda x,y: x*y, qe) == p-1
      print(qe)
      Lg=[pow(g,(p-1)//i,p) for i in qe]
      Lh=[pow(h,(p-1)//i,p) for i in qe]
      length=len(qe)
      La=[]
      for i in range(length):
            if p<1000000000000:#p较小Shanks's算法可以接受就使用Shanks's解决
                  La.append(sDLP(Lg[i],Lh[i],p))
            else:#p-1的最大质因子较小的话暴力枚举法也有很好的表现
                  La.append(BruteForceDLP(Lg[i],Lh[i],p))
            #print(Lg[i],Lh[i],La[i])
      X=CRT(qe,La)
      if pow(g,X,p)==h:
            print("x is Right ! x = ",X)
      else:print("Wrang Answer")

print("g^x = h (mod p)")
p=int(input("p= "))
g=int(input("g= "))
h=int(input("h= "))
start_time=time.time()
PohligHellman(g,h,p)
print("it takes ",time.time()-start_time," seconds",)