有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?
(资料图)
该问题出自《孙子算经》,具体问题的解答口诀由明朝数学家程大位在《算法统宗》中给出:
三人同行七十希,五树梅花廿一支,七子团圆正半月,除百零五便得知。
\(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\) 两两互质):
\[\begin{cases}x \bmod p_1 = a_1\\x \bmod p_2 = a_2\\\cdots\\x \bmod p_n = a_n\\\end{cases}\]过程对于方程
\[\begin{cases}x \bmod p_1 = a_1\\x \bmod p_2 = a_2\\\end{cases}\]我们可以得知 \(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 solve(int p1, int a1, int p2, int a2, int &p, int &a) {// x % p1 = a1// x % p2 = a2// x % p = a// x = k1 * p1 + a1 = k2 * p2 + a2// k1 * p1 - k2 * p2 = a2 - a1 int g, k1, k2; g = exgcd(p1, p2, k1, k2);// k1 * p1 + k2 * p2 = g k2 = -k2;// k1 * p1 - k2 * p2 = g if ((a2 - a1) % g != 0) { p = -1, a = -1; } else { int k = (a2 - a1) / g; k1 = k1 * k, k2 = k2 * k; // k1 * p1 - k2 * p2 = a2 - a1 int x = k1 * p1 + a1; p = p1 / g * p2; a = (x % p + p) % p; }}
通过观察,你会发现,这种解法不需要 \(p_1 \perp p_2\),因此扩展中国剩余定理也是这个解法,这是一种通解。
题目【模板】扩展中国剩余定理(EXCRT)卡 __int128
。
#include using namespace std;typedef long long ll;const int N = 1e5 + 5;int n;ll a[N], b[N];ll qtimes(ll x, ll y, ll p) {if (y == 0) return 0;ll z = qtimes(x, y / 2, p);z = (z + z) % p;if (y & 1) {z = (1ll * z + x) % p;}return z;}ll exgcd(ll a, ll b, ll &x, ll &y) {if (b == 0) {x = 1, y = 0;return a;}ll xp, yp;ll g = exgcd(b, a % b, xp, yp);x = yp, y = xp - yp * (a / b);return g;}void CRT(ll p1, ll a1, ll p2, ll a2, ll &p, ll &a) {ll g, k1, k2;g = exgcd(p1, p2, k1, k2);k2 = -k2;if ((a2 - a1) % g != 0) {p = -1, a = -1;}else {ll k = (a2 - a1) / g;__int128 K1 = k1 * k, K2 = k2 * k;__int128 x = K1 * p1 + a1;p = p1 / g * p2;a = (x % p + p) % p;}}int main() {ios::sync_with_stdio(false);cin.tie(0), cout.tie(0);cin >> n;for (int i = 1; i <= n; ++ i) {cin >> a[i] >> b[i];}ll mod = a[1], rest = b[1];for (int i = 2; i <= n; ++ i) {CRT(mod, rest, a[i], b[i], mod, rest);}cout << rest % mod << "\n";return 0;}
【模板】中国剩余定理(CRT)/ 曹冲养猪
#include using namespace std;typedef long long ll;const int N = 1e5 + 5;int n;ll a[N], b[N];ll qtimes(ll x, ll y, ll p) {if (y == 0) return 0;ll z = qtimes(x, y / 2, p);z = (z + z) % p;if (y & 1) {z = (1ll * z + x) % p;}return z;}ll exgcd(ll a, ll b, ll &x, ll &y) {if (b == 0) {x = 1, y = 0;return a;}ll xp, yp;ll g = exgcd(b, a % b, xp, yp);x = yp, y = xp - yp * (a / b);return g;}void CRT(ll p1, ll a1, ll p2, ll a2, ll &p, ll &a) {ll g, k1, k2;g = exgcd(p1, p2, k1, k2);k2 = -k2;if ((a2 - a1) % g != 0) {p = -1, a = -1;}else {ll k = (a2 - a1) / g;__int128 K1 = k1 * k, K2 = k2 * k;__int128 x = K1 * p1 + a1;p = p1 / g * p2;a = (x % p + p) % p;}}int main() {ios::sync_with_stdio(false);cin.tie(0), cout.tie(0);cin >> n;for (int i = 1; i <= n; ++ i) {cin >> a[i] >> b[i];}ll mod = a[1], rest = b[1];for (int i = 2; i <= n; ++ i) {CRT(mod, rest, a[i], b[i], mod, rest);}cout << rest % mod << "\n";return 0;}
标签:
X 关闭
X 关闭