Pohlig-Hellman算法解决DLP问题

Pohlig-Hellman算法解决DLP问题及Python实现

$Instruction$

Pohlig-Hellman算法适用于求解特定条件下的DLP问题（离散对数问题），以数学表达该问题即是
$$a^x\equiv b(\mathrm{mod}p)$$
其中$x$即是所需要求得的数，已知$a,b,p$；该算法应用的特定条件是指$p-1$是光滑数（光滑数即是该数的质因数都很小，其每个质因数大小在该算法中对应着算法复杂度）

该算法的原理及运行过程如下，需要假设其模数$p-1$是光滑数

假设模$p$的最小本原元为$g$，那么可以用本原元来表示在模$p$意义下的$a$和$b$
$$\begin{gathered} a\equiv g^u(\mathrm{mod}p) \\ b\equiv g^v(\mathrm{mod}p) \end{gathered}$$
代回到初始问题上得到
$$\begin{gathered} a^x\equiv b(\mathrm{mod}p) \\ \Rightarrow g^{ux}\equiv g^v(\mathrm{mod}p) \end{gathered}$$
从上式根据数论定理得到
$$ux\equiv v(\mathrm{mod}p-1)$$
这样我们可以用$u,v$来表示出$x$
$$x\equiv v\cdot u^{-1}(\mathrm{mod}p-1)$$
到这里我们将问题转换到了如何得到$v,u$的大小

---

回到用本原元表示$a,b$的同余式上，

• 以解$g^u\equiv a(\mathrm{mod}p)$为例，

我们已知$a,g,p$；根据算法做出的假设，$p-1$可以被多个小质因数分解，那么
$$p-1=p_1^{q_1}\cdot p_2^{q_2}\cdot p_3^{q_3}\cdots p_n^{q_n}$$
其中$q_i$表示的是质因数的指数；我们将目标$u$用每个素因子的大小作为进制表示，如以$p_1$进制表示（进制表示如果不熟悉的话可以参照二进制的表示方式）
$$u=c_0+c_1{p}_1+c_2{p}_1^2+c_3{p}_1^3+\cdots +c_{q_1-1}{p}_1^{q_1-1}+\cdots +c_k{p}_1^k$$
其中$c_i$均小于$p_1$，变量$k$意味着可以到无限大；同理我们也将$u$以$p_2,p_3,\cdots ,p_n$进制表示，这样就得到$n$个的等式；但是这样的进制表示并没有实质作用，因为我们很可能为了表示$u$，而将$k$拉到很大的数值，具体大小我们无从得知；所以我们用取模来解决这个问题，将等式两侧的大小变到我们能掌握的大小

我们在每个等式中取模$p_i^{q_i}$，将其转换为同余式，以$p_1$进制为例
$$u\equiv c_0+c_1{p}_1+c_2{p}_1^2+c_3{p}_1^3+\cdots +c_{q_1-1}{p}_1^{q_1-1}(\mathrm{mod}{p}_1^{q_1})$$
那么将原来的$n$个等式转换为$n$个同余式后，联立起来得到同余式组
$$\{\begin{array}{rl}& u\equiv c_0+c_1{p}_1+c_2{p}_1^2+c_3{p}_1^3+\cdots +c_{q_1-1}{p}_1^{q_1-1}(\mathrm{mod}{p}_1^{q_1})\\ & u\equiv c_0^{\prime }+c_1^{\prime }{p}_2+c_2^{\prime }{p}_2^2+c_3^{\prime }{p}_2^3+\cdots +c_{q_2-1}^{\prime }{p}_2^{q_2-1}(\mathrm{mod}{p}_2^{q_2})\\ & u\equiv c_0^{\prime \prime }+c_1^{\prime \prime }{p}_3+c_2^{\prime \prime }{p}_3^2+c_3^{\prime \prime }{p}_3^3+\cdots +c_{q_3-1}^{\prime \prime }{p}_3^{q_3-1}(\mathrm{mod}{p}_3^{q_3})\\ & \cdots \\ & u\equiv c_0^{\prime \prime \prime }+c_1^{\prime \prime \prime }{p}_n+c_2^{\prime \prime \prime }{p}_n^2+c_3^{\prime \prime \prime }{p}_n^3+\cdots +c_{q_n-1}^{\prime \prime \prime }{p}_n^{q_n-1}(\mathrm{mod}{p}_n^{q_n})\end{array}$$
分别对每个同余式求解$c_i$得到满足单个同余式的$u$，再对整个同余式组使用中国剩余定理得到最终满足所有同余式的解$u$

• 以解$u\equiv c_0+c_1{p}_1+c_2{p}_1^2+c_3{p}_1^3+\cdots +c_{q_1-1}{p}_1^{q_1-1}(\mathrm{mod}{p}_1^{q_1})$为例，

接下来为了将$p_1$进制表示的$u$应用于求解过程，我们构造
$$(g^u{)}^{\frac{p-1}{p_1^r}}\equiv a^{\frac{p-1}{p_1^r}}(\mathrm{mod}p)$$
其中$r$是我们设置的变量，后面会对其进行赋值；该式的模数应该也可以换作这里的$p_1^{q_1}$或者是其他满足条件的数

将$u$替换成$p_1$进制表示，得到
$$(g^{c_0+c_1{p}_1+c_2{p}_1^2+c_3{p}_1^3+\cdots +c_{q_1-1}{p}_{p_1}^{q_1-1}}{)}^{\frac{p-1}{p_1^r}}\equiv a^{\frac{p-1}{p_1^r}}(\mathrm{mod}p)$$
我们令$r$为$1$，得到
$$g^{c_0\frac{p-1}{p_1}}\cdot g^{c_1(p-1)}\cdot g^{c_2(p-1)p_1}\cdots g^{c_{q_1-1}(p-1)p_1^{q_1-2}}\equiv a^{\frac{p-1}{p_1}}(\mathrm{mod}p)$$
由费马小定理$g^{p-1}\equiv 1(\mathrm{mod}p)$，上式等价于
$$g^{c_0\frac{p-1}{p_1}}\equiv a^{\frac{p-1}{p_1}}(\mathrm{mod}p)$$
此时$c_0$的大小范围已知为$(0,p_1-1)$之间，其余变量大小已知，我们爆破$c_0$大小直到满足该同余式；这里就对应了为什么该算法应用条件需要$p-1$有较小的质因数

这样就可以得到$c_0$，但是我们要知道所有$c_i$的值来表示出在模$p_1^{q_1}$意义下的$u$

我们通过更改$r$的值，并应用费马小定理来达到求解$c_i$的目的，如我们再令$r$为$2$，得到
$$g^{c_0\frac{p-1}{p_1^2}}\cdot g^{c_1\frac{p-1}{p_1}}\cdot g^{c_2(p-1)}\cdots g^{c_{q_1-1}(p-1)p_1^{q_1-3}}\equiv a^{\frac{p-1}{p_1^2}}(\mathrm{mod}p)$$
由费马小定理$g^{p-1}\equiv 1(\mathrm{mod}p)$，上式等价于
$$g^{c_0\frac{p-1}{p_1^2}}\cdot g^{c_1\frac{p-1}{p_1}}\equiv a^{\frac{p-1}{p_1^2}}(\mathrm{mod}p)$$
其中同余式左侧第一项实际上我们已知，那么只用爆破$c_1$的大小即可；其他$c_i$以此类推

当我们求出整个同余式组里的满足各个同余式模意义下的$u$，使用孙子定理即可得到满足原始问题$g^u\equiv a(\mathrm{mod}p)$的$u$

那么用相同的方法再求解出$g^v\equiv b(\mathrm{mod}p)$的$v$即可，代回同余式$x\equiv v\cdot u^{-1}(\mathrm{mod}p-1)$求得$x$

看完整个流程，可能会想为什么不直接对$a^x\equiv b(\mathrm{mod}p)$使用之后将$x$转化为$p-1$的素因子的进制等等操作，这样就可以直接得到$x$，而不用拐弯抹角地去求$u,v$，再来求$x$；这是因为该算法只针对原根起效，也就是说，如果$a$是原根，那么这样操作是可以的；但如果$a$不是原根，就无法求得正确结果，其具体原因可以见$Reference$部分的文章

$PythonImplementation$

from Crypto.Util.number import *
from sympy.polys.galoistools import gf_crt
from sympy.polys.domains import ZZ
import gmpy2
import random

def Pohlig_Hellman(g, h, p):
    # 模数分解
    p_1 = p - 1
    d, factors = 2, []
    while d*d <= p_1:
        while (p_1 % d) == 0:
            factors.append(d)
            p_1 //= d
        d += 1
    if p_1 > 1:
        factors.append(p)
    factors = [[i, factors.count(i)] for i in set(factors)]

    # 求每个素因子进制下的c_i
    x = []
    for factor in factors:
        c_i_list = []
        for i in range(factor[1]):
            if i != 0:
                beta = (beta * pow(g, -(c_i_list[-1] * (factor[0] ** (i - 1))), p)) % p
            else:
                beta = h
            e1 = pow(beta, (p-1) // (factor[0] ** (i + 1)), p)
            e2 = pow(g, (p-1) // factor[0], p)
            for c_i in (range(factor[0])):
                if pow(e2, c_i, p) == e1:
                    c_i_list.append(c_i)
                    break
        x.append(c_i_list)

    # 将模p_i意义下的p_i进制表示转换为模p_i意义下的十进制
    system = []
    for i, factor in enumerate(factors):
        res = 0
        for j, x_j in enumerate(x[i]):
            res += x_j * (factor[0] ** j)
        res = res % (factor[0] ** factor[1])
        system.append(res)

    # 中国剩余定理
    factors = [factors[i][0] ** factors[i][1] for i in range(len(factors))]
    result = gf_crt(system, factors, ZZ)
    return result


if __name__ == "__main__":
    p = 7863166752583943287208453249445887802885958578827520225154826621191353388988908983484279021978114049838254701703424499688950361788140197906625796305008451719
    a = random.randint(0, 2 ** 512)
    x = random.randint(0, 2 ** 256)
    b = pow(a, x, p)
    print("real_x = {}".format(x))

    res = Pohlig_Hellman(a, b, p)
    print("x1 = {}".format(res))

# sage: p = 7863166752583943287208453249445887802885958578827520225154826621191353388988908983484279021978114049838254701703424499688950361788140197906625796305008451719
# sage: primitive_root(p)
# 13
    g = 13 # p的本原元为13
    u = Pohlig_Hellman(g, a, p)
    v = Pohlig_Hellman(g, b, p)
    try:
        x = gmpy2.invert(u, p - 1) * v % (p - 1)
        print("x2 = {}".format(x))
    except: # u 和 p-1 可能不互素，导致之间没有逆元；该情况将u除以最大公因数使得两者之间互素即可，最后结果除以相应的数即可
        i = 0 
        gcd = gmpy2.gcd(u, p - 1)
        while True:
            if gmpy2.gcd(u, p - 1) != 1:
                u = u // gmpy2.gcd(u, p - 1)      
                i += 1
            else:
                break
        assert gmpy2.gcd(u, p - 1) == 1
        x = (gmpy2.invert(u // gmpy2.gcd(u, p - 1), p - 1) * v) % (p - 1) // (gcd ** i)
        print("x = {}".format(x))

     
"""
a ^ x \equiv b (mod p) 之中a不是原根
real_x = 17475167858715014509948693871106677723065342839264165381680037187400995034209
x1 = 3931583376291971643604226624722943901442979289413760112577413310595676694494454509217307369704071534867821221958389972909818020158235480633350085553499260068
x = 17475167858715014509948693871106677723065342839264165381680037187400995034209
"""

"""
a ^ x \equiv b (mod p) 之中a是原根
real_x = 59538927048508916825466804038603833711975305674466217174924215808245351055236
x1 = 59538927048508916825466804038603833711975305674466217174924215808245351055236
x2 = 59538927048508916825466804038603833711975305674466217174924215808245351055236
"""

$Reference$

(10条消息) Pohlig-Hellman算法求解离散对数问题_国科大网安二班的博客-CSDN博客_pohlig-hellman

Pohlig_hellman/Pohlig_hellman.py at master · Sarapuce/Pohlig_hellman (github.com)

【应用代数】本原元和极小多项式 | Hope (leohope.com)