数论

数学能不能去死啊！！！

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

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) ，所以将这个式子拓展出去得到：

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

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} ：

其他分析

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

类欧几里德算法

是一类用于求类似欧几里得算法过程的题目，一般可以在题目中发现 ​\lfloor \rfloor 之类的符号求和，或者一条直线截取一个平面之类的，当然也有可能是很多类似的构造，比较灵活，应题而变。

求解过程

接下来我们基于模板题研究欧几里得算法的过程。

求解 ​f(a, b, c, n) = \sum_{i = 0}^n \lfloor \dfrac {ai + b}{c} \rfloor

首先，我们可以把欧几里得的过程拉出来看一眼。

就可以分情况讨论了：

对于 ​a = 0 的情况，我们可以得到 ​f(a, b, c, n) = (n + 1)\lfloor \dfrac b c \rfloor 。

对于 ​a \geq c 或者 ​b \geq c 的情况，

接下来来到相对复杂的情况，对于 ​a < c 且 ​b < c 的情况，可以考虑如下变换，

至于剩下那两个就随她去吧。

几何意义

其实类欧的过程就类似于一个直线 ​l : ax - cy + b = 0 下方（包含直线上）整点计数 的问题。

​a \geq c 或者 ​b \geq c 的情况，相当于再直线 ​l 下方找一条 最接近 ​l 的 整系数直线 ​q （也就是 ​y = kx + b 中的 ​k \in \mathbb{Z} ） ，然后将 ​q 减掉得到 ​l' = l - q ，再对 ​l' 求解与 ​q 的贡献相加。

​q 由于是一个整系数的直线，所以截得的贡献等价于一个 梯形面积 ，套用梯形的面积公式，可以得到与上文相同的式子。

​a < c 且 ​b < c 的情况，也是很好理解的，就是将直线下方的贡献转变为总点数减去直线上方贡献，操作就是 翻转坐标轴 ，依然可以得到上文相同的式子。

所以其实 ​a < c 且 ​b < c 的情况，就是对原来的式子 取倒数 的过程（翻转坐标轴会使得斜率 ​k 取倒数）。

可以通过几何意义，比较显然的认为复杂度是 ​O(\log \min\{a, c, n\}) 。

欧拉函数

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

定义

表示小于等于 ​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;
}

欧拉反演

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

内容：

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

应用实例

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

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

也可以 ​O(N) 求解。

欧拉定理 & 费马小定理

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

费马小定理

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

证明：

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

因为 ​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 ，得到费马小定理。

扩展欧拉定理

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

线性同余方程

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

定义

求形如

​[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} 。

所以，就有 ​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 。

求线性同余方程组

最小的一个答案。

解法

令 ​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} ，所以

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 。

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

已知两同余方程：

可以转化成为

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

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

卢卡斯定理

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

定理内容

对于质数 ​p ，则满足：

其中若 ​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 （用二项式定理可证），借此进行生成函数。

对于 ​x ^ b 求系数：

稍微变换一下，得证。

实用案例

对于较小的质数 ​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) 使得：

其中 ​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 那么重复以上操作，就会最终变成这个形式：

把形式转化成 BSGS 的，

注意特判答案小于 ​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(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 ，易证故此处不证。

就可以以这种方式求解：

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) 是一个容斥。

性质

积性函数

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

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

常用性质

狄利克雷卷积的形式为 ​\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 可以优雅证明，下给出一般的证明方式：

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

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

应用实例

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

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

杜教筛

只要构造出打破偏见的桥梁，就能到达彼岸。

构造思想

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

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

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

如果能做到快速算 ​\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 上下的两部分。

所以有，

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

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

实用案例

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

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

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

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

Min_25 筛

我虽然只能看见你最朴素的一面，却可以合成出你的全貌。

定义

与杜教筛相同，都是一种解决数论函数求和的亚线性筛法，因为由 Min_25 提出，所以叫 Min_25 筛。

一般来讲， Min_25 筛受限制较小，可以打过 Powerful Number 筛，等其他看起来很厉害的筛法。

主要思想

首先， Min_25 筛要求能够对于所有 质数 ​p ，快速求出 ​f(p ^ c) ，并且可以将数论函数 ​f 拆分成多个能够快速求出 ​f(p) 的 完全积性函数 的和（ 注意此处的完全积性函数是对于质数处的拆分 ，合数不需要考虑，因为之后需要利用该性质的只是用于计算质数 ​p 的 ​f(p) 之和）。

一个很典型的例子就是 多项式函数 。

主要的思想是根据 质因数 分讨，由于质数个数较多，所以把原来的求和分成 质数和合数 讨论，保证合数 最小质因子 的个数在 ​O(\dfrac {\sqrt N}{\log N}) 级别。

另一种理解方式就是考虑 埃氏筛 的过程。

记 ​\operatorname{lpf}(n) 表示 ​n 的最小质因数， ​p_i 表示第 ​i 个质数，并且钦定 ​p_0 = 1 。

首先考虑一个式子 ​S(n, j) = \sum_{i = 2}^n f(i)[\operatorname{lpf}(i) > p_j]

接下来，基于质数和合数分类讨论。

对于质数部分的贡献就是 ​\big( \sum_{i = 1}^n f(i)[i \in \operatorname{Prime}] \big) - \sum_{i = 1}^j f(p_i) 。

对于合数部分的贡献稍微有点复杂但是也很简单，直接枚举最小质因子及其次数，那么贡献就是 ​\sum _ {k > j} \sum _ {c = 1} ^ {p_k^c \leq n} f(p_k^c)(S(\dfrac{n}{p_k^c}, k) + [c > 1]) 。

因为 ​f(p) 的值之前已经算过了，所以这次 ​c = 1 时，是不需要算上 ​1 的贡献的，至于为什么要加 ​1 是因为 ​S 函数的 ​i 是从 ​2 开始算的。

为了方便表示，我们记 ​h(n) = \sum_{i = 1}^n f(i)[i \in \operatorname{Prime}], S_p(j) = \sum_{i = 1}^j f(p_i) ，于是就有。

那么答案就是 ​S(n, 0) + f(1) ，至于为什么加 ​f(1) ，原因同上。

第一部分

现在可以考虑先算 ​h(n) ，然后 ​S_p(j) 可以线筛求出。

还是考虑构造一下，让 ​f'(n) 表示 ​f(n) 的拆分的一个 完全积性函数 ，令 ​g(n, j) = \sum_{i = 2}^n f'(i) [i \in \operatorname{Prime} \operatorname{or} \operatorname{lpf}(i) > p_j] 。

很显然的是 ​h(n) = g(n, |p|) 。

接着考虑递推计算 ​g 函数的值。

首先考虑 ​p_j^2 > n 时，由于 ​p_j 是最小质因子，所以无法造成新的贡献，因此 ​g(n, j) = g(n, j - 1) 。

接下来考虑 ​p_j^2 \leq n ，我们就需要在 ​g(n, j - 1) 的基础上去除以 ​p_j 为最小质因子的贡献并加上多余的减去的质数的贡献。

具体来讲，就是 ​g(n, j - 1) 中包括了 ​\lfloor \dfrac {n}{p_j} \rfloor 个含有质因子的数，由于 ​\operatorname{lpf}(i) > p_{j - 1} 的限制，所以剩余的质因子一定大于等于 ​p_j ，那么这一部分的贡献应该是 ​g(\lfloor \dfrac {n}{p_j} \rfloor, j - 1) ，但是这也包括了一部分符合条件的质数，这部分的答案是 ​\sum_{i = 1}^{j - 1} f'(p_i) 。

综上所述，我们可以出来一个关于 ​g 的递推式：

然后经过了注意力的考察，我们发现所有的要求求的 ​g(x, j) 一定满足 ​x = \lfloor \dfrac n k \rfloor 的形式。

很自然地，会想到数论分块，对于 ​k 与 ​\sqrt n 进行分讨：​k \leq \sqrt n 的情况，用一个 ​\operatorname{id1} 记录 ​k ； ​k > \sqrt n 的情况，用一个 ​\operatorname{id2} 记录 ​\lfloor \dfrac nk \rfloor 。

由于我们只要求 ​g(n, |p|) 的值，所以 ​j 可以滚掉。

其中初值 ​g(n, 0) = \sum_{i = 2}^n f'(i) 。

接下来，我们只要对于所有的 ​f' 对应的 ​g 求和就可以得到 ​h 函数需要的值。

关于时间复杂度

复杂度是比较好算的

第二部分

对于 ​S 函数的计算，其实已经非常简单了，就是十分暴力的枚举计算就好了。

这个复杂度依然是亚线性的，无法爆炸。

这个部分的复杂度被证明为 ​O(N^{1 - \varepsilon}) ，这个我不会证。

其实还有第二种方法，本质上是洲阁筛的第二部分，理论复杂度也为 ​O(\dfrac {N^{\frac 34}}{\log N})，但是常数巨大，不如暴力。

实用案例

首先还是考虑求 ​\sum_{i = 1}^n \mu(i) 。

首先 ​f(p) = -1 可以拆分为 ​f'(p) = -1 ，则 ​g(n, 0) = 1 - n ，直接套板子就行。

依然是建议自己手推 ​\sum_{i = 1}^n \varphi(i) 。

接下来是模板题 ​f(p ^ k) = p^k (p ^ k - 1) ，显然可以拆为 ​p^2 和 ​-p 的和，于是就有了两个 ​g ，其中 ​g(n, 0) 分别为 ​\dfrac {n(n + 1)}{2} - 1 和 ​\dfrac {n(n + 1)(2n + 1)}{6} - 1 。

放一下模板题的代码

#include <bits/stdc++.h>
#define endl "\n"
typedef long long i64;

const int _ = 2e5 + 10;
const int P = 1e9 + 7;

bool np[_];
i64 prs[_], sum1[_], sum2[_];
int tot = 0;
void sieve(int n){
    for(int i = 2;i <= n;i++){
        if(!np[i]) prs[++ tot] = i;
        for(int j = 1;j <= tot && prs[j] * i <= n;j++){
            np[prs[j] * i] = true;
            if(i % prs[j] == 0) break;
        }
    }
    for(int i = 1;i <= tot;i++){
        sum1[i] = (sum1[i - 1] + prs[i]) % P;
        sum2[i] = (sum2[i - 1] + prs[i] * prs[i] % P) % P;
    }
    // 线筛求 p 和 p^2 的和
    return ;
}

i64 n, sq;
int id1[_], id2[_], cnt = 0;
i64 g1[_], g2[_], w[_];
i64 S1(i64 x){ return x * (x + 1) / 2 % P; }
i64 S2(i64 x){ return x * (x + 1) % P * (2 * x % P + 1) % P * 166666668 % P; }
//    1/6 = 166666668 (mod 10^9 + 7)
i64 id(i64 x){
    if(x <= sq) return id1[x];
    return id2[n / x];
} // 对于两类 id 进行分讨
i64 S(i64 x, int j){
    if(prs[j] > x) return 0;
    i64 res = ((g2[id(x)] - g1[id(x)] + P) % P - (sum2[j] - sum1[j] + P) % P + P) % P;
    for(int i = j + 1;i <= tot && prs[i] * prs[i] <= x;i++){
        for(i64 c = 1, sp = prs[i];sp <= x;sp *= prs[i], c ++){
            res = (res + sp % P * (sp % P - 1) % P * (S(x / sp, i) + (c > 1) % P)) % P;
        } // 枚举次数
    } // 枚举最小质因数
    return res;
}

void min_25(){
    sq = std::sqrt(n);
    sieve(sq);
    // 最小质因数最大只有 \sqrt n

    for(i64 l = 1, r;l <= n;l = r + 1){
        r = n / (n / l);
        w[++ cnt] = n / l;
        g1[cnt] = S1(n / l % P) - 1;
        g2[cnt] = S2(n / l % P) - 1;
        // 求 g(x, 0)
        if(n / l <= sq) id1[n / l] = cnt;
        else id2[n / (n / l)] = cnt;
        // 对于两类 id 进行分讨
    }

    for(int i = 1;i <= tot;i++){
        for(int j = 1;j <= cnt && prs[i] * prs[i] <= w[j];j ++){
            int k = id(w[j] / prs[i]);
            g1[j] = (g1[j] - prs[i] * (g1[k] - sum1[i - 1] + P) % P + P) % P;
            g2[j] = (g2[j] - prs[i] * prs[i] % P * (g2[k] - sum2[i - 1] + P) % P + P) % P;
        } // 递推 + 滚动 求 g(n, |p|)
    }

    std::cout << (S(n, 0) + 1) % P << endl;
    // 一般来讲 f(1) = 1 ，具体情况具体对待
}

signed main(){
    std::cin.tie(nullptr) -> sync_with_stdio(false);
    std::cin >> n;
    min_25();
    return 0;
}