2021年9月29日 星期三

「学习笔记」(扩展)中国剩余定理 世界热消息

时间:2023-05-28 20:58:32来源 : 博客园

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


(资料图片)

该问题出自《孙子算经》,具体问题的解答口诀由明朝数学家程大位在《算法统宗》中给出:

三人同行七十希,五树梅花廿一支,七子团圆正半月,除百零五便得知。

\(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;}

关键词:

(责任编辑:黄俊飞)

推荐内容

Back to Top