Loading

「学习笔记」(扩展)中国剩余定理

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

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

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

\(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 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)

//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)/ 曹冲养猪

//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;
}
posted @ 2023-05-28 21:11  yi_fan0305  阅读(37)  评论(0编辑  收藏  举报