基础数论算法详解：最大公约数、质数筛与同余方程

秦九韶算法

秦九韶算法是一种将一元 n 次多项式的求值问题转化为 n 个一次式的算法。其大大简化了计算过程，即使在现代，利用计算机解决多项式的求值问题时，秦九韶算法依然是最优的算法。

利用唯一分解定理

2. 质数的判定

质数和合数

• 一个大于 1 的自然数，除了 1 和它自身外，不能被其他自然数整除的数叫做质数；
• 否则称为合数。
• 其中，质数又称素数。
• 规定 1 既不是质数也不是合数。

试除法判断质数

• 对于一个数 x，根据定义，可以从 [2,x−1] 一个一个尝试，判断 x 能否被整除。
• 优化思路：若 a 是 x 的约数，那么 x/a 也是 x 的约数。因此仅需判断较小的 a，枚举到 sqrt(x) 即可。

代码实现

bool isprime(int x) { 
    if(x <= 1) return false; // 小于等于 1 的数不考虑
    // 试除法判断是否是质数 - 只需枚举到 sqrt(x)
    for(int i = 2; i <= x / i; i++) // 防溢出的写法
    {
        if(x % i == 0) return false;
    }
    return true;
}

时间复杂度

枚举到 sqrt(n)，因此时间复杂度为 O(sqrt(N))。

模板题

洛谷 P5736 质数筛

3. 筛质数

【引入】

上面学习了如何判断一个数是否是质数，如果此时想知道 [1, n] 中有多少个素数呢？或者是 [1, n] 中的素数里面，第 k 个素数是多少？

• 一个自然的想法就是从 2 开始，依次向后对每一个自然数进行一次质数检验。
• 但是这种解法相对暴力，我们这里介绍两种方法，能够快速地将 [1, n] 中的素数全部记录下来。

第一种筛法：埃氏筛法

算法思想

• 对于任意一个大于 1 的正整数，那么它的 k(k > 1) 倍就是合数。
• 因此，如果我们从小到大考虑每个数，然后同时把当前这个数的所有倍数记为合数，没有被标记的数就是素数。

小优化

• 找到一个质数 x 之后，可以从该数的 x 倍向后筛，因此小于 x 的倍数一定被之前筛过了。

代码实现

const int N = 1e6 + 10;
bool st[N]; // 当前这个数有没有被筛掉
int p[N]; // 记录质数
int cnt; // 统计质数个数
// 埃氏筛
void get_prime(int n) {
    for(int i = 2; i <= n; i++) {
        if(!st[i]) // 没有被标记，说明是质数
        {
            p[++cnt] = i; // 记录这个质数
            // 从 i*i 开始，因为小于 i 的倍数已经被划掉了
            for(int j = i * i; j <= n; j += i) // 筛掉这个质数的倍数
            {
                st[j] = true;
            }
        }
    }
}

时间复杂度

埃氏筛的时间复杂度为：O(n log log n)。

第二种筛法：线性筛法

线性筛法，又称欧拉筛法。

算法思想

• 在埃氏筛法中，它会将一个合数重复多次标记。如果能让每个合数都只被标记一次，那么时间复杂度就可以降到 O(n) 了。
• 我们的做法是，让每一个合数被它的最小质因数筛掉。
• 从前往后遍历每一个树 i；
• 1.用 i 的质数倍去筛（因为要被最小质因数筛掉）
• 2.如果 i 是质数的倍数时，停止

代码实现：

const int N = 1e6 + 10;
int n, q;
bool st[N];
int p[N], cnt;
// 欧拉筛
void get_prime(int n) {
    for(int i = 2; i <= n; i++) {
        if(!st[i]) p[++cnt] = i; // 如果没标记过，就是质数
        // 枚举所有的质数
        for(int j = 1; 1LL * i * p[j] <= n; j++) {
            st[i * p[j]] = true;
            if(i % p[j] == 0) break;
            /* 这个判定条件能让每一个合数被自己的最小质因数筛掉。
             1. 如果 i 是合数，枚举到最小质因数的时候跳出循环
             2. 如果 i 是质数，枚举到自身时跳出循环
             注意，在筛的过程中，我们还能知道 p[j] 是 i 的最小质因数 */
        }
    }
}

时间复杂度

每个数只会被自身最小的质因数筛掉一次，时间复杂度为 O(n)。

4. 算术基本定理

算术基本定理

分解质因数

• 分解质因数就是将一个合数用质因数相乘的形式表示出来，例如：360=2^3×3^2×5^1

试除法分解质因数

• n 的所有因数中，不会有两个大于 sqrt(n)。

因此，枚举 [2, sqrt(N)] 中所有的数，如果能整除 n 就一直除下去。如果最后剩下的数大于 1，那就是大于 sqrt(n) 的因子。

小优化

• 如果我们提前打出来一个素数表，枚举所有素数的话，时间复杂度能进一步降低。

代码实现

const int N = 1e6 + 10;
int c[N]; // c[i] 表示 i 这个质数出现的次数
void deprime(int x) {
    for(int i = 2; i <= x / i; i++) {
        int cnt = 0;
        while(x % i == 0) // 只要有这个因子，就除尽，并且计数
        {
            x /= i;
            cnt++;
        }
        c[i] += cnt;
    }
    if(x > 1) c[x]++; // 不要忘记判断最后一个质数
}

时间复杂度

枚举到 sqrt(N)，因此时间复杂度为 O(sqrt(N))。但是最优情况下会达到 O(logn)。

模板题：

洛谷 P2043 质因子分解

洛谷 P10495 阶乘分解

P10495 题解：

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 1e6 + 10;
const int M = 2e4 + 10;
const int INF = 0x3f3f3f3f;
using PII = pair<int, int>;
using ull = unsigned long long;
#define endl '\n'
int dx[] = {0, 0, 1, -1};
int dy[] = {1, -1, 0, 0};
int dxx[] = {0, 0, 1, -1, 1, -1, 1, -1};
int dyy[] = {1, -1, 0, 0, 1, 1, -1, -1};
ull base = 131;
int n;
int c[N];
int p[N];
int cnt;
bool st[N];

void get_prime(int n) {
    for(int i = 2; i <= n; i++) {
        if(!st[i]) p[++cnt] = i;
        for(int j = 1; i * p[j] <= n; j++) {
            st[i * p[j]] = true;
            if(i % p[j] == 0) break;
        }
    }
}

void solve() {
    cin >> n;
    get_prime(n); // 用质数来筛
    for(int i = 1; i <= cnt; i++) {
        int t = 0;
        for(int j = p[i]; j <= n; j *= p[i]) {
            c[p[i]] += n / j;
        }
        cout << p[i] << " " << c[p[i]] << endl;
    }
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int _ = 1; // cin >> _;
    while(_--) {
        solve();
    }
    return 0;
}

5. 约数

求一个整数的所有约数

对于一个整数 x，若 d 是 x 的约数，那么 x/d 也是 x 的约数。也就是说，除了完全平方数，约数都是成对出现的。因此，可以用试除法求一个整数的所有约数。

枚举 1∼sqrt(x) 之间的整数，判断是否能整除 x。

试除法也能求出一个整数的约数个数以及约数总和。

代码实现

const int N = 1e6 + 10;
int d[N], cnt;
void get_d(int x) {
    // 注意从 1 开始循环
    for(int i = 1; i <= x / i; i++) {
        if(x % i == 0) {
            d[++cnt] = i;
            if(i != x / i) d[++cnt] = x / i;
        }
    }
}

时间复杂度

枚举到 sqrt(N)，因此时间复杂度为 O(sqrt(N))。

因此，一个整数 n 的约数个数的上限为 2*sqrt(N)。

求 [1,n] 每个数的约数集合

如果用试除法分别求每一个数的约数，时间复杂度过高。

可以反过来想，对于每个数 d，[1,n] 中以 d 为约数的数就是 d 的倍数，也就是 d,2d,3d... n/d *d。因此可以用倍数法求出 [1,n] 每个数的约数集合。

代码实现

const int N = 1e6 + 10;
int n;
vector<int> d[N];
void get_d() {
    for(int i = 1; i <= n; i++) // 枚举所有约数
        for(int j = 1; j <= n / i; j++) // 约数的倍数
            d[i * j].push_back(i);
}

时间复杂度

• n 倍的调和数，约为 O(nlogn)。

约数个数的和

#include <iostream>
using namespace std;
#define int long long
signed main() {
    int n;
    cin >> n;
    int sum = 0;
    for(int i = 1; i <= n/2; i++) {
        sum += (n/i);
    }
    sum += n - n/2;
    cout << sum << endl;
    return 0;
}

6. 欧拉函数

互质

• a 与 b 互质：说明 a 与 b 的最大公约数是 1。

欧拉函数

• 对于一个数 n，在 1∼n 中，与 n 互质的数的个数，就是欧拉函数，用 φ(n) 表示。
• 比如：φ(1)=1，φ(2)=1，φ(3)=2，φ(4)=2，φ(5)=4……

积性函数

• 若函数 f(n) 满足：f(1)=1，且 f(xy)=f(x)f(y)，x,y∈N∗ 都成立。则 f(n) 为积性函数。

欧拉函数的性质

1. 若 p 为质数，则 φ(p)=p−1；
2. 若 p 为质数，则 φ(p^k)=(p−1)p^k−1；
3. 欧拉函数为积性函数：φ(xy)=φ(x)φ(y)，其中 gcd(x,y)=1。

试除法求单个数的欧拉函数

int phi(int n) {
    int ret = n;
    for(int i = 2; i <= n / i; i++) {
        if(n % i == 0) {
            ret = ret / i * (i - 1); // 先除后乘，保证不会溢出
            while(n % i == 0) n /= i;
        }
    }
    // 别忘记判断最后一个数
    if(n > 1) ret = ret / n * (n - 1);
    return ret;
}

时间复杂度：枚举到 sqrt(N )，因此时间复杂度为 O(sqrt(N ))。

欧拉筛打表欧拉函数

问题背景：需要知道 [1,n] 中，每一个数的欧拉函数。

代码实现

const int N = 1e6 + 10;
int n;
bool st[N];
int p[N], cnt;
int phi[N];
void get_phi() {
    phi[1] = 1;
    for(int i = 2; i <= n; i++) {
        if(!st[i]) {
            p[cnt++] = i;
            phi[i] = i - 1;
        }
        for(int j = 1; 1LL * i * p[j] <= n; j++) {
            int x = i * p[j]; // 当前要筛的数
            st[x] = true;
            if(i % p[j] == 0) {
                phi[x] = p[j] * phi[i];
                break;
            } else {
                phi[x] = (p[j] - 1) * phi[i];
            }
        }
    }
}

时间复杂度：与线性筛时间复杂度一致，为 O(n)。

7. 费马小定理

同余式

• 如果两个整数 a,b 模 m 的余数相同，则称 a, b 模 m 同余，记作：a ≡ b(mod m)。
• 例如：7 ≡ 4(mod 3)

性质

• 同加性：若 a ≡ b(mod m)，则 a + c ≡ b + c (mod m)
• 同乘性：若 a ≡ b(mod m)，则 ac ≡ bc(mod m)
• 也就是说，对于同余的两个数，执行完相同的加减乘操作后，它们两者还是同余的。
• 但是同余不满足同除性，例如 20≡2(mod3)，左右都除以 10，结果就不是同余的。

这就解释了为什么之前遇到需要取模的题目，如果只涉及加减乘运算，我们可以边计算边取模。而存在除法的时候，不能随意取模。

费马小定理

• 如果 p 为质数，且 a,p 互质，则 a^(p−1) ≡ 1(mod p)。
• 例如：2^(5−1)≡1(mod5)

乘法逆元

• 例如：8x≡1(mod5)，解得 x=2,7,…

若 a,m 互质，且满足同余方程 ax ≡ 1(mod m)，则称 x 是 a 模 m 的乘法逆元，记作 a^-1

应用

• 如果涉及除法，是不满足分配律的。利用乘法逆元就可以解决这一问题。

对于算式 (b÷a)modp，可以写成 (b×a^-1) mod p，其中 a×a^-1 ≡1(modp)。因此，对于这个除数 a，如果能够找到一个数 x，使的两者的乘积在模 m 时为 1，即 ax≡1(modp)，那么就可以用 x 替代 a^-1 去取模。

例如：(25÷5)mod3，对于 5 而言，(5×2)≡1(mod3)，因此用 2 代替 5^-1，原式变为 (25×2)mod3=50mod3=2。

求乘法逆元（一）--- 费马小定理 + 快速幂

我们会学习三种求乘法逆元的方式，一直要注意每种方式的使用条件。

• 例如：(25÷5)mod3，其中 5 在模 3 下的一个乘法逆元为 5^(3-2)=5。因此用 5 代替 5^-1，原式变为 (25×5)mod3=125mod3=2。

那么 a 在模 p 意义下的其中一个乘法逆元就是 a^(p-2)，用快速幂就可以快速求出。

由费马小定理得，当 a,p 互质，且 p 为质数时，a^(p-1) ≡ 1(mod p)，也就是 a×a^(p-2) ≡1(modp)。

一定要注意，利用费马小定理求逆元的时候，模数 p 一定要是质数，且 a,p 互质。不过在我们一般做题的时候，取模的值一般都是质数，比如 1e9+7,1e9+9,998244353，并且 a 一般都比这些模数小，所以可以用费马小定理求逆元。如果 p 不是质数，我们会用后面学到的扩展欧几里得算法求逆元。

代码实现

#include <iostream>
using namespace std;
typedef long long LL;
// 必须要保证 a,p 互质，且 p 为质数。
LL qpow(LL a, LL b, LL p) {
    LL ret = 1;
    while(b) {
        if(b & 1) ret = ret * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return ret;
}
int main() {
    LL x, p;
    cin >> x >> p;
    cout << qpow(x, p - 2, p) << endl;
    return 0;
}

时间复杂度

• 与快速幂时间复杂度一致，为 O(logn)。

序列求和

序列求和

8. 欧拉定理与扩展欧拉定理

模板题：

洛谷 P5091 扩展欧拉定理

#include <iostream>
using namespace std;
#define int long long
int a, m;
string s;

int get_phi(int m) {
    int ret = m;
    for(int i = 2; i <= m/i; i++) {
        if(m % i == 0) {
            ret = ret / i * (i - 1);
            while(m % i == 0) m /= i;
        }
    }
    if(m > 1) ret = ret / m * (m - 1);
    return ret;
}

int get_b(string& s, int p) {
    int t = 0;
    bool flag = false;
    for(int i = 0; i < s.size(); i++) {
        t = t * 10 + s[i] - '0';
        if(t >= p) {
            flag = true;
            t %= p;
        }
    }
    if(flag) t += p;
    return t;
}

int qpow(int a, int b, int p) {
    int ret = 1;
    while(b) {
        if(b & 1) ret = (ret * a) % p;
        a = a % p;
        a = a * a % p;
        b >>= 1;
    }
    return ret;
}

signed main() {
    cin >> a >> m >> s;
    int phi = get_phi(m);
    int b = get_b(s, phi);
    cout << qpow(a,b,m);
    return 0;
}

9. 裴蜀定理

裴蜀定理

裴蜀定理又称贝祖定理：

• 一定存在整数 x,y，满足 ax+by=gcd(a,b)。
• 可以用来判定不定方程是否存在整数解。

例如：6x+8y=gcd(6,8)=2，一定有整数解。其中 x=−1,y=1 就是一组解。

裴蜀定理的推论

1. 一定存在整数 x,y，满足 ax+by=gcd(a,b)×n。因此，对于一个二元一次不定方程 ax+by=c，当 gcd(a,b)∣c 时，有解。
2. 一定存在整数 x1,x2,x3,…，满足 a1x1+a2x2+a3x3+⋯+akxk=gcd(a1,a2,a3,…,ak)×n。

注意

• a,b 的正负是不影响结果的。因为 a,b 如果存在解，那么 a,−b 也一定存在解，只不过在 a,b 解的基础上添上一个负号。
• 因此，在用欧几里得算法求 gcd(a,b) 时，为了避免出现负数，可以带入 a,b 的绝对值。

洛谷 P4549 裴蜀定理

#include <iostream>
#include <cmath>
using namespace std;
int gcd(int a, int b) {
    return b == 0 ? a : gcd(b, a % b);
}
int main() {
    int n;
    cin >> n;
    int ret;
    cin >> ret;
    ret = abs(ret);
    for(int i = 2; i <= n; i++) {
        int x;
        cin >> x;
        ret = gcd(ret, abs(x));
    }
    cout <<ret << endl;
    return 0;
}

10. 扩展欧几里得算法

扩展欧几里得算法

裴蜀定理遗留问题：如何求 ax+by=gcd(a,b) 的一组整数解？利用扩展欧几里得算法即可求解。

• 当 b=0 时，ax+by=ax，因此有一组解为：x=1,y=0；

当 b != 0 时：

代码实现

typedef long long LL;
LL exgcd(LL a, LL b, LL& x, LL& y) {
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    LL x1, y1, d; // 先求出下一层的 x1 和 y1
    d = exgcd(b, a % b, x1, y1); // 求当前层
    x = y1, y = x1 - a / b * y1;
    return d;
}
int main() {
    int a, b, x, y;
    int d = exgcd(a, b, x, y);
    cout << x << " " << y << endl;
    return 0;
}

时间复杂度

与欧几里得算法时间复杂度一致，为 O(logn)。

扩欧求二元一次不定方程的解

求解二元一次不定方程 ax+by=c 的解:

模板题：

洛谷 P5656 二元一次不定方程 (exgcd)

代码：

#include <iostream>
using namespace std;
#define int long long
int exgcd(int a, int b, int& x, int& y) {
    if(b == 0) {
        x = 1,y = 0;
        return a;
    }
    int x1, y1, d;
    d = exgcd(b, a%b, x1, y1);
    x = y1;
    y = x1 - a/b * y1;
    return d;
}
signed main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout.tie(nullptr);
    int t;
    cin >> t;
    while(t--) {
        int a, b, c;
        cin >> a >> b >> c;
        int x0, y0; // 特解
        int d = exgcd(a, b, x0, y0);
        if(c % d != 0) {
            cout << -1 << endl;
            continue;
        }
        // c 的特解：
        int x1 = x0* c / d;
        int y1 = y0* c / d;
        int k1 = b/d;
        int k2 = a/d;
        int minx = (x1 % k1 + k1) % k1;
        minx = minx == 0 ? k1 : minx;
        int maxy = (c - a * minx) / b;
        int miny = (y1 % k2 + k2) % k2;
        miny = miny == 0 ? k2 : miny;
        int maxx = (c - b*miny) / a;
        if( maxy <= 0) {
            cout << minx << " " << miny << endl;
        } else {
            int t = (maxx - minx) / k1 + 1;
            cout << t << " " << minx << " " << miny << " " << maxx << " " << maxy << "\n";
        }
    }
    return 0;
}

11. 同余方程

同余方程

求同余方程 ax ≡ b(mod m) 的解。

例如：求 4x ≡ 2(mod 6)，其中一个解为 x = 2。

扩欧算法求同余方程的解

同余方程可以转换为二元一次不定方程：

• 由 ax ≡ b (mod m) 可得：ax = ym + b，移项得 ax − my=b；
• 令 y = −y，就变成二元一次不定方程：ax + my=b；
• 由裴蜀定理得，当 gcd(a,m)∣b 时，有解。

那么，就可以用扩展欧几里得算法求出 x 的值。

扩欧算法求乘法逆元

• 特殊的，如果 b 等于 1，原方程变为 ax≡1(modm)，解出的 x 实际上就是 a 在模 m 意义下的乘法逆元。

因此，利用扩欧算法也可以求出乘法逆元。并且不需要模数 m 是质数，仅需 gcd(a,m)=1，也就是 a,m 互质即可。

模板题：

同余方程

12. 逆元

乘法逆元

乘法逆元一般用于求 a/b (modp) 的值。先算出 b 在模 p 意义下的乘法逆元 b^−1，然后计算 a×b^−1(modp) 的值，就可以将除法转化为乘法。

下面总结一下求逆元的各种方法。

方法一：费马小定理 + 快速幂

• 问题：求 a 在模 m 意义下的逆元。
• 前提条件：a,m 互质，且 m 是质数。
• 方法：求出 a^(m−2)，就是 a 在模 m 意义下的逆元。
• 时间复杂度：O(logm)。

代码实现

// 必须要保证 a,m 互质，且 m 为质数。
LL qpow(LL a, LL b, LL m) {
    LL ret = 1;
    while(b) {
        if(b & 1) ret = ret * a % m;
        b >>= 1;
        a = a * a % m;
    }
    return ret;
}
int main() {
    LL x, m;
    cin >> x >> m;
    cout << qpow(x, m - 2, m) << endl;
    return 0;
}

方法二：扩欧算法求逆元

• 问题：求 a 在模 m 意义下的逆元。
• 前提条件：a,m 互质。
• 方法：解同余方程 ax≡1(modm)。
• 时间复杂度：O(logm)。

代码实现

LL exgcd(LL a, LL b, LL& x, LL& y) {
    if(b == 0) {
        x = 1, y = 0;
        return a;
    }
    LL x1, y1, d;
    d = exgcd(b, a % b, x1, y1);
    x = y1, y = x1 - a / b * y1;
    return d;
}
int main() {
    LL a, m;
    cin >> a >> m;
    LL x, y, d;
    d = exgcd(a, m, x, y);
    if(d == 1) cout << (x % m + m) % m << endl;
    return 0;
}

方法三：递推法

• 问题：求 1∼n 中所有的数在模 p 意义下的逆元。
• 前提条件：1∼n 中所有的数都与 p 互质。
• 方法：线性递推：
  • 当 i=1：inv[1]=1
  • 当 i>1：inv[i]=p−(p/i)×inv[pmodi]modp
• 时间复杂度：O(n)。

代码实现

LL n, p;
LL inv[N];
void get_inv() {
    inv[1] = 1;
    for(int i = 2; i <= n; i++) {
        inv[i] = p - p / i * inv[p % i] % p;
    }
}

递推原理：

模板题：

洛谷 P3811 模意义下的乘法逆元

方法四：打表 i 的阶乘的逆元（求组合数）

打表所有的阶乘

• f[i] 表示 i 的阶乘 mod p 的值，f[i]=f[i−1]×i。
• 这样，从左往右递推一遍，就可以把所有的阶乘全部维护出来。

模数 p 一般是一个较大的质数，比如 1e9+7。此时 p > n > m 时，并且 p 是质数，那么对于 [1,n] 里的所有数，p 都与它们互质，因此逆元一定存在。对于阶乘求逆元，也可以递推。

代码实现

int n, p;
LL f[N], g[N];
LL qpow(LL a, LL b, LL p) {
    LL ret = 1;
    while(b) {
        if(b & 1) ret = ret * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return ret;
}
void init() {
    f[0] = 1;
    for(int i = 1; i <= n; i++) {
        f[i] = f[i - 1] * i % p;
    }
    g[n] = qpow(f[n], p - 2, p);
    for(int i = n - 1; i >= 0; i--) {
        g[i] = (i + 1) * g[i + 1] % p;
    }
}
LL C(int n, int m) {
    if(n < m) return 0;
    return f[n] * g[n - m] % p * g[m] % p;
}

时间复杂度

• 打表：O(n)
• 查询：O(1)

13. 中国剩余定理

引入 - 线性同余方程组

南北朝时期的《孙子算经》中有个问题：有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

用数学公式翻译过来就是一个线性同余方程组：

线性同余方程组：求

的最小非负整数解。

中国剩余定理 (CRT)

• 前提：所有的模数 m1, m2 ,…, mn 两两互质。因此中国剩余定理的使用比较局限。

原理

中国剩余定理是基于'构造法'得出的结果。以下给出 n=3 的构造过程，n 等于任意数的构造过程是一样的。

对于方程组，记 M = m1 × m2 × m3，构造解：

只要能构造出这样的 x，那么 x 就是我们要的解。

中国剩余定理的构造方式

当 x 加减 M 的若干倍时，依旧是满足方程，因此最小非负整数解就是

n 取任意数的时候，构造方式也是同理：

例如

计算结果：

计算每一个方程的 Ci：

算法流程

代码实现

const int N = 1e5 + 10;
int n;
LL r[N], m[N];
LL qmul(LL a, LL b, LL p) {
    LL ret = 0;
    while(b) {
        if(b & 1) ret = (ret + a) % p;
        b >>= 1;
        a = (a + a) % p;
    }
    return ret;
}
void exgcd(LL a, LL b, LL &x, LL &y) {
    if(b == 0) {
        x = 1, y = 0;
        return;
    }
    LL x1, y1;
    exgcd(b, a % b, x1, y1);
    x = y1, y = x1 - a / b * y1;
}
LL crt() {
    LL M = 1;
    for(int i = 1; i <= n; i++) M *= m[i];
    LL ret = 0;
    for(int i = 1; i <= n; i++) {
        LL c = M / m[i];
        LL x, y;
        exgcd(c, m[i], x, y);
        x = (x % m[i] + m[i]) % m[i]; // 补成最小非负整数
        ret = (ret + qmul(qmul(c, x, M), r[i], M)) % M;
    }
    return ret;
}

模板题：

洛谷 P1495 中国剩余定理（CRT）/ 曹冲养猪

15. 扩展中国剩余定理

中国剩余定理的缺陷

中国剩余定理（CRT）只能处理所有模数两两互质的情况，如果模数不互质，CRT 就不适用了。

扩展中国剩余定理 - EXCRT

求解线性同余方程组：

算法思想：将任意整数 X 不断迭代成最终解。

算法步骤

例如

补上一个方程 x ≡ 0 (mod1)，最终解 ret = x × 1+0。

模板题：

洛谷 P4777 扩展中国剩余定理（EXCRT）

代码：

#include <iostream>
using namespace std;
typedef long long LL;
const int N = 1e5 + 10;
int n;
LL m[N], r[N];
LL qmul(LL a, LL b, LL p) {
    LL sum = 0;
    while(b) {
        if(b & 1) sum = (sum + a) % p;
        b >>= 1;
        a = (a + a) % p;
    }
    return sum;
}
LL exgcd(LL a, LL b, LL& x, LL& y) {
    if(b == 0) {
        x = 1, y = 0;
        return a;
    }
    LL x1, y1;
    LL d = exgcd(b, a % b, x1, y1);
    x = y1, y = x1 - a / b * y1;
    return d;
}
LL excrt() {
    LL M = 1, ret = 0;
    for(int i = 1; i <= n; i++) {
        LL a = M, b = m[i], c = r[i] - ret; // 防溢出：把 c 补成最小非负整数
        c = (c % b + b) % b;
        // ax + by = c
        LL x, y, d;
        d = exgcd(a, b, x, y);
        if(c % d) return -1; // x -> ax + by = 1
        // x -> ax + by = c -> x * (c / d)
        LL k1 = b / d;
        x = qmul(x, c / d, k1);
        x = (x % k1 + k1) % k1;
        ret = ret + x * M;
        M = k1 * M;
        ret = (ret % M + M) % M;
    }
    return ret;
}
int main() {
    cin >> n;
    for(int i = 1; i <= n; i++) cin >> m[i] >> r[i];
    cout << excrt() << endl;
    return 0;
}

16. 求两点的连线上对多少个顶点

可以求两点横纵坐标的差值的 gcd