【学习笔记 / 长期更新】OI 中的数论

· · 算法·理论

目录

- Preface

0.1 前言

本文意为作者从 0 开始学习数论，同时也对 OI Wiki 的某些内容做补充说明。

如果你看到有一些小标题没有内容，很正常，作者\color{white}\small\textbf{不}会填坑的。

同时本文用了大量 \LaTeX，加载公式较慢属正常现象。

都看到这了不点个赞吗 /kel

0.2 常用符号

  1. 整除符号：x \mid y，表示 x 整除 y。

  2. 取模符号：x \bmod{y}，表示 x 除以 y 的余数。

  3. 互质符号：x \bot y，表示 x 与 y 互质。

  4. 最大公约数：\gcd(x,y) 或 (x,y)。

  5. 最小公倍数：\operatorname{lcm}(x,y) 或 [x,y]。

  6. 阶乘符号：n!，表示 1 \times 2 \times 3 \times \cdots \times n，特别的，0!=1。

  7. 求和符号：\sum，表示特定条件下几个数的和，例如：

  8. 求积符号：\prod，表示特定条件下几个数的积，例如：

  9. 向上取整符号：\lceil x \rceil，表示大于等于 x 的最大的整数。

  10. 向下取整符号：\lfloor x \rfloor，表示小于等于 x 的最大的整数。

  11. 组合数：\binom{x}{y}。

  12. 整数集：\mathbf{Z}，表示所有的整数。

  13. 自然数集：\mathbf{N}，表示所有的自然数。

  14. 实数集：\mathbf{R}，表示所有的实数。

  15. 存在：\exists x:P(x)，表示至少存在一个 x 使得 P(x) 的值为真。

0.3 快速幂 \Theta(\log n)

从熟悉的快速幂开始学起。

题意：给定 a,b，求 a^b。

朴素算法是 \Theta(b) 的，较慢。

设 b 的二进制为 (n_tn_{t-1}\cdots n_2n_1)_2，考虑把 b 二进制分解，即 n_t2^{p}+n_{t-1}2^{p-1}+ \cdots + n_12^1+n_02^0，其中 n \in\{0,1\}，p 为非负整数，易得 p= \lfloor \log_2n\rfloor。

又由幂的性质可得，a^{x+y} = a^x \times a^y。

于是 a^b 就可写为如下形式：

\begin{aligned} a^b &= a^{n_t2^{p}+n_{t-1}2^{p-1}+ \cdots + n_12^1+n_02^0}\\ &= a^{n_t2^{p}} \times a^{n_{t-1}2^{p-1}} \times \dots a^{n_02^0} \end{aligned}

又有 x^{2y}=x^y \cdot x^y=(x^y)^2，我们就可以在 \log_2 n 的时间算完上式了。

\begin{aligned} 13&=(1101)_2\\ &=1 \times 2^3+ 1 \times2^2 + 0 \times 2^1 + 1 \times 2^0\\ \end{aligned} \downarrow \begin{aligned} 5^{13}&=5^{2^3}\times 5^{2^2} \times 5^{2^0}\\ &=5^8 \times 5^4 \times 5^1\\ &=390625 \times 625 \times 5\\ &=1220703625 \end{aligned}

由于取模并不影响乘法的运算，所以这种方法也可以求 a^b \bmod k。

Sample Code

// 0.3 qpow
template<class T = long long> T qpow(T a, T b,const T mod = b - 2){
    T ans = 1, base = a;
    while(b){
        if(b & 1) ans = ans * base % mod;
        base = base * base % mod;
        b >>= 1;
    }
    return ans;
}

- 数论

1 数论基础

1.1 整除

设 a,b \in \mathbf{Z},a \neq 0，且如果 \exists q \in \mathbf{Z}，使得 b = aq，那么我们就说 a 整除 b，或者 b 能被 a 整除，记作 a \mid b，其中 b 是 a 的倍数，a 是 b 的因数。

否则我们就称 b 不能被 a 整除，记作 a \nmid b。

性质：

以上性质的证明。

对于整数 b \ne 0，b 的约束只有有限个，0 是任何不等于 0 的整数的倍数。

对于整数 b \ne 0,1，b 有四个显然因数为 \pm 1,\pm b，特别的，1 只有两个显然因数。

对于整数 b \ne 0 的真因数称为 b 非显然因数的因数。

1.2 带余除法

设 a,b \in \mathbf{Z}，a \neq 0，则有 b=qa+r，其中 q=\left\lfloor \dfrac{b}{a} \right\rfloor,0 \le r \le |a|。

这里的 r 被称为最小非负余数，也可以写成 r = b\bmod{a}。

其中 r=b-aq=b-a\left\lfloor\dfrac{b}{a}\right\rfloor。

性质

1.3 同余

1.4 数论函数

1.4.1 积性函数

指的是如果 \forall x,y\in\mathbf{N}^*,\gcd(x,y)=1 且 f(1)=1，那么 f(x\times y)=f(x) \times f(y)。

特殊的，如果 \forall x,y\in\mathbf{N}^* 且 f(1)=1，那么 f(x\times y)=f(x) \times f(y)，则称 f(n) 为完全积性函数。

2 素数

定义：若整数 b \ne 0, \pm1 只有四个显然因数的时候，b 是素数，否则 b 是合数。

### 2.1 唯一分解定理 - 合数总能分解成几个素数的乘积。 形式化的，设 $a$ 为合数，$p$ 为素数集，则 $$a=p_1^{\alpha_1}p_2^{\alpha_2}\cdots p_n^{\alpha_n}$$ - 例如 $$ \begin{aligned} 24&=2 \times 2 \times 2 \times 3\\ &=2^3 \times 3^1 \end{aligned} $$ - 应用：<https://atcoder.jp/contests/abc280/tasks/abc280_d> ### 2.2 判定素数 $\Theta(\sqrt{n})

通过素数的性质可知，对于素数 n，\forall k\in \mathbb{Z},k \in[2,n)，有 k \nmid n，于是不难想到一个 \Theta(n) 的朴素程序，用 2 \sim n-1 的所有正整数对 n 进行试除，如果 n \bmod{k}=0，说明 p 除了显然因数外，还有 k 和 \frac{n}{k} 的因数，因此判定为合数。

Sample Code

// 1.2.2 prime check O(n)
bool check_linear(int x){
    for(int i = 2; i < x; i++){
        if(x % i == 0) return false;
    }
    return true;
}

显然，当 n 到 10^{9} 级别时，此种方法便会超时。

注意到 1.1 的一句话：设 b> 0，那么当 d 遍历 b 的正因数的时候，\frac{b}{d} 也在遍历 b 的正因数。

也就是说，当 k 访问到 \sqrt{n} 的时候，如果此时还存在 n \bmod{k}=0，那么它就会在访问到 \frac{n}{k} 的时候被判定为合数，因此，\sqrt{n} 之后的判定完全是多余的，因此时间复杂度便可降到 \Theta(\sqrt{n}) 级别。

Sample Code

// 1.2.2 prime check O(sqrt(n))
bool check_sqrt(int x){
    for(int i = 2; i * i <= x; i++){ //注意这里是 <=x，不然可能把一些完全平方数放过去
        if(x % i == 0) return false;
    }
    return true;
}

2.3 Miller-Rabin 判定素数 \Theta(k \log^3n)

有时候，要判定的数的级别是 10^{18} 级别的，这时候 \Theta(\sqrt{n}) 的 check_sqrt 就会 TLE，这时我们可以对其进行 Miller-Rabin 素数测试，来判断其是否为素数。

2.3.1 费马小定理（Fermat's Little Theorem）

有质数 p，则对于正整数 a \nmid p，有 a^{p-1} \equiv 1\pmod{p}。

考虑它的逆否命题，我们就有了一种判断素数的方法：

如果 a^{p-1} \not \equiv 1 \pmod{p}，则 p 不是素数。

2.3.2 二次探测定理

对于上面的 \text{Carmichael} 数，FLT 不能判断其素性。此时就可以用到二次探测定理。

如果 p 是一个素数，0 < x < p，那么 x^2 \equiv 1 \pmod{p} 的解只有两个：x_1=1,x_2=p-1。

证明：

\begin{aligned} x^2 \equiv 1 \pmod{p}\\ x^2-1 \equiv 0 \pmod{p}\\ (x+1)(x-1) \equiv 0\pmod{p} \end{aligned}

于是有 (x+1)(x-1) \mid p，故得出 1 \sim p-1 的正整数只有 1 和 p-1 两个解，故原命题得证。

2.3.3 Miller-Rabin

于是我们有了一种判定素数的方法：Miller-Rabin。

具体步骤：

  1. 如果 n=1 或 n \bmod{2}=0，直接返回 false;
  2. 把 a^{n-1} 化成 (a^{\frac{n-1}{2}})^2,(a^{\frac{n-1}{4}})^{2^2} 的形式，直到不能化为止，设最后的底数为 a^u。
  3. 选择素数 a，计算 a^u \bmod n。
  4. 不断运用二次探测定理，判断当 x^2 \equiv 1 \pmod{p} 成立时，是否 x=1 或 x=n-1，如果不是即探测失败，返回 false。
  5. 如果满足条件，就判断 a^{n-1} \equiv1\pmod{p} 是否成立，如成立，则这一轮 Miller-Rabin 检测结束。

一般来说素数 a 的选择在 long long 范围内选择 2,3,5,7,11,13,17,19,23 这九个数即可 100\% 正确。

时间复杂度：\Theta(k \log^3n)，k 即为选择了几个素数进行检测。

Sample Code

// 1.2.3.3 Miller-Rabin

typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;

ll mul(ll a,ll b,ll mod){
// 快速乘，原理是用 ull 的溢出解决 ll 的溢出
    ll c=(ld)a/mod*b;
    ll res=(ull)a*b-(ull)c*mod;
    return (res+mod)%mod;
}

ll qpow(ll a,ll b,const ll p){
// 快速幂
    ll ans=1ll,base=a;
    while(b){
        if(b&1) ans=mul(ans,base,p);
        base=mul(base,base,p);
        b>>=1;
    }
    return ans;
}

ll ud[]={2,3,5,7,11,13,17,19,23}; // 素数集

bool MRtest(ll n){
    if(n%2==0 || n<3) return n==2; // 特判
    ll u=n-1,t=0;
    while(u%2==0) u>>=1,++t; // 分解为 a^u
    ll x,pre;
    for(auto p:ud){ // 选择素数
        if(p==n) return true; // 特判
        x=qpow(p,u,n);
        pre=x;
        for(int i=1;i<=t;i++){
            x=mul(x,x,n);
            if(x==1 && pre!=1 && pre!=n-1) return false;
            pre=x; // 二次探测
        }
        if(x!=1) return false; // 费马小定理
    }
    return true;
}

2.4 素数筛

通过 1.2.1 的素数判定，我们知道了如何判定一个数是否是素数，但是，如果要判断 1 \sim n 中所有的数是否是素数呢？

一个显然的做法是扫一遍区间，用根号做法或者 Miller-Rabin 判断每一个数是否是素数，但是 \Theta(n \sqrt{n}) 或 \Theta(n k\log^3 n) 的复杂度显然太慢，于是我们引进了素数筛的思想。

2.4.1 埃氏筛 \Theta(n \log \log n)

由唯一分解定理（2.1），每个合数都可以被分解成若干个质数的乘积，换句话说，每个合数都是它的质因数的倍数。

于是我们可以得到这样一种筛法（埃拉托斯特尼筛法，简称埃氏筛）：建立一个标记数组 vis，每次从小到大选择未被标记的数，把它的倍数都标记为合数，这样在结束后未被标记的数就是素数。

时间复杂度为 \Theta(n \log \log n)，近乎线性。

Sample Code

// 1.2.4.1 Eratosthenes sieve
bool vis[1000005]; //标记此数是不是素数
vector<int> prime; //素数集
void Eratosthenes_sieve(int n){
    vis[1] = 1;
    for(int i = 2; i <= n; i++){
        if(vis[i] == 0){
            prime.push_back(i);
            for(int j = i + i; j <= n; j += i){
                vis[j] = 1;
            }
        }
    }
}

2.4.2 欧拉筛 / 线性筛 \Theta(n)

Eratosthenes sieve 已经足够快了，但是碰到一些 n=10^7 的筛就显然会 TLE。继续改进算法。

我们注意到，在埃氏筛中，42 被筛了三遍（分别为 2,3,7），比较耗费时间，事实上，如果一个数只被被其最小的素因子筛一次，那么复杂度就会降到 \Theta(n)，此种方法被称作欧拉筛或线性筛，见代码：

Sample Code

// 1.2.4.2 Euler sieve
bool vis[1000005];
vector<int> prime;
void Euler_sieve(int n){
    vis[1] = 1;
    for(int i = 2; i <= n; i++){
        if(!vis[i]) prime.push_back(i);
        for(int j = 0; (size_t)j < prime.size(); j++){
            int it = prime.at(j);
            if(1ll * it * i > n) break;
            vis[1ll * it * i] = 1;
            if(i % it == 0){
                // i % pri[j] == 0
                // 换言之，i 之前被 pri[j] 筛过了
                // 由于 pri 里面质数是从小到大的，所以 i 乘上其他的质数的结果一定也是
                // pri[j] 的倍数 它们都被筛过了，就不需要再筛了，所以这里直接 break
                // 掉就好了
                // From OI-Wiki
                break;
            }
        }
    }
}

2.4.3 筛法的一些优化

2.5 分解质因数

给定一个数 n \in \mathbf{N}^+，分解出它的质因数。

2.5.1 朴素算法 \Theta(\sqrt{n})

根据唯一分解定理，一个合数总能被分解成几个素数的乘积，于是我们就有了一个类似根号时间内判定素数的程序，其时间复杂度也为 \Theta(\sqrt n)。

Sample Code

vector<pair<int,int> > factor(ll x){ // 返回一个 pair<int,int> 类型的 vector
    vector<pair<int,int> > v; // first 存储素因子，second 存储素因子的次数
    for(ll i=2;i*i<=x;i++){
        if(x%i==0){
            int t=0;
            while(x%i==0) t++,x/=i;
            v.push_back(make_pair(i,t));
        }
    }
    if(x>1) v.push_back(make_pair(x,1)); // 如果 x 本身就是个素数
    return v;
}

2.5.2 素数筛优化 \Theta(\sqrt{\frac{n}{\ln n}})

我们知道，n 之内的素数个数大约为 \dfrac{n}{\ln n} 个，我们如果事先把这些素数筛出来再加上上面朴素算法就能做到以上的复杂度。

代码和上文几乎一样，只是把 i 放到素数集里面循环，代码就不放了。

2.5.3 Pollard-Rho \Theta(n^{\frac{1}{4}})

2.6 欧拉函数 \varphi(n)

2.6.1 性质

首先有欧拉函数是积性函数，也就是说，如果 \gcd(a,b)=1，有 \varphi(a \times b)=\varphi(a) \times \varphi(b)。

由素数定义可知，如果 p 为素数，\varphi(p)=p-1。

2.6.2 欧拉函数（Euler's totient function）

凭借 2.6.1 和 2.1 唯一分解定理 于是我们就有了欧拉函数的计算方法，设合数 n=\prod^s_{i=1} p_i^{k_i}，则有

\varphi(n)=n \times\prod_{i=1}^{s} \dfrac{p_i-1}{p_i}

我们来分析一下这个公式是怎么来的。

2.6.3 引理

设 p 为素数集，则 \varphi(p^k)=p^k-p^{k-1}。

这很好证明，p^k 个数中，只有 p^{k-1} 个是 p 的倍数，其他数都不含 p 这个因子，故 \varphi(p^k)=p^k-p^{k-1}。

于是，根据欧拉函数的积性和唯一分解定理可知

\begin{aligned} \varphi(n)&=\prod^s_{i=1}\varphi(p_i^{k_i})\\ &=\prod^s_{i=1} p_i^{k_i}-p_i^{k_i-1}\\ &=\prod^s_{i=1} p_i^{k_i-1}\cdot (p_i-1)\\ &=\prod^s_{i=1} p_i^{k_i}\cdot \dfrac{p_i-1}{p_i}\\ &=n \times \prod^s_{i=1} \dfrac{p_i-1}{p_i} \end{aligned}

2.6.4 求一个数的欧拉函数

2.6.4.1 朴素求法 \Theta(\sqrt n)

2.6.4.2 Pollard-Rho \Theta(n ^{\frac{1}{4}})

2.6.5 筛法求欧拉函数 \Theta(n)

我们注意到，欧拉筛每次都会以它最小的素因子筛到这个数。设 n=pq，其中 p 为它最小的素因子，分类讨论：

\begin{aligned} \varphi(n)&=n \times\prod_{i=1}^{s} \dfrac{p_i-1}{p_i}\\ &=p\times q \times\prod_{i=1}^{s} \dfrac{p_i-1}{p_i}\\ &= p \times \varphi(q) \end{aligned}
\varphi(n)=\varphi(p)\times \varphi(q)

Sample Code

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

const int MAXN = 1000000 + 5;

int pri[MAXN],phi[MAXN],cnt;
bool is_prime[MAXN];

void get_euler(const int N){
    for(int i=1;i<=N;i++) is_prime[i]=1;
    is_prime[1]=0;
    phi[1]=1;
    cnt=0;
    for(int i=1;i<=N;i++){
        if(is_prime[i]){
            pri[++cnt]=i;
            phi[i]=i-1; //i为素数，phi[i]=i-1
        }
        for(int j=1;j<=cnt;j++){
            if(i*pri[j]>N) break;
            is_prime[i*pri[j]]=0;
            if(i%pri[j]) // 注意分辨 phi 和 pri
                phi[i*pri[j]]=phi[i]*phi[pri[j]]; // 情况1
            else{
                phi[i*pri[j]]=phi[i]*pri[j]; // 情况2
                break;
            }
        }
    }
}

int n,T;

int main(){
    get_euler(MAXN-5);
    cin>>T;
    while(T--){
        cin>>n;
        cout<<phi[n]<<'\n';
    }
}

3 最大公约数

3.1 最大公约数 \gcd(a,b)

定义正整数 a,b 的最大公约数为他们共同的约数中最大的一个（Greatest Common Divisor），简写为 \gcd。

对于如何快速地求出 \gcd(a,b)，有欧几里得算法。

下面给出这个算法的证明过程（个人觉得 OI-Wiki 上的有些不严谨）：

Proof

接下来我们来推出 \gcd(a,b) \mid \gcd(b,a \bmod{b})：

设 d=\gcd(a,b),r=a \bmod{b}。

由 等式 3 得 r=a-b\left\lfloor\dfrac{a}{b}\right\rfloor。

于是根据 等式 2 ，由 d \mid a \land d \mid b 且 r=a-b\left\lfloor\dfrac{a}{b}\right\rfloor，推出 d \mid r，即 d \mid a \bmod{b}。

由 推论 1 ，因为 d \mid b \land d \mid a \bmod{b}，于是 d \mid \gcd(a,a \bmod{b})。

即 \gcd(a,b) \mid \gcd(b,a \bmod{b})。

接下来我们来推出 \gcd(b,a \bmod{b}) \mid \gcd(a,b)：

设 c=\gcd(b,a \bmod{b})，则有 c\mid b \land c \mid(a \bmod{b})。

由带余除法，a=bq+r，其中 q=\left\lfloor\dfrac{a}{b}\right\rfloor,r=a \bmod b。

于是根据 等式 2 ，由 c \mid b \land c \mid (a\bmod{b}) 且 a=b \left\lfloor\dfrac{a}{b}\right\rfloor+a \bmod{b}，得出 c \mid a。

于是有 c \mid a \land c \mid b 又因为 推论 1 ，得出 c \mid \gcd(a,b)。

即 \gcd(b,a \bmod{b}) \mid \gcd(a,b)。

于是根据 等式 1 ，\gcd(b,a \bmod{b}) \mid \gcd(a,b) \land \gcd(a,b) \mid \gcd(b,a \bmod{b}) \implies \gcd(a,b)=\gcd(b,a\bmod{b}) \square。

至此，欧几里得定理（或辗转相除法、\gcd 递归定理） 证明完毕。

因为每次 a \bmod{b} 的值至多变为 a 的一半，所以该算法的时间复杂度为 \Theta(\log n)。

常见的几种写法：

// 以下默认 typedef long long ll 
ll gcd(ll a,ll b){ // 非递归版
    int tmp;
    while(b!=0){
        tmp=b;
        b=a%b;
        a=tmp;
    }
    return a;
}

ll gcd(ll a,ll b){ // 递归版
    if(b==0) return a;
    return gcd(b,a%b);
}

ll gcd(ll a,ll b){ // 递归版简化
    return b==0?a:gcd(b,a%b);
}

ll gcd(ll a,ll b){ // 位运算版
    while(b){
        a%=b;
        b^=a;
        a^=b;
        b^=a;
    }
    return a;
}

ll gcd(ll a,ll b){ // 位运算简化
    while(b^=a^=b^=a%=b);
    return a;
}

3.2 最小公倍数 \operatorname{lcm}(a,b)

定义正整数 a,b 的最小公倍数为他们共同的倍数中最小的一个（Least Common Multiple），简写为 \operatorname{lcm}。

我们接下来来推一推 \operatorname{lcm} 要如何快速地求。

由 唯一分解定理，可得

a=p_1^{ka_1}p_2^{ka_2}\cdots p_s^{ka_s},b=p_1^{kb_1}p_2^{kb_2}\cdots p_s^{kb_s}

因为 \gcd 为二者的最大公因数，所以有

\gcd(a,b)=p_1^{\min(ka_1,kb_1)}p_2^{\min(ka_2,kb_2)}\cdots p_s^{\min(ka_s,kb_s)}

同理，\operatorname{lcm} 即为

\operatorname{lcm}(a,b)=p_1^{\max(ka_1,kb_1)}p_2^{\max(ka_2,kb_2)}\cdots p_s^{\max(ka_s,kb_s)}

因为 ka_1+kb_1=\min(ka_1+kb_1)+\max(ka_2+kb_2)，所以 \gcd(a,b) \times \operatorname{lcm}(a,b)=a \times b。

也就是说，只要求出了两个数的 \gcd，两个数的 \operatorname{lcm} 就能 \Theta(1) 求出来，容易看出，求解 \operatorname{lcm} 的时间复杂度也为 \Theta(\log n)。

ll lcm(ll a,ll b){
    return a/gcd(a,b)*b; // 防止 a*b 爆 ll，所以先 /gcd 再 *b
}

3.3 扩展欧几里得算法 \text{exgcd}

扩展欧几里得算法（Extended Euclidean algorithm），简写为 \text{exgcd}，是用来求解形如 ax+by=\gcd(a,b) 的一组可行解。

接下来我们来分析一下这个算法：

首先，当 b=0 时，显然有 \gcd(a,0)=a，此时 \begin{cases}x=1\\y=0\end{cases} 是方程的一组可行解。

设 ax_1+by_1=\gcd(a,b)，bx_2+(a \bmod b)y_2=\gcd(b,a \bmod{b})。

于是根据 欧几里得定理，显然有 ax_1+by_1=bx_2+(a \bmod{b})y_2

于是

\begin{aligned} ax_1+by_1&=bx_2+(a\bmod b)y_2\\ ax_1+by_1&=bx_2+\left( a-b \left\lfloor \dfrac{a}{b}\right\rfloor\right) y_2\\ ax_1+by_1&=bx_2+ay_2-b\left\lfloor \dfrac{a}{b}\right\rfloor y_2\\ ax_1+by_1&=ay_2+b\left(x_2-\left\lfloor \dfrac{a}{b}\right\rfloor y_2\right)\\ x_1=y_2&,y_1=x_2-\left\lfloor \dfrac{a}{b}\right\rfloor y_2 \end{aligned}

用 \gcd 递归求解即可。当 b=0 时，带入上文的特解返回。

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

3.3.1 用 \text{exgcd} 求解的一个问题

题意：给定 a,b，求出 ax \equiv 1 \pmod{b} 中 x 的最小正整数解，保证有解。

问题转化：ax \equiv 1 \pmod{b} \implies ax+by=1，其中 y 是我们引进的一个负整数。

明显地，\text{exgcd} 可以用来求解形如 ax+by=\gcd(a,b) 的解，现在的问题是如何把 1 和 \gcd(a,b) 联系起来。

引理：ax+by=m 有整数解的必要条件是 \gcd(a,b)\mid m。

Proof：因为 \gcd(a,b) \mid a \land \gcd(a,b) \mid b，所以 \gcd(a,b) \mid (ax+by)，即 \gcd(a,b) \mid m。

于是再联系到题目上，就有 \gcd(a,b)\mid 1，因为 \gcd(a,b) 显然是个正整数，所以 \gcd(a,b) 只能等于 1 了。

这也就是说，a,b 互质，于是 \gcd(a,b)=1，正好符合 \text{exgcd} 的原始形式，所以 \text{exgcd} 就可以在 \Theta(\log n) 的时间来求解。

最小正整数解

因为 ax+by=1，所以 ax+by+kab-kab=1，所以 a(x\pm kb)+b(y \mp ka)=1。

这也就是说，如果 x>0，那么就把 x 一直减 b，直到不能减为止，反之亦然。

这在 C++ 的运算中，用 x = (x % b + b) % b 即可解决，后面的 + b) % b 是为了让 x 从负数变为正数。

Sample Code

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

ll a,b,x,y;

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

int main(){
    cin>>a>>b;
    exgcd(a,b,x,y);
    x=(x%b+b)%b;
    cout<<x<<'\n';
    return 0;
}

3.3.2 有理数取余：用 \text{exgcd} 求解的另一个问题

这个值被定义为 bx \equiv a \pmod p 的解。

我们来分析一下。

因为同余有性质： 如果 b \equiv a \pmod p，那么两边乘上一个数之后，还是不变： b \times d \equiv a \times d \pmod p，带进原式，就变成了 bx \equiv a \pmod p。

慢着，我们先看上一个问题，他求解的是 bx_1 \equiv 1 \pmod{p}，然后两边同乘 a，变成了 b (ax_1) \equiv a \pmod p。

于是有 x=ax_1，问题解决。

等等，那无解情况呢？我们分情况讨论一下：

  1. - $a \bmod p = 0$，即 $a$ 也为 $p$ 的倍数，于是原同余式有无数解。 - $a \bmod p \ne 0$，即 $a$ 不为 $p$ 的倍数，于是原同余式无解（因为等式左边 $\mod p$ 永远为 $0$）。
  2. - 即 $b,p$ 互质，$\gcd(b,p)=1$，于是根据 **3.3.1** 中的一个结论，原方程必定有解。

总结：如果 b \bmod p=0，则原方程无解。

Sample Code

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

inline ll read(){
    char ch=getchar(); ll x(0),f(0);
    for(; !isdigit(ch); ch=getchar()) f|=(ch=='-');
    for(;  isdigit(ch); ch=getchar()) x=((x<<1)+(x<<3)+(ch^48))%p;
    return f?-x:x;
}

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

ll a,b,x,y;

int main(){
    a=read(); b=read();
    if(b==0)
        return puts("Angry!")&0;
    exgcd(b,p,x,y);
    ll ans=(a*x%p+p)%p;
    printf("%lld\n",ans);
}

3.3.3 \text{exgcd} 例题

4 数论分块

4.1 算法

给定 n\in \mathbf N^{+}, 求 \sum^n_{i=1} \left\lfloor\dfrac{n}{i}\right\rfloor 的值。

朴素算法是 \Theta(n) 的，考虑优化。

有结论：对于 \left\lfloor\dfrac n i\right\rfloor =\left\lfloor\dfrac n j\right\rfloor 且 1 \le i \le j \le n，有 j 的最大值为 \left\lfloor\dfrac n i\right\rfloor。

证明：设 k = \left\lfloor\dfrac n i\right\rfloor，显然有 k \le \dfrac n i。

\therefore \left\lfloor\dfrac n k\right\rfloor \ge \left\lfloor\frac n {\frac n i}\right\rfloor = \left\lfloor i \right\rfloor=i \therefore j=\max_{\left\lfloor\frac n i\right\rfloor =\left\lfloor\frac n j\right\rfloor} i=\left\lfloor\dfrac n k\right\rfloor=\left\lfloor\frac n {\frac n i}\right\rfloor \qquad \square

容易证得，符合要求的区间约有 \sqrt n 种。

时间复杂度 \Theta(\sqrt n) 。

Sample Code

ll H(int n){
    ll res = 0;
    int l = 1, r;
    while(l <= n){
        r = n / (n / l);
        res += 1ll * (r - l + 1) * (n / l);
        l = r + 1;
    }
    return res;
}

4.2 例题

4.3 用数论分块解决的一个问题

以 P2261 [CQOI2007]余数求和 为例。

给出 n,k，求 \sum^n_{i=1} k \bmod i。

题目转化：由带余除法一章，我们知道 k \bmod i = k - i \times \left \lfloor \dfrac k i \right \rfloor。

于是问题就变成了 \sum^n_{i=1} k - i \times \left \lfloor \dfrac k i \right \rfloor。

\begin{aligned} G(n,k)&=\sum^n_{i=1} k - i \times \left \lfloor \dfrac k i \right \rfloor\\ &=n\times k-\sum^n_{i=1} i \times \left \lfloor \dfrac k i \right \rfloor \end{aligned}

其中 \sum^n_{i=1} i \times \left \lfloor \dfrac k i \right \rfloor 部分可以数论分块。因为 \left \lfloor \dfrac k i \right \rfloor 相等，由乘法分配律可得 \sum^r_{l=1} l \times \left \lfloor \dfrac k l \right \rfloor=\left \lfloor \dfrac k l \right \rfloor \times \sum^r_{l=1} l。后面的部分可以数论分块顺带解决掉。于是就做完了。

注意这题和上面那题不同。这题 n,k 是分开的，所以当 \left \lfloor \dfrac k l\right \rfloor=0 时要特判，不然会 RE。

以及 r 要和 n 取 \min。

复杂度：\Theta(\sqrt n)。

ll H(ll n, ll k){
    ll ans = k * n;
    int l = 1, r;
    while(l <= n){
        if(k / l != 0) r = k / (k / l);
        else r = n;
        chkmin(r, n);
        ans -= 1ll * (l + r) * (r - l + 1) / 2 * (k / l);
        l = r + 1;
    }
    return ans;
}

- Reference