模p^k的同余方程和离散对数求解

modular Equation

考虑求解多项式同余方程
f ( x ) = 0   m o d   p k , f ( x ) = 0   m o d   2 k f(x)=0\bmod{p^k},f(x)=0\bmod{2^k} f(x)=0modpk,f(x)=0mod2k其中 0 ≤ x < p k 0\le x\lt p^k 0≤x<pk，考虑将 x x x写成 p p p进制的形式
x = a 0 + a 1 p + . . . + a k − 1 p k − 1 , 0 ≤ a i < p x=a_0+a_1p+...+a_{k-1}p^{k-1},0\le a_i\lt p x=a0​+a1​p+...+ak−1​pk−1,0≤ai​<p首先两边同余 p p p，把 a 0 a_0 a0​视作未知数，得到 f ( a 0 ) = 0   m o d   p f(a_0)=0\bmod{p} f(a0​)=0modp，可以穷举得到 a 0 a_0 a0​；

下一次两边同余 p 2 p^2 p2，那么 f ( a 0 + a 1 p ) = 0   m o d   p 2 f(a_0+a_1p)=0\bmod{p^2} f(a0​+a1​p)=0modp2，将多项式展开得到 g ( a 1 ) = 0   m o d   p 2 g(a_1)=0\bmod{p^2} g(a1​)=0modp2，可以穷举得到 a 1 a_1 a1​；

同理，第 i i i次两边同余 p i p^i pi，此时 a 0 , . . . , a i − 2 a_0,...,a_{i-2} a0​,...,ai−2​均已知，那么
f ( a 0 + a 1 p + . . . + a i − 2 p i − 2 + a i − 1 p i − 1 ) = 0   m o d   p i f(a_0+a_1p+...+a_{i-2}p^{i-2}+a_{i-1}p^{i-1})=0\bmod{p^i} f(a0​+a1​p+...+ai−2​pi−2+ai−1​pi−1)=0modpi将多项式展开得到 g ( a i − 1 ) = 0   m o d   p i g(a_{i-1})=0\bmod{p^i} g(ai−1​)=0modpi，穷举得到解 a i − 1 a_{i-1} ai−1​；

以此类推得到所有 a 0 , a 1 , . . . , a k − 1 a_0,a_1,...,a_{k-1} a0​,a1​,...,ak−1​即可得到方程的解 x x x；

显然，方程的解不唯一；也就是说存在多条 ( a 0 , a 1 , . . . , a i − 1 ) (a_0,a_1,...,a_{i-1}) (a0​,a1​,...,ai−1​)路径，每次穷举得到一个 a i − 1 a_{i-1} ai−1​时，不可直接结束循环，而应该把所有这些 a i − 1 a_{i-1} ai−1​存起来；下一次分别从这些 a i − 1 a_{i-1} ai−1​出发寻找 a i a_i ai​，以此类推，得到所有满足满足方程 f ( x ) = 0   m o d   p k f(x)=0\bmod{p^k} f(x)=0modpk的解 x = a 0 + a 1 p + . . . + a i − 1 p k − 1 x=a_0+a_1p+...+a_{i-1}p^{k-1} x=a0​+a1​p+...+ai−1​pk−1；初步实现如下

def solve_equation_over_pk(fx, p, k): # f(x) = 0 mod p^k
    pre_possible = [0]  # 初值为0
    pi = p
    pi_ = 1
    for i in range(1, k + 1):  # 进行k次
        new_possible = []
        for x_ in pre_possible:  # 分别从上一次得到的所有a_{i-1}出发寻找a_i
            temp = deepcopy(x_)
            PR_ = PolynomialRing(Zmod(pi), 'y')
            y = PR_.gen()
            gx = fx.change_ring(Zmod(pi))
            gx = gx(temp + y * pi_)  # g(x+yp^{i-1})
            for an in range(p):  # 穷举
                if gx(an) == 0:  # 满足
                    temp += an * pi_  # 更新并添加
                    new_possible.append(temp)
                    temp = deepcopy(x_)  # 继续找
        pre_possible = new_possible
        pi *= p
        pi_ *= p
    return pre_possible # 返回所有可能的解

实际上，SageMath也有相关的方法解该方程（真方便啊），例如

p = 2
k = 100
m = pow(p, k)

x = var('x')
fx = 31635333913961551790218176796 * x ^ 2 + 343355432375781642567844278672 * x + 569420970791634503307822359156

# solve the equation f(x)=0 over m=p^k
solutions = solve_mod([fx == 0], m)  
for solution in solutions:
    print(int(fx(int(solution[0]))) % m == 0) # True

但是理解一下背后的原理并自己实现还是很有意义的。

Application

可以用于解决一类已知低位的问题，转化成多项式同余方程 f ( x ) = 0   m o d   p k f(x)=0\bmod{p^k} f(x)=0modpk，再解出低位；进而一般通过Coppersmith恢复完整目标数据，以下是两个例子。

LSB of p+q

考虑已知RSA模数 N = p × q N=p\times q N=p×q和 p + q p+q p+q的低位 s = ( p + q ) m o d    2 k s=(p+q)\mod 2^{k} s=(p+q)mod2k；求解 p p p和 q q q的低位，进一步通过一元Coppersmith分解模数 N N N
s = p + q   m o d   2 k    ⟹    s p = p 2 + N   m o d   2 k    ⟹    f ( x ) = x 2 − s x + N   m o d   2 k \begin{align*} &s=p+q\bmod{2^k}\\ \implies &sp=p^2+N\bmod{2^k}\\ \implies &f(x)=x^2-sx+N\bmod{2^k} \end{align*} ⟹⟹​s=p+qmod2ksp=p2+Nmod2kf(x)=x2−sx+Nmod2k​通过解该同余方程 f ( x ) f(x) f(x)，即可得到解 ( p   m o d   2 k ) (p\bmod{2^k}) (pmod2k)和 ( q   m o d   2 k ) (q\bmod{2^k}) (qmod2k)

def solve_equation_over_pk(fx, p, k):
    pre_possible = [0]  # 初值为0
    pi = p
    pi_ = 1
    for i in range(1, k + 1):  # 进行k次
        new_possible = []
        for x_ in pre_possible:  # 分别从上一次得到的所有a_{i-1}出发寻找a_i
            temp = deepcopy(x_)
            PR_ = PolynomialRing(Zmod(pi), 'y')
            y = PR_.gen()
            gx = fx.change_ring(Zmod(pi))
            gx = gx(temp + y * pi_)  # g(x+yp^{i-1})
            for an in range(p):  # 穷举
                if gx(an) == 0:  # 满足
                    temp += an * pi_  # 更新并添加
                    new_possible.append(temp)
                    temp = deepcopy(x_)  # 继续找
        pre_possible = new_possible
        pi *= p
        pi_ *= p
    return pre_possible

from Crypto.Util.number import getPrime

p = getPrime(512)
q = getPrime(512)

N = p * q
s = (p + q) % pow(2, 230)

PR = PolynomialRing(ZZ, 'x')
x = PR.gen()

fx = x * x - s * x + N
all_solutions = solve_equation_over_pk(fx, 2, 230)
print(p % pow(2, 230) in all_solutions)
print(q % pow(2, 230) in all_solutions)
print(len(all_solutions))
'''
True
True
8
'''

能够得到目标解，且解数并不多，那么分别尝试Coppersmith即可完成对模数 N N N的分解。当然如果是模 2 k 2^k 2k下，也可以直接剪枝求解，这里给出EXP；每轮筛选满足情况的节点，并记录进位；下一轮从这些节点出发直到找到所有节点，恢复目标解。

from Crypto.Util.number import *

p = getPrime(512)
q = getPrime(512)
n = p * q
k = 230
Lsb = (p + q) % pow(2, k)
p_Lsb = bin(p % pow(2, k))[2:].zfill(k)
q_Lsb = bin(q % pow(2, k))[2:].zfill(k)
# Lsb = (p+q)%pow(2,230)

pqs = [(0, '', '')]  # 进位up,p,q
Lsb_bin = bin(Lsb)[2:].zfill(k)
for b in range(k-1, -1, -1):
    pre_pqs = []
    for up, p, q in pqs:
        for i in range(2):
            for j in range(2):
                if (i + j + up) % 2 == int(Lsb_bin[b]):  # 满足1
                    p_ = str(i) + p
                    q_ = str(j) + q
                    if int(p_, 2) * int(q_, 2) % pow(2, k - b) == n % pow(2, k - b):  # 满足2
                        pre_pqs.append(((i + j + up) // 2, p_, q_))
    pqs = pre_pqs
pqs = set([i[1] for i in pqs])
print(p_Lsb in pqs)
print(q_Lsb in pqs)
print(len(pqs))
'''
True
True
8
'''

LSB of d

对于RSA，假设 N = p × q N=p\times q N=p×q， e d = 1   m o d   ( p − 1 ) ( q − 1 ) ed=1\bmod{(p-1)(q-1)} ed=1mod(p−1)(q−1)；通常来说针对于 e e e很小的情况
e d = 1 + k ( p − 1 ) ( q − 1 ) ed=1+k(p-1)(q-1) ed=1+k(p−1)(q−1)其中 d < ( p − 1 ) ( q − 1 ) , k < e d<(p-1)(q-1),k<e d<(p−1)(q−1),k<e；由于 e e e很小所以可以穷举得到 k k k；

假设已知 d 0 = ( d   m o d   2 a ) , s = ( p + q )   m o d   2 a d_0=(d\bmod{2^a}),s=(p+q)\bmod{2^a} d0​=(dmod2a),s=(p+q)mod2a；那么也对上式两边同余 2 a 2^a 2a；实际上就是得到 ( p + q ) (p+q) (p+q)的低位，继而得到 p p p和 q q q的低位。如下
e d 0 = 1 + k ( N − s + 1 )   m o d   2 a k ( N + p 2 ) = ( k N + k + 1 − e d 0 ) p   m o d   2 a f ( p ) = k p 2 − ( k N + k + 1 − e d 0 ) p + k N = 0   m o d   2 a ed_0=1+k(N-s+1)\bmod{2^a}\\ k(N+p^2)=(kN+k+1-ed_0)p\bmod{2^a}\\ f(p)=kp^2-(kN+k+1-ed_0)p+kN=0\bmod{2^a} ed0​=1+k(N−s+1)mod2ak(N+p2)=(kN+k+1−ed0​)pmod2af(p)=kp2−(kN+k+1−ed0​)p+kN=0mod2a

from Crypto.Util.number import getPrime, inverse

def solve_equation_over_pk(fx, p, k):
    pre_possible = [0]  # 初值为0
    pi = p
    pi_ = 1
    for i in range(1, k + 1):  # 进行k次
        new_possible = []
        for x_ in pre_possible:  # 分别从上一次得到的所有a_{i-1}出发寻找a_i
            temp = deepcopy(x_)
            PR_ = PolynomialRing(Zmod(pi), 'y')
            y = PR_.gen()
            gx = fx.change_ring(Zmod(pi))
            gx = gx(temp + y * pi_)  # g(x+yp^{i-1})
            for an in range(p):  # 穷举
                if gx(an) == 0:  # 满足
                    temp += an * pi_  # 更新并添加
                    new_possible.append(temp)
                    temp = deepcopy(x_)  # 继续找
        pre_possible = new_possible
        pi *= p
        pi_ *= p
    return pre_possible

p = getPrime(512)
q = getPrime(512)
N = p * q
e = 47
d = inverse(e, (p - 1) * (q - 1))
a = 230
d0 = d % pow(2, a)

PR = PolynomialRing(ZZ, 'x')
x = PR.gen()
all_sol = []
for k in range(1, e):
    fx = k * x ^ 2 - (k * N + k + 1 - e * d0) * x + k * N
    sol = solve_equation_over_pk(fx, 2, 230)
    all_sol += sol
print(p % pow(2, a) in all_sol)
print(q % pow(2, a) in all_sol)
'''
True
True
'''

当然纯爆破 k k k再求解会出现大量的解，如果有办法直接确定 k k k，那么就会快很多。

conclusion

如果仅仅是 f ( x ) = 0   m o d   m f(x)=0\bmod{m} f(x)=0modm，其中 m m m是一个素数的话，可以直接利用Sagemath的roots方法解该同余方程。

如果 m m m可以分解为 p 1 e 1 p 2 e 2 . . . p i e i {p_1}^{e_1}{p_2}^{e_2}...p_i^{e_i} p1​e1​p2​e2​...piei​​；则可以分别得到 x = x i   m o d   p i e i x=x_i\bmod{{p_i}^{e_i}} x=xi​modpi​ei​，再通过中国剩余定理恢复完整 ( x   m o d   m ) (x\bmod{m}) (xmodm)。
当然如果 m m m是一个未知分解的大合数，一般就只有考虑Coppersmith求小根，或者两个或多个多项式求最大公因式（存在同根，如RSA的相关明文攻击）

DLP

考虑离散对数问题
y = g x   m o d   p k , y = g x   m o d   2 k y=g^x\bmod{p^k},y=g^x\bmod{2^k} y=gxmodpk,y=gxmod2k首先 ϕ ( p k ) = p k − p k − 1 = ( p − 1 ) p k − 1 \phi(p^k)=p^k-p^{k-1}=(p-1)p^{k-1} ϕ(pk)=pk−pk−1=(p−1)pk−1；

假设 g g g是模 p k p^k pk的一个原根，那么 0 ≤ x < ϕ ( p k ) 0\le x\lt\phi(p^k) 0≤x<ϕ(pk)；同时注意欧拉定理
y = g x = g x    m o d     ϕ ( p k )   m o d   p k y=g^x=g^{x\,\bmod{\,\phi(p^k)}}\bmod{p^k} y=gx=gxmodϕ(pk)modpk将 x x x写为 p p p进制的形式
x = a 0 + a 1 p + . . . + a k − 2 p k − 2 + a k − 1 p k − 1 , 0 ≤ a i < p x=a_0+a_1p+...+a_{k-2}p^{k-2}+a_{k-1}p^{k-1},0\le a_i\lt p x=a0​+a1​p+...+ak−2​pk−2+ak−1​pk−1,0≤ai​<p首先，两边同时 ϕ ( p k ) / p \phi(p^k)/p ϕ(pk)/p次幂；（ 注意指数模 ϕ ( p k ) \phi(p^k) ϕ(pk) ）
y ϕ ( p k ) / p = g x ϕ ( p k ) / p = g a 0 ϕ ( p k ) / p   m o d   p k y ϕ ( p k ) / p = ( g ϕ ( p k ) / p ) a 0   m o d   p k y^{\phi(p^k)/p}=g^{x\phi(p^k)/p}=g^{a_0\phi(p^k)/p}\bmod{p^k{}}\\ y^{\phi(p^k)/p}=(g^{\phi(p^k)/p})^{a_0}\bmod{p^k} yϕ(pk)/p=gxϕ(pk)/p=ga0​ϕ(pk)/pmodpkyϕ(pk)/p=(gϕ(pk)/p)a0​modpk其中 0 ≤ a 0 < p 0\le a_0\lt p 0≤a0​<p，此时就可以通过穷举或者大步小步快速求出 a 0 a_0 a0​；

同理，下次两边同时 ϕ ( p k ) / p 2 \phi(p^k)/p^2 ϕ(pk)/p2次幂，那么
y ϕ ( p k ) / p 2 = g x ϕ ( p k ) / p 2 = g a 0 ϕ ( p k ) / p 2 + a 1 ϕ ( p k ) / p   m o d   p k y ϕ ( p k ) / p 2 ( g ϕ ( p k ) / p 2 ) − a 0 = ( g ϕ ( p k ) / p ) a 1   m o d   p k y^{\phi(p^k)/p^2}=g^{x\phi(p^k)/p^2}=g^{a_0\phi(p^k)/p^2+a_1\phi(p^k)/p}\bmod{p^k}\\ y^{\phi(p^k)/p^2}(g^{\phi(p^k)/p^2})^{-a_0}=(g^{\phi(p^k)/p})^{a_1}\bmod{p^k} yϕ(pk)/p2=gxϕ(pk)/p2=ga0​ϕ(pk)/p2+a1​ϕ(pk)/pmodpkyϕ(pk)/p2(gϕ(pk)/p2)−a0​=(gϕ(pk)/p)a1​modpk同理可以穷举或者大步小步得到 a 1 a_1 a1​；

第 i i i次时，两边同时 ϕ ( p k ) / p i \phi(p^k)/p^i ϕ(pk)/pi次幂，那么
y ϕ ( p k ) / p i = g x ϕ ( p k ) / p i = g a 0 ϕ ( p k ) / p i + a 1 ϕ ( p k ) / p i − 1 + . . . + a i − 1 ϕ ( p k ) / p   m o d   p k y ϕ ( p k ) / p i ( g ϕ ( p k ) / p i ) − ( a 0 + a 1 p + . . . + a i − 2 p i − 2 ) = ( g ϕ ( p k ) / p ) a i − 1   m o d   p k y^{\phi(p^k)/p^i}=g^{x\phi(p^k)/p^i}=g^{a_0\phi(p^k)/p^i+a_1\phi(p^k)/p^{i-1}+...+a_{i-1}\phi(p^k)/p} \bmod{p^k}\\ y^{\phi(p^k)/p^i}(g^{\phi(p^k)/p^i})^{-(a_0+a_1p+...+a_{i-2}p^{i-2})}=(g^{\phi(p^k)/p})^{a_{i-1}}\bmod{p^k} yϕ(pk)/pi=gxϕ(pk)/pi=ga0​ϕ(pk)/pi+a1​ϕ(pk)/pi−1+...+ai−1​ϕ(pk)/pmodpkyϕ(pk)/pi(gϕ(pk)/pi)−(a0​+a1​p+...+ai−2​pi−2)=(gϕ(pk)/p)ai−1​modpk其中 a 0 , a 1 , . . . , a i − 2 a_0,a_1,...,a_{i-2} a0​,a1​,...,ai−2​已知，那么就可以穷举得到 a i − 1 a_{i-1} ai−1​；

以此类推得到 a 0 , a 1 , . . . , a k − 2 a_0,a_1,...,a_{k-2} a0​,a1​,...,ak−2​；但由于 ϕ ( p k ) = ( p − 1 ) p k − 1 \phi(p^k)=(p-1)p^{k-1} ϕ(pk)=(p−1)pk−1，所以上述操作只能进行 k − 1 k-1 k−1次，第 k k k次时， ϕ ( p k ) / p k \phi(p^k)/p^k ϕ(pk)/pk是无意义的；也就是说 a k − 1 a_{k-1} ak−1​得单独求；没有太大的关系。写成
y = g a 0 + a 1 p + . . . + a k − 2 p k − 2 + a k − 1 p k − 1   m o d   p k y g − ( a 0 + a 1 p + . . . + a k − 2 p k − 2 ) = ( g p k − 1 ) a k − 1   m o d   p k y=g^{a_0+a_1p+...+a_{k-2}p^{k-2}+a_{k-1}p^{k-1}}\bmod{p^k}\\ yg^{-(a_0+a_1p+...+a_{k-2}p^{k-2})}=(g^{p^{k-1}})^{a_{k-1}}\bmod{p^k} y=ga0​+a1​p+...+ak−2​pk−2+ak−1​pk−1modpkyg−(a0​+a1​p+...+ak−2​pk−2)=(gpk−1)ak−1​modpk同理可求得 a k − 1 a_{k-1} ak−1​；才真正恢复 x = a 0 + a 1 p + . . . + a k − 1 p k − 1 x=a_0+a_1p+...+a_{k-1}p^{k-1} x=a0​+a1​p+...+ak−1​pk−1；完成DLP求解。

上述推导假设 g g g是模 p k p^k pk的一个原根，那么 ϕ ( p k ) \phi(p^k) ϕ(pk)是满足 g x ≡ 1   m o d   p k g^x\equiv 1\bmod{p^k} gx≡1modpk的最小正整数；当 g g g不是模 p k p^k pk的一个原根时，根据性质一定有阶 o r d p k ( g ) ∣ ϕ ( p k ) ord_{p^k}(g)|\phi(p^k) ordpk​(g)∣ϕ(pk)，那么阶也形如
o r d p k ( g ) = a p b ord_{p^k}(g)=ap^b ordpk​(g)=apb其中 a ∈ { 1 , p − 1 } , b ∈ { 1 , . . . , k − 1 } a\in \{1,p-1\},b\in \{1,...,k-1\} a∈{1,p−1},b∈{1,...,k−1}；同理可以按照上面的方式来进行求解。实际上如果不论 g g g的阶是多少，都采用 ϕ ( p k ) \phi(p^k) ϕ(pk)来算会得到多个解，这些解在模 o r d p k ( g ) ord_{p^k}(g) ordpk​(g)是同余的。因此进行统一化处理：采用 ϕ ( p k ) \phi(p^k) ϕ(pk)求出所有在模 ϕ ( p k ) \phi(p^k) ϕ(pk)意义下的解，那么目标解一定在其中，如果题目解 x < o r d p k ( g ) x\lt ord_{p^k}(g) x<ordpk​(g)；那么目标解一定是最小的那个解。

对于多解的处理同前，记录当前满足的所有节点，下一次从这些节点出发寻找下一个节点；一条完整的路径组成一个解。

from random import getrandbits, randrange
from gmpy2 import powmod, invert

def solve_DLP_over_pk(g, y, p, k):
    N = (p - 1) * pow(p, k - 1)
    m = pow(p,k)
    g_ = powmod(g, N // p, m)
    pre_possible = [0]
    pi = 1
    Ni = N // p
    for i in range(k - 1):
        tg = powmod(g, Ni, m)
        now_possible = []
        for x_ in pre_possible:
            temp = x_
            y_ = powmod(y, Ni, m) * invert(powmod(tg, temp, m), m) % m
            for ai in range(p):
                if powmod(g_, ai, m) == y_:
                    now_possible.append(temp + ai * pi)
        pi *= p
        Ni //= p
        pre_possible = now_possible

    solutions = []
    for temp in pre_possible:
        y_ = y * invert(powmod(g, temp, m), m) % m
        g_ = powmod(g, pi, m)
        for an in range(p):
            if y_ == powmod(g_, an, m):
                x_ = temp + an * pi
                assert powmod(g, x_, m) == y
                solutions.append(x_ % N)
    print(set(solutions))


p = 2
g = 3
k = 250
m = pow(p, k)
N = (p - 1) * pow(p, k - 1)
x = randrange(1, N)
print(x)
y = powmod(g, x, m)
solve_DLP_over_pk(g, y, p, k)

Pohlig Hellman

上述问题和Pohlig Hellman算法是一致的；假设 g g g是模素数 p p p的原根，那么 g g g在模 p p p下的阶 o r d p ( g ) = p − 1 ord_p(g)=p-1 ordp​(g)=p−1；如果 p − 1 p-1 p−1光滑，设最大因子为 l l l，则可以在 O ( l ) O(\sqrt{l}) O(l ​)下求解离散对数问题
y = g x   m o d   p , 0 ≤ x ≤ p − 1 p − 1 = p 1 e 1 p 2 e 2 . . . p k e k y=g^x\bmod{p},0\le x\le p-1\\ p-1={p_1}^{e_1}{p_2}^{e_2}...p_{k}^{e_k} y=gxmodp,0≤x≤p−1p−1=p1​e1​p2​e2​...pkek​​具体来说先求出
{ x = x 1   m o d   p 1 e 1 x = x 2   m o d   p 2 e 2 ⋯ x = x k   m o d   p k e k \begin{cases} x=x_1\bmod{p_1}^{e_1}\\ x=x_2\bmod{p_2}^{e_2}\\ \cdots\\ x=x_k\bmod{p_k}^{e_k} \end{cases} ⎩ ⎨ ⎧​x=x1​modp1​e1​x=x2​modp2​e2​⋯x=xk​modpk​ek​​再通过中国剩余定理恢复 ( x   m o d   p − 1 ) (x\bmod{p-1}) (xmodp−1)；以求解 x = x k   m o d   p k e k x=x_k\bmod{{p_k}^{e_k}} x=xk​modpk​ek​为例
y ( p − 1 ) / p k e k = g x ( p − 1 ) / p k e k = ( g ( p − 1 ) / p k e k ) x k   m o d   p y^{(p-1)/{p_k}^{e_k}}=g^{x(p-1)/p_k^{e_k}}=(g^{(p-1)/{p_k}^{e_k}})^{x_k}\bmod{p}\\ y(p−1)/pk​ek​=gx(p−1)/pkek​​=(g(p−1)/pk​ek​)xk​modp其中 g ′ = g ( p − 1 ) / p k e k g'=g^{(p-1)/{p_k}^{e_k}} g′=g(p−1)/pk​ek​的阶为 p k e k {p_k}^{e_k} pk​ek​，所以相当于指数 x k = x   m o d   p k e k x_k=x\bmod{{p_k}^{e_k}} xk​=xmodpk​ek​；即
y ′ = ( g ′ ) x k   m o d   p y'=(g')^{x_k}\bmod{p} y′=(g′)xk​modp将 x k x_k xk​记为 p k p_k pk​进制的形式
x k = a 0 + a 1 p k + . . . + a e k − 1 p k e k − 1 , a i ∈ [ 0 , p k − 1 ] x_k=a_0+a_1p_k+...+a_{e_k-1}{p_k}^{e_k-1},a_i\in [0,p_k-1] xk​=a0​+a1​pk​+...+aek​−1​pk​ek​−1,ai​∈[0,pk​−1]方法和前面类似，先两边同时 p k e k − 1 {p_k}^{e_k-1} pk​ek​−1次幂，那么
( y ′ ) p k e k − 1 = ( g ′ ) a 0   m o d   p (y')^{{p_k}^{e_k-1}}=(g')^{a_0}\bmod{p} (y′)pk​ek​−1=(g′)a0​modp可以通过穷举或者大步小步快速求出 a 0 a_0 a0​；以此类推恢复 x k x_k xk​。椭圆曲线上的Pohlig Hellman同理。总之都是在子群上求解小规模离散对数，再通过中国剩余定理得到目标解。

例题

DASCTF2024十二月赛【数论的香氛】

from sympy import isprime
from sympy.ntheory import legendre_symbol
import random
from Crypto.Util.number import bytes_to_long

k=79    #<-- i couldn't stress more

def get_p():
    global k
    while True:
        r=random.randint(2**69,2**70)
        p=2**k*r+1
        if isprime(p):
            return p
        else:
            continue

def get_q():
    while True:
        r=random.randint(2**147,2**148)
        q=4*r+3
        if isprime(q):
            return q
        else:
            continue


def get_y():
    global n,p,q
    while True:
        y=random.randint(0,n-1)
        if legendre_symbol(y,p)==1:
            continue
        elif legendre_symbol(y,q)==1:
            continue
        else:
            return y


flag=b'DASCTF{redacted:)}'
flag_pieces=[flag[0:10],flag[11:21],flag[22:32],flag[33:43],flag[44:]]
#assert int(bytes_to_long((flag_pieces[i] for i in range(5)))).bit_length()==k

p=get_p()
q=get_q()
n=p*q
print(f'{n=}')

y=get_y()
print(f'{y=}')


def encode(m):
    global y,n,k
    x = random.randint(1, n - 1)
    c=(pow(y,m,n)*pow(x,pow(2,k),n))%n
    return c

cs=[]
for i in range(len(flag_pieces)):
    ci=encode(bytes_to_long(flag_pieces[i]))
    cs.append(ci)

print(f'{cs=}')

'''
n=542799179636839492268900255776759322356188435185061417388485378278779491236741777034539347
y=304439269593920283890993394966761083993573819485737741439790516965458877720153847056020690
cs=[302425991290493703631236053387822220993687940015503176763298104925896002167283079926671604, 439984254328026142169547867847928383533091717170996198781543139431283836994276036750935235, 373508223748617252014658136131733110314734961216630099592116517373981480752966721942060039, 246328010831179104162852852153964748882971698116964204135222670606477985487691371234998588, 351248523787623958259846173184063420603640595997008994436503031978051069643436052471484545]
'''

参数 k = 79 k=79 k=79， 69 \text{69} 69比特随机数 r r r，生成素数 p = r 2 k + 1 p=r2^k+1 p=r2k+1； 147 \text{147} 147比特随机数 r r r，生成素数 q = 4 r + 3 q=4r+3 q=4r+3；生成随机数 y y y，确保勒让德符号 ( y p ) = ( y q ) = − 1 \binom{y}{p}=\binom{y}{q}=-1 (py​)=(qy​)=−1；根据欧拉判别式
( y p ) = y p − 1 2 = − 1   m o d   p y p − 1 = 1   m o d   p \binom{y}{p}=y^{\frac{p-1}{2}}=-1\bmod{p}\\ y^{p-1}=1\bmod{p} (py​)=y2p−1​=−1modpyp−1=1modp即 y y y是模 p p p下的原根。模数 n = p × q n=p\times q n=p×q，对于每一个flag分片 m i m_i mi​， m i < 2 k m_i\lt 2^{k} mi​<2k；生成随机数 x x x，加密为
c i = y m i × x 2 k   m o d   n c_i=y^{m_i}\times x^{2^k}\bmod{n} ci​=ymi​×x2kmodn首先对于 p = r 2 k + 1 p=r2^k+1 p=r2k+1， r r r只有64比特比较小，所以一元Coppersmith求解
f ( r ) = r 2 k + 1 = 0   m o d   p f(r)=r2^k+1=0\bmod{p} f(r)=r2k+1=0modp对于加密式，两边同余 p p p
c i = y m i × x 2 k   m o d   p c_i=y^{m_i}\times x^{2^k}\bmod{p} ci​=ymi​×x2kmodp其中 x x x随机数未知，所以设法消去， ϕ ( p ) = p − 1 = r 2 k \phi(p)=p-1=r2^k ϕ(p)=p−1=r2k，两边同时 r r r次幂
c i r = y r m i x 2 k r = ( y r ) m i   m o d   p {c_i}^r=y^{rm_i}x^{2^kr}=(y^r)^{m_i}\bmod{p} ci​r=yrmi​x2kr=(yr)mi​modp此时 y r y^r yr的阶为 2 k 2^k 2k，而 m i < 2 k m_i\lt 2^k mi​<2k；所以即为求解离散对数
y = g x   m o d   p y = c i r   m o d   p , g = y r   m o d   p , x < 2 k , o r d g ( p ) = 2 k y=g^x\bmod{p}\\ y={c_i}^r\bmod{p},g=y^r\bmod{p},x\lt 2^k,ord_g(p)=2^k y=gxmodpy=ci​rmodp,g=yrmodp,x<2k,ordg​(p)=2k将 x x x写成二进制形式
x = a 0 + a 1 2 + . . . + a k − 1 2 k − 1 x=a_0+a_12+...+a_{k-1}2^{k-1} x=a0​+a1​2+...+ak−1​2k−1然后方法同前，不再赘述。当然也可以直接用sagemath的discrete_log方法直接求解。

n = 
y = 
cs = 
k = 79
R.< x > = PolynomialRing(Zmod(n))

f = x * 2 ^ k + 1
f = f.monic()
r = f.small_roots(X=2 ^ 70, beta=0.45, epsilon=0.02)[0]
r = int(r)
p = gcd(int(f(r)), n)  # 分解得到n

from gmpy2 import invert,powmod
from Crypto.Util.number import long_to_bytes

g = powmod(y, r, p)
for ci in cs:
    yi = powmod(ci, r, p)  # pow(g,xi,p) == yi where ord(g) = 2^k
    Ni = pow(2, k - 1)
    g_ = powmod(g, Ni, p)
    pi = 1

    x_ = 0
    for _ in range(k):
        tg = powmod(g, Ni, p)
        temp_y = powmod(yi, Ni, p) * powmod(tg, -x_, p) % p

        # for an in range(2): # an=0 or 1
        #    if powmod(g_, an, p) == temp_y:
        #        x_ += an * pi

        if temp_y != 1:
            x_ += pi
        pi *= 2
        Ni //= 2
    print(long_to_bytes(x_))