中国剩余定理

中国剩余定理

求解模线性方程组。

{x1≡a1(modr1)x2≡a2(modr2)⋯xk≡ak(modrk)" role="presentation" style="text-align: center; position: relative;">⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪x1≡a1(modr1)x2≡a2(modr2)⋯xk≡ak(modrk)$$\begin{array}{r}\{\begin{array}{l}{x}_1\equiv a_1(\mathrm{mod}{r}_1)\\ x_2\equiv a_2(\mathrm{mod}{r}_2)\\ \cdots \\ x_k\equiv a_k(\mathrm{mod}{r}_k)\end{array}\end{array}$$

crt

当 r1,r2,⋯,rk" role="presentation" style="position: relative;">r1,r2,⋯,rk$r_1,r_2,\cdots ,r_k$ 为互质的，设 M=∏i=1kri,Ai=Mri" role="presentation" style="position: relative;">M=∏ki=1ri,Ai=Mri$M=\prod _{i=1}^k{r}_i,A_i={\displaystyle \frac{M}{r_i}}$ ，可的 M" role="presentation" style="position: relative;">M$M$ 是 r1,⋯,rk" role="presentation" style="position: relative;">r1,⋯,rk$r_1,\cdots ,r_k$ 的 lcm

可得 Ai,ri" role="presentation" style="position: relative;">Ai,ri$A_i,r_i$ 是互质的，则必存在 Ai" role="presentation" style="position: relative;">Ai$A_i$ 模 ri" role="presentation" style="position: relative;">ri$r_i$ 意义下的逆元 ti" role="presentation" style="position: relative;">ti$t_i$ 使得 Aiti≡1(modri)" role="presentation" style="position: relative;">Aiti≡1(modri)$A_i{t}_i\equiv 1(\mathrm{mod}{r}_i)$

解 xi=aiAiti" role="presentation" style="position: relative;">xi=aiAiti$x_i=a_i{A}_i{t}_i$ 是满足第 i" role="presentation" style="position: relative;">i$i$ 个方程的。同时 ∀j≠i" role="presentation" style="position: relative;">∀j≠i$\mathrm{\forall }j\ne i$ ，都有 xi≡0(modrj)" role="presentation" style="position: relative;">xi≡0(modrj)$x_i\equiv 0(\mathrm{mod}{r}_j)$ ，因为 Ai" role="presentation" style="position: relative;">Ai$A_i$ 中有一个因子 rj" role="presentation" style="position: relative;">rj$r_j$

可以得到通解 x=∑i=1kaiAiti+pM(p∈Z)" role="presentation" style="position: relative;">x=∑ki=1aiAiti+pM(p∈Z)$x=\sum _{i=1}^k{a}_i{A}_i{t}_i+pM(p\in \mathbb{Z})$ 。显然满足条件

#include 
using namespace std;
typedef unsigned long long uLL;
typedef long double LD;
typedef long long LL;
typedef double db;
const int N = 55;
LL exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) return x = 1, y = 0, a;
    register LL gcd = exgcd(b, a % b, y, x);
    return y -= a / b * x, gcd;
}
inline LL inv(LL a, LL b) {
    register LL x, y;
    exgcd(a, b, x, y);
    return (x % b + b) % b;
}
int n;
LL a[N], r[N], mul = 1, A[N], ans;
int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) scanf("%lld%lld", &r[i], &a[i]), mul *= r[i];
    for (int i = 1; i <= n; i++) {
        A[i] = mul / r[i];
        ans += a[i] * A[i] * inv(A[i], r[i]);
    }
    printf("%lld", ans % mul);
}

excrt

当 r1,⋯,rk" role="presentation" style="position: relative;">r1,⋯,rk$r_1,\cdots ,r_k$ 不互质，此时使用扩展 crt

{x1≡a1(modr1)x2≡a2(modr2)⋯xk≡ak(modrk)" role="presentation" style="text-align: center; position: relative;">⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪x1≡a1(modr1)x2≡a2(modr2)⋯xk≡ak(modrk)$$\begin{array}{r}\{\begin{array}{l}{x}_1\equiv a_1(\mathrm{mod}{r}_1)\\ x_2\equiv a_2(\mathrm{mod}{r}_2)\\ \cdots \\ x_k\equiv a_k(\mathrm{mod}{r}_k)\end{array}\end{array}$$

看前两个方程。

x=k⋅r1+a1=p⋅r2+a2" role="presentation" style="position: relative;">x=k⋅r1+a1=p⋅r2+a2$x=k\cdot r_1+a_1=p\cdot r_2+a_2$ ，其中 k,p∈Z" role="presentation" style="position: relative;">k,p∈Z$k,p\in \mathbb{Z}$

k⋅r1−p⋅r2=a2−a1" role="presentation" style="position: relative;">k⋅r1−p⋅r2=a2−a1$k\cdot r_1-p\cdot r_2=a_2-a_1$ ，这是形如 ax+by=c" role="presentation" style="position: relative;">ax+by=c$ax+by=c$ 的形式，用 exgcd 求解。

此处可得 gcd(r1,r2)∤(a2−a1)" role="presentation" style="position: relative;">gcd(r1,r2)∤(a2−a1)$gcd(r_1,r_2)\nmid (a_2-a_1)$ 时无解。

求出一个解 k0" role="presentation" style="position: relative;">k0$k_0$ ，则通解 k=k0+o∗r2gcd(r1,r2)(o∈Z)" role="presentation" style="position: relative;">k=k0+o∗r2gcd(r1,r2)(o∈Z)$k=k_0+o\ast {\displaystyle \frac{r_2}{gcd(r_1,r_2)}}(o\in \mathbb{Z})$

将 k" role="presentation" style="position: relative;">k$k$ 带入 k⋅r1+a1" role="presentation" style="position: relative;">k⋅r1+a1$k\cdot r_1+a_1$ 可得

(k0+o∗r2gcd(r1,r2))⋅r1+a1=k0∗r1+a1+o∗lcm(r1,r2)" role="presentation" style="text-align: center; position: relative;">(k0+o∗r2gcd(r1,r2))⋅r1+a1=k0∗r1+a1+o∗lcm(r1,r2)$$(k_0+o\ast {\displaystyle \frac{r_2}{gcd(r_1,r_2)}})\cdot r_1+a_1=k_0\ast r_1+a_1+o\ast \text{lcm}(r_1,r_2)$$

即 x≡k0∗r1+a1(modlcm(r1,r2))" role="presentation" style="position: relative;">x≡k0∗r1+a1(modlcm(r1,r2))$x\equiv k_0\ast r_1+a_1(\mathrm{mod}\text{lcm}(r_1,r_2))$

相当于合并了两个方程，合并 n−1" role="presentation" style="position: relative;">n−1$n-1$ 次即可得到答案。

注意 k0" role="presentation" style="position: relative;">k0$k_0$ 应乘上 a2−a1gcd(r1,r2)" role="presentation" style="position: relative;">a2−a1gcd(r1,r2)${\displaystyle \frac{a_2-a_1}{gcd(r_1,r_2)}}$

模板题需要用龟速乘。

#include 
using namespace std;
typedef unsigned long long uLL;
typedef long double LD;
typedef long long LL;
typedef double db;
const int N = 100005;
LL exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) return x = 1, y = 0, a;
    register LL gcd = exgcd(b, a % b, y, x);
    return y -= a / b * x, gcd;
}
LL G;
inline LL mul(LL x, LL y, LL p) {
    return (x * y - (LL)((LD)x / p * y) * p + p) % p;
}
inline bool sol(LL a, LL b, LL c, LL &x, LL &y) {
    G = exgcd(a, b, x, y);
    register LL t;
    if (c % G) return false;
    t = b / G;
    x = mul(x, c / G, t);
    x = (x + t) % t;
    y = (c - a * x) / b;
    return true;
}
int n, isL;
LL a[N], r[N];
int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) scanf("%lld%lld", &r[i], &a[i]);
    isL = 1;
    for (int i = 1; i < n; i++) {
        register LL x, y;
        if (!sol(r[i], r[i + 1], a[i + 1] - a[i], x, y)) {
            isL = 0; break;
        }
        a[i + 1] = x * r[i] + a[i];
        r[i + 1] = r[i] / G * r[i + 1];
    }
    printf("%lld", a[n]);
}