中国剩余定理 与 扩展中国剩余定理
中国剩余定理 与 扩展中国剩余定理
一、问题引入
先看一道\(poj\)上的题目:【\(POJ1006\)】 \(Biorhythms\)
题意:
人自出生起就有体力,情感和智力三个生理周期,分别为\(23\),\(28\)和\(33\)天。一个周期内有一天为峰值,在这一天,人在对应的方面(体力,情感或智力)表现最好。通常这三个周期的峰值不会是同一天。现在给出三个日期,分别对应于体力,情感,智力出现峰值的日期。然后再给出一个 起始日期\(d\),要求从这一天开始,算出最少再过多少天后三个峰值同时出现。
分析:
首先我们知道,任意两个峰值之间一定相距整数倍的周期。假设一年的第\(N\)天达到峰值,则 下次 达到峰值的时间为\(N+kT\)(\(T\)是周期,\(k\)是任意正整数)。
所以,三个峰值同时出现的那一天(\(S\))应满足
\(N_1,N_2,N_3\)分别为为体力,情感,智力出现峰值的日期, \(T_1,T_2,T_3\) 分别为体力,情感,智力周期。 需要求出\(k_1,k_2,k_3\)三个非负整数使上面的等式成立。
想直接求出\(k_1,k_2,k_3\)很难,但是目的是求出\(S\), 可以考虑 从结果逆推。根据上面的等式,\(S\) 满足三个要求:
- 除以\(T_1\)余数为\(N_1\)
- 除以\(T_2\)余数为\(N_2\)
- 除以\(T_3\)余数为\(N_3\)
找到最小的一个这样的数就可以了。
这就是著名的 中国剩余定理,老祖宗在几千年前已经对这个问题想出了一个精妙的解法。依据此解法的算法,时间复杂度可达到\(O(1)\)。
二、中国剩余定理
在《孙子算经》中有这样一个问题:“今有物不知其数,三三数之剩二(除以\(3\)余\(2\)),五五数之剩三(除以\(5\)余\(3\)),七七数之剩二(除以\(7\)余\(2\)),问物几何?”这个问题称为 孙子问题 ,该问题的一般解法国际上称为 中国剩余定理。
解法共三步
-
① 找三个数:
- 从\(3\)和\(5\)的公倍数中找出被\(7\)除余\(1\)的最小数\(15\)
- 从\(3\)和\(7\)的公倍数中找出被\(5\)除余\(1\)的最小数\(21\)
- 从\(5\)和\(7\)的公倍数中找出被\(3\)除余\(1\)的最小数\(70\)
-
② 找符合条件的一个数:
- 用\(15\)乘以\(2\)(\(2\)为最终结果除以\(7\)的余数)
- 用\(21\)乘以\(3\)(\(3\)为最终结果除以\(5\)的余数)
- 用\(70\)乘以\(2\)(\(2\)为最终结果除以\(3\)的余数)
然后把三个乘积相加\(15∗2+21∗3+70∗2=233\)
-
③ 用\(233\)除以\(3,5,7\)三个数的最小公倍数\(105\),得到余数,即\(233\%105=23\)
答:\(23\)就是符合条件的最小数
\(Q\):这有数学依据吗?
现在将 孙子问题 拆分成几个简单的小问题,从零开始,揣测古人是如何推导出这个解法的:
首先,我们假设:
- \(n_1\)是满足除以\(3\)余\(2\)的一个数,比如\(2,5,8\)等等,也就是满足\(3∗k_1+2(k_1>=0)\) 的一个任意数。
- \(n_2\)是满足除以\(5\)余\(3\)的一个数,比如\(3,8,13\)等等,也就是满足\(5∗k_2+3(k_2>=0)\) 的一个任意数。
- \(n_3\)是满足除以\(7\)余\(2\)的一个数,比如\(2,9,16\)等等,也就是满足\(7∗k_3+2(k_3>=0)\) 的一个任意数。
有了前面的假设,我们先从\(n_1\)这个角度出发,已知\(n_1\)满足除以\(3\)余\(2\),能不能使得\(n_1+n_2\)的和仍然满足除以\(3\)余\(2\),进而使得\(n_1+n_2+n_3\)的和仍然满足除以\(3\)余\(2\)?
这就牵涉到一个 最基本数学定理,如果有\(a\%b=c\),则有\((a+k∗b)\%b=c\)(\(k\)为非零整数),换句话说,如果一个除法运算的余数为\(c\),那么被除数与\(k\)倍的除数相加(或相减)的和(差)再与除数相除,余数不变。
以此定理为依据,如果\(n_2\)是\(3\)的倍数,\(n_1+n_2\)就依然满足除以\(3\)余\(2\)。同理,如果\(n_3\)也是\(3\)的倍数,那么\(n_1+n_2+n_3\)的和就满足除以\(3\)余\(2\)。这是从\(n_1\)的角度考虑的,再从\(n_2\),\(n_3\)的角度出发,我们可推导出以下三点:
-
为使\(n_1+n_2+n_3\)的和满足除以\(3\)余\(2\),\(n_2\)和\(n_3\)必须是\(3\)的倍数
-
为使\(n_1+n_2+n_3\)的和满足除以\(5\)余\(3\),\(n_1\)和\(n_3\)必须是\(5\)的倍数
-
为使\(n_1+n_2+n_3\)的和满足除以\(7\)余\(2\),\(n_1\)和\(n_2\)必须是\(7\)的倍数
因此,为使\(n_1+n_2+n_3\)的和作为 孙子问题 的一个最终解,需满足:
- \(n_1\)除以\(3\)余\(2\),且是\(5\)和\(7\)的公倍数
- \(n_2\)除以\(5\)余\(3\),且是\(3\)和\(7\)的公倍数
- \(n_3\)除以\(7\)余\(2\),且是\(3\)和\(5\)的公倍数
所以,孙子问题解法的 本质 是:
- 从\(5\)和\(7\)的公倍数中找一个除以\(3\)余\(2\)的数\(n_1\)
- 从\(3\)和\(7\)的公倍数中找一个除以\(5\)余\(3\)的数\(n_2\)
- 从\(3\)和\(5\)的公倍数中找一个除以\(7\)余\(2\)的数\(n_3\)
再将三个数相加得到解。在求\(n_1,n_2,n_3\)时又用了一个 小技巧,以\(n_1\)为例,并非从\(5\)和\(7\)的公倍数中直接找一个除以\(3\)余\(2\)的数,而是
先找一个 除以\(3\)余\(1\)的数,再乘以\(2\)。也就是先求出\(5\)和\(7\)的公倍数模\(3\)下的逆元 ,再用逆元去乘余数!
\(Q\):为啥这样干呢?
答:不这样干不好整啊,\(5\)和\(7\)的公倍数可不少,一个个暴力去找太累了。不这样暴力还能咋办?
(a / b) % p = (a%p / b%p) %p (错)
因为除法没法进行取模运算,需要找到乘法逆元。
\(Q\):哪里来的除法,哪里来的取模?
答:\(M_i=M/m_i\)就是除法,需要对\(m_i\)取模。
这里又有一个数学公式,如果\(a\%b=c\),那么\((a∗k)\%b=a\%b+a\%b+…+a\%b=c+c+…+c=k∗c(k>0)\),也就是说,如果一个除法的余数为\(c\),那么被除数的\(k\)倍与除数相除的余数为\(k∗c\),展开式中已证明。
最后,我们还要清楚一点,\(n_1+n_2+n_3\)只是问题的一个解,并不是最小的解。如何得到最小解?我们只需要从中最大限度的减掉\(3,5,7\)的公倍数\(105\)即可。道理就是前面讲过的定理 如果\(a\%b=c\),则有\((a−k∗b)\%b=c\) ,所以\((n_1+n_2+n_3)\%105\)就是最终的最小解。
【中国剩余定理】
设\(x \in Z\),有
令$ \large \left{\begin{matrix}
A=a_1a_2a_3...a_n & \
A_i=A/a_i& \
t_i 是A_i关于a_i 的逆元,即A_it_i \equiv 1 (mod \ a_i)
\end{matrix}\right.$
则
【扩展欧几里得求逆元】
扩展欧几里得算法,顾名思义,是欧几里得算法的扩展,这意味着该算法不仅保留了欧几里得算法的核心——求最大公约数,还具备欧几里得算法没有的功能——求解贝祖等式,即可以求出方程\(ax+by=c\)的一组解\((x_0,y_0)\)
特别的,当 \(gcd(a,b)=1\) ,即\(a,b\)互素时,\(ax+by=1\)中的解
- \(x\) 是\(a\)模\(b\)的乘法逆元,即\(a * x \equiv 1(mod \ b)\)
- \(y\) 是\(b\)模\(a\)的乘法逆元,即\(b * x \equiv 1(mod \ a)\)
解释:方程\(ax+by=1\)写成同余方程,就是:\(a * x \equiv 1(mod \ b)\),也就是逆元的定义式
【r = (r % a + a) % a
】
作用:将余数加上若干个循环节,使得结果落在\(0 \sim a-1\)之间。
举栗子:
按\(5\)个一组划分,最后缺少\(2\)个
本质:就是多了\(3\)个,\(r\)其实是\(3\)
\(Q\):那就直接\(r+a=3\)不就行了吗,\((r\%a+a)\%a\)是什么意思?
\(A:(1)\)如果\(r>=0, r+a>=a\)了,这时应该再多分出去一组才对,不能直接\(r+a\) ,但 \((r\%a+a)\%a\)就没有问题了
$ (2)\(如果\)r<0,(r%a+a)%a$也一样正确
结论:公式 \((r\%a+a)\%a\)可以把负的余数改为正的余数
#include <bits/stdc++.h>
using namespace std;
int main() {
int r = -2, a = 5;
r = (r % a + a) % a;
cout << r << endl;
// 输出3
return 0;
}
三、扩展中国剩余定理
\(Q\):中国剩余定理要求所有的\(m_i\)互质,那么如果不互质,怎么求解同余方程组?
答:对于模数不互质的情况,我们可以通过合并相邻两组方程的方式来逐步得到解。
先计算两个线性同余方程组的解,之后依此类推
① 可以写成\(x=k_1*a_1+m_1\) ③
② 可以写成\(x=k_2*a_2+m_2\) ④
③=④ $ \rightarrow $ \(k_1*a_1+m_1=k_2*a_2+m_2\) $ \rightarrow$ \(k_1*a_1-k_2*a_2=m_2-m_1\) ⑤
根据 裴蜀定理,当
\(m_2-m_1\)为 \(gcd(a_1,-a_2)\)的倍数时,方程 ⑤ 有无穷多组整数解
而\(k_1*a_1-k_2*a_2=gcd(a_1,-a_2)=d\),可以用 拓展欧几里得算法 来求解出一组特解(\(k_1',k_2'\)),即\(exgcd(a_1,-a_2,k_1',k_2')\)
如何求解\(k_1*a_1-k_2*a_2=gcd(a_1,-a_2)=d\)中的通解\(k_1,k_2\)呢?
结论:\(k_1,k_2\)的通解
证明:
\(k_1\)的通解为\(s_1\),\(k_2\)的通解为\(s_2\),\(k_1\)的特解为\(k_1'\),\(k_2\)的特解为\(k_2'\)
则
\(\because\) \(gcd(a_1,-a_2)=d\),\(\therefore \frac{a_1}{d}\)与\(\frac{a_2}{d}\)互质。
解释:最大公约数嘛,如果除完了还有约数,那还是什么最大公约数~
\(\therefore \frac{a_2}{d} \mid s_1-k_1', \frac{a_1}{d} \mid s_2-k_2'\)
\(\therefore s_1-k_1'=k*\frac{a_2}{d} \rightarrow k_1=k_1'+k*\frac{a_2}{d}\)整除式变形
\(\therefore s_2-k_2'=k*\frac{a_1}{d} \rightarrow k_2=k_2'+k*\frac{a_1}{d}\)整除式变形
将\(k_1\)通解代入 ③ 得
对照原方程\(x=a_1*k_1+m_1,x=a_2*k_2+m_2\),我们发现,两个方程合并后得到的新方程,其实和原来的两个方程都是差不多的样子,只不过系数发生了变化,而合并后的方程肯定是即符合 ① 双符合 ②, 如果继续迭代下去,迭代\(n-1\)次后,可以将\(n\)个线性同余方程合并为一个方程,最终方程就能满足所有的约束要求。为了能够实现代码循环迭代,那么就可以
这样反复循环迭代就可以找出最后一全兼容方程了!
\(Q\):此处为什么要取绝对值呢
\(A\):因为不知道\(\frac{a_2}{d}\)的正负性,我们在原基础上要尽量减多个\(abs(\frac{a_2}{d})\),使其为正整数且最小。
四、练习题
模数互质
\(P1495\) 【模板】中国剩余定理(\(CRT\))/曹冲养猪
\(Code\) 中国剩余定理解法
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
const int N = 20;
typedef __int128 INT128;
int n; // 同余方程个数
int a[N]; // 除数数组
int m[N]; // 余数数组
int A = 1; // 除数累乘积
int x, y; // (x,y)是方程 Mi * x + m[i] * y = 1的一个解,x是Mi关于m[i]的逆元
int res; // 最小整数解
// 扩展欧几里得模板
int exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1, y = 0;
return a;
}
int d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
// 中国剩余定理模板
void CRT() {
// 一、预处理
for (int i = 1; i <= n; i++) {
cin >> a[i] >> m[i]; // 读入除数数组和余数数组
m[i] = (m[i] % a[i] + a[i]) % a[i]; // ① 预算余数为正数: r[i]可能为负数,预处理成正数,本题没有这个要求,但考虑到通用模板,加上了这块代码
A *= a[i]; // ② 预处理除数连乘积
}
// 二、计算
// 功能: 求 x ≡ r_i mod m_i 方程组的解
for (int i = 1; i <= n; i++) {
int Ai = A / a[i]; // ①除数连乘积除去当前的除数,得到Mi
exgcd(Ai, a[i], x, y); // ②扩欧求逆元
x = (x % a[i] + a[i]) % a[i]; // ③逆元防负数
res = (res + INT128(m[i] * A / a[i] * x) % A) % A; // ④累加res(看公式) 快速乘,防止LL也在乘法过程中爆炸
}
}
signed main() {
cin >> n;
CRT();
cout << res << endl;
}
\(Code\) 扩展中国剩余定理解法
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
const int N = 100010;
int a[N], m[N];
int n;
int exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1, y = 0;
return a;
}
int d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
int exCRT() { // m代表合并多少次,有n个方程,其实m=n-1
int a1 = a[1], m1 = m[1], a2, m2, k1, k2;
for (int i = 2; i <= n; i++) { // a1,m1:做为先遣军,准备与后续的方程合并
a2 = a[i], m2 = m[i];
int d = exgcd(a1, -a2, k1, k2); // 扩展欧几里得求方程特解,x有用,y无用
if ((m2 - m1) % d) return -1; // 拼出来的同余方程无解,继续不下去
k1 *= (m2 - m1) / d; // 同比扩大指定倍数
int t = abs(a2 / d); // 最小累计单元
k1 = (k1 % t + t) % t; // 取得最小正整数解,否则可以会在计算过程中爆long long
m1 = k1 * a1 + m1; // 准备下一轮的m[1],a[1]
a1 = abs(a1 / d * a2); // 两个顺序不能反
}
return (m1 % a1 + a1) % a1; // 取得最小正整数解
}
signed main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> a[i] >> m[i];
cout << exCRT() << endl;
}
\(Code\) 中国剩余定理解法
#include <cstdio>
#include <iostream>
using namespace std;
#define endl "\r"
#define int long long
const int N = 10;
int n = 3;
int a[N] = {0, 23, 28, 33};
int m[N];
int A = 1;
int x, y;
int res;
int exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1, y = 0;
return a;
}
int d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
int qmul(int a, int b, int mod) {
int res = 0;
while (b) {
if (b & 1) res = (res + a) % mod;
a = (a + a) % mod;
b >>= 1;
}
return res;
}
void CRT() {
for (int i = 1; i <= n; i++) {
m[i] = (m[i] % a[i] + a[i]) % a[i];
A *= a[i];
}
for (int i = 1; i <= n; i++) {
int Ai = A / a[i];
exgcd(Ai, a[i], x, y);
x = (x % a[i] + a[i]) % a[i];
res = (res + qmul(m[i], A / a[i] * x, A)) % A;
}
}
int T;
int p, e, i, d;
signed main() {
while (cin >> p >> e >> i >> d && ~p) {
T++;
A = 1;
res = 0;
m[1] = p, m[2] = e, m[3] = i;
CRT();
if (res <= d) res += A;
printf("Case %lld: the next triple peak occurs in %lld days.\n", T, res - d);
}
}
\(Code\) 扩展中国剩余定理
#include <cstdio>
#include <iostream>
#include <stdlib.h>
using namespace std;
#define int long long
#define endl "\r"
const int N = 10;
int n = 3;
int a[N] = {0, 23, 28, 33};
int m[N];
int exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1, y = 0;
return a;
}
int d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
int exCRT() { // m代表合并多少次,有n个方程,其实m=n-1
int a1 = a[1], m1 = m[1], a2, m2, k1, k2;
for (int i = 2; i <= n; i++) { // a1,m1:做为先遣军,准备与后续的方程合并
a2 = a[i], m2 = m[i];
int d = exgcd(a1, -a2, k1, k2); // 扩展欧几里得求方程特解,x有用,y无用
if ((m2 - m1) % d) return -1; // 拼出来的同余方程无解,继续不下去
k1 *= (m2 - m1) / d; // 同比扩大指定倍数
int t = abs(a2 / d); // 最小累计单元
k1 = (k1 % t + t) % t; // 取得最小正整数解,否则可以会在计算过程中爆long long
m1 = k1 * a1 + m1; // 准备下一轮的m[1],a[1]
a1 = abs(a1 / d * a2); // 两个顺序不能反
}
return (m1 % a1 + a1) % a1; // 取得最小正整数解
}
int T;
int m1, m2, m3, d;
signed main() {
while (cin >> m[1] >> m[2] >> m[3] >> d && ~m[1]) {
T++;
int res = exCRT();
int A = a[1] * a[2] * a[3];
if (res <= d) res += A;
printf("Case %lld: the next triple peak occurs in %lld days.\n", T, res - d);
}
}
\(Code\) 中国剩余定理实现代码
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
const int N = 20;
int qmul(int a, int b, int mod) {
int res = 0;
while (b) {
if (b & 1) res = (res + a) % mod;
a = (a + a) % mod;
b >>= 1;
}
return res;
}
int n;
int a[N];
int m[N];
int A = 1;
int x, y;
int res;
int exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1, y = 0;
return a;
}
int d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
void CRT() {
for (int i = 1; i <= n; i++) cin >> m[i];
for (int i = 1; i <= n; i++) cin >> a[i];
for (int i = 1; i <= n; i++) {
m[i] = (m[i] % a[i] + a[i]) % a[i];
A *= a[i];
}
for (int i = 1; i <= n; i++) {
int Mi = A / a[i];
exgcd(Mi, a[i], x, y);
x = (x % a[i] + a[i]) % a[i];
res = (res + qmul(m[i], A / a[i] * x, A)) % A;
}
}
signed main() {
cin >> n;
CRT();
cout << res << endl;
}
\(Code\) 扩展中国剩余定理实现代码
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
const int N = 100010;
int a[N], m[N];
int n;
int exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1, y = 0;
return a;
}
int d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
int exCRT() { // m代表合并多少次,有n个方程,其实m=n-1
int a1 = a[1], m1 = m[1], a2, m2, k1, k2;
for (int i = 2; i <= n; i++) { // a1,m1:做为先遣军,准备与后续的方程合并
a2 = a[i], m2 = m[i];
int d = exgcd(a1, -a2, k1, k2); // 扩展欧几里得求方程特解,x有用,y无用
if ((m2 - m1) % d) return -1; // 拼出来的同余方程无解,继续不下去
k1 *= (m2 - m1) / d; // 同比扩大指定倍数
int t = abs(a2 / d); // 最小累计单元
k1 = (k1 % t + t) % t; // 取得最小正整数解,否则可以会在计算过程中爆long long
m1 = k1 * a1 + m1; // 准备下一轮的m[1],a[1]
a1 = abs(a1 / d * a2); // 两个顺序不能反
}
return (m1 % a1 + a1) % a1; // 取得最小正整数解
}
signed main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> m[i];
for (int i = 1; i <= n; i++) cin >> a[i];
cout << exCRT() << endl;
}
模数不互质
洛谷\(P4777\) 【模板】扩展中国剩余定理(\(EXCRT\)).
上面两道题是同一道题,代码是一样的。
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
const int N = 100010;
int a[N], m[N];
int n;
int exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1, y = 0;
return a;
}
int d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
int exCRT() { // m代表合并多少次,有n个方程,其实m=n-1
int a1 = a[1], m1 = m[1], a2, m2, k1, k2;
for (int i = 2; i <= n; i++) { // a1,m1:做为先遣军,准备与后续的方程合并
a2 = a[i], m2 = m[i];
int d = exgcd(a1, -a2, k1, k2); // 扩展欧几里得求方程特解,x有用,y无用
if ((m2 - m1) % d) return -1; // 拼出来的同余方程无解,继续不下去
k1 *= (m2 - m1) / d; // 同比扩大指定倍数
int t = abs(a2 / d); // 最小累计单元
k1 = (k1 % t + t) % t; // 取得最小正整数解,否则可以会在计算过程中爆long long
m1 = k1 * a1 + m1; // 准备下一轮的m[1],a[1]
a1 = abs(a1 / d * a2); // 两个顺序不能反
}
return (m1 % a1 + a1) % a1; // 取得最小正整数解
}
signed main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> a[i] >> m[i];
cout << exCRT() << endl;
}