【学习笔记】初等数论

[学习笔记]初等数论

最大公约数 \(\gcd\)

对于两个整数 \(a, b\)，有 \(\gcd(a, b) \times \operatorname{lcm}(a, b) = a \times b\)。

欧几里得算法（辗转相除法）：

\[\gcd(a, b) = \gcd(b, a \bmod b) \]

代码：

int gcd(int a, int b){
    return b ? gcd(b, a%b) : a;
}

或者直接使用 __gcd(a, b)。

辗转相减法：

\[\gcd(a, b) = \gcd(a, b-a) \]

推广到 \(n\) 项：

\[\gcd(a_1, a_2, \dots, a_n) = \gcd(a_1, a_2-a_1, a_3-a_2,\dots, a_n-a_{n-1}) \]

注意：因为 \(\gcd\) 不可能是负数，所以当权值为负时答案应该取绝对值。

线性筛

欧拉筛。不多讲了。

int n;
int prime[N], cnt;
bool is_prime[N];

void getprime(){
    for(int i=2; i<=n; i++) is_prime[i] = 1;
    for(int i=2; i<=n; i++){
        if(is_prime[i]){
            prime[++cnt] = i;
        }
        for(int j=1; j<=cnt && i*prime[j]<=n; j++){
            is_prime[i*prime[j]] = 0;
            if(i%prime[j] == 0) break;
        }
    }
}

积性函数

数论函数：

数论函数指定义域为正整数的函数。数论函数也可以视作一个数列。

积性函数：

若一个数论函数 \(f\) 满足，对所有互质的 \(n, m\)，都满足 \(f(nm) = f(n)f(m)\)，则称 \(f\) 是一个积性函数。

此外，如果对所有 \(n, m\)，都满足 \(f(nm) = f(n)f(m)\)，则称 \(f\) 为完全积性函数。

性质：都可以通过筛法来快速求解。

常见的积性函数：

欧拉函数 \(\varphi(n)\)：

表示 \(\left[ 1, n \right]\) 中与 \(n\) 互质的数的数量。

欧拉公式：记 \(n = \prod p_i^{c_i}\)，则 \(\varphi(n) = \prod (p_i-1)p_i^{c_i-1}\)。

性质：
\(\sum\limits_{d \mid n} \varphi(d) = n\)

莫比乌斯函数 \(\mu(n)\)：

\(\mu(n) = \begin{cases} (-1)^k & n = p_1 \hspace{0.1em} p_2 \dots p_k \hspace{0.3em} (\hspace{0.1em}p_i\hspace{0.1em}是质数) \\ 0 & n \hspace{0.1em} 的某个质因子指数 \hspace{0.1em} c \ge 2 \end{cases}\)

性质：
\(\sum\limits_{d \mid n} \mu(d) = \left[n = 1\right]\)

单位函数 \(\epsilon(n)\)：

\(\epsilon(n) = \begin{cases} 1 & n = 1 \\ 0 & n \neq 1 \end{cases}\)

筛 \(\varphi(n)\)

从欧拉筛出发考虑。下面的 \(i\) 指当前的数，\(p\) 指质数表中的质数。

• 当 \(i\) 是质数时，显然 \(\varphi(i) = i-1\)。

• 当 \(i \not\equiv 0 \hspace{0.3em} (\bmod \hspace{0.3em} p)\) 时，即 \(i, p\) 互质时，根据积性函数定义得：\(\varphi(ip) = \varphi(i)\varphi(p) = \varphi(i) \times (p-1)\)。

• 当 \(i \equiv 0 \hspace{0.3em} (\bmod \hspace{0.3em} p)\) 时，即 \(p\) 是 \(i\) 的一个因子时。设 \(i\) 唯一分解中 \(p\) 的指数是 \(c\)，因为 \(\varphi(p^c) = (p-1)p^{c-1}\)，因此 \(\varphi(ip) = \varphi(i) \times p\)。

  解释：

  • 对于 \(\varphi(p^c) = (p-1)p^{c-1}\)，\(\left[ 1, p\right]\) 内有 \(p-1\) 个数与其互质。所以 \(\left[ 1, p^c\right]\) 内有 \((p-1)p^{c-1}\) 个数与其互质。

  • 对于 \(\varphi(ip) = \varphi(i) \times p\)。因为 \(\varphi(ip) = \varphi(\frac{i}{p^c}) \cdot \varphi(p^{c+1})\)，又因为 \(\frac{i}{p^c}\) 与 \(p^c\) 互质，所以 \(\varphi(\frac{i}{p^c}) \varphi(p^c) = \varphi(i)\)。又因为 \(\varphi(p^{c+1}) = (p-1)p^{c}\)。所以将 \(\varphi(\frac{i}{p^c})\) 和 \(\varphi(p^{c+1})\) 两项拆开消元后即可得到 \(\varphi(ip) = \varphi(i) \times p\)。

    也可通过欧拉公式（欧拉函数的定义）理解。

代码实现：

int n;
int prime[N], phi[N], cnt;
bool is_prime[N];

void getprime(){
    phi[1] = 1;
    for(int i=2; i<=n; i++) is_prime[i] = 1;
    for(int i=2; i<=n; i++){
        if(is_prime[i]){
            phi[i] = i-1;
            prime[++cnt] = i;
        }
        for(int j=1; j<=cnt && i*prime[j]<=n; j++){
            is_prime[i*prime[j]] = 0;
            if(i%prime[j] == 0){
                phi[i*prime[j]] = prime[j]*phi[i];
                break;
            }
            else phi[i*prime[j]] = (prime[j]-1)*phi[i];
        }
    }
}

筛 \(\mu(n)\)

还是考虑三种情况。

• 当 \(i\) 是质数时，显然 \(\mu(i) = -1\)。
• 当 \(i \not\equiv 0 \hspace{0.3em} (\bmod \hspace{0.3em} p)\) 时，即 \(i, p\) 互质时，根据积性函数定义得：\(\mu(ip) = \mu(i)\mu(p) = -\mu(i)\)。
• 当 \(i \equiv 0 \hspace{0.3em} (\bmod \hspace{0.3em} p)\) 时，即 \(p\) 是 \(i\) 的一个因子时。根据该函数定义得 \(\mu(ip) = 0\)。

代码实现应该差不多。

---

P2568 GCD

首先我们枚举质数：

\[\sum_{p \in \rm prime} \sum_{i=1}^{n} \sum_{j=1}^{n} \left[ \gcd(i, j) = p \right] \]

变形：

\[\sum_{p \in \rm prime} \sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{p} \rfloor} \left[ \gcd(i, j) = 1 \right] \]

接下来改变 \(j\) 的枚举上界（其中 \(-1\) 的原因是 \(i=j=1\) 时的答案会被重复统计，因此这里的 \(-1\) 是在第二层中的）：

\[\sum_{p \in \rm prime} \left( \sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} \left( 2\sum_{j=1}^{i} \left[ \gcd(i, j) = 1 \right] \right) -1 \right) \]

此时发现 \(\sum\limits_{j=1}^{i} \left[ \gcd(i, j) = 1 \right]\) 就是 \(\varphi(i)\)，故原式可化为：

\[\sum_{p \in \rm prime} \left( 2\sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} \varphi(i) -1 \right) \]

可以前缀和求 \(\sum \varphi(i)\) 的值，最后枚举 \(p \in \rm prime\) \({\rm and}\) \(p \le n\) 并统计答案即可。

时间复杂度：\(O(n)\)。

Code

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 1e7 + 5;

int prime[N], phi[N], cnt;
bool is_prime[N];
ll ans, sum[N];

int main(){
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int n; cin>>n;
    phi[1] = 1;
    for(int i=2; i<=n; i++) is_prime[i] = 1;
    for(int i=2; i<=n; i++){
        if(is_prime[i]){
            phi[i] = i-1;
            prime[++cnt] = i;
        }
        for(int j=1; j<=cnt && i*prime[j]<=n; j++){
            is_prime[i*prime[j]] = 0;
            if(i%prime[j] == 0){
                phi[i*prime[j]] = prime[j]*phi[i];
                break;
            }
            else phi[i*prime[j]] = (prime[j]-1)*phi[i];
        }
    }
    for(int i=1; i<=n; i++) sum[i] = sum[i-1]+phi[i];
    for(int i=1; i<=cnt; i++) ans += 2*sum[n/prime[i]]-1;
    cout<<ans;
    return 0;
}

P10463 Interval GCD

辗转相减法（前置知识）：

\[\gcd(a, b) = \gcd(a, b-a) \]

推广到 \(n\) 项：

\[\gcd(a_1, a_2, \dots, a_n) = \gcd(a_1, a_2-a_1, a_3-a_2,\dots, a_n-a_{n-1}) \]

看到题意中把 \(A_l, A_{l+1}, \dots, A_r\) 都加上 \(d\)，自然想到差分。

令差分完之后的数组为 \(B\)，那么 \(\gcd(B_{l+1}, B_{l+2}, \dots, B_r) = \gcd(A_{l+1}, A_{l+2}, \dots A_r)\)，而 \(\gcd(B_l,B_{l+1}) \neq \gcd(A_l, A_{l+1})\)，所以 \(A_l\) 需要单独维护。

以下用线段树维护差分数组的 \(\gcd\)，用树状数组维护差分数组，最后通过前缀和求得 \(A_l\)。

Code

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 5e5+5;
#define lowbit(x) x&(-x)
#define ls(x) ((x)<<1)
#define rs(x) ((x)<<1|1)
#define gcd(x, y) abs(__gcd((x), (y)))

int n, m;
ll a[N], cha[N], c[N];
struct node{
    int l, r;
    ll v;
}tr[N<<2];

void add(int x, ll y){
    while(x <= n){
        c[x] += y;
        x += lowbit(x);
    }
}

ll sum(int x){
    ll ret = 0;
    while(x){
        ret += c[x];
        x -= lowbit(x);
    }
    return ret;
}

void pushup(int x){
    tr[x].v = gcd(tr[ls(x)].v, tr[rs(x)].v);
}

void buildtr(int x, int l, int r){
    tr[x].l = l, tr[x].r = r;
    if(l == r){
        tr[x].v = cha[l];
        return;
    }
    int mid = (l+r) >> 1;
    buildtr(ls(x), l, mid);
    buildtr(rs(x), mid+1, r);
    pushup(x);
}

void upd(int x, int tar, ll d){
    if(tr[x].l > tar || tr[x].r < tar)
        return;
    if(tr[x].l == tr[x].r && tr[x].l == tar){
        tr[x].v += d;
        return;
    }
    int mid = (tr[x].l+tr[x].r) >> 1;
    if(tar<=mid) upd(ls(x), tar, d);
    upd(rs(x), tar, d);
    pushup(x);
}

ll query(int x, int l, int r){
    if(tr[x].l > r || tr[x].r < l)
        return 0;
    if(tr[x].l >= l && tr[x].r <= r)
        return tr[x].v;
    return gcd(query(ls(x), l, r), query(rs(x), l, r));
}

int main(){
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    cin>>n>>m;
    for(int i=1; i<=n; i++){
        cin>>a[i];
        cha[i] = a[i] - a[i-1];
        add(i, cha[i]);
    }
    buildtr(1, 1, n);
    while(m--){
        char op; int l, r; cin>>op>>l>>r;
        if(op == 'C'){
            ll d; cin>>d;
            upd(1, l, d);
            upd(1, r+1, -d);
            add(l, d);
            add(r+1, -d);
        } else if(op == 'Q'){
            cout<<gcd(sum(l), query(1, l+1, r))<<"\n";
        }
    }
    return 0;
}

USACO08DEC Patting Heads S

从枚举倍数的角度思考：对于一个数 \(i\)，若它在原数组中出现了 \(cnt_i\) 次，那么 \(i\) 这个数会对它的每一个倍数产生 \(cnt_i\) 的贡献。这样就可以通过查询这样产生的贡献的总和来计算答案了。

注意最后要减去一个对自己的贡献。

Code

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 1e6+5;

int a[N], w[N], c[N];

int main(){
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int m = 0; // 值域
    int n; cin>>n;
    for(int i=1; i<=n; i++){
        cin>>a[i];
        c[a[i]]++;
        m = max(m, a[i]);
    }
    for(int i=1; i<=m; i++){
        for(int j=i; j<=m; j+=i){
            w[j] += c[i];
        }
    }
    for(int i=1; i<=n; i++)
        cout<<w[a[i]]-1<<"\n"; // 减去自己对自己的贡献
    return 0;
}

Divan and Kostomuksha (easy version)

考虑 \(dp_i\) 表示目前前缀 \(\gcd = i\)，前缀 \(\gcd\) 和的最大值。

首先计算出序列 \(a\) 中是 \(i\) 的倍数的数量，记为 \(cnt_i\)。（这里处理方法和上一题相似）

接下来考虑转移，假设是从 \(j\) 转移，转移到 \(i\)，相当于在 \(j+1\) 至 \(f_j\) 所在位置全部填上 \(a\) 中是 \(i\) 的倍数但不是 \(j\) 的倍数的数，容易发现多的是 \(i\) 的倍数的数量为 \(cnt_i - cnt_j\)。所以转移式为：

\[dp_i = \max_{i \mid j}(dp_j + (cnt_i - cnt_j) \times i) \]

Code

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 5e6;

ll cnt[N+5], a[N+5];
ll dp[N+5];
//cnt 表示 a 中能被 i 整除的数的个数
int main(){
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    ll n, m = 0; cin>>n; // m 是值域
    for(int i=1; i<=n; i++){
        cin>>a[i];
        cnt[a[i]]++;
        m = max(m, a[i]);
    }
    for(int i=1; i<=m; i++){
        for(int j=i*2; j<=m; j+=i){
            cnt[i] += cnt[j];
        }
    }
    ll ans = 0;
    for(int i=m; i>=1; i--){
        dp[i] = cnt[i]*i;
        for(int j=i*2; j<=m; j+=i){
            dp[i] = max(dp[i], dp[j]+1ll*(cnt[i]-cnt[j])*i);
        }
        ans = max(ans, dp[i]);
    }
    cout<<ans;
    return 0;
}

费马小定理

若 \(p\) 为质数，\(\gcd(a, p) = 1\)，则 \(a^{p-1} \equiv 1 \, (\bmod \, p)\)。

所以可以求得 \(a\) 在模意义下的乘法逆元为 \(a^{p-2}\)。快速幂求解即可。

裴蜀定理

设 \(a, b\) 是不全为零的整数，那么对任意整数 \(x, y\)，满足 \(\gcd(a, b) \mid ax+by\)，且存在整数 \(x, y\)，使得 \(ax+by=\gcd(a, b).\)

简单来说，我们设 \(d = \gcd(a, b)\)，那么对于方程 \(ax + by = d\)，一定存在一组整数解。并且对于方程 \(ax + by = z\)，如果满足 \(d \mid z\)，那么方程一定有整数解，否则无整数解。

扩展欧几里得算法（Exgcd）

应用：

1. 求解不定方程。
2. 求解模意义下的逆元。
3. 求解线性同余方程。

求解 \(ax+by = \gcd(a, b).\)

设当前的式子为 \(ax+by = \gcd(a, b)\)，则肯定存在式子 \(bx^{\prime}+(a \bmod b)y^{\prime} = \gcd(b, a \bmod b).\)

根据欧几里得算法，\(\gcd(a, b) = \gcd(b, a \bmod b).\)

\[\therefore ax + by = bx^{\prime} + (a \bmod b)y^{\prime} \]

\[\because a \bmod b = a - \lfloor \frac{a}{b} \rfloor \times b \]

\[\therefore ax + by = bx^{\prime} + \left( a - \lfloor \frac{a}{b} \rfloor \times b \right) y^{\prime} \]

\[\therefore ax + by = bx^{\prime} + ay^{\prime} - \lfloor \frac{a}{b} \rfloor \times b \times y^{\prime} \]

\[\therefore ax + by = ay^{\prime} + b \cdot \left( x^{\prime} - \lfloor \frac{a}{b} \rfloor \times y^{\prime} \right) \]

所以得到 \(\begin{cases}x = y^{\prime} \\ y = x^{\prime} - \lfloor \frac{a}{b} \rfloor \times y^{\prime} \end{cases}\)

当 \(b=0\) 时得到 base case：\(\begin{cases}x = 1 \\ y = 0 \end{cases}\)，所以递归求解即可。

求解二元一次不定方程 \(ax+by = c.\)

根据裴蜀定理得：

若 \(\gcd(a, b) \nmid c\)，则无解，故问题转化为求解 \(ax+by = \gcd(a, b).\)

过程：

1. 用 \(ax+by = \gcd(a, b)\) 的特解来求 \(ax+by = c\) 的特解。

   设 \(d = \gcd(a, b), c = k \cdot d\)，已知 \(ax_0+by_0 = d\)，两边同乘 \(k\)，得 \(ax_0 \cdot k + by_0 \cdot k = d \cdot k = c\)，所以 \(\begin{cases}x = x_0 \cdot k \\ y = y_0 \cdot k \end{cases}\) 就是原方程的一组解。

2. 用这组特解 \(x_1, y_1\) 找出方程所有的解。

   设 \(a(x_1+m) + b(y_1+n) = c\)，展开得 \(ax_1+am+by_1+bn = c\)，所以只需要让 \(am+bn = 0\)。

   继续设 \(d = \gcd(a, b)\)，我们让 \( \begin{cases} m = t \times \frac{b}{d} \\ n = -t \times \frac{a}{d} \end{cases}\)，即可成立。

   所以只要任取一个 \(t\)，就可以算出相应的 \(m, n\)，进而推出 \(x, y\)，即可得到通解。

求解线性同余方程 \(ax \equiv b \pmod n.\)

以下内容摘自OI_wiki。

定理 1：线性同余方程 \(ax\equiv b \pmod n\) 可以改写为如下线性不定方程：

\[ax + nk = b \]

（注：\(ax \div n = -k \dots b\)）

其中 \(x\) 和 \(k\) 是未知数。这两个方程是等价的，有整数解的充要条件为 \(\gcd(a,n) \mid b\)。

应用扩展欧几里德算法可以求解该线性不定方程。根据定理 1，对于线性不定方程 \(ax+nk=b\)，可以先用扩展欧几里得算法求出一组 \(x_0,k_0\)，也就是 \(ax_0+nk_0=\gcd(a,n)\)，然后两边同时除以 \(\gcd(a,n)\)，再乘 \(b\)。就得到了方程

\[a\dfrac{b}{\gcd(a,n)}x_0+n\dfrac{b}{\gcd(a,n)}k_0=b \]

于是找到方程的一个解。

定理 2：若 \(\gcd(a,n)=1\)，且 \(x_0、k_0\) 为方程 \(ax+nk=b\) 的一组解，则该方程的任意解可表示为：

\[x=x_0+nt \]

\[k=k_0-at \]

并且对任意整数 \(t\) 都成立。

根据定理 2，可以从已求出的一个解，求出方程的所有解。实际问题中，往往要求出一个最小整数解，也就是一个特解

\[x=(x \bmod t+t) \bmod t \]

其中有

\[t=\dfrac{n}{\gcd(a,n)} \]

如果仔细考虑，用扩展欧几里得算法求解与用逆元求解，两种方法是等价的。

P5656 【模板】二元一次不定方程 (exgcd)

在上述求解二元一次不定方程的基础上还要增加第三步：考虑解的最大或最小值。

令 \(d = \gcd(a, b)\)，首先：\(\begin{cases} x = x_0 + t \times \frac{b}{d} \\ y = y_0 - t \times \frac{a}{d} \end{cases}\)

可以看出，当 \(t\) 增大的时候，\(x\) 越来越大，\(y\) 越来越小。

假设任选的 \(t > 0\)。接下来设 \(d_x = \frac{b}{d}, d_y = \frac{a}{d}\)，也就是\(\begin{cases} x = x_0 + t \cdot d_x \\ y = y_0 - t \cdot d_y \end{cases}\)。

若限制 \(x > 0\)，则 \(x_0+td_x > 0\)，\(t > -\frac{x_0}{d_x}\)。

若限制 \(y > 0\)，则 \(y_0-td_y > 0\)，\(t < \frac{y_0}{d_y}\)。

因为 \(t\) 一定是整数，所以有 \(\lceil \frac{-x_0+1}{d_x} \rceil \le t \le \lfloor \frac{y_0-1}{d_y} \rfloor\)。

• 对不等式左侧解释，因为 \(> \rightarrow \ge\)，而 \(t\) 一定为整数，所以分式分子要 \(+1\)，排除 \(- \frac{x_0}{d_x}\) 刚好为整数的可能，右侧同理。

若 \(\lceil \frac{-x_0+1}{d_x} \rceil > \lfloor \frac{y_0-1}{d_y} \rfloor\)，那么两个解不可能都是正整数。

• 此时，答案即为没有正整数解但有整数解，\(x\) 为正整数且最小即为 \(t=\lceil \frac{-x_0+1}{d_x} \rceil\) 时 \(x\) 的值，\(y\) 为正整数且最小即为 \(t = \lfloor \frac{y_0-1}{d_y} \rfloor\) 时 \(y\) 的值。

否则就有正整数解，正整数解个数即为 \(t\) 合法取值范围内的整数数量。

• 此时 \(x\) 的最大值和 \(y\) 的最小值都是在 \(t\) 最大时取到，而
  \(x\) 的最小值和 \(y\) 的最大值都是在 \(t\) 最小时取到，直接计算即可。

Code

#include <bits/stdc++.h>
using namespace std;
#define ll long long

ll x, y;

ll exgcd(ll a, ll b){
    if(b == 0){
        x = 1, y = 0;
        return a;
    }
    ll d = exgcd(b, a % b);
    ll t = x;
    x = y;
    y = t - a/b * y;
    return d;
}

int main(){
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int T; cin>>T;
    while(T--){
        ll a, b, c; cin>>a>>b>>c;
        ll d = exgcd(a, b);
        if(c % d != 0){
            cout<<"-1\n";
            continue;
        }
        x = x*c/d, y = y*c/d;
        ll dx = b/d, dy = a/d;
        ll down = ceil((-x+1)*1.0/dx), up = floor((y-1)*1.0/dy);
        if(down > up){
            cout<<x+down*dx<<" "<<y-up*dy<<"\n";
            continue;
        }
        cout<<up-down+1<<" "<<x+down*dx<<" "<<y-up*dy<<" "<<x+up*dx<<" "<<y-down*dy<<"\n";
    }   
    return 0;
}

中国剩余定理（CRT）

以下内容大部分摘自OI_wiki。

中国剩余定理可求解如下形式的一元线性同余方程组（其中 \(n_1, n_2, \cdots, n_k\) 两两互质）：

\[\begin{cases} x \equiv a_1 \pmod {n_1} \\ x \equiv a_2 \pmod {n_2} \\ \vdots \\ x \equiv a_k \pmod {n_k} \\ \end{cases} \]

求解过程：

1. 计算所有模数的积 \(n\)。

2. 对于第 \(i\) 个方程：

   a. 计算 \(m_i=\frac{n}{n_i}\);

   b. 计算 \(m_i\) 在模 \(n_i\) 意义下的逆元 \(m_i^{-1}\);

   c. 计算 \(c_i=m_im_i^{-1}\)（不要对 \(n_i\) 取模）。

3. 方程组在模 \(n\) 意义下的唯一解为：\(x=\sum_{i=1}^k a_ic_i \pmod n\)。

证明：

我们需要证明上面算法计算所得的 \(x\) 对于任意 \(i=1,2,\cdots,k\) 满足 \(x\equiv a_i \pmod {n_i}\)。

枚举 \(j \in [1, k]\)，当 \(i\neq j\) 时，因为 \(n_i \mid \frac{n}{n_j}\)，所以有 \(m_j \equiv 0 \pmod {n_i}，故 c_j \equiv m_j \equiv 0 \pmod {n_i}\)。又有 \(c_i \equiv m_i \cdot (m_i^{-1} \bmod {n_i}) \equiv 1 \pmod {n_i}\)，所以我们有：

\[\begin{aligned} x&\equiv \sum_{j=1}^k a_jc_j &\pmod {n_i} \\ &\equiv a_ic_i &\pmod {n_i} \\ &\equiv a_i \cdot m_i \cdot (m^{-1}_i \bmod n_i) &\pmod {n_i} \\ &\equiv a_i &\pmod {n_i} \end{aligned} \]

即对于任意 \(i=1,2,\cdots,k\)，上面算法得到的 \(x\) 总是满足 \(x\equiv a_i \pmod{n_i}\)，即证明了解同余方程组的算法的正确性。

因为我们没有对输入的 \(a_i\) 作特殊限制，所以任何一组输入 \(\{a_i\}\) 都对应一个解 \(x\)。

另外，若 \(x\neq y\)，则总存在 \(i\) 使得 \(x\) 和 \(y\) 在模 \(n_i\) 下不同余。

故系数列表 \(\{a_i\}\) 与解 \(x\) 之间是一一映射关系，方程组总是有唯一解。

P1495 【模板】中国剩余定理（CRT）

板子。注意求逆元方式。

Code

#include <bits/stdc++.h>
using namespace std;
#define ll long long

int n;
int a[100005], b[100005];
ll sum = 1, x, y;
ll ans;

void exgcd(ll a, ll b){
    if(b == 0){
        x = 1, y = 0;
        return;
    }
    exgcd(b, a%b);
    ll z = x;
    x = y, y = z - a/b * y;
}

int main(){
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    cin>>n;
    for(int i=1; i<=n; i++){
        cin>>a[i]>>b[i];
        sum *= a[i];
    }
    for(int i=1; i<=n; i++){
        ll m = sum/a[i];
        x = 0, y = 0;
        exgcd(m, a[i]);
        ans = (ans + (__int128)b[i]*m*(x<0 ? x+a[i] : x)) % sum;
    }
    cout<<ans;
    return 0;
}

扩展中国剩余定理（ExCRT）

给定 \(n\) 组非负整数 \(a_i, m_i\) ，求解关于 \(x\) 的方程组的最小非负整数解（注意这里 \(m_1, m_2, \cdots, m_n\) 并不一定两两互质）。

\[\begin{cases}x\equiv a_1\pmod{m_1}\\x\equiv a_2\pmod{m_2}\\\vdots\\x\equiv a_n\pmod{m_n}\end{cases} \]

先考虑只有两个方程的情况。设两个方程分别是 \(x\equiv a_1 \pmod {m_1}、x\equiv a_2 \pmod {m_2}\)。

将它们转化为不定方程：\(x=m_1p+a_1=m_2q+a_2\)，其中 \(p, q\) 是整数，则有 \(m_1p-m_2q=a_2-a_1\)。

由裴蜀定理，当 \(\gcd(m_1,m_2) \nmid (a_2-a_1)\) 时，无解。

其他情况下，可以通过扩展欧几里得算法解出来一组可行解 \((p, q)\)。

则原来的两方程组成的模方程组的解为 \(x\equiv b\pmod M\)，其中 \(b=m_1p+a_1，M=\text{lcm}(m_1, m_2)\)。

多个方程用上面的方法两两合并即可，相当于跑 \(n-1\) 次扩欧。

P4777 【模板】扩展中国剩余定理（EXCRT）

Code

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define i128 __int128

i128 x, y;
i128 a1, a2, m1, m2;

i128 exgcd(i128 a, i128 b){
    if(b == 0){
        x = 1, y = 0;
        return a;
    }
    i128 d = exgcd(b, a % b);
    i128 t = x;
    x = y, y = t - a/b * y;
    return d;
}

void merge(){
    i128 c = a2 - a1;
    i128 d = exgcd(m1, m2);
    if(c % d != 0){
        cout<<"-1";
        exit(0);
    }
    i128 M = m1 * m2 / d;
    x = (x * c / d % (m2 / d) + (m2 / d)) % (m2 / d);
    a1 = ((m1 * x + a1) % M + M) % M;
    m1 = M;
}

int main(){
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int n; cin>>n;
    for(int i=1; i<=n; i++){
        ll A_, B_; cin>>A_>>B_; // A 是模数
        m2 = A_, a2 = B_;
        if(i != 1) merge();
        else m1 = A_, a1 = B_;
    }
    cout<<(ll)a1;
    return 0;
}

同余最短路

当出现形如「给定 \(n\) 个整数，求这 \(n\) 个整数能拼凑出多少的其他整数（\(n\) 个整数可以重复取）」，以及「给定 \(n\) 个整数，求这 \(n\) 个整数不能拼凑出的最小（最大）的整数」，或者「至少要拼几次才能拼出模 \(K\) 余 \(p\) 的数」的问题时可以使用同余最短路的方法。

同余最短路利用同余来构造一些状态，可以达到优化空间复杂度的目的。

类比差分约束方法，利用同余构造的这些状态可以看作单源最短路中的点。同余最短路的状态转移通常是这样的 \(f(i+y) = f(i) + y\)，类似单源最短路中 \(f(v) = f(u) +edge(u,v)\)。

P3403 跳楼机

题目大意：给定 \(x, y, z, h\)，对于 \(k \in [1,h]\)，有多少个 \(k\) 能够满足 \(ax+by+cz=k。（0\leq a,b,c，1\le x,y,z\le 10^5，h\le 2^{63}-1）\)

不妨假设 \(x < y < z\)。并让 \(h \rightarrow h-1\)，这样对 \(x\) 取模比较好表示。（通常选取一组 \(a_i\) 中最小的那个数对它取模，也就是此处的 \(x\)，这样可以尽量减小空间复杂度（剩余系最小）。）

首先利用同余的方式分类，即将所有楼层模 \(x\) 的值分成 \(x\) 类。不难发现，若某层可以到达，则所有同类楼层中，高于它的部分均可到达。

因此，我们只需要求出每类楼层中最低可到达的楼层，即可求出这一类中有几层可达。对于每一类楼层都用一个点来代替，比如点 \(a\) 包括了第 \(a, a+x, a+2x, \dots\) 层。显然 \(0 \le a < x\)。

以点 \(0\) 为起点，到点 \(a\) 的最短距离表示其中最低可达的楼层。

接下来在图中连边，点 \(i\) 一共向外连两条边：

\[i \xrightarrow{y} (i+y) \bmod x \]

\[i \xrightarrow{z} (i+z) \bmod x \]

上述两种边分别代表通过向上跳 \(y\) 和 \(z\) 层，到达其他类的操作。

用 Dijkstra 算法（SPFA 也可）即可求出每类楼层中最低可达的楼层（也就是最短距离），最后再统计答案即可。

Code

#include <bits/stdc++.h>
using namespace std;
#define ll unsigned long long

ll h, dis[100005], ans;
int x, y, z;

struct node{
    int fi; ll se;
    friend bool operator<(node a, node b){
        return a.se > b.se;
    }
};
vector<node> g[100005];
priority_queue<node> q;
bool vis[100005];

void dij(){
	for(int i=1; i<x; i++)
		dis[i] = INT64_MAX;
	dis[0] = 0;
	q.push({0, 0});
	while(!q.empty()){
		node t = q.top(); q.pop();
		if(vis[t.fi]) continue;
		vis[t.fi] = true;
		for(auto p : g[t.fi]){
			if(dis[t.fi]+p.se < dis[p.fi]){
				dis[p.fi] = dis[t.fi]+p.se;
				q.push({p.fi, dis[p.fi]});
			}
		}
	}
}

int main(){
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    cin>>h>>x>>y>>z;
    h--;
    for(int i=0; i<x; i++){
        g[i].push_back({(i+y)%x, y});
        g[i].push_back({(i+z)%x, z});
    }
    dij();
    for(int i=0; i<x; i++){
        if(h >= dis[i])
            ans += (h - dis[i]) / x + 1;
    }
    cout<<ans;
    return 0;
}