中国剩余定理 (Chinese remainder theorem),定理中有中国二字,算是我国古代的文化瑰宝。

# 介绍

先引入一道著名的例题,出自《孙子算法》,叫做 “物不知数”。

1
有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?

意思很简单,就是要求出那个特定的数。

中国剩余定理具体的定义是这样:有以下一次线性同余方程组

{xa1(modm1)xa2(modm2)xan(modmn)\begin{cases} x \equiv {a_1} \pmod {m_1} \\ x \equiv {a_2} \pmod {m_2} \\ \quad \vdots \\ x \equiv {a_n} \pmod {m_n} \end{cases}

我的理解就是已知a1a2...an{a_1}、{a_2}、...、{a_n}m1m2...mn{m_1}、{m_2}、...、{m_n}。求出未知的xx

代入例题就是:

{x2(mod3)x3(mod5)x2(mod7)\begin{cases} x \equiv {2} \pmod {3} \\ x \equiv {3} \pmod {5} \\ x \equiv {2} \pmod {7} \end{cases}

定理的证明过程略去主要讲应用,这里直接写结论:

①首先定义N=m1×m2××mn=i=1nN={m_1} \times {m_2} \times \cdots \times {m_n}=\prod_{i=1}^{n}

②对于每一项中的 m 令{t_i} = \frac{N}

③ans=i=1nai×ti×=\sum_{i=1}^{n}{a_i}\times {t_i}\timesinv(ti,mi)+k×N({t_i},{m_i})+k \times N

回到原题,代入公式我们可求得:

ans=i=1nai×ti×inv(ti,N)=2×(5×7)×inv(35,3)+3×(3×7)×inv(21,5)+2×(3×5)×inv(15,7)=233=23+105×2ans=\sum_{i=1}^{n}{a_i}\times {t_i}\times inv({t_i},N) = 2 \times (5 \times 7) \times inv(35,3) + 3 \times (3 \times 7) \times inv(21,5) + 2 \times (3 \times 5) \times inv(15,7) = 233 = 23 + 105 \times 2

得到的答案为 23,代入验证也符合答案。这是中国剩余定理最简单的应用。转化为 python 语言也很方便,直接贴一个板子:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import gmpy2
from functools import reduce
ai = [2,3,2]
mi = [3,5,7]
def basic_CRT(ai,mi):
assert reduce(gmpy2.gcd,mi) == 1
assert len(ai) == len(mi)
N = reduce(lambda x,y:x * y,mi)
ans = 0
for a,m in zip(ai,mi):
t = N // m
ans += a * t * gmpy2.invert(t,m)
return ans % N,N
ans = basic_CRT(ai,mi)[0]
print(ans)

这是基础的 CRT,前提是模数 m 之间需要互质,否则逆元将无法求得。

如果涉及到模数 m 之间不互质,即 reduce (gcd,mi) != 1,则是扩展中国剩余定理的范畴。

以下内容参考 ruanx 的 blog

先从两个式子的例子讲起。设有下一次线性同余方程组:

{xa1(modm1)xa2(modm2)\begin{cases} x \equiv {a_1} \pmod {m_1} \\ x \equiv {a_2} \pmod {m_2} \end{cases}

其中 gcd (m1,m2)!=1

该方程组等价于:
x=a1+k1×m1=a2+k2×m2x={a_1} + {k_1} \times {m_1} = {a_2} + {k_2} \times {m_2}
其中 a,m 是已知量,令 gcd(m1,m2)=dp=m×d(m1,m2)=d,p=m \times d,再移项得到
a1a2=k2×m2k1×m1{a_1}-{a_2}={k_2} \times {m_2} - {k_1} \times {m_1}
a1a2d=k2×p2k1×p1\frac{a_1-a_2}{d}={k_2} \times {p_2} - {k_1} \times {p_1}
联系扩展欧几里得定理:
a,bZ,(x,y)Z\forall a,b \in Z,\quad \exists (x,y) \in Z

ax+by=gcd(a,b)ax+by=gcd(a,b)

换言之就是,对应任意的整数 a,b,都存在一对整数 x,y,使得 ax+by=gcd (a,b) 成立。通俗理解就是无论 a,b 怎么变,总有一组 (x,y) 满足该式子。

易得推论:

若存在:

ax+by=max+by=m

对于任意的 a,b 总存在 (x,y) 满足该等式,则有

ax+by=m=k×gcd(a,b)ax+by=m=k \times gcd(a,b)

(x,y)\exists(x,y),则一定有gcd(a,b)mgcd(a,b)|m

对应扩展欧几里得定理和该例子,对应一下:

ax+by=k×gcd(a,b)ax+by=k \times gcd(a,b)
p2×k2p1×k1=a1a2d×gcd(p1,p2){p_2} \times {k_2} - {p_1} \times {k_1}=\frac{a_1-a_2}{d} \times gcd(p_1,p_2)
gcd(p1,p2)=1gcd(p_1,p_2)=1

设通过扩展欧几里得定理得到了该方程的解 (λ1,λ2{\lambda_1},{\lambda_2}),即:

p2×λ2+p1×λ1=1{p_2} \times {\lambda_2} + {p_1} \times {\lambda_1} = 1

{\lambda}_2=\frac{d\times k_2}{a_1-a_2}\quad{\lambda}_1=-\frac{d\times k_1}

k2=a1a2d×λ2k1=a1a2d×k_2=\frac{a_1-a_2}{d}\times{\lambda_2}\quad k_1=-\frac{a_1-a_2}{d}\times

即我们得到了一组通解:

x=a1+k1×m1=a1a1a2dλ1x=a_1 + k_1 \times m_1 =a_1 -\frac{a_1-a_2}{d}{\lambda_1}

其中,a,m 为已知量,λ1{\lambda_1} 为可求量,即我们得到了一个 x 满足:

{xa1(modm1)xa2(modm2)\begin{cases} x \equiv {a_1} \pmod {m_1} \\ x \equiv {a_2} \pmod {m_2} \end{cases}

而整个解系为:

x=a1a1a2dλ1m1+lcm(m1,m2)x=a_1 -\frac{a_1-a_2}{d}{\lambda_1}{m_1}+lcm(m1,m2)

或者是理解为我们成功将两个模数不互素的一元线性同余方程合并成了一个方程:

x \equiv a_1 -\frac{a_1-a_2}{d}{\lambda}_1{m_1} \pmod

所以如果给了很多个方程组,只需将这些同余方程合并成一个,得到的就是最后的结果了。思路清晰,开始写板子:(看着 ruanx 大佬的 blog 学的所以写出来的板子都差不多)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import gmpy2
from functools import reduce
ai = [2,3,2]
mi = [3,5,7]
def union(x1, x2):
a1, m1 = x1
a2, m2 = x2
d = gmpy2.gcd(m1, m2)
assert (a2 - a1) % d == 0
p1,p2 = m1 // d,m2 // d
_,l1,l2 = gmpy2.gcdext(p1,p2)
k = -((a1 - a2) // d) * l1
lcm = gmpy2.lcm(m1,m2)
ans = (a1 + k * m1) % lcm
return ans,lcm


def excrt(ai,mi):
tmp = zip(ai,mi)
return reduce(union, tmp)
ans = excrt(ai,mi)[0]
print(ans)

# 例题

引入一道 ctf 的例题:

2022D.I.E&Trinity 战队新生赛 dpdqdr

附件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#coding:utf-8
import gmpy2
from Crypto.Util.number import *

p = getPrime(512)
q = getPrime(512)
r = getPrime(512)
e = getPrime(32)
n = p*q*r
phi = (p-1)*(q-1)*(r-1)
d = gmpy2.invert(e,phi)
dp = d%((q-1)*(r-1))
dq = d%((p-1)*(r-1))
dr = d%((p-1)*(q-1))
flag = 'flag{XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX}'+PADDING
m = bytes_to_long(flag.encode())
c = pow(m,e,n)

print(p)
print(q)
print(r)
print(dp)
print(dq)
print(dr)
print(c)

#12922128058767029848676385650461975663483632970994721128398090402671357430399910236576943902580268365115559040908171487273491136108931171215963673857907721
#10395910293559541454979782434227114401257890224810826672485874938639616819909368963527556812339196570118998080877100587760101646884011742783881592586607483
#8104533688439480164488403019957173637520526666352540480766865791142556044817828133446063428255474375204188144310967625626318466189746446739697284656837499
#73360412924315743410612858109886169233122608813546859531995431159702281180116580962235297605024326120716590757069707814371806343766956894408106019058184354279568525768909190843389534908163730972765221403797428735591146943727032277163147380538250142612444372315262195455266292156566943804557623319253942627829
#40011003982913118920477233564329052389422276107266243287367766124357736739027781899850422097218506350119257015460291153483339485727984512959771805645640899525080850525273304988145509506962755664208407488807873672040970416096459662677968243781070751482234692575943914243633982505045357475070019527351586080273
#21504040939112983125383942214187695383459556831904800061168077060846983552476434854825475457749096404504088696171780970907072305495623953811379179449789142049817703543458498244186699984858401903729236362439659600561895931051597248170420055792553353578915848063216831827095100173180270649367917678965552672673
#220428832901130282093087304800127910055992783874826238869471313726515822196746908777026147887315019800546695346099376727742597231512404648514329911088048902389321230640565683145565701498095660019604419213310866468276943241155853029934366950674139215056682438149221374543291202295130547776549069333898123270448986380025937093195496539532193583979030254746589985556996040224572481200667498253900563663950531345601763949337787268884688982469744380006435119997310653

这题直接套板子就行,当时没看出这个是 crt,还在模仿之前学过的 dpdq 泄露在推。

{ddp(mod(q1)×(r1))ddq(mod(p1)×(r1))ddr(mod(p1)×(q1))\begin{cases} d \equiv {dp} \pmod {(q-1)\times(r-1)} \\ d \equiv {dq} \pmod {(p-1)\times(r-1)} \\ d \equiv {dr} \pmod {(p-1)\times(q-1)} \end{cases}

就直接贴 exp 了:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
from Crypto.Util.number import *
import gmpy2
from functools import reduce
p = 12922128058767029848676385650461975663483632970994721128398090402671357430399910236576943902580268365115559040908171487273491136108931171215963673857907721
q = 10395910293559541454979782434227114401257890224810826672485874938639616819909368963527556812339196570118998080877100587760101646884011742783881592586607483
r = 8104533688439480164488403019957173637520526666352540480766865791142556044817828133446063428255474375204188144310967625626318466189746446739697284656837499
dp = 73360412924315743410612858109886169233122608813546859531995431159702281180116580962235297605024326120716590757069707814371806343766956894408106019058184354279568525768909190843389534908163730972765221403797428735591146943727032277163147380538250142612444372315262195455266292156566943804557623319253942627829
dq = 40011003982913118920477233564329052389422276107266243287367766124357736739027781899850422097218506350119257015460291153483339485727984512959771805645640899525080850525273304988145509506962755664208407488807873672040970416096459662677968243781070751482234692575943914243633982505045357475070019527351586080273
dr = 21504040939112983125383942214187695383459556831904800061168077060846983552476434854825475457749096404504088696171780970907072305495623953811379179449789142049817703543458498244186699984858401903729236362439659600561895931051597248170420055792553353578915848063216831827095100173180270649367917678965552672673
c = 220428832901130282093087304800127910055992783874826238869471313726515822196746908777026147887315019800546695346099376727742597231512404648514329911088048902389321230640565683145565701498095660019604419213310866468276943241155853029934366950674139215056682438149221374543291202295130547776549069333898123270448986380025937093195496539532193583979030254746589985556996040224572481200667498253900563663950531345601763949337787268884688982469744380006435119997310653

def union(x1, x2):
a1, m1 = x1
a2, m2 = x2
d = gmpy2.gcd(m1, m2)
assert (a2 - a1) % d == 0
p1,p2 = m1 // d,m2 // d
_,l1,l2 = gmpy2.gcdext(p1,p2)
k = -((a1 - a2) // d) * l1
lcm = gmpy2.lcm(m1,m2)
ans = (a1 + k * m1) % lcm
return ans,lcm


def excrt(ai,mi):
tmp = zip(ai,mi)
return reduce(union, tmp)

mi = [(q - 1) * (r - 1),(p - 1) * (r - 1),(p - 1) * (q - 1)]
ai = [dp,dq,dr]
d,lcm = excrt(ai,mi)
n = p * q * r
m = pow(c,d,n)
print(long_to_bytes(m))

总结一下就是,中国剩余定理应用于已知不同模数通过不同模数求出的不同密文的情况。本文就简单介绍了中国剩余定理及其用法,更深层的证明过程就没有涉及。

編集日 閲覧数

*~( ̄▽ ̄)~[お茶]を一杯ください

tsuppari Alipay

Alipay

tsuppari PayPal

PayPal