数论

数学能不能去死啊！！！

本文缺少一些比较困难证明，且对于部分杂乱的定理缺少记录，主要原因是实力不够。

By Pan_g

筛法

筛到你的时候，你应该感到幸运。

埃拉托斯特尼筛法

简称埃氏筛。

思路

思路就是把每一个素数的倍数标记一遍，非常简单。

优化

因为遍历到 $x$ 时， $2 \sim x - 1$ 都筛过了，所以只用从 $x^2$ 开始标记，这时的复杂度是 $O(N \log\log N)$ （有的题会用上）。

剩下的就是一些常数优化，比如筛到平方根、只筛奇数啥的。

不过有一种分块优化，不过应该不重要。

欧拉筛

也称线性筛，由其复杂度为线性得名。

思路

埃氏筛的每一个合数会标记很多遍，所以考虑如何只让其标记一遍。

考虑反过来想，把每个数的素数倍标记。

接着再思考，当遍历到一个当前数 $x$ 的质因子 $d$ 时，那么接着往后标记就没有意义了。

为什么呢？因为再往后的倍数一定已经被 $d$ 及其他的倍数（包括后面的）标记掉了，所以一定就没有意义了，这个感性理解一下，反正想想也就明白了，不会细说。

std::vector<int> sieve(int n){
    std::vector<int> np(n, false), prime;
    for(int i = 2;i < n;i++){
        if(!np[i]) prime.emplace_back(i);
        for(auto p : prime){
            if(1ll * i * p >= n) break;
            np[i * p] = true;
            // 切记先标记再判断倍数
            if(i % p == 0) break;
        }
    }
    return prime;
} // Eular Sieve

应用

筛法除了筛素数之外还有求积性函数以及类似性质函数的功能。

这个就是按照性质筛就好了，推式子也就分为三类：

• $x$ 为质数的答案；
• $p | x$ 且 $p$ 是质数的转移；
• $p \bot x$ 且 $p$ 是质数的转移；

约数类

两人相约之时，一定有着人格上的相等。

最大公约数 & 最小公倍数

最大公约数一般记为 $\gcd$ ，最小公倍数一般记为 $\operatorname{lcm}$ 。

辗转相除法

都知道辗转相除法，简单证明一下：

令 $\gcd(x, y) = k, x = ak, y = bk, a = pb + q(p, q\in \mathbb{Z}, q < b)$

所以 $x = pbk + qk, y = bk, x \bmod y = qk$

由于 $q < b$ ，所以 $q\bot b$ ，那么就有 $\gcd(x,y) = \gcd(y, x \bmod y)$ 。

实现比较简单。

都知道 $\operatorname{lcm}(x, y) = \dfrac{xy}{\gcd(x, y)}$ 不写了。

更相减损术

$\gcd(a, b) = \gcd(b, a - b)$ 证明比较简单。

有一个优化就是当 $2 | a$ 且 $2 | b$ 时，有 $\gcd(a, b) = 2\gcd(\dfrac a2, \dfrac b2)$ ，然后可以大幅提升效率。

主要适用于大数的相关计算，也就高精度上有用。

扩展欧几里得算法

辗转相除法也叫欧几里得算法，而扩展欧几里得也就叫 exgcd 。

原理

用于求解一个不定方程 $ax + by = \gcd(a, b)$ 的可行解，根据裴蜀定理，该方程一定有解。

注意到 $\gcd(a, b) = \gcd(b, a \bmod b)$ ，所以将这个式子拓展出去得到：

$$ \begin{cases} ax_1 + by_1 = \gcd(a, b) \\ bx_2 + (a \bmod b)y_2 = \gcd(b, a \bmod b) \\ \gcd(a, b) = \gcd(b, a \bmod b) \\ \end{cases} $$

$$ \begin{aligned} ax_1 + by_1 &= bx_2 + (a \bmod b)y_2 \\ ax_1 + by_1 &= bx_2 + (a - b \lfloor \dfrac a b \rfloor)y_2 \\ x_1a + y_1b &= y_2 a + (x_2 - \lfloor \dfrac a b \rfloor y_2)b \end{aligned} $$

$$ \begin{cases} x_1 = y_2 \\ y_1 = x_2 - \lfloor \dfrac a b \rfloor y_2 \end{cases} $$

所以可以在辗转相除的过程中求出一个解。

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

通解为，下面 $t\in \mathbb{Z}$ ：

$$ \begin{cases} x = x_0 + bt \\ y = y_0 - at \end{cases} $$

其他分析

复杂度同欧几里得算法为 $O(\log \max\{a, b\})$ ，且运算过程中 $|x| \leq |b|, |y| \leq |a|$ ，妈妈再也不用担心我爆 long long 了。

欧拉函数

我们融合吧，发现并不是超融合，只是重新成为了我们本身。

定义

表示小于等于 $n$ 的所有与 $n$ 互质的数的个数，一般记作 $\varphi(n)$ ，所以 $\varphi(1) = 1$ ， $\varphi(p) = p - 1$ （ $p$ 是素数）。

重要结论

积性函数

欧拉函数是积性函数，满足当 $a \bot b$ 时， $\varphi(ab) = \varphi(a)\varphi(b)$ ，这个我不会证。

决定了可以用筛法求。

std::vector<int> sieve(int n){
    std::vector<int> np(n, false), prime, phi(n);
    phi[1] = 1;
    for(int i = 2;i < n;i++){
        if(!np[i]){
            prime.emplace_back(i);
            phi[i] = i - 1;
        }
        for(auto p : prime){
            if(1ll * i * p >= n) break;
            np[i * p] = true;
            if(i % p == 0){
                phi[i * p] = phi[i] * p;
                break;
            }
            phi[i * p] = phi[i] * phi[p];
        }
    }
    return phi;
} // Eular Sieve

约数性质

对于 $n = p^k$ ， $p$ 为质数，则 $\varphi(n) = p^k - p^{k - 1}$ ，这是极其显然的。

更显然的，由唯一分解定理，设 $n = \prod_{i = 1}^{c}p_i^{k_i}$ ，其中 $p_i$ 都是质数，则有 $\varphi(n) = n\prod_{i = 1}^c \dfrac{p_i - 1}{p_i}$ 。

决定了可以用 $O(\sqrt N)$ 的复杂度求，可以用 Pollard-Rho 优化。

int phi(int n){
    int res = n;
    for(int i = 2;i * i <= n;i++){
        if(n % i == 0){
            res = res / i * (i - 1);
            while(n % i == 0) n /= i;
        }
    }
    if(n > 1) res = res / n * (n - 1);
    return res;
}

欧拉反演

要注意欧拉反演仅限在国内这么叫。

内容：

$$ n = \sum_{d | n} \varphi(d) $$

证明是简单的，变为狄利克雷卷积的形式就是 $\operatorname{id} = \varphi \circ 1$ 。

应用实例

$$ \sum_{i = 1} ^ n \sum_{j = 1} ^ n \gcd(i, j) $$

两种求法（这个例子之后还会提到）：

$$ \begin{aligned} &\sum_{i = 1} ^ n \sum_{j = 1} ^ n \gcd(i, j) \\ =& \sum_{k = 1} ^ n k \sum_{i = 1} ^ n \sum_{j = 1} ^ n [\gcd(i, j) = k] \\ =& \sum_{k = 1} ^ n k \sum_{i = 1} ^ {\lfloor \frac nk \rfloor} \sum_{j = 1} ^ {\lfloor \frac nk \rfloor} [\gcd(i, j) = 1] \\ =& \sum_{k = 1} ^ n k \Big(\big(2\sum_{i = 1} ^ {\lfloor \frac nk \rfloor} \varphi(i)\big) - 1\Big) \end{aligned} $$

我们对 $\varphi(x)$ 进行前缀和，处理掉后面的一坨，可以 $O(N)$ 求解。

$$ \begin{aligned} &\sum_{i = 1} ^ n \sum_{j = 1} ^ n \gcd(i, j) \\ =& \sum_{i = 1} ^ n \sum_{j = 1} ^ n \sum_{d | \gcd(i, j)} \varphi(d) \\ =& \sum_{d = 1} ^ n \varphi(d) \sum_{i = 1} ^ {\lfloor \frac nk \rfloor} \sum_{j = 1} ^ {\lfloor \frac nk \rfloor}1 \\ =& \sum_{d = 1} ^ n \varphi(d) \lfloor \frac n{d} \rfloor ^ 2 \end{aligned} $$

也可以 $O(N)$ 求解。

欧拉定理 & 费马小定理

两个毫无关联之人，却在最初的起点相遇。

费马小定理

内容：若 $p$ 是素数，且 $a\bot p$ ，则 $a^{p - 1} \equiv 1 \pmod p$ 。

证明：

构造数列 $\{a_{p - 1}\}$ ，其中 $a_i = i$ ，先证：

$$ \prod_{i = 1}^{p - 1} A_i \equiv \prod_{i = 1}^{p - 1} aA_i \pmod p $$

因为 $a\bot p, A_i\bot p$ ，所以 $aA_i \bot p$ 。

那么可以推知， $aA_i \bmod p$ 两两不同，得证。

令 $x = (p - 1)!$ ，有 $a^{p - 1}x \equiv x \pmod p$ 。

消去 $x$ 得证。

欧拉定理

内容：若 $a\bot p$ ，则 $a^{\varphi(p)} \equiv 1 \pmod p$ 。

证明同上，因为 $p$ 有 $\varphi(p)$ 个与之互质的数，将这些数拎出来成一个序列，同上证明即可（一些简单性质可以反证）。

当 $p$ 为素数时， $\varphi(p) = p - 1$ ，得到费马小定理。

扩展欧拉定理

可以用于光速幂 / 幂塔等问题，内容如下：

$$ a^b \equiv \begin{cases} a^{b \bmod \varphi(p)}, &a \bot p \vee b < \varphi(p) \\ a^{b \bmod \varphi(p) + \varphi(p)}, &a \not \bot p \end{cases} \pmod p $$

线性同余方程

你怎么跑的比我还慢？不，我套了你一圈。

定义

求形如

$$ ax \equiv b \pmod p $$

$[0, p - 1]$ 内唯一解或所有解。

逆元解

将原式变为 $x = a^{-1}b \bmod p$ 。

费马小定理

当 $a \bot p$ 且 $p$ 为质数时，可以用费马小定理： $a^{-1} = a^{p - 2} \bmod p$ 。

然后 $x = a^{p - 2}b \bmod p$ ，同理在 $p$ 不为素数时，也可以用欧拉定理。

线性求逆元

假设已经计算掉了 $1 \sim x - 1$ 的所有 $\mod p$ 意义下的逆元，来求 $x$ 的逆元。

令 $a = \lfloor \dfrac p x\rfloor, b = p \bmod x$ ，则有 $p = ax + b$ 。

在 $\mod p$ 意义下，写同余式，然后两边同乘 $x^{-1}b^{-1}$ 。

$$ ab^{-1} + x^{-1} \equiv 0 \pmod p $$

所以，就有 $x^{-1} = \lfloor \dfrac p x\rfloor (p \bmod x)^{-1} \mod p$ ，同时显然有 $1^{-1} \equiv 1 \pmod p$ 。

以上计算全部基于 $p$ 为素数，因为 $p$ 非素时，会存在 $b^{-1}$ 不存在的情况。

可以做到 $O(N)$ 求 $1 \sim N$ 的逆元。

对于同余方程，解法同上。

扩展欧几里得解

变为 $ax + py = b$ 的形式，那么无解情况与多解情况一目了然。

无解，可以用裴蜀定理；多解，跑完扩欧后直接套公式就好了。

中国剩余定理

有物不知其数，问物几何？什么几何？平面几何怎么比立体几何还难？世界上的巧合真多啊！

定义

最开始由中国人发现这个问题，所以称为中国剩余定理，简称 CRT 。

求线性同余方程组

$$ \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} $$

最小的一个答案。

解法

令 $p = \prod_{i = 1}^n m_i$ 。

对于第 $i$ 个方程：

• 计算 $k_i = \dfrac p{m_i}$ ；
• 计算 $k_i$ 在 $\mod m_i$ 意义下的逆元 ；
• 计算 $c_i = k_i(k_i ^{-1} \mod m_i)$ 。

答案 $x = (\sum_{i = 1}^n a_ic_i) \mod p$ 。

证明：

显然，对于 $i \neq j$ ，$c_j \equiv k_j \equiv 0 \pmod {m_i}$ ； $c_i \equiv k_i(k_i^{-1} \mod m_i) \equiv 1 \pmod {m_i}$ ，所以

$$ \begin{aligned} x &\equiv \sum_{j = 1}^n a_jc_j &\pmod {m_i} \\ &\equiv \sum_{j = 1}^n a_jk_j(k_j^{-1} \mod m_j) & \pmod {m_i} \\ &\equiv a_ik_i(k_i^{-1} \mod m_i) &\pmod {m_i} \\ &\equiv a_i &\pmod {m_i} \end{aligned} $$

int exgcd(int a, int b, int &x, int &y){
    if(b == 0){
        x = 1, y = 0;
        return 0;
    }
    int t = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return t;
}
int CRT(int n, std::vector<int> a, std::vector<int> m){
    int prod = 1, ans = 0;
    for(int i = 0;i < n;i++) prod *= m[i];
    for(int i = 0;i < n;i++){
        int k = prod / m[i], inv, y;
        exgcd(k, m[i], inv, y);
        ans = (ans + a[i] * m % prod * inv % prod + prod) % prod;
    }
    return ans;
}

扩展中国剩余定理

简称 exCRT ，我认为比 CRT 的一般解法好理解。

思路

首先，对于 CRT 的思路，仅限模数 $m_i$ 全部互质的情况下才能使用。所以，就有了 exCRT 。

思路很简单，就是将两个同余方程合并，最后得到一个同余方程。

已知两同余方程：

$$ \begin{cases} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \\ \end{cases} $$

可以转化成为

$$ \begin{aligned} x = m_1p + a_1 &= m_2q + a_2 \\ m_1p - m_2q &= a_2 - a_1 \end{aligned} $$

无解情况一目了然，对于可行解 $p, q$ ，可以将两个方程合并为

$$ x \equiv m_1p + a_1 \pmod {\operatorname{lcm}(m_1, m_2)} $$

两两合并即可，两种求解复杂度相同。

卢卡斯定理

“学 $p$ 进制有什么用啊？”“可以写 Lucas 定理”“ Lucas 定理又不考，学它干嘛？”“可以 Lu ”

定理内容

对于质数 $p$ ，则满足：

$$ \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \lfloor \frac a p \rfloor \\ \lfloor \frac bp \rfloor \end{pmatrix} \begin{pmatrix} a \bmod p \\ b \bmod p \end{pmatrix} $$

其中若 $n < m$ ，则 $\begin{pmatrix} n \\ m \end{pmatrix} = 0$ 。

证明

设 $a = \sum m_ip^i, b = \sum n_ip^i$ ，考虑 $(a + b) ^ p \equiv a ^ p + b ^ p \pmod p$ （用二项式定理可证），借此进行生成函数。

$$ \begin{aligned} &(1 + x) ^ a \\ =& (1 + x) ^ {\sum m_i p ^ i} \\ =& \prod (1 + x ^ {p ^ i})^{m_i} \pmod p \end{aligned} $$

对于 $x ^ b$ 求系数：

$$ \begin{aligned} \begin{pmatrix} a \\ b \end{pmatrix} =& [x ^ b](1 + x) ^ a \\ =& [x ^ b]\prod (1 + x ^ {p ^ i})^{m_i} \\ =& \prod \begin{pmatrix} m_i \\ n_i \end{pmatrix} \end{aligned} $$

稍微变换一下，得证。

实用案例

对于较小的质数 $p$ ，可以预处理比较小的阶乘，然后算组合数。

接着利用 Lucas 定理递归求解组合数。

i64 C(int n, int m){
    if(n < m) return 0LL;
    return fac[n] * inv[fac[m]] % p * inv[fac[n - m]] % p; 
}
i64 Lucas(int n, int m){
    if(n == 0) return 1LL % p;
    return Lucas(n / p, m / p) * C(n % p, m % p) % p;
}

威尔逊定理

除了我，就是一切。

内容

对于每个 $p > 1$ ，且 $p$ 为质数，则有 $(p - 1)! \equiv -1 \pmod p$ ，而且这是充要条件。

证明

不想证，不会，考得不多。

原根

原来不是每个人都有属于自己的那个根。

概念

我们考虑对模意义下的乘法域建立一个群，对于最小的 $k$ 满足 $a^k \equiv 1 \pmod p$ ，则称 $k$ 为 $a$ 在 $\mod p$ 意义下的阶，记为 $\delta _p(a)$ 。

对于原根 $g$ 来讲，满足 $\delta _p(g) = \varphi(p)$ 。

阶的性质

一般来讲，常用的就前两个。

• $a^k(k \in [1, \delta _p(a)]\cap\mathbb{Z}) \bmod p$ 两两不同。
• $\forall a^k \equiv 1 \pmod p$ ，都有 $\delta _p(a) | k$
• 当 $a, b\bot p$ 时，则 $\delta _p(ab) = \delta _p(a)\delta _p(b) \Leftrightarrow \delta _p(a) \bot \delta _p(b)$
• 当 $a\bot p$ 时， $\delta _p(a^k) = \dfrac{\delta _p(a)}{\gcd(\delta _p(a), k)}$

原根性质

原根判定

对于 $p\geq 3, g \bot p$ ，则 $g$ 为原根有充要条件：

对于 $\varphi(p)$ 的所有质因数 $k$，满足 $g ^ {\frac {\varphi(p)}k} \not\equiv 1 \pmod p$ 。

原根个数

对于 $p$ 有原根，则个数为 $\varphi(\varphi(p))$ ，且满足 $g \equiv g_0 ^k \pmod p$ ， $g_0, g$ 是其中的原根， $k$ 是正整数。

原根存在定理

当 $p$ 存在原根，当且仅当 $p = 2, 4, a^k, 2a^k$ ， $a$ 时不为 $2$ 的质数， $k$ 是正整数。

最小原根范围

对于最小的原根，王元证明了质数 $p$ 的最小原根为 $O(p^{0.25 + \varepsilon})$

离散对数

我对你一心一意，但是你却……

定义

取有原根的 $p$ ，设其中一个原根为 $g$ ，若 $a \bot p$ ，可以证明有且一个 $0\leq k < \varphi(p)$ 使得：

$$ g ^ k \equiv a \pmod p $$

其中 $k$ 记为 $\operatorname{ind}_g a$ ，称为离散对数。

性质

其性质与对数相似。

若 $a\bot p, b\bot p$ ， $g, g_1, g_2$ 都是 $p$ 的原根。

• $\operatorname{ind}_g (ab) \equiv \operatorname{ind}_g a + \operatorname{ind}_g b \pmod {\varphi(p)}$
• $\operatorname{ind}_{g_1} a \equiv \operatorname{ind}_{g_2} a \operatorname{ind}_{g_1} g_2 \pmod {\varphi(p)}$
• $a \equiv b \pmod p \iff \operatorname{ind}_g a = \operatorname{ind}_g b$

大步小步算法

简称 BSGS(baby-step giant-step) ，可以用于求解 $a^x \equiv b \pmod p$ 的问题，保证 $a\bot p$ 。

思路

考虑将答案分为 $x = jS - i$ 两部分，显然可以猜到 $S = \sqrt p$ 最优。

那么就变为 $a^{jS - i} \equiv b \pmod p$ ，简单变换一下， $a^{jS} \equiv b a^i \pmod p$ （此处要保证 $a^i \bot p$ ）。

于是将 $ba^i \mod p$ 存下来（小步），然后用 $a^{jS}$ 去找相等，得到答案，时间复杂度为 $O(\sqrt p \log p)$ （如果用的是 std::map ）。

i64 BSGS(i64 a, i64 b, i64 p){
    i64 S = std::sqrt(p) + 1;
    std::map<i64, i64> mp;
    i64 G = 1;
    for(i64 i = 0;i < S;i++, G = G * a % p){
        mp[b * G % p] = i;
    }
    for(i64 i = 1, Gj = G;i <= S;i++, Gj = Gj * G % p){
        if(mp.count(Gj)){
            return i * S - mp[Gj];
        }
    }
    return -1;
}

小小拓展

对于求解 $x^a \equiv b \pmod p$ ， $p$ 是质数。

对于任意一个 $x$ ，如果 $p$ 有原根 $g$ 的话，那么可以表示为 $x = g^k$ ，于是就会有 $g ^ {ka} \equiv b \pmod p$ ，也可用 BSGS 求解。

我们再用 $b = g^t$ ，于是 $g^{ka} \equiv g^t \pmod p$ 。

同时取离散对数 $ka \equiv t \pmod p$ ，这是线性同余方程，易于求解。

找到特解后，通解也是容易求的。

exBSGS

对于 $a\not\bot p$ 的情况，无法用 BSGS 求解，于是考虑构造让 $a' \bot p'$ 。

再次之前，可以向讲讲其它做法，其实同样也可以结合扩展欧拉定理，将答案上界设为 $2\varphi(p)$ ，同样的扩展欧拉定理是一个 rho 型（想象一下好了），所以满足 $a^{jS} \equiv ba^i \pmod p$ 的就只有两个可能的答案，然后一通乱搞就好了。

在考虑较简单实现的做法（上面的挺烦的）：

通过取模的性质，去掉 $d = \gcd(a, p)$ ，得到

$\dfrac ad a^{x - 1} \equiv \dfrac bd \pmod {\dfrac pd}$

注意如果 $d \not | b$ ，那么就会无解，如果对于上述情况，依旧不满足 $a\bot \dfrac pd$ 那么重复以上操作，就会最终变成这个形式：

$$ \dfrac {a ^ k} {D} a ^ {x - k} \equiv \dfrac bD \pmod {\dfrac pD} $$

把形式转化成 BSGS 的，

$$ a ^ {x - k} \equiv \dfrac b {a ^ k} \pmod {\dfrac pD} $$

注意特判答案小于 $k$ 的情况，最后的 BSGS 的答案一定要加 $k$ 。

i64 gcd(i64 a, i64 b){ return !b ? a : gcd(b, a % b); }
i64 exgcd(i64 a, i64 b, i64 &x, i64 &y){
    if(b == 0){
        x = 1, y = 0;
        return 0;
    }
    i64 t = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return t;
}
i64 BSGS(i64 a, i64 b, i64 p){
    i64 S = std::sqrt(p) + 1;
    std::map<i64, i64> mp;
    i64 G = 1;
    for(i64 i = 0;i < S;i++, G = G * a % p){
        mp[b * G % p] = i;
    }
    for(i64 i = 1, Gj = G;i <= S;i++, Gj = Gj * G % p){
        if(mp.count(Gj)){
            return i * S - mp[Gj];
        }
    }
    return -1;
}
i64 exBSGS(i64 a, i64 b, i64 p){
    b %= p;
    if(b == 1 || p == 1) return 0;
    i64 k = 0, val = 1;
    while(true){
        i64 g = gcd(a, p);
        if(g == 1) break;
        if(b % g) return -1;
        b /= g, p /= g;
        val = val * a / g % p;
        k ++;
        if(b == val) return k;
    }

    i64 inv, y;
    exgcd(val, p, inv, y);
    inv = (inv % p + p) % p;
    b = b * inv % p;
    a %= p;
    i64 res = BSGS(a, b, p);
    return res == -1 ? -1 : res + k;
}

数论分块

不平均的东西也有如同平均一样的美，就像你一样。

求解问题

用以求解

$$ F(n) = \sum_{i = 1}^n f(i)g(\lfloor \dfrac ni \rfloor) $$

的类似式子，但前提是要能求出 $f(i)$ 的前缀和。

其实也叫整除分块。

操作方式

因为 $1\sim n$ 的区间内，存在很多相同的 $\lfloor \dfrac ni \rfloor$ ，若果将这些一起考虑，就会变得很简单了。

对于 $\lfloor \dfrac nx \rfloor = \lfloor \dfrac nk \rfloor$ 的最大的 $x$ 为 $\lfloor \dfrac n{\lfloor \frac nk \rfloor} \rfloor$ ，易证故此处不证。

$$ \sum_{i = 1}^n f(i)g(\lfloor \dfrac ni \rfloor) $$

就可以以这种方式求解：

int F(int n){
    int res = 0;
    for(int l = 1, r;l <= n;l = r + 1){
        r = (n / (n / l));
        res += (pre_f[r] - pre_f[l - 1]) * g(n / l);
    }
    return res;
}

复杂度分析

考虑分类讨论：

• $i < \sqrt N$ 时，对于 $\lfloor \dfrac ni \rfloor$ 显然只有 $O(\sqrt N)$ 种取值；
• $i > \sqrt N$ 时，对于 $\lfloor \dfrac ni \rfloor$ 一定有 $\lfloor \dfrac ni \rfloor < \sqrt N$ ，所以也只有 $O(\sqrt N)$ 种取值。

所以复杂度是 $O(\sqrt N)$ 的。

莫比乌斯函数

听说莫比乌斯队长发明了莫队算法。

定义

一般记为 $\mu(n)$ ，定义为

$$ \mu(n) = \begin{cases} 1, &n = 1 \\ 0, &\exist d > 1, d^2 | n \\ (-1)^k, &\operatorname{otherwise.} (k = \sum_{d | n, d \in \operatorname{Prime}} 1) \end{cases} $$

本质上来讲 $\mu(n)$ 是一个容斥。

性质

积性函数

跟欧拉函数一样，这也是个积性函数，也可以筛法求。

std::vector<int> sieve(int n){
    std::vector<int> np(n, false), prime, mu(n);
    mu[1] = 1;
    for(int i = 2;i < n;i++){
        if(!np[i]){
            prime.emplace_back(i);
            mu[i] = -1;
        }
        for(auto p : prime){
            if(1ll * i * p >= n) break;
            np[i * p] = true;
            if(i % p == 0){
                mu[i * p] = 0;
                break;
            }
            mu[i * p] = -mu[i];
        }
    }
    return mu;
} // Eular Sieve

常用性质

$$ \sum_{d | n} \mu(d) = [n = 1] $$

狄利克雷卷积的形式为 $\mu \circ 1 = \varepsilon$

这也说明了 $\mu$ 是 $1$ 的逆元。

证明：

由唯一分解定理，设 $n = \prod_{i = 1}^k p_i^{c_i}, m = \prod_{i = 1}^k p_i$ 。

$\sum_{d | n} \mu(d) = \sum_{d | m} \mu(d) = \sum_{i = 0}^k \begin{pmatrix} k \\ i \end{pmatrix} (-1) ^i = \big(1 + (-1)\big) ^ k$

对于 $k$ 仅在 $n = 1$ 时为 $0$ ，得证。

莫比乌斯变换

对于两个数论函数 $f(n), g(n)$ ，满足 $f = g \circ 1$ ，则有 $g = f \circ \mu$ 。

用先前 $\mu \circ 1 = \varepsilon$ 可以优雅证明，下给出一般的证明方式：

$$ \begin{aligned} g(n) &= \sum_{d | n} \mu(d)f(\lfloor \dfrac n d \rfloor) \\ &= \sum_{d | n} \mu(d) \sum_{k | \lfloor \frac n d \rfloor} g(\lfloor \dfrac n {kd} \rfloor) \\ &= \sum_{k | n} g(k) \sum_{d | \lfloor \frac nk \rfloor} \mu(d) \\ &= g(n) \end{aligned} $$

最后一步是因为后面依托只有在 $\lfloor \dfrac nk \rfloor = 1$ 时才为 $1$ ，否则为 $0$ ，此时 $k = n$ ，所以就为 $g(n)$ 。

所以就会有很多有意思的结论，比如先前 $\operatorname{id} = \varphi \circ 1$ ，由上变换可得 $\varphi = \operatorname{id} \circ \mu$ ，有的题目会让你求一些抽象的 $\varphi$ 可能会用到。

应用实例

$$ \sum_{i = 1} ^ n \sum_{j = 1} ^ n \gcd(i, j) $$

还是这个老朋友，但是换一个形式，我们求这个（下面的推算，保证 $n < m$ ）：

$$ \begin{aligned} &\sum_{i = 1} ^ n \sum_{j = 1} ^ m \operatorname{lcm}(i, j) \\ =& \sum_{i = 1} ^ n \sum_{j = 1} ^ m \dfrac {ij}{\gcd(i, j)} \\ =& \sum_{d = 1}^n \sum_{i = 1} ^ n \sum_{j = 1} ^ m \dfrac {ij}{d}[\gcd(i, j) = d] \\ =& \sum_{d = 1}^n \sum_{i = 1} ^ {\lfloor \frac nd \rfloor} \sum_{j = 1} ^ {\lfloor \frac md \rfloor} dij[\gcd(i, j) = 1] \\ =& \sum_{d = 1}^n d\sum_{i = 1} ^ {\lfloor \frac nd \rfloor} \sum_{j = 1} ^ {\lfloor \frac md \rfloor} ij \sum_{k | \gcd(i, j) = 1} \mu(k) \\ =& \sum_{d = 1}^n d\sum_{k = 1}^n k^2\mu(k) \sum_{i = 1} ^ {\lfloor \frac n{kd} \rfloor} \sum_{j = 1} ^ {\lfloor \frac m{kd} \rfloor} ij \\ =& \sum_{d = 1}^n d\sum_{k = 1}^n k^2\mu(k)f(\lfloor \frac n{kd} \rfloor, \lfloor \frac m{kd} \rfloor) & f(x, y) &= \dfrac{x(x + 1)y(y + 1)}4\\ =& \sum_{d = 1}^n g(\lfloor \frac nd \rfloor, \lfloor \frac md \rfloor)d & g(x, y) &= \sum_{k = 1}^n k^2\mu(k)f(\lfloor \frac n{kd} \rfloor, \lfloor \frac m{kd} \rfloor) \end{aligned} $$

其中 $g(x, y)$ 可以用数论分块在 $O(\sqrt N)$ 的时间内完成，然后再观察最后一个式子发现其实也是一个数论分块，所以也可以在 $O(\sqrt N)$ 的时间内解决，所以时间复杂度为 $O(N)$ 。

杜教筛

只要构造出打破偏见的桥梁，你们就能在一起了。

构造思想

对于数论函数 $f(n)$ ，求 $S(n) = \sum_{i = 1}^n f(i)$ 。

对于一个数论函数 $g(n)$ ，一定有：

$$ \begin{aligned} &\sum_{i = 1}^n (f\circ g)(i)\\ =& \sum_{i = 1}^n\sum_{d | i} g(d)f(\dfrac i d) \\ =& \sum_{d = 1}^n g(d) \sum_{i = 1}^{\lfloor\frac nd\rfloor} f(i) \\ =& \sum_{i = 1}^n g(i)S(\lfloor\dfrac nd\rfloor) \end{aligned} $$

然后我们知道，随手相减，就可以得到：

$$ \begin{aligned} g(1)S(n) &= \sum_{i = 1}^n g(i)S(\lfloor\dfrac nd\rfloor) - \sum_{i = 2}^n g(i)S(\lfloor\dfrac nd\rfloor) \\ &= \sum_{i = 1}^n (f\circ g)(i) - \sum_{i = 2}^n g(i)S(\lfloor\dfrac nd\rfloor) \end{aligned} $$

如果能做到快速算 $\sum_{i = 1}^n (f\circ g)(i)$ ，那么就可以数论分块算 $\sum_{i = 2}^n g(i)S(\lfloor\dfrac nd\rfloor)$ 轻松解决这个问题，最后除以 $g(1)$ 就是 $S(n)$ 。

时间复杂度分析

对于数论分块我们知道，可以分为 $\sqrt N$ 上下的两部分。

所以有，

$$ \begin{aligned} T(n) &= \sum_{i = 1}^{\lfloor\sqrt N\rfloor} O(\sqrt i) + \sum_{i = 2}^{\lfloor\sqrt N\rfloor} O(\sqrt {\dfrac n i}) \\ &= O( \int_{0}^{\sqrt n} (\sqrt x + \sqrt{ \dfrac n x}) \operatorname{d}x) \\ &= O(\dfrac 23 (\sqrt n) ^ {\frac 32} + 2\sqrt n \sqrt{\sqrt N}) \\ &= O(n ^ {\frac 3 4}) \end{aligned} $$

但是我们可以预处理其中一部分，所以还有，（下设阈值为 $B(B \geq \sqrt N)$ ）

$$ \begin{aligned} T(n) &= T'(B) + \sum_{i = 1}^{\lfloor \frac n B \rfloor} O(\sqrt{ \dfrac n i}) \\ &= T'(B) + O( \int_{0}^{\frac n B} \sqrt {\dfrac n x }\operatorname{d}x) \\ &= T'(B) + O(2 \sqrt n \sqrt {\dfrac n B}) \\ &= O(B + 2 \dfrac n {\sqrt B}) \\ &\geq O(n ^ {\frac 2 3}) \end{aligned} $$

此处默认预处理复杂度为 $O(n)$ ，即线性筛，此时取等条件是 $B = n ^ {\frac 23}$

实用案例

求 $S(n) = \sum_{i = 1}^n \mu(i)$

已知 $\mu \circ 1 = \varepsilon$ ，所以可以：

$$ S(n) = 1 - \sum_{i = 2}^n S(\lfloor\dfrac nd\rfloor) $$

对于较大的 $n$ 所对应的 $S(n)$ 要用 std::map 或 std::unordered_map 存一下，其它正常数论分块。

对于初学者，可以尝试手推一下用杜教筛求 $S(n) = \sum_{i = 1}^n \varphi(i)$ 。