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=loggβ。x = {log}_g{\beta}。x=loggβ。

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

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

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

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

（二）、实现思路

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

（三）、代码实现

## n为g的阶
import gmpy2
def little_big_step_method(g,beta,N,p):# 求解x = (log_g)betam = int(math.sqrt(N))+1gj_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 = betafor 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使得β=gxmodp\beta = g^x mod pβ=gxmodp，g为素数p的原根，阶为p-1 ，如果p−1=Πi=1tqieip-1= \Pi_{i=1}^{t}q_i^{e^i}p−1=Πi=1tqiei，且qiq^{i}qi互素，那么若能高效求出xmodq1e1,xmodq2e2，...，xmodqieix\;mod\;q_1^{e_1},x\;mod\;q_2^{e_2} ，...，x\;mod\;q_i^{e_i}xmodq1e1,xmodq2e2，...，xmodqiei，就可以使用中国剩余定理快速求出xmod(p−1)x\;mod\;(p-1)xmod(p−1)

令xi=xmodqieix_i = x\;mod\;q_i^{e_i}xi=xmodqiei，则会有x=ki⋅qiei+ri，β=gxmodx=k_i·q_i^{e_i}+r_i，\beta = g^x\;mod\;x=ki⋅qiei+ri，β=gxmod，p记g的阶为N=p-1则有：
βNqieimodp=(gx)Nqiei=gki⋅N⋅gxi⋅Nqiei\beta^{\cfrac{N}{q_i^{e_i}}} mod\;p = (g^x)^{\cfrac{N}{q_i^{e_i}}} = g^{ki·N}·g^{x_i · {\cfrac{N}{q_i^{e_i}}}}βqieiNmodp=(gx)qieiN=gki⋅N⋅gxi⋅qieiN
由于gki⋅Nmodp=gki⋅(p−1)modp=1modpg^{ki·N}\;mod\; p = g^{ki·(p-1)}\;mod\; p = 1\;mod\;pgki⋅Nmodp=gki⋅(p−1)modp=1modp，所以：
βNqieimodp=gxi⋅Nqiβei=(gNqiei)ximodp\beta^{\cfrac{N}{q_i^{e_i}}} mod\;p = g^{x_i · {\cfrac{N}{q_i^\beta{e_i}}}} = {(g^{{\cfrac{N}{q_i^{e_i}}}})}^{\;x_i}\;mod \;pβqieiNmodp=gxi⋅qiβeiN=(gqieiN)ximodp
令β′=βNqiei，g′=gNqiei\beta' =\beta^{\cfrac{N}{q_i^{e_i}}}，g' = g^{{\cfrac{N}{q_i^{e_i}}}}β′=βqieiN，g′=gqieiN，所以到这步，就需求解问题β′=(g′)ximodp\beta' = (g')^{x_i}\;mod\;pβ′=(g′)ximodp，即：xi=logg′β′x_i = {log}_{g'}^{\beta'}xi=logg′β′，就用小步大步法或是穷举法求出

（二）、实现思路

输入：已知大素数P的元根为g ，β\betaβ 是一个小于P的数，g模P的阶是N=p-1, 以及N的素数分解N=Πi=1tqieiN = \Pi_{i=1}^tq_i^{e_i}N=Πi=1tqiei（t为N分解得到的素数的个数，相同的素数只记为一个）
输出：整数x∈ZNx\in Z_Nx∈ZN ，且gx=βmodPg^x= \beta \;mod\;Pgx=βmodP
1、对于每一个i∈1,...,ti \in {1,...,t}i∈1,...,t
(1)计算gi=gNqieig_i = g^{\cfrac{N}{q_i^{e_i}}}gi=gqieiN，易知g_i的阶是qieiq_i^{e_i}qiei(详见下图定理10.5)
(2)计算hi=βNqieih_i = \beta^{\cfrac{N}{q_i^{e_i}}}hi=βqieiN，hi∈<gi>h_i\in<g_i>hi∈<gi>
(3)使用穷举法或者大步小步法得到xix_ixi，使得gixi=hi{g_i}^{x_i}= h_igixi=hi
2.利用中国剩余定理求解x，其中x=ximodqieix=x_i \;mod\;q_i^{e_i}x=ximodqiei

（三）、代码实现

from hashlib import new
import math
import gmpy2def little_big_step_method(g,beta,N,p):# 求解x = (log_g)betam = int(math.sqrt(N))+1gj_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 = betafor 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 = 14632691854639937953996750549254161821338360p=26622572818608571599593915643850055101138771g=65539N = p-1L = [[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)#             breakfor 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) % Nprint("私钥x为：%d"%pow(g,x,p))