Pohlig-Hellman算法求Elgamal算法私钥+大步小步法解决DLP（离散对数问题）

问题描述：在利用Elgamal算法公钥中，已知参数（本题中数字均为十进制整数）

p=26622572818608571599593915643850055101138771
g=65539,
小明的公钥为：
14632691854639937953996750549254161821338360
如果你知道模数p-1可以分解为
{{2, 1}, {3, 1}, {5, 1}, {7, 1}, {11, 1}, {13, 1}, {17, 1}, {19, 1}, {29, 1}, {31, 1}, {37, 1}, {41, 1}, {47, 1}, {53, 1}, {61, 1}, {73, 1}, {97, 1}, {101, 1}, {103, 1}, {107, 1}, {113, 1}, {137, 1}, {139, 1}, {151, 1}, {167, 1}, {173, 1}, {179, 1}}
请利用Pohlig–Hellman算法求私钥x;

一、大步小步法

（一）、原理：

g为p的原根，p-1是g的阶，求解问题为$x={log}_g\beta 。$

若使用穷举方法对x尝试0到p-1之间的数，验证$\beta =g^x$，则算法复杂度为O(p-1)

小步大步法基于以下事实：

$\beta =g^x$，设$x=i\cdot m+j$，其中$0\le i,j<m$，如果令$m=\sqrt{p-1}$，此时建立一个条目为$（j,g^j）$的表(j表的长度为$\sqrt{m})$，则$i\cdot m+j$可以遍历1到p-1的全部内容，并且遍历的复杂度取决于i，为$O(\sqrt{p-1})$，（m为定值，j已经存放在表中，建立j表的复杂度为$O(\sqrt{p-1}）$，j表排序的复杂度为$O(\sqrt{p-1}\cdot log\sqrt{p-1})$,所以总复杂度为$O(\sqrt{p-1})$
$ps:在x=i\cdot m+j中$，对于每个i,j，其变化较小每次加1，这是小步，而导致x每次变化的间隔较大为m（因为i·m），这是大步，即为小步大步法命名

大步小步法分析参考文章：https://blog.csdn.net/m0_66201040/article/details/124453140
复杂度分析：

（二）、实现思路

输入：生成元g的阶p-1和元$\beta$
输出：离散对数$x={log}_g\beta$
– 设置$m=\lceil \sqrt{p-1}\rceil$
– 建立一个条目为 的表(j,$g^j(modp)$)，其中 。以条目中的第二项对表排序
– 计算 和设置g^m
– 从0到m-1进行如下循环
• 检查$\gamma$是否为表中某个第2项。
• 如果$\gamma \equiv g^i(modp)$,则返回(x=i·m+j)
• 设置 $\gamma =\gamma \cdot g^{-m}$ (其中$g^{-m}为g^m$的模逆元)

（三）、代码实现

## n为g的阶
import gmpy2
def little_big_step_method(g,beta,N,p):# 求解x = (log_g)beta
    m = int(math.sqrt(N))+1
    gj_j_dict= {}
	# 建立（g^j，j）表，与上述原理（j，g^j表内容相反）
    for j in range(0,m):
        gj_j_dict[pow(g,j,p)] = j 
    gj_sorted_list = sorted(gj_j_dict)
    gamma = beta
    for i in range(0,m):
        if gamma in gj_sorted_list:
            return i*m + gj_j_dict[gamma]
        gamma = (gamma * gmpy2.invert(pow(g,m,p),p) )%p

二、Pohlig-Hellman

（一）、 原理

如果要求解x使得$\beta =g^{x}modp$，g为素数p的原根，阶为p-1 ，如果$p-1={\Pi }_{i=1}^t{q}_i^{e^i}$，且$q^i$互素，那么若能高效求出$xmod{q}_1^{e_1},xmod{q}_2^{e_2}，...，xmod{q}_i^{e_i}$，就可以使用中国剩余定理快速求出$xmod(p-1)$

令$x_i=xmod{q}_i^{e_i}$，则会有$x=k_i\cdot q_i^{e_i}+r_i，\beta =g^{x}mod$，p记g的阶为N=p-1则有：
${\beta }^{\frac{N}{q_i^{e_i}}}modp=(g^x{)}^{\frac{N}{q_i^{e_i}}}=g^{ki\cdot N}\cdot g^{x_i\cdot \frac{N}{q_i^{e_i}}}$
由于$g^{ki\cdot N}modp=g^{ki\cdot (p-1)}modp=1modp$，所以：
${\beta }^{\frac{N}{q_i^{e_i}}}modp=g^{x_i\cdot \frac{N}{q_i^{\beta }{e}_i}}={(g^{\frac{N}{q_i^{e_i}}})}^{x_i}modp$
令${\beta }^{\prime }={\beta }^{\frac{N}{q_i^{e_i}}}，g^{\prime }=g^{\frac{N}{q_i^{e_i}}}$，所以到这步，就需求解问题${\beta }^{\prime }=(g^{\prime }{)}^{x_i}modp$，即：$x_i={log}_{g^{\prime }}^{{\beta }^{\prime }}$，就用小步大步法或是穷举法求出

（二）、实现思路

输入：已知大素数P的元根为g ，$\beta$ 是一个小于P的数，g模P的阶是N=p-1, 以及N的素数分解$N={\Pi }_{i=1}^t{q}_i^{e_i}$（t为N分解得到的素数的个数，相同的素数只记为一个）
输出：整数$x\in Z_N$ ，且$g^x=\beta modP$
1、对于每一个$i\in 1,...,t$
(1)计算$g_i=g^{\frac{N}{q_i^{e_i}}}$，易知g_i的阶是$q_i^{e_i}$(详见下图定理10.5)
(2)计算$h_i={\beta }^{\frac{N}{q_i^{e_i}}}$，$h_i\in <g_i>$
(3)使用穷举法或者大步小步法得到$x_i$，使得${g_i}^{x_i}=h_i$
2.利用中国剩余定理求解x，其中$x=x_{i}mod{q}_i^{e_i}$

（三）、代码实现

from hashlib import new
import math
import gmpy2


def little_big_step_method(g,beta,N,p):# 求解x = (log_g)beta
    m = int(math.sqrt(N))+1
    gj_j_dict= {}
	# 建立（g^j，j）表，与上述原理（j，g^j表内容相反）
    for j in range(0,m):
        gj_j_dict[pow(g,j,p)] = j 
    gj_sorted_list = sorted(gj_j_dict)
    gamma = beta
    for i in range(0,m):
        if gamma in gj_sorted_list:
            return i*m + gj_j_dict[gamma]
        gamma = (gamma * gmpy2.invert(pow(g,m,p),p) )%p
def solve_g_i(g,q_i,N):
    g_i = []
    for q in q_i:
        g_i.append(pow(g,(N//q),N+1))
    return g_i
def solve_h_i(beta,q_i,N):
    h_i = []
    for q in q_i:
        h_i.append(pow(beta,(N//q),N+1))
    return h_i
if __name__ == '__main__':
    beta = 14632691854639937953996750549254161821338360
    p=26622572818608571599593915643850055101138771
    g=65539
    N = p-1
    L = [[2, 1], [3, 1], [5, 1], [7, 1], [11, 1], [13, 1], [17, 1], [19, 1], [29, 1], [31, 1], [37, 1], [41, 1], [47, 1], [53, 1], [61, 1], [73, 1], [97, 1], [101, 1], [103, 1], [107, 1], [113, 1], [137, 1], [139, 1], [151, 1], [167, 1], [173, 1], [179, 1]]
    q_i = []

    for i in range(len(L)):
        q_i.append(L[i][0])
    g_i = [pow(g,N//q_i[i],p) for i in range(len(q_i))]
    h_i = [pow(beta,N//q_i[i],p) for i in range(len(q_i))]
    x_i = []
    # searchAll method
    # for i in range(len(q_i)):
    #     for j in range(q_i[i]):
    #         if pow(g_i[i],j,p) == h_i[i]:
    #             x_i.append(j)
    #             break
    for i in range(len(q_i)):
        x_i.append(little_big_step_method(g_i[i],h_i[i],q_i[i],p))
x = 0
for i in range(len(x_i)):
    Mi = N // q_i[i]
    Mii = gmpy2.invert(Mi, q_i[i])
    x = (x + x_i[i] * Mi * Mii) % N
    
print("私钥x为：%d"%pow(g,x,p))