中国剩余定理

模数不互素的情况

当中国剩余定理在应用中模数不互素时，在求解MiM_iMi对模数mim_imi的逆元Mi′M_i'Mi′时可能会报错说：两个数不是互素的，无法进行求解逆元的操作（我这里用的封装函数gmpy2.invert()求解逆元）

那么假定现在有一组同余式需要用到中国剩余定理求解，但是模数不互素；

将式1：x≡c1(modm1)x\equiv c_1(mod~m_1)x≡c1(mod m1) 式2：x≡c2(modm2)x\equiv c_2(mod~m_2)x≡c2(mod m2) 转换为

x=c1+x1∗m1x = c_1 +x_1*m_1x=c1+x1∗m1，x=c2+x2∗m2x=c_2+x_2*m_2x=c2+x2∗m2（为什么这里乘数要写作xix_ixi呢，因为在之后的推导公式中xix_ixi会作为新的一轮值再次计算）

联立得

c1+x1∗m1=c2+x2∗m2c_1+x_1*m_1=c_2+x_2*m_2c1+x1∗m1=c2+x2∗m2

移项

x1∗m1=c2−c1+x2∗m2x_1*m_1=c_2-c_1+x_2*m_2x1∗m1=c2−c1+x2∗m2

转换为同余式形式

x1∗m1≡c2−c1(modm2)x_1*m_1\equiv c_2-c_1(mod~m_2)x1∗m1≡c2−c1(mod m2)

这里就类似于a∗x≡b(modm)a*x\equiv b(mod~m)a∗x≡b(mod m)（这里x1x_1x1是xxx，m1m_1m1是aaa）

所以该同余式要有解需要满足gcd(a,m)=1gcd(a,m)=1gcd(a,m)=1

所以求d=gcd(m1,m2)d=gcd(m_1,m_2)d=gcd(m1,m2)

所以同余式两边同时除以ddd

得到x1∗m1d≡c2−c1d(modm2d)x_1*\frac{m_1}{d}\equiv \frac{c_2-c_1}{d}(mod~\frac{m_2}{d})x1∗dm1≡dc2−c1(mod dm2)

那么求解x1≡(m1d)−1∗c2−c1d(modm2d)x_1 \equiv (\frac{m_1}{d})^{-1}*\frac{c_2-c_1}{d}(mod~\frac{m_2}{d})x1≡(dm1)−1∗dc2−c1(mod dm2)

这里有一点就是求m1d\frac{m_1}{d}dm1对模数m2d\frac{m_2}{d}dm2的逆元

可以这样来推导：

m1d∗(m1d)−1≡1(modm2d)\frac{m_1}{d}*(\frac{m_1}{d})^{-1}\equiv 1 (mod~\frac{m_2}{d})dm1∗(dm1)−1≡1(mod dm2)

相当于：

m1d∗(m1d)−1=1+k∗m2d\frac{m_1}{d}*(\frac{m_1}{d})^{-1}=1+k*\frac{m_2}{d}dm1∗(dm1)−1=1+k∗dm2

写到关于m1,m2m_1,m_2m1,m2的扩展欧几里得算法

k1∗m1+k2∗m2=dk_1*m_1+k_2*m_2=dk1∗m1+k2∗m2=d

同时除以ddd并移项

k1∗m1d=1−k2∗m2dk_1*\frac{m_1}{d}=1-k_2*\frac{m_2}{d}k1∗dm1=1−k2∗dm2

对比一下这两个等式，是不是项都一样呢？那么m1d\frac{m_1}{d}dm1的系数k1k_1k1实际上就是等于(m1d)−1(\frac{m_1}{d})^{-1}(dm1)−1

另外如果直接求解ta的逆元，可能会得到不正确的值~~（我试过但是不知道为什么）~~

那么x1≡k1∗c2−c1d(modm2d)x_1\equiv k_1*\frac{c_2-c_1}{d}(mod~\frac{m_2}{d})x1≡k1∗dc2−c1(mod dm2)

设X=k1∗c2−c1dX=k_1*\frac{c_2-c_1}{d}X=k1∗dc2−c1

则x1=m2d∗y+Xx_1=\frac{m_2}{d}*y+Xx1=dm2∗y+X

又因为x=c1+m1∗x1x=c_1+m_1*x_1x=c1+m1∗x1

代入后，

x=c1+m1∗m2d∗y+X∗m1x=c_1+\frac{m_1*m_2}{d}*y+X*m_1x=c1+dm1∗m2∗y+X∗m1

也就相当于同余式，

x≡c1+X∗m1(modm1∗m2d)x\equiv c_1+X*m_1(mod~\frac{m_1*m_2}{d})x≡c1+X∗m1(mod dm1∗m2)

那么这就是成功合并完了两个同余式的结果，成为一个新的同余式，其中c1+X∗m1c_1+X*m_1c1+X∗m1重新作为新的同余式里的ccc，m1∗m2d\frac{m_1*m_2}{d}dm1∗m2也是新的同余式里的mmm

另外这里m1∗m2d\frac{m_1*m_2}{d}dm1∗m2实际上也就是lcm[m1,m2]=m1∗m2gcd(m1,m2)lcm[m_1,m_2]=\frac{m_1*m_2}{gcd(m_1,m_2)}lcm[m1,m2]=gcd(m1,m2)m1∗m2

所以可以直接求m1,m2m_1,m_2m1,m2的最小公倍数去代替m1∗m2d\frac{m_1*m_2}{d}dm1∗m2

将这个新的同余式和剩下没有合并的同余式依次代入相同的计算中，就可以得到最后的结果

Python代码实现

i = 0  # 这里我使用了一个i作为全局变量
# 我定义的函数需要重复执行且在不同同余式组进行求解的时候i要重新赋值
def CRT(cl,ml):global iclist = clmoudlelist = mlif i == 0:m1,m2 = ml[i],ml[i+1]c1,c2 = cl[i],cl[i+1]else:m1,m2 = ml[i],ml[len(ml)-1]c1,c2 = cl[i],cl[len(cl)-1]d = gmpy2.gcd(m1,m2)print(d)c = c2 - c1assert c % d == 0# 这里代码的注释部分是最初的写法，但是结果会出错，我不知道为什么# X = c // d * gmpy2.invert(m1 // d,m2 // d)_,k1,k2 = gmpy2.gcdext(m1,m2) X = c // d * k1 * m1if i == 0:clist.append((X + c1) % gmpy2.lcm(m1,m2))# moudlelist.append(m1*m2 // d)moudlelist.append(gmpy2.lcm(m1,m2))i += 2else:clist[len(clist)-1] = (X + c1) % gmpy2.lcm(m1,m2)# moudlelist[len(moudlelist)-1] = m1*m2 // dmoudlelist[len(moudlelist)-1] = gmpy2.lcm(m1,m2)i += 1return clist,moudlelist

例题

e = 7
cl1 = [int(powmod(bytes_to_long(m1), e, x)) for x in ul]
cl2 = [int(powmod(bytes_to_long(m2), e, y)) for y in vl]
# ul,vl已知

很显然这里使用低加密指数广播攻击

但是在使用中国剩余定理的过程出现了报错：两数不互素

所以应用模数不互素的情况来求解

代码实现

import gmpy2
from Crypto.Util.number import *ml1 = [10537190383977432819948602717449313819513015810464463348450662860435011008001132238851729268032889296600248226221086420035262540732157097949791756421026015741477785995033447663038515248071740991264311479066137102975721041822067496462240009190564238288281272874966280,121723653124334943327337351369224143389428692536182586690052931548156177466437320964701609590004825981378294358781446032392886186351422728173975231719924841105480990927174913175897972732532233,1440176324831562539183617425199117363244429114385437232965257039323873256269894716229817484088631407074328498896710966713912857642565350306252498754145253802734893404773499918668829576304890397994277568525506501428687843547083479356423917301477033624346211335450]
ml2 = [168450500310972930707208583777353845862723614274337696968629340838437927919365973736431467737825931894403582133125917579196621697175572833671789075169621831768398654909584273636143519940165648838850012943578686057625415421266321405275952938776845012046586285747,1921455776649552079281304558665818887261070948261008212148121820969448652705855804423423681848341600084863078530401518931263150887409200101780191600802601105030806253998955929263882382004,25220695816897075916217095856631009012504127590059436393692101250418226097323331193222730091563032067314889286051745468263446649323295355350101318199942950223572194027189199046045156046295274639977052585768365501640340023356756783359924935106074017605019787]
cl1 = [2852589223779928796266540600421678790889067284911682578924216186052590393595645322161563386615512475256726384365091711034449682791268994623758937752874750918200961888997082477100811025721898720783666868623498246219677221106227660895519058631965055790709130207760704, 21115849906180139656310664607458425637670520081983248258984166026222898753505008904136688820075720411004158264138659762101873588583686473388951744733936769732617279649797085152057880233721961, 301899179092185964785847705166950181255677272294377823045011205035318463496682788289651177635341894308537787449148199583490117059526971759804426977947952721266880757177055335088777693134693713345640206540670123872210178680306100865355059146219281124303460105424]
cl2 = [148052450029409767056623510365366602228778431569288407577131980435074529632715014971133452626021226944632282479312378667353792117133452069972334169386837227285924011187035671874758901028719505163887789382835770664218045743465222788859258272826217869877607314144, 1643631850318055151946938381389671039738824953272816402371095118047179758846703070931850238668262625444826564833452294807110544441537830199752050040697440948146092723713661125309994275256, 10949587016016795940445976198460149258144635366996455598605244743540728764635947061037779912661207322820180541114179612916018317600403816027703391110922112311910900034442340387304006761589708943814396303183085858356961537279163175384848010568152485779372842]i = 0
def CRT(cl,ml):global iclist = clmoudlelist = mlif i == 0:m1,m2 = ml[i],ml[i+1]c1,c2 = cl[i],cl[i+1]else:m1,m2 = ml[i],ml[len(ml)-1]c1,c2 = cl[i],cl[len(cl)-1]d = gmpy2.gcd(m1,m2)print(d)c = c2 - c1assert c % d == 0# X = c // d * gmpy2.invert(m1 // d,m2 // d)_,k1,k2 = gmpy2.gcdext(m1,m2) X = c // d * k1 * m1if i == 0:clist.append((X + c1) % gmpy2.lcm(m1,m2))# moudlelist.append(m1*m2 // d)moudlelist.append(gmpy2.lcm(m1,m2))i += 2else:clist[len(clist)-1] = (X + c1) % gmpy2.lcm(m1,m2)# moudlelist[len(moudlelist)-1] = m1*m2 // dmoudlelist[len(moudlelist)-1] = gmpy2.lcm(m1,m2)i += 1return clist,moudlelist
ii = 3
temp_cl,temp_ml = CRT(cl1,ml1)
while len(cl1) - ii > 0:temp_cl,temp_ml = CRT(temp_cl,temp_ml)ii += 1
final_c = temp_cl[len(temp_cl)-1]
print(long_to_bytes(gmpy2.iroot(final_c,7)[0]))
ii = 3
i = 0
temp_cl,temp_ml = CRT(cl2,ml2)
while len(cl2) - ii > 0:temp_cl,temp_ml = CRT(temp_cl,temp_ml)ii += 1
final_c = temp_cl[len(temp_cl)-1]
print(long_to_bytes(gmpy2.iroot(final_c,7)[0]))

参考文章

(12条消息) 中国剩余定理算法详解（互质与不互质情况）_codeswarrior的博客-CSDN博客_中国剩余算法