不超過 x 的素數個數

\(\pi(x)\) 表示不超過 \(x\) 的素數個數。容易看出可以在 \(O(N)\) 時間復雜度，\(O(N)\) 空間復雜度離線預處理求出小於 \(N\) 的素數全體。但是如果 \(N=10^{14}\) 或者更大，這種做法必然是不現實的。因此下面給出高效的求解方法...

理論基礎: 參考潘承洞《數論基礎》以及論文包.zip

\(\psi(x,s)\)

\(\psi(x,s)\) 表示不超過 \(x\) 且能不能被前 \(s\) 個素數整除的正整數個數。即

\[\psi(x,s) = \sum_{n \leq x} \sum_{d|(n,m_s)} \mu(d) = \sum_{d|m_s} u(d)\lfloor \frac{x}{d} \rfloor \]

其中 \(m_s = p_1 \cdots p_s\) 為前 \(s\) 個素數的積。

另一方面，顯然我們有

\[\psi(x,s) = \psi(x,s-1) - \psi(\frac{x}{p_s},s-1) \]

\(\pi(x)\)

我們知道一個數 \(n>1\) 是素數當且僅當不存在素數 \(p \leq \sqrt{n}\) 使得 \(p \mid n\)。因此當 \(s \geq \pi(\sqrt{x})\) 時，

\[\psi(x,s) = \pi(x) - s + 1 \]

\(P_k(x,s)\)

設 \(P_k(x,s)\) 為 不超過 \(x\) 且每個素因子都大於 \(p_s\) 且素因子(按重根計)個數為 \(k\) 的整數個數（方法屬於 Lehmer）。
進一步設 \(P_0(x,s)=1\)。則

\[\psi(x,s) = \sum_{k=0} ^{\infty} P_k(x,s) \]

顯然 \(P_1(x,s) = \pi(x)-s\)。

若 \(\pi(\sqrt[3]{x}) \leq s \leq \pi(\sqrt{x})\) 則 \(P_k(x,s)=0,k \geq 3\) 此時

\[\psi(x,s) = 1 + \pi(x)-s + P_2(x,s) \]

其中

\[P_2(x,s) = \sum_{k=s+1}^{\pi(\sqrt{x})} \left( \pi(\frac{x}{p_k}) - k + 1 \right) \]

注意到上式中 \(\frac{x}{p_k} < x^{\frac{2}{3}}\)

\(\pi(x)\) 的計算公式

\[\pi(x) = \psi(x,s) + \frac{(\pi(\sqrt{x})+s-2)(\pi(\sqrt{x})- s+1)}{2} - \sum_{k=s+1}^{\pi(\sqrt{x})} \pi(\frac{x}{p_k}) \]

取上面 $s = \pi(\sqrt[3]{x}) $
因此問題最終轉化成求 \(\psi(x,\pi(\sqrt[3]{x}))\)。它可以利用

1. \(\psi(x,0) = \lfloor x \rfloor\)
2. \(\psi(x,s) = \psi(x,s-1) - \psi(\frac{x}{p_s},s-1)\)

至此問題貌似就這么解決了。但是由於這個遞歸會使得程序效率大大降低，因此需要一些預處理操作。

1. 若 \(x<p_s\) 則 \(\psi(x,s) ＝ 1\)
2. 給定一個小整數 M，預處理出 \(\psi(x,s)\),其中 \(x < q=p_1 \cdots p_s,\quad s<=M\)
   則 \(\psi(x,s) = \psi(x \mod q,s) + \lfloor \frac{x}{q} \rfloor \psi(q,s)\)

lehmer 計算公式

我自己寫的代碼沒有上面的快，兩種計算各有優勢

令 \(s = \pi(\sqrt[4]{x}), t= \pi(\sqrt[3]{x})\)。則，對任意 \(i>3, P_i(x,s) = 0\),

\[\begin{array}{rl} \psi(x,s) &= 1 + \pi(x) - s + P_2(x,s) + P_3(x,s) \\ &= 1+ \pi(x) - s + P_2(x,s) + \sum_{k=s+1}^{t} P_2(\frac{x}{p_k},k-1) \\ \end{array} \]

即：

\[\pi(x) = \psi(x,s)-1+s-P_2(x,s) - \sum_{k=s+1}^{t} P_2(\frac{x}{p_k},k-1) \]

注意到 \(\frac{x}{p_k} < \sqrt{x}\) ，所以最后一個式子可以用下式求，最后計算復雜度在於 \(P_2(x,s)\)

\[\sum_{k=s+1}^{t} P_2(\frac{x}{p_k},k-1) = \sum_{k=s+1}^{t} \sum_{j=k}^{\pi(\sqrt{\frac{x}{p_k}})} \pi(\frac{x}{p_k p_j}) - j+1 \]

穩定簡潔的 DP 做法

我們令 \(dp(x,s) = \psi(x,s)+s-1\) 它的意義是，\(2～x\) 中被前 \(s\) 個素數篩完后的偽素數個數。因此我們有 \(dp(0,0)=0,dp(x,0)=x-1,x>1,dp(x,\pi(\sqrt{x})) = \pi(x)\) 且有狀態轉移

\[dp(x,s) = dp(x,s-1)-dp(\frac{x}{p_s},s-1)+s-1 \]

因為 \(dp(p_{s-1},s-1) = s-1\)，最后一項可以寫成 \(dp(p_{s-1},s-1)\)。雖然上面需要二維數組，但是實際上我們可以優化成一維數組的情況。因為

\[dp(x,s) = dp(x,s-1)-dp(\frac{x}{p_s},s-1)+ dp(p_{s-1},s-1) \]

另外我們也不可能開 \(O(n)\) 的數組，但是可以利用一種黑科技開 \(O(\sqrt{n})\) 的數組即可達到我們的目的。
即我們用 \(L[x]\) 表示 \(dp(x,s)\) 用 \(R[x]\) 表示 \(dp(\frac{n}{x},s)\)。

若\(L[m]!=L[m-1]\) ，則說明 \(m\) 不能被前 \(s\) 個素數整除是第 \(s+1\) 個素數。

我們需要的是 \(R[1]\)

上述做法的時間復雜度為 \(O(\frac{n}{\log n})\) 且常數特別小，代碼十分簡潔。

這個騷方法還目前還沒有找到其它的應用 0.0

主要還沒法對它用樹狀數組

求第 n 個素數的方法

• 根據概率分布找到大致界限
• 再牛頓梯度法（或者二分查找）使得 \(\pi(m)= n\)
• 素數判斷遞減 \(m\) 直到 \(m\) 為素數
• 參考這里

\(\psi(x,s)\) 計算太慢了，需要優化

我們知道，若 \(\pi(\sqrt[3]{x}) \leq s \leq \pi(\sqrt{x})\) 則

\[\psi(x,s) = 1 + \pi(x)-s + P_2(x,s) \]

其中

\[P_2(x,s) = \sum_{k=s+1}^{\pi(\sqrt{x})} \left( \pi(\frac{x}{p_k}) - k + 1 \right) = \sum_{k=s+1}^{\pi(\sqrt{x})} \left( \psi(\frac{x}{p_k},s)+ s - k \right) \]

注意到上式中 $\frac{x}{p_k} < x^{\frac{2}{3}} $

我們之前的操作本質是遞歸讓 \(x,s\) 變小，通過打表預處理來加速遞歸使得滿足一定的效率需要。

現在我們來直接計算得到我們的答案。

符號約定

取定整數 \(\sqrt[3]{x} \leq y = x^{\frac{1}{3}+\epsilon }< \sqrt{x}, z = \frac{x}{y}, s = \pi(y)\) ，約定 \(p, q\) 是素數。 預處理 \(y\) 以內的數組：isp[], p[], mu[], pi[] ，對 \([1,z]\) 內的 \(\psi(x,s)\) 用樹狀數組（可以在我的博客網站搜索：樹狀數組）維護（注意到這需要 \(O(z)\) 的內存，單次維護不現實，所以我們可以一段段的維護，保證每一段的長度為 \(2^{\lfloor \log_2(y) \rfloor +1}\) 來提高效率）

做完上面的預處理后，我們發現 \(P_2(x,s)\) 可以直接計算了。

\(\psi(x,s)\) 直接計算

在本博文的最開始有：

\[\psi(x,s) = \sum_{n \leq x} \sum_{d|(n,p)} \mu(d) = \sum_{d|p} u(d)\lfloor \frac{x}{d} \rfloor \]

其中 \(p\) 為前 \(s\) 個素數的積。

但是最右邊本質上有很多項，所以直接算，其實復雜度特別高。

我們還有遞歸公式：

\[\begin{array}{cl} \psi(x,s) &= \psi(x,s-1) - \psi(\frac{x}{p_s},s-1) \\ &= \psi(x,s-2) - \psi(\frac{x}{p_{s-1}},s-2) - \psi(\frac{x}{p_s},s-2) + \psi(\frac{x}{p_sp_{s-1}},s-2) \end{array} \]

可以一直分解下去，如果一直分解下去就可以得到最上面的公式了。

所以我們約定：對於節點 \(\mu(n) \psi(\frac{x}{n},b)\) ，如果滿足

• 原始節點： \(b = c, n \leq y\)
• 特殊節點：\(n > y\)

就不再分解。這也等價於說 如果 \(n < y\) 且 \(b>c\) 就分解。因為一開始 \(n=1,b=a>c\)。而且 \(n\) 會增大，\(b\) 會減小，所以節點一定會有限步內，落入上述兩種框架中的一種。並且 特殊節點的父節點 \(-\mu(n) \psi(\frac{x}{\frac{n}{p_{d+1}}},b+1)\) 必然滿足 \(\frac{n}{ p_{d+1} } \leq y <n\) 且 \(b+1>c\)。綜上：

以前設置 \(c = 7\)，但是后來發現沒必要，\(c=0\) 就挺好。

\[\psi(x,s) = \sum_{n=1} ^y \mu(n) \lfloor \frac{x}{n} \rfloor + \sum_{\frac{n}{\delta(n)} \leq y < n} \mu(n) \psi(\frac{x}{n}, \pi(\delta(n))-1) = S_0 + S \]

\(S_0\) 很好處理，計算 \(S\)： 對 \(p = \delta(n)\) 一起求：

\[S = - \sum_{p \leq y} \sum_{\begin{array}{c} \delta(m) >p \\ m \leq y < mp \end{array}} \mu(m) \psi(\frac{x}{mp},\pi(p)-1) \]

注意到：\(\frac{x}{mp} < z\) ，所以我們已經可以把 \(\psi(x,s)\) 直接計算出來了。

但是我們可以避免很多計算來提高效率。於是我們有下列一系列的操作

• \(p \geq \sqrt{y}\)，則 \(m\) 為素數 ，且此時 \(mp > p^2 \geq y\) (若 \(m\) 為合數，則 \(m \geq \delta(m) ^2 >p^2 \geq y\) 矛盾)
• \(\frac{x}{mp} < p\) 時，\(\psi(\frac{x}{mp},\pi(p)-1) = 1\)
• \(\sqrt{\frac{x}{mp}} < p\) 時，\(\psi(\frac{x}{mp},\pi(p)-1) = \pi(\frac{x}{mp}) - \pi(p) +2\)
• \(\sqrt{y} < \sqrt{z} < x^{\frac{1}{4}} < x^{\frac{1}{3}} < y\)

由此對\(S\)分段計算

\[S_1 = - \sum_{\sqrt[3]{x} < p \leq y} \sum_{\begin{array}{c} \delta(m) >p \\ m \leq y < mp \end{array}} \mu(m) \psi(\frac{x}{mp},\pi(p)-1) \]

\[S_2 = - \sum_{\sqrt{z} < p \leq \sqrt[3]{x}} \sum_{\begin{array}{c} \delta(m) >p \\ m \leq y < mp \end{array}} \mu(m) \psi(\frac{x}{mp},\pi(p)-1) \]

\[S_3 = - \sum_{\sqrt{y} < p \leq \sqrt{z}} \sum_{\begin{array}{c} \delta(m) >p \\ m \leq y < mp \end{array}} \mu(m) \psi(\frac{x}{mp},\pi(p)-1) \]

\[S_4 = - \sum_{p \leq \sqrt{y}} \sum_{\begin{array}{c} \delta(m) >p \\ m \leq y < mp \end{array}} \mu(m) \psi(\frac{x}{mp},\pi(p)-1) \]

由限制關系式，我們化簡 \(S_1, S_2, S_3\)

\[\begin{array}{cl} S_1 &= - \sum_{\sqrt[3]{x} < p \leq y} \sum_{p<q \leq y} \mu(q) \psi(\frac{x}{pq},\pi(p)-1) \\ &= \sum_{\sqrt[3]{x} < p < q \leq y} 1 \\ &= {\pi(y)-\pi(\sqrt[3]{x}) \choose 2} \end{array} \]

\[\begin{array}{cl} S_2 &= - \sum_{\sqrt{z} < p \leq \sqrt[3]{x}} \sum_{p<q \leq y} \mu(q) \psi(\frac{x}{pq},\pi(p)-1) \\ &= \sum_{\sqrt{z} < p \leq \sqrt[3]{x}} \sum_{p<q \leq y} \pi(\frac{x}{pq}) - \pi(p) +2 \\ &= \sum_{\sqrt{z} < p \leq \sqrt[3]{x}} \sum_{p<q \leq y} \pi(\frac{x}{pq}) + \sum_{\sqrt{z} < p \leq \sqrt[3]{x}} (\pi(p)-\pi(y))(\pi(p)-2) \end{array} \]

\[\begin{array}{cl} S_3 &= - \sum_{\sqrt{y} < p \leq \sqrt{z}} \sum_{p < q \leq y} \mu(q) \psi(\frac{x}{pq},\pi(p)-1) \\ &= \sum_{\sqrt{y} < p \leq \sqrt{z}} \sum_{p < q \leq y} \psi(\frac{x}{pq},\pi(p)-1) \end{array} \]

\(S_2\) 也可以用 \(S_3\) 的式子求，只是效率不高。

\(S_2\) 中 \(\frac{x}{pq}, p<y\)，即可以直接求。

當然了還可以繼續細化，但是我嫌麻煩就不想再細化了！

也就是現在的核心就是樹狀數組分段維護數據，然后每一段的總值要用數組存起來就好了。然后用這個數據結果計算 \(S_2,S_3,S_4,P_2(x,s)\)，然后就大功告成了 0.0

用 \(\psi(x,s)\) 計算 \(\pi(x)\)，還是用 \(\pi(x)\) 計算 \(\psi(x,s)\) ，這是一個問題。

用樹狀數組維護的時候會有一個很大的問題就是：求和式中 每此動 \(p\) 整個維護就要從 \(1 \to p\) 重新維護一次很麻煩。這個問題沒解決所以我不想寫代碼。

想把 \(\frac{x}{pq}\) 的所有可能的值單調遞增排列但是又不現實。

第 n 個素數

做法依賴於 \(\pi(x)\) 的計算

素數定理( 這里就不證了)

\[\lim _{x \to \infty} \frac{\pi(x)}{x/\ln x} = 1 \]

從而我們知道：

\[\lim _{x \to \infty} \frac{p_n}{n \ln n} = 1 \]

其中，\(p_n\) 為第 \(n\) 個素數，顯然 \(p_n\) 是 \(\pi(x) = n\) 最小的解。

\(p_n\) 求解

• 預處理小於 \(N\) 的素數
• 初始值 \(n\ln n\)
• 牛頓梯度法

第 \(n\) 個素數和 \(\pi(x)\) 的網站

世界紀錄保持者的求法

\(\sum_{p \leq n, p \text{ is prime}} p\)

我們可以在不求出 不超過 \(x\) 的所有素數 的情況下，求出最終結果。

\[f_k(x) \doteq \sum_{p \leq x} p^k \]

• \(f_0(x) = \pi(x)\)
• \(f_1(x) = \sum_{p \leq x} p\) 這是我們關心的結果
• 對於一般的 \(k\) 借助自然數方冪和快速算法也可以求

\(f_k(x,s)\)：最小素因子大於 \(p_s\) 且不超過 \(x\) 的數的 \(k\) 次方和

\[f_k(x,s) = \sum_{m \leq x, \delta(m)>p_s} m^k \]

其中，\(\delta(m)\) 表示 \(m\) 的最小素因子（約定 \(\delta(1) = + \infty\)）。

\(f_k(x,s)\) 的遞推公式

\[\begin{aligned} f_k(x,s) &= \sum_{m \leq x, \delta(m)>p_s} m^k \\ &= f(x,s-1) - \sum_{m \leq x, \delta(m) = p_s} m^k \\ &= f(x,s-1) - p_s ^k f(\frac{x}{p_s},s-1) \end{aligned} \]

\(\displaystyle f_k(x,0) = \sum_{i=1} ^ {\lfloor x \rfloor} i^k\)

\(f_k(x,s)\) 和 \(f_k(x)\) 的關系

若 \(s >= \pi(\sqrt{x})\)，則

\[f_k(x) = f_k(x,s) - 1 + f_k(p_s) \]

其實我們還可以，類似於 \(\pi(x)\) 的計算一樣, 取 \(\pi(\sqrt[3]{x}) \leq s \leq \pi(\sqrt{x})\)

\[f_k(x) = f_k(x,s) - 1 + f_k(p_s) - \sum_{i=s+1} ^ {\pi(\sqrt{x})} (f(\frac{x}{p_i})-f(p_{i-1})) p_i ^k \]