《一些特殊的数论函数求和问题》阅读随笔

好至少它教会了我如何把质数求和转化成积分的渐进 对着 $\pi(x)$ 微就行了 然后直接 $u\text dv = uv- v\text du$

18.3k……阿巴阿巴

目录

目录

• 引言
• $1$ 基础知识
• $2$ 积性函数求和
  • $2.1$ $2.2$ 什么是 $\text{Min_25}$ 筛？
  • $2.3$ 时间复杂度
• $3$ 素数的 $k$ 次幂前缀和
  • $3.1$ 关于 $P_{2,k}(n, a)$
  • $3.2$ 关于 $\phi_k(n, a)$
    • $3.2.1$ $p\le \sqrt B$
    • $3.2.2$ $p> n^{\frac 13}$
    • $3.2.3$ $\sqrt B < p \le n^{\frac 13}$
  • $3.3$ 优化
• $4$ 约数函数前缀和
  • $4.1$ 题面
  • $4.2$ 朴素转化
  • $4.3$ $\text{Stern-Brocot Tree}$ 和 $\text{Farey}$ 序列
  • $4.4$ 利用 $\text{Stern-Brocot Tree}$ 切割曲线
    • $4.4.1$ 坐标系的转换
    • $4.4.2$ 计算整点数
  • $4.5$ 时间复杂度
  • $4.6$ 优化
  • $4.7$ 拓展

引言

这玩意挺常见的。而且你会发现有些做法能套在很多题上。咱选几个做法聊聊。
首先我讲讲拓展埃式筛求前缀和，然后再聊聊素数和，最后给你讲讲怎么用面积拟合函数。
反正不要钱，多少看一点~

$1$ 基础知识

定义 $\textbf{1.1}$（数论函数）：一个函数 $f$ 被称为数论函数，当且仅当其定义域为正整数集，陪域为复数集。

定义 $\textbf{1.2}$（积性函数）：一个函数 $f$ 被称为积性函数，当且仅当 $\forall a,b\text{ s.t. } (a,b) = 1,\ f(ab) = f(a)f(b)$。

引理 $\textbf{1.1}$（数论分块）：给定正整数 $n$，对于任意正整数 $m$，$\left\lfloor\frac nm \right\rfloor$ 的取值只有 $O(\sqrt n)$ 种。

证明：当 $m\le \sqrt n$ 时 $m$ 的取值有 $O(\sqrt n)$ 个，而当 $m > \sqrt n$ 时 $0 \le \left\lfloor\frac nm \right\rfloor \le \sqrt n$，因此 $\left\lfloor\frac nm \right\rfloor$ 的取值有 $O(\sqrt n)$ 个。

定理 $\textbf{1.1}$（素数定理）：令 $\pi(x)$ 为不大于 $x$ 的素数个数。有 $\pi(x) = O(\frac {x}{\ln x})$。

证明略。

$2$ 积性函数求和

你发现这种题很常见。16 年任之洲同样写了积性函数求和的几种方法，并得到了一种适用范围更广的方法，又称“洲阁筛法”。事实上其更概括性的名称为“扩展埃式筛”。
但是这种筛法也有不足。首先是这个做法从 $O(\frac {n}{\log n})$ 优化到 $O(\frac{n^{\frac 34}}{\log n})$ 的过程中很多部分较难理解，初学者很难上手，其次是这个做法常数比较大，代码实现较复杂。但扩展埃式筛并不只有这一种做法。实际上，存在一种实现难度更低且在小范围数据下表现更好的做法。它与洲阁筛的思想类似，但理解较为容易。由于该算法由 $\text{Min_25}$ 使用后才开始普及，其也被称为 $\text{Min_25}$ 筛。
本文中讨论的扩展埃式筛即 $\text{Min_25}$ 筛。

$2.1$ $2.2$ 什么是 $\text{Min_25}$ 筛？

$2.3$ 时间复杂度

对于质数部分的计算，复杂度证明和洲阁筛完全相同。考虑每个 $i = \left\lfloor\frac nm \right\rfloor$ 只会在枚举不超过 $\sqrt i$ 的质数时产生贡献。由定理 1.1，复杂度渐进地是

$$\begin{aligned} &\sum_{i=1}^{\sqrt n}\mathcal O \left( \frac{\sqrt i}{\log \sqrt i} \right) + \sum_{i=1}^{\sqrt n}\mathcal O \left( \frac{\sqrt { \left\lfloor\frac ni \right\rfloor } }{\log \sqrt { \left\lfloor\frac ni \right\rfloor } } \right) \\ = \ & \mathcal O \left( \sum_{i=1}^{\sqrt n}\frac{\sqrt { \left\lfloor\frac ni \right\rfloor } }{\log \sqrt { \left\lfloor\frac ni \right\rfloor } } \right) \\ = \ & \mathcal O \left( \int_1^{\sqrt n} \frac{\sqrt { \frac nx } }{\log \sqrt { \frac nx } } \text dx \right) \\ = \ & \mathcal O \left(\frac{n^{\frac 34}}{\log n}\right) \end{aligned}$$

$3\to 4$ 考虑首先将分母变成渐进的 $\log n$，然后只需要对分子积分，而这个处理是平凡的。可以写出

$$\mathcal O \left( \int_1^{\sqrt n} \frac{\sqrt { \frac nx } }{\log \sqrt { \frac nx } } \text dx \right) = \mathcal O \left( \frac{ \int_1^{\sqrt n} \sqrt { \frac nx } \text dx }{\log n } \right) = \mathcal O \left( \frac {n^{\frac 12} n^{\frac 14} }{\log n } \right) = \mathcal O \left( \frac { n^{\frac 34} }{\log n } \right)$$

然后考虑统计答案部分。

定义 $\textbf{2.1}$：对正整数 $n$，我们定义 $big(n)$ 为 $n$ 最大的质因子，$small(n)$ 为 $n$ 最小的质因子，$cnt(n)$ 为 $n$ 的可重质因子个数。并令 $big(1) = small(1) = \infty, \ cnt(1) = 0$。

如果我们直接按 $s(n,k)$ 的转移式 dp，空间复杂度 $O(\sqrt n)$，时间复杂度按质数部分的分析能得到 $\mathcal O \left( \frac { n^{3/4} }{\log n } \right)$。
如果直接暴力计算呢？在小数据的测试下，该算法似乎比 $\mathcal O \left( \frac { n^{3/4} }{\log n } \right)$ 更优秀一些。我们尝试分析该算法的复杂度。

考虑计算 $s(n,k)$ 时的递归。我们的时间复杂度就是进入递归的次数，即 $g(n) - g(p_k)$ 被累加入答案的次数。对于从一个值开始的递归，我们不妨假设累加入的所有质数 $p$ 处的贡献，实际上产生于 $l \times p$ 处，那我们将这次递归对复杂度的贡献计入 $l\times p_k$。$p_k$ 显然是 $big(l)$。由于每个函数值只统计一次更新，因此每次的 $l\times p_k$ 不同。同时 $\forall l \times p_k \le n$ 都会出现，因此该算法的复杂度就是 $l\times big(l) \le n$ 的 $l$ 的个数。

引理 $\textbf{2.1}$：对于给定的实常数 $0 < \alpha < 1$，我们令 $\mathcal Q(n) = \left\{k \le n : \ big(k) \le k^\alpha\right\}$，则 $|\mathcal Q(n)|\sim n\rho\left(\frac 1\alpha\right)$，其中 $\rho$ 为 Dickman 函数。

引理 $\textbf{2.2}$：对于任意 $x > 1$，$\rho(x) > 0$ 随 $x$ 增大而迅速衰减，$\rho(x) \approx x^{-x}$。

查阅此处内容可以获得关于这两个引理的信息。

定理 $\textbf{2.1}$：令 $\mathcal M(n) = \left\{i : i \times big(i) \le n \right\}$，则对于任意 $0 < \alpha < 1$，$|\mathcal M(n)| = \Omega(n^\alpha)$。

证明：
不妨取 $\mathcal P = \left\{i \le n^{\alpha} :big(i) \le i^{1 - \alpha} \right\}$。由引理知 $|\mathcal Q(n)| = \Omega(n^\alpha)$。而 $\forall x\in \mathcal P,\ x \times big(x) \le x\times x^{1-\alpha} \le n^{2\alpha - \alpha^2} \le n$，且 $\mathcal P \in \mathcal M(n)$，于是有 $|\mathcal M(n)| = \Omega(n^{\alpha})$。

定理 $\textbf{2.2}$：$|\mathcal M(n)| = \Theta(n^{1 - \epsilon})$。

证明：
我们取前 $\log \log n$ 个质数。令 $p_i$ 表示第 $i$ 个质数，最大质因子不超过 $p_i$ 的（小于 $n$ 的）数的个数不超过 $(\log n)^{\log log n} \le \mathcal O(\frac {n}{\log \log n})$。而其他 $k \text{ s.t. }k \times big(k) \le n$ 不超过 $\frac{n}{\log \log n}$，因此 $|\mathcal M(n)| = \mathcal O(\frac{n}{\log \log n})$。
由定理 $2.1$，$|\mathcal M(n)| = \Theta(n^{1 - \epsilon})$。

所以这玩意是 $\Theta(n^{1 - \epsilon})$ 的，但是实际跑出来的效果比较快？为啥？

引理 $\textbf{2.3}$：$\sum_{p\ge n,\ p \text{是质数}} \frac {1}{p^2} = \Theta\left(\frac {1}{n\ln n}\right)$

证明：
进行初等积分。

$$\begin{aligned} & \sum_{p\ge n,\ p \text{是质数}} \frac {1}{p^2} \\ = \ & \int_n^\infty \frac{\text d \pi(x)}{x^2} \\ = \ & \frac{\pi(x)}{x^2}\bigg\lvert_n^{\infty} - \int_n^{\infty} \pi(x) \text d\left(\frac{1}{x^2}\right) \\ = \ & \Theta\left( \frac{1}{x\ln x}\bigg\lvert_n^{\infty} + 2\int_n^{\infty} \frac{\text dx}{x^2 \ln x} \right) \\ = \ & \Theta\left( -\frac{1}{n\ln n} + 2\int_{\frac 1n}^{0} \frac{x^2\text d\frac 1x}{ \ln \frac 1x} \right) \\ = \ & \Theta\left( -\frac{1}{n\ln n} - 2\int_{0}^{\frac 1n} \frac{x^2 (-x^{-2})\text dx}{ -\ln x} \right) \\ = \ & \Theta\left( -\frac{1}{n\ln n} - 2\int_{0}^{\frac 1n} \frac{\text dx}{ \ln x} \right) \\ = \ & \Theta\left( -\frac{1}{n\ln n} - 2\ \text{li}\ \frac 1n \right) \end{aligned}$$

其中 $\text{li}$ 为对数积分函数。由经典结论可知当 $0 < n < 1$ 时 $\text{li}\ n \sim \frac{n}{\ln n}$。然后证完。

$1\to 2$ 考虑 $u\text dv = uv- v\text du$。任何发现 $\pi()$ 没了的情况考虑 $\pi(n)\sim \frac {n}{\ln n}$。
为容易理解，这里的推导和论文里有些许不同。

定理 $\textbf{2.3}$：当 $big(i) \ge n^{\frac 14}$ 时，满足 $i\times big(i) \le n$ 的 $i$ 至多有 $\mathcal O\left(\frac{n^{3/4}}{\log n}\right)$ 个。

证明：
当 $big(i) \ge n^{\frac 14}$ 时，我们枚举每个 $p = big(i)$，至多有 $\left\lfloor\frac{n}{p^2}\right\rfloor$ 个满足条件的数，因此这部分贡献了 $\mathcal O\left(\frac{n^{3/4}}{\log n}\right)$ 的复杂度。
打表可得，在 $n\le 10^{13}$ 时，对于 $\le n^{\frac 14}$ 的质数 $p$，满足 $big(i) = p$ 的数 $i$ 的数量 $=O(\sqrt n)$。只有 $\mathcal O\left(\frac{n^{1/4}}{\log n}\right)$ 个质数，因此这部分的复杂度也是 $\mathcal O\left(\frac{n^{3/4}}{\log n}\right)$。
这就证明了定理。

因此我们能发现，这样执行的复杂度是 $\mathcal O\left(\frac{n^{3/4}}{\log n}\right)$ 的。

$3$ 素数的 $k$ 次幂前缀和

先前的算法涉及到对素数处取值求前缀和的操作。复杂度是 $\mathcal O\left(\frac{n^{3/4}}{\log n}\right)$。然而有时我们需要的只是 $\sum_{p\le n,\ p\text{是质数}}p^k$。这时我们存在更优的做法。下面复杂度的推导中不包括求 $k$ 次幂的影响。

定义 $\textbf{3.1}$（前缀和函数）：对正整数 $n$，定义 $\pi_k(n)$ 为素数的 $k$ 次幂前缀和函数，即 $\pi_k(n) = \sum_{p\le n,\ p\text{是质数}}p^k$。特别的，$\pi_0(x)$ 即为素数计数函数 $\pi(x)$。

定义 $\textbf{3.2}$（部分筛函数）：对正整数 $n,a$，定义部分筛函数 $\phi_k(n, a) = \sum_{i=1}^n[small(i) > p_a]i^k$。当 $a=0$ 时 $\phi_k(n, a) = \sum_{i=1}^n i^k$。

定义 $\textbf{3.3}$（特定筛函数）：对正整数 $n,a$，定义特定筛函数 $P_{s,k} (n, a) = \sum_{i=1}^n[small(i) > p_a, cnt(i) = s]i^k$。

当 $n^{\frac 13} \le p_a = B \le x^{\frac 12}$ 时，可以枚举质因子得到 $\phi_k(n, a) = P_{0,k}(n, a) + P_{1,k}(n, a) + P_{2,k}(n, a) = 1 + \pi_k(n) - \pi_k(p_a) + P_{2,k}(n, a)$。

则 $\pi_k(n) = \phi_k(n, a) + \pi_k(p_a) - P_{2,k}(n, a) - 1$。
由于 $p_a \le \sqrt n$，$\pi_k(p_a)$ 可以直接筛出来。我们现在只需要考虑 $\phi_k(n, a)$ 和 $P_{2,k}(n, a)$ 的计算。

$3.1$ 关于 $P_{2,k}(n, a)$

考虑 $s = 2$。$i$ 能贡献 $i^k$，当且仅当 $i = pq$，且 $p_a < p \le q \le \frac {n}{p},\ B < q < \frac nB$。枚举 $p$ 后计算 $q$ 的贡献即可。
复杂度 $O(\frac{n}{B})$。

$3.2$ 关于 $\phi_k(n, a)$

引理 $\textbf{3.1}$：$\phi_k(n, a) = \phi_k(n, a - 1) - p_a^k \phi_k(\left\lfloor\frac{n}{p_a}\right\rfloor, a - 1)$

显然。

定理 $\textbf{3.1}$：$\phi_k(n, a) = \sum_{big(i) \le p_a, small(i) > p_b}\mu(i)i^k \phi_k(\left\lfloor\frac{n}{i}\right\rfloor, b)$

就是把引理 $3.1$ 多用了几次。
不难(?)导出 $\phi_k(n, a) = \sum_{big(i) \le p_a}\mu(i)i^k \left\lfloor\frac{n}{i}\right\rfloor^k$。

首先朴素计算 $n \le p_a = B$ 的贡献，然后只需要考虑 $\sum_{big(i) \le p_a < i} \mu(i)i^k \left\lfloor\frac{n}{i}\right\rfloor^k$ 的贡献。奇妙结论：

$$\sum_{big(i) \le p_a < i} \mu(i)i^k \left\lfloor\frac{n}{i}\right\rfloor^k = \sum_{B < i \le B \times small(i)} \mu(i)i^k \phi_k(\left\lfloor\frac{n}{i}\right\rfloor, \pi_0(small(i)) - 1)$$

这式子可以建出递归对应的二叉树得到。

换元。令 $p = small(i), m = i / p$。$\mu(i) = -\mu(m)$，也就导出了

$$-\sum_{p \le B}\sum_{small(m) > p, m \le B < mp} \mu(m) (mp)^k \phi_k(\left\lfloor\frac{n}{mp}\right\rfloor, \pi_0(p) - 1)$$

阈值法。

$3.2.1$ $p\le \sqrt B$

可以朴素求解。我们计算出所有 $\phi$，随后枚举 $p,m$ 统计。
由引理 $1.1$，有效状态是 $\frac{\sqrt{nB}}{\log B}$ 的。

$\text{otherwise } p > \sqrt B$。
如果 $m$ 不是素数，则由于 $small(m) > p$，$m > p^2 > B$，这样的 $m$ 不存在。
因此我们只需要考虑 $m$ 为质数的情况。也就是说下面的讨论中 $\mu(m) = -1$。
改变求和式为

$$\sum_{p \le B}\sum_{m > p, m \le B < mp} (mp)^k \phi_k(\left\lfloor\frac{n}{mp}\right\rfloor, \pi_0(p) - 1)$$

$3.2.2$ $p> n^{\frac 13}$

$\left\lfloor\frac{n}{mp}\right\rfloor < n^{\frac 13} < p$，则根据定义有 $\phi_k(\left\lfloor\frac{n}{mp}\right\rfloor, \pi_0(p) - 1) = 1$。
改变求和式为

$$\sum_{n^{1/3} < p \le B} p^k \sum_{m > p, m \le B < mp} m^k$$

我们可以枚举 $p$，求 $m$ 的贡献。筛出质数 $k$ 次方和计算。复杂度 $\mathcal O(B)$。

$3.2.3$ $\sqrt B < p \le n^{\frac 13}$

直接拿定义冲

$$\begin{aligned} & \sum_{\sqrt B < p \le n^{1/3}}\sum_{m > p, m \le B < mp} (mp)^k \phi_k(\left\lfloor\frac{n}{mp}\right\rfloor, \pi_0(p) - 1) \\ = \ & \sum_{\sqrt B < p \le n^{1/3}}\sum_{m > p, m \le B < mp} (mp)^k \sum_{mpr \le n, small(r) \ge p} r^k \\ = \ & \sum_{r = 1}^{\left\lfloor\frac nB\right\rfloor} r^k\sum_{\sqrt B < p \le n^{1/3}}\sum_{m > p, m \le B < mp} (mp)^k [mpr \le n, small(r) \ge p] \\ = \ & \sum_{r = 1}^{\left\lfloor\frac nB\right\rfloor} r^k\sum_{\sqrt B < p \le \min(n^{1/3}, small(r))}p^k\sum_{p < m \le \min(\lfloor\frac n{pr}\rfloor, b)} m^k \end{aligned}$$

于是问题来到了维护这边。
枚举 $p$，他能贡献的 $r$ 满足 $small(r) \ge p$，这里使用树状数组维护。数论分块维护 $\left\lfloor\frac n{pr}\right\rfloor$ 得到 $m$ 的贡献。在树状数组上查询的次数是 $\mathcal O(\sqrt {\frac np})$ 的。

于是问题来到了复杂度这边。

首先是修改。

定理 $\textbf{3.2}$：对于任意 $0 < \alpha < 1$，$n$ 以内恰有 $k \ge 1$ 个质因子且最小质因子不小于 $n^{\alpha}$ 的数的个数是 $\mathcal O\left(\frac{n}{\log^k n}\right)$ 的。

证明：
$k=1$ 即为定理 $1.1$。
归纳法。假设 $k < k_0$ 的所有情况都成立，现在证明 $k = k_0$ 时同样成立。
进行初等积分：

$$\begin{aligned} & \mathcal O\left(\int_{n^{\alpha}}^{\sqrt n}\frac{ \frac np }{\log^{k-1}\frac np} \text d \pi(p))\right) \\ = \ & \mathcal O\left( \frac{ \frac np }{\log^{k-1}\frac np} \frac {p}{\log p}\bigg\lvert_{n^{\alpha}}^{\sqrt n} \right) - \mathcal O\left( \int_{n^{\alpha}}^{\sqrt n}\frac {p}{\log p} \text d\frac{ \frac np }{\log^{k-1}\frac np} \right ) \\ = \ & \mathcal O\left( \frac{n}{\log^k n} \right) - \mathcal O\left( \int_{n^{\alpha}}^{\sqrt n}\frac {np}{\log p} \left(\frac{1}{p\log^{k-1}\frac np}\right)' \text dp \right ) \\ = \ & \mathcal O\left( \frac{n}{\log^k n} \right) - \mathcal O\left( \int_{n^{\alpha}}^{\sqrt n}\frac {np}{\log p} \left(\frac{-\log\frac np + k - 1}{p^2 \log^{k} \frac np}\right) \text dp \right ) \\ = \ & \mathcal O\left( \frac{n}{\log^k n} \right) + \mathcal O\left( \int_{n^{\alpha}}^{\sqrt n}\frac {n}{p\log p\log^{k-1} \frac np}\text dp \right ) \\ = \ & \mathcal O\left( \frac{n}{\log^k n} \right) \end{aligned}$$

感想：积分是用来分析的，不是用来创人的。
$3\to 4$ 用渐进擦除了 $k-1$，把 $\log$ 前的负号提了出去。

归纳成立。

upd: 这个似乎是错的，因为积分出了问题，实际上应当趋于 $\mathcal O(n/\log n)$。

我们在树状数组上进行操作的时候需要保证 $r$ 的最小质因子不小于 $\sqrt B$。有用的 $r$ 只有 $\mathcal O\left(\frac{n}{B\log n}\right)$ 种。乘入树状数组的复杂度后为 $\mathcal O\left(\frac{n}{B}\right)$。

然后是查询。

定理 $\textbf{3.3}$：$\sum_{p\le m}\sqrt{\frac np} = \mathcal O\left(\frac{\sqrt {mn} }{\log m}\right)$。

证明：
进行初等积分：

$$\begin{aligned} & \sum_{p\le m}\sqrt{\frac np} \\ = \ & \sqrt n \sum_{p\le m}\frac 1{\sqrt p} \\ = \ & \mathcal O\left( \sqrt n \int_{1}^m \frac{\text d\pi(p)}{\sqrt p} \right) \\ = \ & \mathcal O\left(\sqrt n \left(\frac{\pi(p)}{\sqrt p}\bigg\lvert_1^m - \int_{1}^m \pi(p)\text d\frac{1}{\sqrt p}\right) \right) \\ = \ & \mathcal O\left( \sqrt n \left(\frac{\pi(m)}{\sqrt m} +\frac 12 \int_{1}^m \frac{1}{\sqrt p\ln p}\text d p\right) \right) \\ = \ & \mathcal O\left( \sqrt n\frac{\pi(m)}{\sqrt m} \right) \\ = \ & \mathcal O\left(\frac{\sqrt {mn} }{\log m}\right) \end{aligned}$$

这个挺好理解。

然后带入 $m = n^{\frac 13}$ 可知询问次数是 $\mathcal O\left(\frac{n^{2/3}}{\log n}\right)$ 的。带上树状数组就是 $\mathcal O\left(n^{2/3}\right)$。

到这里分析完了计算 $\phi$ 的所有情况，总时间复杂度总结一下就是 $\mathcal O\left(n^{\frac 23} + \frac nB + \frac{\sqrt{nB}}{\log B}\right)$，取 $B = n^{\frac 13}\log^{\frac 23}n$ 得到后半部分复杂度是 $\mathcal O\left(\left(\frac n{\log n}\right)^{2/3}\right)$，总复杂度还是 $\mathcal O\left(n^{2/3}\right)$。

$3.3$ 优化

可以发现复杂度瓶颈在 $\sqrt B < p \le n^{\frac 13}$ 的查询部分。
于是问题又来到了维护这边。

抄式子：

$$\sum_{r = 1}^{\left\lfloor\frac nB\right\rfloor} r^k\sum_{\sqrt B < p \le \min(n^{1/3}, small(r))}p^k\sum_{p < m \le \min(\lfloor\frac n{pr}\rfloor, b)} m^k$$

我们进一步利用那个定理。分别考虑 $r$ 是否为质数：
如果 $r$ 是质数，若 $r \ge n^{\frac 13}$，则对所有满足条件的 $p$ 来说 $r$ 都有贡献。因此枚举 $p$，将质数部分统计完即可。复杂度 $\mathcal O\left( \frac{n^{2/3}}{\log n} \right)$。
$\text{otherwise }$有 $r \ge p^2$，且存在 $q$ 需要使得 $\frac xr > p^2$，则 $p \le x^{\frac 14}$。此时由定理 $3.3$ 可知询问次数为 $\mathcal O\left(\frac{x^{5/8}}{\log x}\right)$，这部分不会影响复杂度。

因此喜提 $\mathcal O\left(\left(\frac n{\log n}\right)^{2/3}\right)$ 的总复杂度。

指路【模板】Meissel–Lehmer 算法。
指路 OI Wiki。
最大的收获是重学积分（

$4$ 约数函数前缀和

$4.1$ 题面

定义 $\sigma(i)$ 为 $i$ 的约数个数，给定正整数 $n < 2^{63}$，求 $\sum_{i=1}^n \sigma(i)$。

$4.2$ 朴素转化

$$\begin{aligned} & \sum_{i=1}^n \sigma(i) \\ = \ & \sum_{x=1}^n \sum_{d|x} 1 \\ = \ & \sum_{d=1}^n \left\lfloor\frac nd\right\rfloor \end{aligned}$$

考虑转化

$$\begin{aligned} & \\ \sum_{i=1}^n \lfloor \frac n i \rfloor = \ & \sum_{i=1}^n \sum_{j=1}^n \ [ij \le n] \qquad (几何意义，变为\ xy = n\ 线下整点求和) \\ =\ & \sum_{i=1}^{\sqrt n}\sum_{j=1}^{\frac n i} 1 + \sum_{i=1}^{\sqrt n}\sum_{j=1}^{\frac n i} 1 - \sum_{i=1}^{\sqrt n}\sum_{j=1}^{\sqrt n} 1 \\ = \ &(xy=n\ 在\ x\in[1,\sqrt n]\ 范围内的线下整点) + (同上，y\in [1,\sqrt n] \ 范围) \\ & - (两次求和的面积覆盖重复的部分) \\ = \ &2\sum_{i=1}^{\sqrt n}\lfloor\frac n i\rfloor - (\lfloor \sqrt n\rfloor)^2 \end{aligned}$$

问题转化为拟合 $xy = n$ 下方的整点数。咋拟合呢？想到

$4.3$ $\text{Stern-Brocot Tree}$ 和 $\text{Farey}$ 序列

《具体数学里面应该有。》
不不不我还会说几句的

定义 $\textbf{4.1}$（$\text{Farey}$ 序列）：将分母不超过 $n$ 的 $0\sim 1$ 的最简分数从小到大排成一个序列，我们称它为 $n$ 阶 $\text{Farey}$ 序列。如果存在一个 $\text{Farey}$ 序列使得 $p,q$ 连续出现，则我们称这两个数字为 $\text{Farey neighbour}$。

引理 $\textbf{4.1}$：任意 $\frac ab > \frac cd$ 是 $\text{Farey neighbour}$，当且仅当 $ad - bc = 1$。

证明：自己查。

在本文中，若 $\frac ab > \frac cd$ 是 $\text{Farey neighbour}$，则称 $\frac ba > \frac dc$ 也是 $\text{Farey neighbour}$。容易发现引理 $4.1$ 在此情况下仍然成立。

定义 $\textbf{4.2}$（$\text{Stern-Brocot Tree}$）：见此处。

$4.4$ 利用 $\text{Stern-Brocot Tree}$ 切割曲线

每次用一条直线切割，统计直线下方的整点数。初始化就是经过 $(\sqrt n, \sqrt n)$ 且斜率为 $-1$ 的直线。

得到一条直线后计算整点数可以直接用皮克定理，不展开。

$4.4.1$ 坐标系的转换

考虑计算 $S(x_0, y_0, a, b, c, d)$。由于曲线在 $PR$ 和 $PQ$ 上方，考虑令 $PR,PQ$ 分别为 $u,v$ 轴。某个方向 $+1$ 相当于是向这个方向走到下一个接触到的整点。

定义有 $(a,b) = (c,d) = 1$，因此对于 $(u,v)$，其在原坐标系下就是

$$x = x_0 + ud - vb$$

$$y = y_0 + uc - va$$

由上，有 $ax + by = ax_0 + by_0 + (ad - bc)u$。由于 $\frac ab$ 与 $\frac cd$ 是 $\text{Farey neighbour}$，有 $ad - bc = 1$。因此 $u$ 是整数。施于 $v$ 得到 $v$ 也是整数。

自然想到找 $xy = n$ 的对应。不妨令 $w_1 = ax_0 + by_0, w_2 = cx_0 + dy_0$，则 $x_0 = dw_1 - bw_2, y_0 = aw_2 - cw_1$。于是 $x = d(u + w_1) - b(v + w_2), y = a(v + w_2) - c(u + w_1)$。

曲线即为 $(u + w_1)(v + w_2) - cd(u + w1)^2 - ab(v + w_2)^2 = n$。由于 $w_1, w_2, cd, ab \ge 0$，该曲线上 $u > 0$ 与 $v > 0$ 一一对应。设 $u = U(v), v = V(u)$。我们能发现 $U,V$ 在第一象限内都是形似 $xy = n$ 的凹函数。可以求导证明。

$4.4.2$ 计算整点数

我们仍然找一条切凹函数且斜率为 $-1$ 的直线，统计它下方的整点数。我们设切点为 $G$，原系坐标是 $(G_x, G_y)$，新系坐标是 $(G_u, G_v)$。随后找到两个整点 $S,T$ 使得这两点在 $u$ 坐标上极小地夹住 $G$。形式化地，$S_u \le G_u < S_u + 1 = T_u$，$S_v = \left\lfloor V(S_u)\right\rfloor, T_v = \left\lfloor V(T_u)\right\rfloor$。这两点显然存在且不同。

过 $S,T$ 作切线的平行线 $SA, TB$ 分别交 $v$ 轴于 $A$、交 $u$ 轴与 $B$。并设 $M(A_x, A_y + 1)$，$C(B_x + 1, B_y)$，作 $MN$ 平行 $SA$ 交曲线于 $N$。然后把折线 $ASTB$ 左下方（被染上橙色的区域，包括边界）的整点统计入答案。这段图形整体上斜率为 $-1$，带换到 $xOy$ 坐标系下斜率为 $-\frac{a + c}{b + d}$，而通过 $\text{Stern-Brocot Tree}$ 二分可以得到 $\frac ab, \frac{a + c}{b + d}, \frac cd$ 是 $\text{Farey neighbour}$，这样可以接着切割了。

我们断言，任何未被统计的点，要么出现在 $MN$ 上方，要么出现在 $AC$ 右侧。因此可以递归 $S(B_x, B_y, a, b, a + c, b + d) + S(C_x, C_y, a + c, b + d, c, d)$ 求解。

为啥这里是 $(B_x, B_y)$ 不是 $(M_x, M_y)$？

$4.5$ 时间复杂度

据说这个算法可以被证到 $\mathcal O(n^{\frac 13}\log n)$，但是论文中只证到 $\mathcal O(n^{\frac 13}\log^3 n)$。但总而言之，现在能够证明的是该算法的复杂度为 $\mathcal{\tilde{O}}(n^{\frac 13})$。

一篇论文中给出了一种 $\mathcal O(n^{\frac 13})$ 的做法：A Successive Approximation Algorithm for Computing the Divisor Summatory Function, Richard Sladkey, 2012。

论文里代码的实现

不保证正确性，代码效率低下

namespace calc {
	int n;
	double C1 = 4000, C2 = 10;

	int delta(int i) { return i * (i + 1) / 2; }
	
	int SR(int w, int h, int a1, int b1, int c1, int a2, int b2, int c2) {
		int s = 0, a3 = a1 + a2, b3 = b1 + b2;
			
		auto H = [&](int u, int v) {
			return (b2 * (u + c1) - b1 * (v + c2)) * (a1 * (v + c2) - a2 * (u + c1));
		};
		auto U_tan = [&] { 
			return (int)(sqrtl( (a1 * b2 + b1 * a2 + 2 * a1 * b1) * (a1 * b2 + b1 * a2 + 2 * a1 * b1) * n / (a3 * b3) ) ) - c1;
		};
		auto V_floor = [&](int u) {
			return (int)(((a1 * b2 + a2 * b1) * (u + c1) - ceil(sqrtl((u + c1) * (u + c1) - 4 * a1 * b1 * n))) / (2 * a1 * b1) - c2);
		};
		auto U_floor = [&](int v) {
			return (int)(((a1 * b2 + a2 * b1) * (v + c2) - ceil(sqrtl((v + c2) * (v + c2) - 4 * a2 * b2 * n))) / (2 * a2 * b2) - c2);
		};

		auto SI = [&](int i_max, int p1, int p2, int q, int r) {
			int s = 0, A = p1 * p1 - 2 * r * n, B = p1 * q, C = 2 * p1 - 1;
			rep(i, 1, i_max - 1) {
				C += 2, A += C, B += q;
				s += (B - floor(sqrtl(A))) / r;
			} return s - (i_max - 1) * p2;
		};
		auto SW = [&]() {
			return SI(w, c1, c2, a1 * b2 + a2 * b1, 2 * a1 * b1);
		};
		auto SH = [&]() {
			return SI(h, c2, c1, a1 * b2 + a2 * b1, 2 * a2 * b2);
		};

		if (h > 0 and H(w, 1) <= n) {
			s = s + w, c2 ++, h --;
		} 
		if (w > 0 and H(1, h) <= n) {
			s = s + h, c1 ++, w --;
		}
		if (w <= C2) return s + SW();
		if (h <= C2) return s + SH();
		int u4 = U_tan(), v4 = V_floor(u4), u5 = u4 + 1, v5 = V_floor(u5);
		int v6 = u4 + v4, u7 = u5 + v5;

		auto SN = [&]() {
			return delta(v6 - 1) - delta(v6 - u5) + delta(u7 - u5);
		};
		
		s += SN();
		s += SR(u4, h - v6, a1, b1, c1, a3, b3, c1 + c2 + v6);
		s += SR(w - u7, v5, a3, b3, c1 + c2 + u7, a2, b2, c2);
		return s;
	}

	int calc() {
		int x_min, x_max, y_min;
		x_max = floor(sqrtl(n));
		y_min = n / x_max;
		x_min = min(x_max, ceil(C1 * pow(2 * n, 1. / 3)));
		lll s = 0;
		int a2 = 1, x2 = x_max, y2 = y_min, c2 = a2 * x2 + y2;

		while (true) {
			int a1 = a2 + 1;
			int x4 = floor(sqrtl(n / a1)), y4 = n / x4, c4 = a1 * x4 + y4;
			int x5 = x4 + 1, y5 = n / x5, c5 = a1 * x5 + y5;
			if (x4 < x_min) break;

			auto SM = [&]() {
				return delta(c4 - c2 - x_min) - delta(c4 - c2 - x5) + delta(c5 - c2 - x5);
			};
			s = s + SM();
			s = s + SR(a1 * x2 + y2 - c5, a2 * x5 + y5 - c2, a1, 1, c5, a2, 1, c2);
			a2 = a1, x2 = x4, y2 = y4, c2 = c4;
		} 

		auto SQ = [&](int x1, int x2) {
			int s = 0, x = x2, bet = n / (x + 1), eps = n % (x + 1), delt = n / x - bet, gam = bet - x * delt;
			while (x >= x1) {
				eps += gam;
				if (eps >= x) {
					delt ++, gam -= x, eps -= x;
					if (eps >= x) {
						delt ++, gam -= x, eps -= x;
						if (eps >= x) break;
					}
				} else if (eps < 0) {
					delt --, gam += x, eps += x;
				} 
				gam += 2 * delt, bet += delt, s += bet, x --;
			} 
			eps = n % (x + 1), delt = n / x - bet, gam = bet - x * delt;
			while (x >= x1) {
				eps = eps + gam;
				int delt2 = eps / x;
				delt += delt2, eps -= x * delt2;
				gam = gam + 2 * delt - x * delt2, bet += delt, s += bet, x --;
			} 
			while (x >= x1) {
				s += n / x, x --;
			} 
			return s;
		};

		auto S1 = [&]() {
			return SQ(1, x_min - 1);
		};
		auto S2 = [&]() { 
			return (x_max - x_min + 1) * y_min + delta(x_max - x_min); 
		};
		auto S3 = [&]() {
			int ret = 0;
			rep(x, x_min, x2 - 1) 
				ret += (n / x) - (a2 * (x2 - x) + y2);
			return ret;
		};

		s = s + S1() + S2() + S3();
		return 2 * s - x_max * x_max;
	}
}

int paper(int n) {
	calc::n = n;
	return calc::calc();
}

$4.6$ 优化

我们发现，我们其实没有必要费时费力地切割原面积。其实问题可以转化为从 $(\sqrt n,\sqrt n)$ 这个点开始拟合 $xy = n$ 的下半部分。我们需要拟合的面积是下图中红色虚线下方的阴影部分：

假设我们现在拟合到了 $(x, y)$，我们要选择一个不在阴影内且在当前点右侧的点 $(p,q)$，满足 $\frac {q - y}{p - x}$ 最大且 $p$ 最小的点，并将 $(x,y)$ 到 $(p,q)$ 的线段加入折线，在这过程中加入折线下整点数。
我们可以维护一棵 $\frac 01, \frac 11$ 为初始值的 $\text{Stern-Brocot Tree}$。每次拿出栈顶分数 $\frac pq$，从 $(x,y)$ 移动到 $(x + a, y - b)$（我们也可以将分数看作一个向量）。就这么一直移动直到再走一步会进入阴影中。然后把栈首弹出。接着对新的栈首执行上述操作直到栈首 $L$ 不行，第二个元素 $R$ 行的时候。
这时我们将栈首弹出，开始二分最接近斜率。
我们在 $\text{Stern-Brocot Tree}$ 上走，设当前 $\text{mid} = (L_x + R_x, L_y + R_y)$，考虑当 $(x,y)$ 加入 $\text{mid}$ 时是否会进入阴影。
如果不会进入则 $\text{mid}$ 入栈，$R = \text{mid}$，接着二分。如果会进入，讨论一下。如果加入后斜率 $\le R$，那再怎么加都没用了，可以直接退出了。反之就令 $L= \text{mid}$。
栈的用处就是回溯到第一个合法位置。

zzt 的实现。

然后复杂度和先前那个做法的证法类似，是 $\mathcal O(n^{\frac 13} \log n)$。
$\text{Min_25}$ 在 SPOJ 说他无法给出一个证明，灵感来源 PE540 的讨论区也无法给出。
但是他推测这个曲线的凸包上有 $\Theta(n^{\frac 13}\log n)$ 种斜率（每个形如 $\left[\frac{\sqrt n}{2^{k+1}}, \frac{\sqrt n}{2^{k}} \right]$ 的区间上有 $\mathcal O(n)$ 种。
所以他推测这种算法的时间复杂度有界 $\Omega(n^{\frac 13}\log^3 n)$。
为了保证复杂度，当前点的 $y< n^{\frac 13}$ 时停止迭代，开始朴素做法。

$4.7$ 拓展

这里除了凸性外没有用到 $xy = n$ 的任何性质，也就是说这个做法可以自然推广到其他奇怪的凹曲线内部整点计数。