《一些特殊的数论函数求和问题》阅读随笔
好至少它教会了我如何把质数求和转化成积分的渐进 对着 \(\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，复杂度渐进地是
\(3\to 4\) 考虑首先将分母变成渐进的 \(\log n\)，然后只需要对分子积分，而这个处理是平凡的。可以写出
然后考虑统计答案部分。
定义 \(\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)\)
证明：
进行初等积分。
其中 \(\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\) 的贡献。奇妙结论：
这式子可以建出递归对应的二叉树得到。
换元。令 \(p = small(i), m = i / p\)。\(\mu(i) = -\mu(m)\)，也就导出了
阈值法。
\(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\)。
改变求和式为
\(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\)。
改变求和式为
我们可以枚举 \(p\)，求 \(m\) 的贡献。筛出质数 \(k\) 次方和计算。复杂度 \(\mathcal O(B)\)。
\(3.2.3\) \(\sqrt B < p \le n^{\frac 13}\)
直接拿定义冲
于是问题来到了维护这边。
枚举 \(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\) 时同样成立。
进行初等积分：
感想：积分是用来分析的，不是用来创人的。
\(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)\)。
证明：
进行初等积分：
这个挺好理解。
然后带入 \(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}\) 的查询部分。
于是问题又来到了维护这边。
抄式子：
我们进一步利用那个定理。分别考虑 \(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\) 朴素转化
考虑转化
问题转化为拟合 \(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)\)，其在原坐标系下就是
由上，有 \(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}\)。
栈的用处就是回溯到第一个合法位置。
然后复杂度和先前那个做法的证法类似，是 \(\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\) 的任何性质，也就是说这个做法可以自然推广到其他奇怪的凹曲线内部整点计数。
以下是博客签名，与正文无关。
请按如下方式引用此页：
本文作者 joke3579，原文链接：https://www.cnblogs.com/joke3579/p/paperessay221204.html。
遵循 CC BY-NC-SA 4.0 协议。
请读者尽量不要在评论区发布与博客内文完全无关的评论，视情况可能删除。
