算法学习笔记(10): 中国剩余定理
中国剩余定理,也叫孙子定理,之所以叫这个名字,是因为《孙子算经》中有这样一个问题:
有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?
这个被叫做“物不知数”的问题本质上是解下面的同余方程组:
\begin{cases}x\equiv b_1 \pmod{a_1}\\x\equiv b_2 \pmod{a_2}\\...\\x\equiv b_n \pmod{a_n}\end{cases}
后来的数学家在研究中发现,这一方程组有解的一个充分条件是 a_1, \ a_2, \ ... \ a_n 两两互质,并用构造法给出了这种情况下方程的通解。而这种方法在算法竞赛中也常常会用到,例如这道模板题:
(洛谷P1495 曹冲养猪)
题目描述
自从曹冲搞定了大象以后,曹操就开始捉摸让儿子干些事业,于是派他到中原养猪场养猪,可是曹冲满不高兴,于是在工作中马马虎虎,有一次曹操想知道母猪的数量,于是曹冲想狠狠耍曹操一把。举个例子,假如有16头母猪,如果建了3个猪圈,剩下1头猪就没有地方安家了。如果建造了5个猪圈,但是仍然有1头猪没有地方去,然后如果建造了7个猪圈,还有2头没有地方去。你作为曹总的私人秘书理所当然要将准确的猪数报给曹总,你该怎么办?
输入格式
第一行包含一个整数n (n <= 10) – 建立猪圈的次数,解下来n行,每行两个整数ai, bi( bi <= ai <= 1000), 表示建立了ai个猪圈,有bi头猪没有去处。你可以假定ai,aj互质.
输出格式
输出包含一个正整数,即为曹冲至少养母猪的数目。
(本质上就是给“物不知数”套了个背景。)
我们从“物不知数”这个具体问题出发。(接下来一大波数学公式,请做好心理再看)
要想直接找到一个 x 使得方程组 \begin{cases}x\equiv 2\pmod{3}\\x\equiv 3\pmod{5}\\x\equiv 2\pmod{7}\end{cases} 成立当然是不容易的,但是要找到 n_1, \ n_2,\ n_3 使得 \begin{cases}n_1\equiv 2\pmod{3}\\n_2\equiv 3\pmod{5}\\n_3\equiv 2\pmod{7}\end{cases} 是相对容易的。
那么令 x=n_1+n_2+n_3 可以吗?那恐怕未必。在什么情况下 n_1\equiv 2\pmod{3} 可以推出 n_1+n_2\equiv 2\pmod{3} 呢?显然,那只有当 n_2 是3的倍数时成立。同理,要使 n_1+n_2+n_3 也符合前式,需要 n_2 和 n_3 都是3的倍数。
这样推下去,x=n_1+n_2+n_3符合方程组的条件是 n_1 是35的倍数, n_2 是21的倍数, n_3 是15的倍数。也就是说,现在我们只需要解三个同余方程 :
\begin{cases}35m_1\equiv 2\pmod{3}\\21m_2\equiv 3\pmod{5}\\15m_3\equiv 2\pmod{7}\end{cases} 。
注意到模数两两互质,则 \gcd(35,3)=\gcd(21,5)=\gcd(15,7)=1 ,所以我们可以用拓展欧几里得的方法解
\begin{cases}35w_1\equiv 1\pmod{3}\\21w_2\equiv 1\pmod{5}\\15w_3\equiv 1\pmod{7}\end{cases} 。(其实相当于求逆元)
解得 w_1=2,\ w_2=1,\ w_3=1 ,然后可得\begin{cases}m_1=2w_1=4\\ m_2=3w_2=3\\ m_3=2w_3=2\end{cases} ,于是 \begin{cases}n_1=35m_1=140\\ n_2=21m_2=63\\ n_3=15m_3=30\end{cases} 。
三者相加,即得一特解233(这里的233不是网络意义下的233,但我算出来不禁233了)。所有与233在模105意义下同余的数都是这个方程组的解,要求最小正数解只需对105取模即可,这里得出来是23。
现在我们把刚刚这个过程一般化。我们设 p=\prod_{i=1}^{n}a_i (即所有模数的乘积),并设 r_i=\frac{p}{a_i} (在“物不知数”中即为35、21和15)。于是 w_i=\left.\text{inv}\left(r_i\right)\right|_{a_i} (表示 r_i 在模 a_i 意义下的逆元), m_i=b_iw_i , 而n_i=r_im_i ,所有 n_i 相加即得 x 。
我们把以上这些综合成一个(看起来可能有点劝退的)公式即是:
x \equiv \sum_{i=1}^{n}b_i{r_i\left.\left[r_i\right]^{-1}\right|_{a_i}}\pmod{p}
现在我们来看(可能相对没那么劝退的)代码吧:
inline ll CRT(ll a[], ll b[], ll n) // a是模数数组,b是余数数组,n是数组长度
{
ll p = 1, x = 0;
for (int i = 0; i < n; ++i)
p *= a[i];
for (int i = 0; i < n; ++i)
{
ll r = p / a[i];
x += (b[i] * r * inv(r, a[i])) % p; // 逆元的求法参见上篇文章,或者下面有完整代码
}
return x % p;
}
这个函数返回的是符合方程组的最小正数解,一般要求的正是这个。
再附上曹冲养猪的完整AC代码:
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll exgcd(ll a, ll b, ll &x, ll &y)
{
if (b == 0)
{
x = 1;
y = 0;
return a;
}
ll d = exgcd(b, a % b, y, x);
y -= (a / b) * x;
return d;
}
inline ll inv(ll a, ll p)
{
ll x, y;
exgcd(a, p, x, y);
return (x % p + p) % p;
}
inline ll CRT(ll a[], ll b[], ll n)
{
ll p = 1, x = 0;
for (int i = 0; i < n; ++i)
p *= a[i];
for (int i = 0; i < n; ++i)
{
ll r = p / a[i];
x += (b[i] * r * inv(r, a[i])) % p;
}
return x % p;
}
int main()
{
ll n, a[10], b[10];
scanf("%lld", &n);
for (int i = 0; i < n; ++i)
scanf("%lld%lld", a + i, b + i);
printf("%lld\n", CRT(a, b, n));
return 0;
}