有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?
该问题出自《孙子算经》,具体问题的解答口诀由明朝数学家程大位在《算法统宗》中给出:
三人同行七十希,五树梅花廿一支,七子团圆正半月,除百零五便得知。
\(2 \times 70 + 3 \times 21 + 2 \times 15 = 233 = 2 \times 105 + 23\),故答案为 \(23\)。
定义
中国剩余定理 (Chinese Remainder Theorem, CRT) 可求解如下形式的一元线性同余方程组(其中 \(n_1, n_2, \cdots, n_k\) 两两互质):
过程
对于方程
我们可以得知 \(x = k_1 p_1 + a_1 = k_2 p_2 + a_2\),进而推出 \(k_1 p_1 - k_2 p_2 = a_2 - a_1\)
我们观察这个式子,是不是跟 \(xa + yg = g\) 很像?
我们可以用扩展欧几里得算法来做。
扩展欧几里得算法代码
int exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1, y = 0;
return a;
}
int xp, yp;
int g = exgcd(b, a % b, xp, yp);
x = yp, y = xp - yp * (a / b);
return g;
}
中国剩余定理代码
void calc(ll a1, ll b1, ll a2, ll b2, ll& tmpa, ll& tmpb) {
ll ansx, ansy, g;
g = exgcd(a1, a2, ansx, ansy);
ansy = -ansy;
if ((b2 - b1) % g) {
tmpa = tmpb = -1;
return ;
}
tmpa = a1 / g * a2;
ll k = (b2 - b1) / g;
ansx = qtimes(ansx, k, tmpa);
ansy = qtimes(ansy, k, tmpa);
ll ans = (qtimes(a1, ansx, tmpa) + b1) % tmpa;
tmpb = (ans % tmpa + tmpa) % tmpa;
}
通过观察,你会发现,这种解法不需要 \(p_1 \perp p_2\),因此扩展中国剩余定理也是这个解法,这是一种通解。
题目
【模板】扩展中国剩余定理(EXCRT)文章来源:https://www.toymoban.com/news/detail-462515.html
//The code was written by yifan, and yifan is neutral!!!
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define bug puts("NOIP rp ++!");
#define rep(i, a, b, c) for (int i = (a); i <= (b); i += (c))
#define per(i, a, b, c) for (int i = (a); i >= (b); i -= (c))
template<typename T>
inline T read() {
T x = 0;
bool fg = 0;
char ch = getchar();
while (ch < '0' || ch > '9') {
fg |= (ch == '-');
ch = getchar();
}
while (ch >= '0' && ch <= '9') {
x = (x << 3) + (x << 1) + (ch ^ 48);
ch = getchar();
}
return fg ? ~x + 1 : x;
}
const int N = 1e5 + 5;
int n;
ll a[N], b[N];
ll exgcd(ll a, ll b, ll& x, ll& y) {
if (!b) {
x = 1, y = 0;
return a;
}
ll g = exgcd(b, a % b, y, x);
y -= (a / b) * x;
return g;
}
ll qtimes(ll x, ll y, ll mod) {
ll ans = 0;
if (y < 0) {
x = -x;
y = -y;
}
while (y) {
if (y & 1) {
ans = (ans + x) % mod;
}
y >>= 1;
x = (x + x) % mod;
}
return ans;
}
void calc(ll a1, ll b1, ll a2, ll b2, ll& tmpa, ll& tmpb) {
ll ansx, ansy, g;
g = exgcd(a1, a2, ansx, ansy);
ansy = -ansy;
if ((b2 - b1) % g) {
tmpa = tmpb = -1;
return ;
}
tmpa = a1 / g * a2;
ll k = (b2 - b1) / g;
ansx = qtimes(ansx, k, tmpa);
ansy = qtimes(ansy, k, tmpa);
ll ans = (qtimes(a1, ansx, tmpa) + b1) % tmpa;
tmpb = (ans % tmpa + tmpa) % tmpa;
}
int main() {
n = read<int>();
rep (i, 1, n, 1) {
a[i] = read<ll>(), b[i] = read<ll>();
}
ll lasa = a[1], lasb = b[1];
rep (i, 2, n, 1) {
calc(lasa, lasb, a[i], b[i], lasa, lasb);
}
cout << lasb << '\n';
return 0;
}
【模板】中国剩余定理(CRT)/ 曹冲养猪文章来源地址https://www.toymoban.com/news/detail-462515.html
//The code was written by yifan, and yifan is neutral!!!
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define bug puts("NOIP rp ++!");
#define rep(i, a, b, c) for (int i = (a); i <= (b); i += (c))
#define per(i, a, b, c) for (int i = (a); i >= (b); i -= (c))
template<typename T>
inline T read() {
T x = 0;
bool fg = 0;
char ch = getchar();
while (ch < '0' || ch > '9') {
fg |= (ch == '-');
ch = getchar();
}
while (ch >= '0' && ch <= '9') {
x = (x << 3) + (x << 1) + (ch ^ 48);
ch = getchar();
}
return fg ? ~x + 1 : x;
}
const int N = 1e5 + 5;
int n;
ll a[N], b[N];
ll exgcd(ll a, ll b, ll& x, ll& y) {
if (!b) {
x = 1, y = 0;
return a;
}
ll g = exgcd(b, a % b, y, x);
y -= (a / b) * x;
return g;
}
ll qtimes(ll x, ll y, ll mod) {
ll ans = 0;
if (y < 0) {
x = -x;
y = -y;
}
while (y) {
if (y & 1) {
ans = (ans + x) % mod;
}
y >>= 1;
x = (x + x) % mod;
}
return ans;
}
void calc(ll a1, ll b1, ll a2, ll b2, ll& tmpa, ll& tmpb) {
ll ansx, ansy, g;
g = exgcd(a1, a2, ansx, ansy);
ansy = -ansy;
if ((b2 - b1) % g) {
tmpa = tmpb = -1;
return ;
}
tmpa = a1 / g * a2;
ll k = (b2 - b1) / g;
ansx = qtimes(ansx, k, tmpa);
ansy = qtimes(ansy, k, tmpa);
ll ans = (qtimes(a1, ansx, tmpa) + b1) % tmpa;
tmpb = (ans % tmpa + tmpa) % tmpa;
}
int main() {
n = read<int>();
rep (i, 1, n, 1) {
a[i] = read<ll>(), b[i] = read<ll>();
}
ll lasa = a[1], lasb = b[1];
rep (i, 2, n, 1) {
calc(lasa, lasb, a[i], b[i], lasa, lasb);
}
cout << lasb << '\n';
return 0;
}
到了这里,关于「学习笔记」(扩展)中国剩余定理的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!