好至少它教會了我如何把質數求和轉化成積分的漸進 對著 \(\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^{\frac 34} }{\log n } \right)\)。
如果直接暴力計算呢？在小資料的測試下，該演算法似乎比 \(\mathcal O \left( \frac { n^{\frac 34} }{\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(n)}{n^2} \\ = \ & \frac{\pi(n)}{n^2}\Big\lvert_n^{\infty} - \int_n^{\infty} \pi(n) d\left(\frac{1}{n^2}\right) \\ = \ & \Theta\left( \frac{1}{n\ln n}\bigg\lvert_n^{\infty} + 2\int_n^{\infty} \frac{\text dn}{n^2 \ln n} \right) \\ = \ & \Theta\left( -\frac{1}{n\ln n} + 2\int_{\frac 1n}^{0} \frac{n^2\text d\frac 1n}{ \ln \frac 1n} \right) \\ = \ & \Theta\left( -\frac{1}{n\ln n} - 2\int_{0}^{\frac 1n} \frac{n^2 (-n^{-2})\text dn}{ -\ln n} \right) \\ = \ & \Theta\left( -\frac{1}{n\ln n} - 2\int_{0}^{\frac 1n} \frac{\text dn}{ \ln n} \right) \\ = \ & \Theta\left( -\frac{1}{n\ln n} - 2\ \text{li}\ \frac 1n \right) \end{aligned}\]

其中 \(\text{li}\) 為對數積分函式。由經典結論可知 \(\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 xB\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 xB\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 xp})\) 的。

於是問題來到了複雜度這邊。

首先是修改。

定理 \(\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\) 前的負號提了出去。

歸納成立。

我們在樹狀陣列上進行操作的時候需要保證 \(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})\)。

實現細節可以參照這裡：Richard Sladkey, A Successive Approximation Algorithm for Computing the Divisor Summatory Function, 2012。

\(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\) 的任何性質，也就是說這個做法可以自然推廣到其他奇怪的凹曲線內部整點計數。