RSA--e与phi不互素

前言

我们知道，在RSA系统中，解密依赖于欧拉定理

$$a^{\phi(p)}\equiv1\;mod\;p\;,\;(a,p)=1\\$$

以及同余等式

$$ed\equiv1\;mod\;\phi(n)\\$$

,那么在已知

$$c\equiv\;m^{e}\;mod\;n$$

下，容易得到

$$c^{d}\equiv\;(m^{e})^{d}\;mod\;n\\$$

$$ed=1+k\phi(n)\\$$

$$c^{d}\equiv\;m*m^{k\phi(n)}\;\equiv\;m*1\;mod\;n\\$$

但是当$e,phi$不互素的情况下，私钥$d$就无法正常计算，在此之下想要顺利求解明文$m$,衍生了一系列问题，本文基于笔者日常中遇到的相关题型，给出对应的处理方法

从简单情况出发

我们不妨令

$$gcd(e,\phi)=t\\$$

$$te^{'}=e\\$$

那么考虑转化为公钥指数为$e^{'}$的RSA，有

$$e^{'}d^{'}\equiv1mod\phi(n)\\$$

$$c\equiv\;m^{e}\equiv\;m^{te^{'}}\;mod\;n\\$$

$$c^{d^{'}}\equiv\;m^{te^{'}d^{'}}\;\equiv\;m^{t}\;mod\;n\\$$

显然此时我们能够计算

$$c^{d^{'}}\;mod\;n\\$$

如果$t$不是很大（个位数），当$m^{t}<n$的时候，我们直接能得到

$$c^{d^{'}}\;mod\;n=m^{t}\\$$

，直接开$t$次方就能得到$m$

如果$m^{t}$不那么大的话，我们姑且还能尝试爆破

$$c^{d^{'}}=kn+m^{t}\\$$

，但是如果数据不支持，大到爆破不了，我们就得进一步观察已有数据中的代数结构了

如果$e$同时与$p-1$，$q-1$互质（CRT加速RSA问题）

这里也就是在RSA中使用CRT算法加速运算的例子，我们都来讲讲如何操作

为什么要用CRT？

在RSA的计算中，我们有

$$ed\equiv1\;mod\;\phi(n)\\$$

$$c\equiv\;m^d\;modn\\$$

，我们知道$n=pq$，$n$往往非常大，那么计算上面两个式子就非常的消耗时间，而且也需要$(e,\phi)=1$,我们注意到$p,q$也是两个大质数，很有可能$e$与即$p-1\;\;(\phi(p))$ ; 即$q-1\;\;(\phi(q))$也互素，那么我能就能把问题转化到模$p,q$ 下来讨论，下面是推导过程

$$由c\equiv\;m^{e}\;mod\;n\\$$

$$c\equiv\;m^e\;mod\;p\\$$

$$c\equiv\;m^e\;mod\;q\\$$

这是由于$q|n\;;p|n$也就是模的传递性，那么我们转过来计算在模$p,q$下的私钥$d_p;d_q$，即

$$edp\equiv1\;mod\;\phi(p)\\$$

$$ed_q\equiv1\;mod\;\phi(q)\\$$

那么在解密的时候就有

$$c^{d_p}\equiv\;m^{ed_p}\equiv\;m^{1+k\phi(p)}\;\equiv\;m\;mod\;p\\$$

$$c^{d_q}\equiv\;m^{ed_q}\equiv\;m^{1+k\phi(q)}\;\equiv\;m\;mod\;q\\$$

$$记c^{d_p}\equiv\;m\;mod\;p\;\;=m_1\;\;c^{d_q}\equiv\;m\;mod\;q\;\;=m_2\\$$

这里是因为$m$可能比$n$要大，所以我们无法直接解出来，然后我们根据中国剩余定理可以构造出一个对$m$的解，如下

由第二个同余式不妨设

$$m=m_2+hq\;\\$$

我们进而只需要解出来$h$,把这个式子带入另一个同余式

$$m_2+hq\equiv\;m_1mod\;p\\$$

$$hq\equiv(m_1-m_2)\;mod\;p\\$$

我们只需计算

$$I_q\equiv\;q^{-1}\;mod\;p\\$$

那么有

$$h\equiv\;I_q(m_1-m_2)\;mod\;p\\$$

得到$h$之后回代即可成功求出$m$

可见这里我们的计算都是在模$p,q$的式子中进行的，乘法运算的次数少了不少，实现了加速

说回正题，从上面不难看出，我们自始至终都没有尝试计算$ed\equiv1\;mod\;\phi(n)$，成功回避了这个问题

来个例子–黑盾杯 2020-Factor

n = 3454083680130687060405946528826790951695785465926614724373
e = 3
c = 1347530713288996422676156069761604101177635382955634367208
gcd(m, n) = 1

这里n不是很大，简单分解后能得到三个素因子，我们记为pqr，如下

p=11761833764528579549
q=17100682436035561357
r=17172929050033177661
print(n-p*q*r) # 0
print(gcd(e,phi)) # 3
print(gcd(e,p-1)) # 1
print(gcd(e,q-1)) # 3
print(gcd(e,r-1)) # 1
print(gcd(e,(p-1)*(r-1))) # 1

简单检查了关系，如上，这题方法很多，在上面的解法下即利用$(e,p-1)=(e,r-1)=1$来使用CRT求解，代码如下

#sage
dp=inverse(e,p-1)
dr=inverse(e,r-1)
m1=pow(c,dp,p)
m2=pow(c,dr,r)
Ir=pow(r,-1,p)
h=pow(Ir*(m1-m2),1,p)
m=m2+h*r
print(long_to_bytes(m))
# CMISCCTF{3_RSA}

这里也能转化到在模$pr$下尝试求解，但是那对明文的大小有一定要求，我们还是使用通法，这里只要能找到两个满足要求的因子就行，但是这不是什么时候都能生效，请看下题

强网杯2022–[ASR]

from Crypto.Util.number import getPrime
from secret import falg
pad = lambda s:s + bytes([(len(s)-1)%16+1]*((len(s)-1)%16+1))

n = getPrime(128)**2 * getPrime(128)**2 * getPrime(128)**2 * getPrime(128)**2
e = 3
flag = pad(flag)
print(flag)
assert(len(flag) >= 48)
m = int.from_bytes(flag,'big')
c = pow(m,e,n)
print(f'n = {n}')
print(f'e = {e}')
print(f'c = {c}')
n = 8250871280281573979365095715711359115372504458973444367083195431861307534563246537364248104106494598081988216584432003199198805753721448450911308558041115465900179230798939615583517756265557814710419157462721793864532239042758808298575522666358352726060578194045804198551989679722201244547561044646931280001
e = 3
c = 945272793717722090962030960824180726576357481511799904903841312265308706852971155205003971821843069272938250385935597609059700446530436381124650731751982419593070224310399320617914955227288662661442416421725698368791013785074809691867988444306279231013360024747585261790352627234450209996422862329513284149

注意到这里特别对m进行了填充，就是保证m比n的因子大，无法转化为在小因子的域下求解，看看下面这段失败的exp

from Crypto.Util.number import *
e = 3
n = 8250871280281573979365095715711359115372504458973444367083195431861307534563246537364248104106494598081988216584432003199198805753721448450911308558041115465900179230798939615583517756265557814710419157462721793864532239042758808298575522666358352726060578194045804198551989679722201244547561044646931280001
c = 945272793717722090962030960824180726576357481511799904903841312265308706852971155205003971821843069272938250385935597609059700446530436381124650731751982419593070224310399320617914955227288662661442416421725698368791013785074809691867988444306279231013360024747585261790352627234450209996422862329513284149

p = 225933944608558304529179430753170813347
q = 260594583349478633632570848336184053653
r = 218566259296037866647273372633238739089
t = 223213222467584072959434495118689164399

print(gcd(e,p-1)) # 3
print(gcd(e,q-1)) # 1
print(gcd(e,r-1)) # 3
print(gcd(e,t-1)) # 1
dt=inverse(e,t-1)
dq=inverse(e,q-1)
m1=pow(c,dt,t)
m2=pow(c,dq,q)
Iq=pow(q,-1,t)
h=pow(Iq*(m1-m2),1,t)
m=m2+h*q
print(long_to_bytes(m))

这里的确找到了满足$(a_i-1,e)=1$的两个因子，但是我们注意到，若$m=hq+m_2$又$m_2\in(0,q)$,所以$m\in(hq,(h+1)q)$我们在上述代码中解出来的$h$满足$h\in(0,t)$，也就是说，如果$m>(1+t)q$的时候，这个方法就暂时失效了，这是在多因子RSA下特别的问题，注意到我们还有两个因子没有使用，为了也能使用上它们，我们给出更加一般的方法

有限域开根

想想我们的核心目标是什么？–求解$m$ 那么我们为什么要找到这个$(a_i-1,e)=1$的$a_i$?重新捋一遍过程，我们就发现，是为了去计算$d_{a_i}$来求解满足$m^e\equiv\;c\;mod\;a_i$的$m\equiv\;c^{d_{a_i}}\;mod\;a_i$,换言之，我们还是在求模不同素因子下的$m$,我们不妨想一下，$d_{a_i}$真的重要吗？,我们看这个方程

$$m^e-c\equiv0\;mod\;a_i\\$$

这不就是在$GF(a_i)$下的方程吗？这里$e$还非常小，完全可以在这个有限域$GF(a_i)$下来开根求解所有符合题意的$m$的取值集合$M_i$，显然可能有多解，但是我们可以依次枚举$M_1,M_2,M_3,M_4$中的值，由上述分析，若$(a_k-1,e)=1$，那么$M_k=\left\{m\;mod\;a_k\right\}$否则也有$m\;mod\;a_j\;\in\;M_j$

那么我们枚举这四个集合，一定能找到满足我们需要的

$$\begin{cases}m_1=m\;mod\;a_1\\m_2=m\;mod\;a_2\\m_3=m\;mod\;a_3\\m_4=m\;mod\;a_4\\\end{cases}$$

这就是中国剩余定理了，我们直接打CRT就能解出来m，而且我们用到了所有的因子，几乎可以保证$m$这里是唯一的，如果不是唯一的，我们确定一下$n=p^{a}q^{b}r^{c}t^{d}$中的$a,b,c,d$，转化到$m$模$p^{a},\;q^{b},\;r^{c},\;t^{d}$下就行，下面是修改后的exp

# sage
from Crypto.Util.number import *
from itertools import product

# Given parameters
e = 3
n = 8250871280281573979365095715711359115372504458973444367083195431861307534563246537364248104106494598081988216584432003199198805753721448450911308558041115465900179230798939615583517756265557814710419157462721793864532239042758808298575522666358352726060578194045804198551989679722201244547561044646931280001
c = 945272793717722090962030960824180726576357481511799904903841312265308706852971155205003971821843069272938250385935597609059700446530436381124650731751982419593070224310399320617914955227288662661442416421725698368791013785074809691867988444306279231013360024747585261790352627234450209996422862329513284149

# Prime factors
p = 225933944608558304529179430753170813347
q = 260594583349478633632570848336184053653
r = 218566259296037866647273372633238739089
t = 223213222467584072959434495118689164399

# Create polynomial rings and find roots in each field
R.<x> = Zmod(p)[]
f = x^e - c
res_p = f.monic().roots()

R.<x> = Zmod(q)[]
f = x^e - c
res_q = f.monic().roots()

R.<x> = Zmod(r)[]
f = x^e - c
res_r = f.monic().roots()

R.<x> = Zmod(t)[]
f = x^e - c
res_t = f.monic().roots()

# Iterate through all possible combinations of roots
for (m_p, _) in res_p:
    for (m_q, _) in res_q:
        for (m_r, _) in res_r:
            for (m_t, _) in res_t:
                # Combine using CRT in stages
                # First combine p and q
                m_pq = crt(int(m_p), int(m_q), p, q)
                # Then combine r and t
                m_rt = crt(int(m_r), int(m_t), r, t)
                # Finally combine the two results
                m = crt(m_pq, m_rt, p*q, r*t)
                
                mes = long_to_bytes(m)
                print(mes)
        
# b'flag{Fear_can_hold_you_prisoner_Hope_can_set_you_free}\x06\x06\x06\x06\x06\x06'

在这个脚本下我们能解决大部分问题，但是当$e$太大的时候我们就得考虑其他的方法来处理了

根据题型，我们有几种应付的思路

UCSCCTF2025–EzCalculate

题面

p = 9586253455468582613875015189854230646329578628731744411408644831684238720919107792959420247980417763684885397749546095133107188260274536708721056484419031
q = 8998523072192453101232205847855618180700579235012899613083663121402246420191771909612939404791268078655630846054784775118256720627970477420936836352759291
n = p * q
e = 65536
c = 74962027356320017542746842438347279031419999636985213695851878703229715143667648659071242394028952959096683055640906478244974899784491598741415530787571499313545501736858104610426804890565497123850685161829628373760791083545457573498600656412030353579510452843445377415943924958414311373173951242344875240776

phi = (p - 1) * (q - 1)
print(gcd(e, phi)) # 4
print(gcd(e, p - 1)) # 2
print(gcd(e, q - 1)) # 2

print(gcd(e,(p-1)//2)) # 1
print(gcd(e,(q-1)//2)) # 1

在这里，我们注意到$e$同时与$\frac{p-1}{2}$和$\frac{q-1}{2}$互素，于是我们能有一个取巧的方法，先转化到对因子的同余式，即

$$c\equiv\;m^e\;mod\;p\\$$

$$c\equiv\;m^e\;mod\;q\\$$

再不妨记

$$s_p=\frac{p-1}{2}\\$$

$$s_q=\frac{q-1}{2}\\$$

容易有

$$(e,s_p)=(e,s_q)=1\\$$

那么必然存在

$$ed_{s_p}\equiv1\;mod\;s_p\\$$

$$ed_{s_q}\equiv1\;mod\;s_q\\$$

写成等式，即

$$ed_{s_q}=1+ks_q\\$$

$$ed_{s_p}=1+ks_p\\$$

我们考虑同余式

$$m^{\phi(p)}\equiv\;m^{p-1}\equiv\;(m^{s_p})^2\equiv1\;mod\;p$$

这里就初见端倪了，我们做一个换元

$$x_p=m^{s_p}$$

那不就是二次剩余的判别式吗

$$是的二次剩余不是的二次剩余m^{s_p}\;mod\;p\;\equiv\;\begin{cases}1\;\;;m^{s_p}\;p\\-1\;\;;m^{s_p}p\end{cases}$$

注意到这里我们得到了重要的模$p$式子，回看解密过程

$$c^{d_{s_p}}\equiv\;m^{ed_{s_p}}\equiv\;m^{1+k_{s_p}}\equiv\;m*(m^{s_p})^k$$

也就是

$$c^{d_{s_p}}\equiv\;\pm\;m\;mod\;p$$

也就是说我们有

$$m_1\equiv\;m\;\equiv\;\pm\;c^{d_{s_p}}\;mod\;p$$

回到了熟悉的地方，再构造CRT就可以解出$m$了

那么这里如果再改成

$$(e,p-1)=(e,q-1)=4,8,16...$$

呢，思路还是考虑这个式子

$$(m^{\frac{p-1}{2^i}})^{2^i}\equiv1\;mod\;p$$

可以考虑去计算$2^i$次的本原单位根，但考虑到虚数的引用带来的一系列不确定性，不再过度深入。

那么如果$e=65537$这种大素数呢？往往会伴随另外一个条件，是$e\;|\;\phi(N)$,在这种情况下，我们可以使用AMM算法来找到一个满足条件的$m$

AMM算法

我们先从上面的RSA题目中抽离出来，来看看整数域开根下的一般问题，我们面对的是同余式

$$m^e\;\equiv\;b\;mod\;p$$

我们往往已知$b,p,e$且$p$是一个素数，我们想要找出一个满足方程的$m$，从简单出发，我们先讨论$e=2$的情况

e=2时

也就是处理方程

$$m^2\;\equiv\;b\;mod\;p$$

因为$p$是素数，我们可以这样来设

$$p-1=2^t*s$$

其中为奇数$t>0,s$

那么我们知道，对于所有的模$p$的二次剩余数$x$ （这里$b$就是一个模$p$的二次剩余）有

$$x^{\frac{p-1}{2}}\equiv1\;mod\;p$$

代入上面的式子，可以得到

$$x^{s*2^{t-1}}\equiv1\;mod\;p$$

再次从简单出发，如果$t=1$

那么直接有

$$x^s\equiv1\;mod\;p$$

两边同乘一个$x$ ，又$(x,p)=1$，直接开方，有

$$x^{\frac{s+1}{2}}\equiv\;x^{\frac{1}{2}}\;mod\;p$$

前面强调了$b$就是这样的一个二次剩余，也就是说，我们代入$x=b$，这个方程是成立的，所以说

$$b^{\frac{s+1}{2}}\equiv\;b^{\frac{1}{2}}\;mod\;p$$

这个方程两边平方后，就是我们非常熟悉的

$$b^{s+1}\;\equiv\;b\;mod\;p$$

所以说，我们就找到了一个符合要求的$m=b^{\frac{s+1}{2}}$

那么，如果$t>1$呢？

我们知道对于所有的模$p$的非二次剩余$y$,有

$$y^{\frac{p-1}{2}}\equiv\;y^{s*2^{t-1}}\equiv-1\;mod\;p$$

我们不难得到

$$(x^{s*2^{t-2}})^2\equiv1\;mod\;p$$

开方，能得到两种结果，即

$$x^{s*2^{t-2}}\equiv1\;mod\;p\\$$

$$x^{s*2^{t-2}}\equiv-1\;mod\;p$$

为了避免负数的出现导致开根开出虚数，我们乘上一个非二次剩余，即

$$x^{s*2^{t-2}}\;\;*\;\;y^{s*2^{t-1}\;*\;k}\;\;\;\;\equiv1\;mod\;p$$

这里用$k$来记录是否为负，$y$是自己生成的

注意到我们这里$x$实际上就是$b$，而$s,t$都是已知的，在代码中我们可以直接算这个$x^{s*2^{t-2}}\;mod\;p$，如果是1，就令$k=0$，否则为1

同余式右边为1后继续开方，即

$$x^{s*2^{t-3}}\;\;*\;\;y^{s*2^{t-2}*k_1}\;\;*\;\;y^{s*2^{t-1}*k_2}\;\equiv\;1\;mod\;p$$

确定$k_2$的方法和上面一样，我们一直开方，直到$x$的指数幂中不含$t$,会得到形似如下的式子

$$x^{s}\;\;*\;\;y^{s*(k_1+2k_2+4k_3+...+2^{t-2}k_{t-1})}\;\equiv\;1\;mod\;p$$

两边同乘一个$x$，再开方就得到了我们需要的$m$

$$m\equiv\;b^{\frac{1}{2}}\;\equiv\;b^{\frac{s+1}{2}}\;\;*\;\;y^{\frac{1}{2}s*(k_1+2k_2+4k_3+...+2^{t-2}k_{t-1})}\;mod\;p$$

这里的核心思路就是消去$t$，构造

$$vb^{u}\equiv\;1\;mod\;p$$

的代数式

这里AMM找到的只是一个满足要求的解，我们知道在$e=2$的情况下，符合条件的$m$有两个，为了找到另一个解，也就是描述在模$p$下的$-1$,我们使用如下方法，任意找到一个$h$，使得它不是模$p$的二次剩余，那么有

$$h^{\frac{p-1}{2}}\equiv\;-1\;mod\;p\\$$

$$m*h^{\frac{p-1}{2}}\equiv\;-m\;mod\;p\\$$

$$(m*h^{\frac{p-1}{2}})^2\equiv\;(-m)^2\equiv\;b\;mod\;p\\$$

也就是说，$m$,$m*h^{\frac{p-1}{2}}$都是满足要求的解，在$GF(p)$下

e更大的时候

我们现在要考虑的是

$$m^e\equiv\;b\;mod\;p$$

这里的问题是从二次剩余转到$e$次剩余，这里我们需要条件$e\;|\;p-1$,设

$$p-1=e^t*s$$

那么对我们要求的$x$，有

$$(x^e)^{\frac{p-1}{e}}\equiv\;c^{\frac{p-1}{e}}\equiv\;c^{s*e^{t-1}}\;\equiv1\;mod\;p\\$$

还是从简单出发，如果$t=1$

即

$$c^{s}\equiv\;1\;mod\;p$$

此时如果$(s,e)=1$，那么可以参考RSA的解法

$$ed\equiv1\;mod\;s\;\;;\;\;ed-1=ks$$

代入，

$$c^{ed}\equiv\;(c^d)^e\;\equiv\;c\;mod\;p$$

可见这里$m\equiv\;c^d\;mod\;p$就是一个满足要求的解

如果$(s,e)\neq1$可以考虑离散对数的方法

如果$t>1$

我们站在高一点的视角下一个结论

AMM算法就是在模p整数域中构造一个e阶子群来开根

先看要讨论的式子

$$m^e\equiv\;c\;mod\;p$$

其中$e\;|\;p$

我们可以令

$$p-1\;=\;se^{t}$$

由欧拉定理,当$(x,p)=1$

$$x^{p-1}\equiv\;1\;mod\;p$$

这里可以视为一个生成元为$x$，阶为$p$的模$p$子群

那么稍微变形一下

$$(x^{\frac{p-1}{e}})^e\equiv\;1\;mod\;p$$

令

$$g=x^{\frac{p-1}{e}}$$

那么我们就构造了一个模$p$下，生成元为$g$的$e$阶子群$G$

$$G=\left\{1,g,g^2,...,g^{e-1}\right\}$$

特别注意$g$不能是模$p$的$e$次剩余，也就是

$$g^{\frac{p-1}{e}}\;\not\equiv\;\;1\;mod\;p$$

否则$g$就是1开$e$根开出来的一个结果，你拿$g$生成的群就无法包含所有可能的根

接下来我们说明$G$中元素进行$e$次幂后模$p$为1

$$g_i^{e}\equiv\;(g^{e})^i\equiv1^i\;\equiv\;1\;mod\;p$$

那么也就是说，在模$p$群下对1开$e$次方根的结果都会落在群$G$中

知道这个结论，我们对上式做一个小变形

$$(x^{s*e^{t-1}})^e\equiv1\;mod\;p$$

再

$$(x^e)^{s*e^{t-1}}\equiv(c)^{s*e^{t-1}}\equiv1\;mod\;p$$

$c,s,t,e,p$我们都知道，仿照上面的思路，先对

$$(c)^{s*e^{t-1}}$$

开$e$次方根，就能得到

$$(c)^{s*e^{t-2}}$$

我们再直接计算这个元素的值，然后在这个群$G$中找到逆元

假设

$$(c)^{s*e^{t-2}}\equiv\;g^i\;mod\;p$$

那么逆元就是

$$k_j\equiv\;g^{e-i}\;mod\;p$$

所以说上面的同余式开根后就直接是

$$(c)^{s*e^{t-2}}\;\;*\;\;k_j\;\;\equiv\;\;1\;\;mod\;p$$

如此往下计算，直到同余式左边只剩下$s$次幂，最后解出来$m$之后还要记得根据这个群$G$的$e$个单位根对应出$e$组解，当$(e,s)=1$的时候，我们来操作一下，就是

分解p-1

$$p-1=se^{t}$$

计算逆元d，找到初始解

$$d\equiv\;e^{-1}\;mod\;s$$

$$m_0\;\equiv\;c^{d}\;mod\;p$$

逐步上升

先找到生成元$g$

$$g^{\frac{p-1}{e}}\;\not\equiv\;\;1\;mod\;p$$

再定义一些初始变量

$$a\equiv\;g^{e^{t-1}}\;mod\;s$$

$$b\;\equiv\;c^{ed-1}\;\equiv\;\frac{(c^d)^e}{c}\;\equiv\;\frac{m_0^{e}}{c}\;mod\;p$$

$$z\;\equiv\;g^s\;mod\;p$$

$$h=1$$

可知

$$z^{e^{t}}\equiv\;g^{s*e^{t}}\;\equiv\;g^{p-1}\;\equiv\;1\;mod\;p$$

再在每一次迭代的过程中，计算

$$d_i\equiv\;b^{e^{t-1-i}}\;mod\;p$$

如果

$$d_i=1$$

说明这次开方的结果开出来就是$e$次剩余根，无需调整，否则要进行修改，并同时对$b,h,z$进行修正

最后对得到的根乘上$e$个单位根，得到解

# sage
def AM1M(residue, e, p):
    """
    返回 x^e ≡ residue mod p 的所有解
    """
    if (p-1) % e != 0:
        raise ValueError(f"e={e} 不能整除 p-1={p-1}，AMM 不适用")

    # 分解 p-1 = e^t * s
    t = valuation(p - 1, e)
    s = (p - 1) // (e^t)
    
    # 计算 d = e^{-1} mod s
    d = inverse_mod(e, s)
    
    # 简单情况：t=1
    if t == 1:
        x = power_mod(residue, d, p)
        F = GF(p)
        zeta = F.zeta(e)
        return [ (x * zeta**i) % p for i in range(e) ]
    
    
    # 寻找 e 次非剩余（生成元）
    F = GF(p)
    g = F.multiplicative_generator()
    noresidue = g
    while (noresidue^((p-1)//e)) == 1:
        noresidue *= g
    
    # 计算单位根和辅助变量
    a = power_mod(noresidue, e^(t-1) * s, p)
    b = power_mod(residue, e * d - 1, p)
    c = power_mod(noresidue, s, p)
    h = 1
    
    # Hensel 提升
    for i in range(1, t):
        d_i = power_mod(b, e^(t-1-i), p)
        if d_i == 1:
            j = 0
        else:
            j = -discrete_log(d_i, a)  
        b = (power_mod(c, e, p)^j * b) % p
        h = (power_mod(c, j, p) * h) % p
        c = power_mod(c, e, p)
    
    x0 = (power_mod(residue, d, p) * h) % p
    ω = F.zeta(e)  # e 阶单位根
    return [ (x0 * ω^i) % p for i in range(e) ]  # 返回所有解

这里能得到所有模p下的解，再整合起来打CRT就解决了

来个实战例子

NCTF2019–EzRSA

e = 0x1337
p = 199138677823743837339927520157607820029746574557746549094921488292877226509198315016018919385259781238148402833316033634968163276198999279327827901879426429664674358844084491830543271625147280950273934405879341438429171453002453838897458102128836690385604150324972907981960626767679153125735677417397078196059
q = 112213695905472142415221444515326532320352429478341683352811183503269676555434601229013679319423878238944956830244386653674413411658696751173844443394608246716053086226910581400528167848306119179879115809778793093611381764939789057524575349501163689452810148280625226541609383166347879832134495444706697124741
n = p * q
c = 10562302690541901187975815594605242014385201583329309191736952454310803387032252007244962585846519762051885640856082157060593829013572592812958261432327975138581784360302599265408134332094134880789013207382277849503344042487389850373487656200657856862096900860792273206447552132458430989534820256156021128891296387414689693952047302604774923411425863612316726417214819110981605912408620996068520823370069362751149060142640529571400977787330956486849449005402750224992048562898004309319577192693315658275912449198365737965570035264841782399978307388920681068646219895287752359564029778568376881425070363592696751183359

phi = (p - 1) * (q - 1)
print(gcd(e, phi)) #4919
print(gcd(e, p - 1)) #4919
print(gcd(e, q - 1)) # 4919
print(int(e)) #4919

exp 这里使用python多进程加快运算

from Crypto.Util.number import *
import itertools
from multiprocessing import Pool, cpu_count
from tqdm import tqdm
from sage.all import *

e = 0x1337
p = 199138677823743837339927520157607820029746574557746549094921488292877226509198315016018919385259781238148402833316033634968163276198999279327827901879426429664674358844084491830543271625147280950273934405879341438429171453002453838897458102128836690385604150324972907981960626767679153125735677417397078196059
q = 112213695905472142415221444515326532320352429478341683352811183503269676555434601229013679319423878238944956830244386653674413411658696751173844443394608246716053086226910581400528167848306119179879115809778793093611381764939789057524575349501163689452810148280625226541609383166347879832134495444706697124741
n = p * q
c = 10562302690541901187975815594605242014385201583329309191736952454310803387032252007244962585846519762051885640856082157060593829013572592812958261432327975138581784360302599265408134332094134880789013207382277849503344042487389850373487656200657856862096900860792273206447552132458430989534820256156021128891296387414689693952047302604774923411425863612316726417214819110981605912408620996068520823370069362751149060142640529571400977787330956486849449005402750224992048562898004309319577192693315658275912449198365737965570035264841782399978307388920681068646219895287752359564029778568376881425070363592696751183359

cp = c % p
cq = c % q

inv_p = pow(p, -1, q)
inv_q = pow(q, -1, p)

def AMM(residue, e, p):
    if (p-1) % e != 0:
        raise ValueError(f"e={e} does not divide p-1={p-1}")

    t = valuation(p - 1, e)
    s = (p - 1) // (e**t)
    
    d = inverse_mod(e, s)
    
    if t == 1:
        x = power_mod(residue, d, p)
        F = GF(p)
        zeta = F.zeta(e)
        return [ (x * zeta**i) % p for i in range(e) ]
    
    F = GF(p)
    g = F.multiplicative_generator()
    noresidue = g
    while (noresidue**((p-1)//e)) == 1:
        noresidue *= g
    
    a = power_mod(noresidue, e**(t-1) * s, p)
    b = power_mod(residue, e * d - 1, p)
    c = power_mod(noresidue, s, p)
    h = 1
    
    for i in range(1, t):
        d_i = power_mod(b, e**(t-1-i), p)
        if d_i == 1:
            j = 0
        else:
            j = -discrete_log(d_i, a)
        b = (power_mod(c, e, p)**j * b) % p
        h = (power_mod(c, j, p) * h) % p
        c = power_mod(c, e, p)
    
    x0 = (power_mod(residue, d, p) * h) % p
    ω = F.zeta(e)
    return [ (x0 * ω**i) % p for i in range(e) ]

def compute_crt(args):
    m1, m2 = args
    m1_int = int(m1)  # Convert to Python int
    m2_int = int(m2)  # Convert to Python int
    x = (m1_int * q * inv_q + m2_int * p * inv_p) % n
    return x

def find_flag():
    resp = AMM(cp, e, p)
    resq = AMM(cq, e, q)

    print(f"[+] Solutions in GF(p): {len(resp)}")
    print(f"[+] Solutions in GF(q): {len(resq)}")
    print(f"[+] Total combinations: {len(resp) * len(resq)}")

    with Pool(cpu_count()) as pool:
        combinations = itertools.product(resp, resq)
        for x in tqdm(
            pool.imap_unordered(compute_crt, combinations, chunksize=10000),
            total=len(resp) * len(resq),
            desc="Brute-forcing CRT"
        ):
            flag = long_to_bytes(x)
            if b'NCTF{' in flag:
                print(f"\n[+] Flag: {flag.decode()}")
                pool.terminate()
                return

    print("[-] Flag not found")

if __name__ == "__main__":
    find_flag()
    
    # 大约跑三分钟