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

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

18.3k……阿巴阿巴

引言

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

\(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\) 的任何性质，也就是说这个做法可以自然推广到其他奇怪的凹曲线内部整点计数。

posted @ 2022-12-04 17:08 joke3579 阅读(335) 评论(5) 编辑收藏举报

刷新页面返回顶部