CRT 和 EXCRT

· · 个人记录

CRT 和 EXCRT

CRT

题目问题：解方程

\begin{cases} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \\ \dots\\ x \equiv a_i \pmod {m_i} \\ \forall i, j, i\neq j , m_i \perp m_j \end{cases}

解法：

M = \sum_{i = 1} ^ n m_i$, $M_i = \frac{M}{ m_i}

可以发现 M_i M_i^{-1} \equiv 1 \pmod {m_i}

诶，我们发现 a \times b\mod c = a \mod c \times b \mod c

所以我们可以将 LHS 和 RHS 同时乘以 a_i 这样就可以构造出一个合法的解

所以 ans = \sum_{i = 1} ^ n {M_iM_i^{-1}a_i}

因为取最小正整数解，所以只要 \mod M 就可以了

之后我们考虑逆元应该怎么求，可以用费马小定理，但是觉得应该学学 exgcd

我们来看看 exgcd，可以解 ax + by = \gcd(a, b) 这个方程

$a \to b, b \to a - b\lfloor \frac{a}{b}\rfloor

所以我们应该变换 x = y', y = x - y'\lfloor \frac{a}{b}\rfloor 就可以了

第一次自己推出来， Happy!

\gcd(M_i, m_i^{-1}) = 1 a = M_i, b = m_i^{-1}, ax = 1 - by \Rightarrow ax \equiv 1 \pmod b

也就是 M_iM_i^{-1} \equiv 1 \pmod {m_i}

P1495 中国剩余定理(CRT)

Code 
#include <bits/stdc++.h>
using namespace std;
template <typename T>
void r1(T &x) {
    x = 0;
    char c(getchar());
    int f(1);
    for(; !isdigit(c); c = getchar()) if(c == '-') f = -1;
    for(; isdigit(c);c = getchar()) x = (x << 1) + (x << 3) + (c ^ 48);
    x *= f;
}
#define int long long
const int maxn = 2e5 + 5;
const int maxm = maxn << 1;

typedef int room[maxn];

int n, m, M(1);
int mi[maxn], a[maxn], Mi[maxn];
int ans(0);

void exgcd(int a,int b,int &x,int &y) {
    if(b == 0) return x = 1, y = 0, void();
    exgcd(b, a % b, x, y);
    int z = x; x = y, y = z - y * (a / b);
    // 注意 y 的处理顺序，需要 floor(a / b) 所以不能让 y 的因子被除
}

signed main() {
    int i, j;
    r1(n);
    for(i = 1;i <= n; ++ i) {
        r1(mi[i]), r1(a[i]);
        M *= mi[i];
    }
    for(i = 1;i <= n; ++ i) {
        Mi[i] = M / mi[i];
        int x(0), y(0);
        exgcd(Mi[i], mi[i], x, y);
        x = (x + mi[i]) % mi[i];
        ans += Mi[i] * x * a[i] % M;
    }
    printf("%lld\n", ans % M);
    return 0;
}

EXCRT

题目意思：解方程

\begin{cases} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \\ \dots \\ x \equiv a_i \pmod {m_i} \\ \forall i, j, m_i \not\perp m_j \end{cases}

我们考虑对于一个方程组，假设已经求出了前 k 个方程的解。

M = \sum_{i = 1} ^ k m_i$ 代码实现的时候让 $M =Lcm(m_1,m_2\dots m_k)

这样不容易炸

那么我们存在一个通解 x + i \times M (i \in Z)

之后我们考虑加入一个元素之后的方程组

x + t \times M \equiv a_k \pmod {m_k},t > 0 t \times M \equiv a_k - x \pmod {m_k}, t > 0

之后再让 x' = x + t \times M

既然我们已知 M, a_k 就可以使用 exgcd 求解

众所周知我们需要求的方程是

M x\equiv a_k - ans_{last} \pmod {m_i}

也就是 ax + cy = b

而 exgcd 求解的方程是 ax + by = gcd(a, b)

诶，我们就可以直接用 exgcd 求出

之后考虑 $x\times \frac{a_k - ans_{last}}{gcd(M, m_i)}$ 就可以得到答案 当然 $nowans \text{ 要加上 } M \times x

我们的 M 也要乘上 \frac{m_i}{gcd} 原因在开头已经说过了

\color{red}\text{重点:}$ 如何判定无解，就是不能 $\gcd(M, m_i) \not| \ a_k - ans_{last}

也可以理解成 (a_k - ans_{last}) 是 \gcd(M, m_i) 的倍数

至于为什么要自己写，因为没有题解看得懂

Code 
#include <bits/stdc++.h>
using namespace std;
template <typename T>
void r1(T &x) {
    x = 0;
    char c(getchar());
    int f(1);
    for(; !isdigit(c); c = getchar()) if(c == '-') f = -1;
    for(; isdigit(c);c = getchar()) x = (x << 1) + (x << 3) + (c ^ 48);
    x *= f;
}
#define int long long
const int maxn = 1e5 + 5;
const int maxm = maxn << 1;

typedef int room[maxn];
int exgcd(int a,int b,int &x,int &y) {
    if(b == 0) return x = 1, y = 0, a;
    int gcd = exgcd(b, a % b, x, y);
    int z = x; x = y, y = z - y * (a / b);
    return gcd;
}
int ksc(int x,int mi,int mod) {
    int ans(0);
    while(mi) {
        if(mi & 1) ans = (ans + x) % mod;
        x = (x + x) % mod;
        mi >>= 1;
    }
    return ans;
}
int n, m[maxn], a[maxn];
int ans, M;

int excrt() {
    M = m[1], ans = a[1];
    for(int i = 2;i <= n; ++ i) {
        int x, y;
        int gcd = exgcd(M, m[i], x, y);
        int Up = ((a[i] - ans) % m[i] + m[i]) % m[i];
        if(Up % gcd) return -1;

        x = ksc(x, Up / gcd, m[i]);
//        x = (x + m[i]) % m[i];
        ans += M * x;
        M = M * (m[i] / gcd);
        ans = (ans + M) % M;
    }
    return ans;
}

signed main() {
    int i, j;
    r1(n);
    for(i = 1; i <= n; ++ i) r1(m[i]), r1(a[i]);
    printf("%lld\n", excrt());
    return 0;
}