【学习笔记】初等数论

[学习笔记]初等数论

最大公约数 $gcd$

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

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

$$gcd(a,b)=gcd(b,amodb)$$

代码：

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$ 为完全积性函数。

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

常见的积性函数：

欧拉函数 $\phi (n)$：

表示 $[1,n]$ 中与 $n$ 互质的数的数量。

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

性质：
$\sum _{d\mid n}\phi (d)=n$

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

$\mu (n)=\{\begin{array}{ll}(-1{)}^k& n=p_1{p}_2\dots p_k(p_i\mathrm{是}\mathrm{质}\mathrm{数})\\ 0& n\mathrm{的}\mathrm{某}\mathrm{个}\mathrm{质}\mathrm{因}\mathrm{子}\mathrm{指}\mathrm{数}c\ge 2\end{array}$

性质：
$\sum _{d\mid n}\mu (d)=[n=1]$

单位函数 $ϵ(n)$：

$ϵ(n)=\{\begin{array}{ll}1& n=1\\ 0& n\ne 1\end{array}$

筛 $\phi (n)$

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

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

• 当 $i\not\equiv 0(modp)$ 时，即 $i,p$ 互质时，根据积性函数定义得：$\phi (ip)=\phi (i)\phi (p)=\phi (i)\times (p-1)$。

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

  解释：

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

  • 对于 $\phi (ip)=\phi (i)\times p$。因为 $\phi (ip)=\phi (\frac{i}{p^c})\cdot \phi (p^{c+1})$，又因为 $\frac{i}{p^c}$ 与 $p^c$ 互质，所以 $\phi (\frac{i}{p^c})\phi (p^c)=\phi (i)$。又因为 $\phi (p^{c+1})=(p-1)p^c$。所以将 $\phi (\frac{i}{p^c})$ 和 $\phi (p^{c+1})$ 两项拆开消元后即可得到 $\phi (ip)=\phi (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(modp)$ 时，即 $i,p$ 互质时，根据积性函数定义得：$\mu (ip)=\mu (i)\mu (p)=-\mu (i)$。
• 当 $i\equiv 0(modp)$ 时，即 $p$ 是 $i$ 的一个因子时。根据该函数定义得 $\mu (ip)=0$。

代码实现应该差不多。

---

P2568 GCD

首先我们枚举质数：

$$\sum _{p\in \mathrm{p}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}}\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)=p]$$

变形：

$$\sum _{p\in \mathrm{p}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}}\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{p}\rfloor }[gcd(i,j)=1]$$

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

$$\sum _{p\in \mathrm{p}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}}(\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\left(2\sum _{j=1}^i[gcd(i,j)=1]\right)-1)$$

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

$$\sum _{p\in \mathrm{p}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}}(2\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\phi (i)-1)$$

可以前缀和求 $\sum \phi (i)$ 的值，最后枚举 $p\in \mathrm{p}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}$ $\mathrm{a}\mathrm{n}\mathrm{d}$ $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})\ne 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$，若它在原数组中出现了 $cn{t}_i$ 次，那么 $i$ 这个数会对它的每一个倍数产生 $cn{t}_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)

考虑 $d{p}_i$ 表示目前前缀 $gcd=i$，前缀 $gcd$ 和的最大值。

首先计算出序列 $a$ 中是 $i$ 的倍数的数量，记为 $cn{t}_i$。（这里处理方法和上一题相似）

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

$$d{p}_i=\underset{i\mid j}{max}(d{p}_j+(cn{t}_i-cn{t}_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(modp)$。

所以可以求得 $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)$，则肯定存在式子 $b{x}^{\prime }+(amodb)y^{\prime }=gcd(b,amodb).$

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

$$\therefore ax+by=b{x}^{\prime }+(amodb)y^{\prime }$$

$$\because amodb=a-\lfloor \frac{a}{b}\rfloor \times b$$

$$\therefore ax+by=b{x}^{\prime }+(a-\lfloor \frac{a}{b}\rfloor \times b)y^{\prime }$$

$$\therefore ax+by=b{x}^{\prime }+a{y}^{\prime }-\lfloor \frac{a}{b}\rfloor \times b\times y^{\prime }$$

$$\therefore ax+by=a{y}^{\prime }+b\cdot (x^{\prime }-\lfloor \frac{a}{b}\rfloor \times y^{\prime })$$

所以得到 $\{\begin{array}{l}x=y^{\prime }\\ y=x^{\prime }-\lfloor \frac{a}{b}\rfloor \times y^{\prime }\end{array}$

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

求解二元一次不定方程 $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$，已知 $a{x}_0+b{y}_0=d$，两边同乘 $k$，得 $a{x}_0\cdot k+b{y}_0\cdot k=d\cdot k=c$，所以 $\{\begin{array}{l}x=x_0\cdot k\\ y=y_0\cdot k\end{array}$ 就是原方程的一组解。

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

   设 $a(x_1+m)+b(y_1+n)=c$，展开得 $a{x}_1+am+b{y}_1+bn=c$，所以只需要让 $am+bn=0$。

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

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

求解线性同余方程 $ax\equiv b(\mathrm{mod}n).$

以下内容摘自OI_wiki。

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

$$ax+nk=b$$

（注：$ax÷n=-k\dots b$）

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

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

$$a{\displaystyle \frac{b}{gcd(a,n)}}{x}_0+n{\displaystyle \frac{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=(xmodt+t)modt$$

其中有

$$t={\displaystyle \frac{n}{gcd(a,n)}}$$

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

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

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

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

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

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

若限制 $x>0$，则 $x_0+t{d}_x>0$，$t>-\frac{x_0}{d_x}$。

若限制 $y>0$，则 $y_0-t{d}_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$。

• 对不等式左侧解释，因为 $>\to \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{array}{l}x\equiv a_1(\mathrm{mod}{n}_1)\\ x\equiv a_2(\mathrm{mod}{n}_2)\\ ⋮\\ x\equiv a_k(\mathrm{mod}{n}_k)\end{array}$$

求解过程：

1. 计算所有模数的积 $n$。

2. 对于第 $i$ 个方程：

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

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

   c. 计算 $c_i=m_i{m}_i^{-1}$（不要对 $n_i$ 取模）。

3. 方程组在模 $n$ 意义下的唯一解为：$x=\sum _{i=1}^k{a}_i{c}_i(\mathrm{mod}n)$。

证明：

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

枚举 $j\in [1,k]$，当 $i\ne j$ 时，因为 $n_i\mid \frac{n}{n_j}$，所以有 $m_j\equiv 0(\mathrm{mod}{n}_i)，\mathrm{故}{c}_j\equiv m_j\equiv 0(\mathrm{mod}{n}_i)$。又有 $c_i\equiv m_i\cdot (m_i^{-1}mod{n}_i)\equiv 1(\mathrm{mod}{n}_i)$，所以我们有：

$$\begin{array}{rlr}x& \equiv \sum _{j=1}^k{a}_j{c}_j& (\mathrm{mod}{n}_i)\\ & \equiv a_i{c}_i& (\mathrm{mod}{n}_i)\\ & \equiv a_i\cdot m_i\cdot (m_i^{-1}mod{n}_i)& (\mathrm{mod}{n}_i)\\ & \equiv a_i& (\mathrm{mod}{n}_i)\end{array}$$

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

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

另外，若 $x\ne 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{array}{l}x\equiv a_1(\mathrm{mod}{m}_1)\\ x\equiv a_2(\mathrm{mod}{m}_2)\\ ⋮\\ x\equiv a_n(\mathrm{mod}{m}_n)\end{array}$$

先考虑只有两个方程的情况。设两个方程分别是 $x\equiv a_1(\mathrm{mod}{m}_1)、x\equiv a_2(\mathrm{mod}{m}_2)$。

将它们转化为不定方程：$x=m_{1}p+a_1=m_{2}q+a_2$，其中 $p,q$ 是整数，则有 $m_{1}p-m_{2}q=a_2-a_1$。

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

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

则原来的两方程组成的模方程组的解为 $x\equiv b(\mathrm{mod}M)$，其中 $b=m_{1}p+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\le a,b,c，1\le x,y,z\le {10}^5，h\le 2^{63}-1）$

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

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

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

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

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

$$i\stackrel{y}{\to }(i+y)modx$$

$$i\stackrel{z}{\to }(i+z)modx$$

上述两种边分别代表通过向上跳 $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;
}