「Meissel-Lehmer 演算法」是一種能在亞線性時間複雜度內求出 \(1\sim n\) 內質數個數的一種演算法。
在看素數相關論文時發現了這個演算法，論文連結：Here。
演算法的細節來自 OI wiki，轉載僅作為學習使用。
目前先 mark 一下這個演算法，等有空的時候再來研究一下，演算法的時間複雜度為 \(\mathcal{O}(n^{\frac23})\) ，所以 \(n\) 的範圍可以擴大至 \(10^{12}\) 的級別；
程式碼實現
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
//通過知道前面的 n^1/3 的質數可以推斷後面n^2/3的質數所以可以適當減小
const int N = 9e3;
const int M = 2; //為了減小記憶體可以不過是質數
const int PM = 2 * 3 * 5; //為了減小記憶體可以不過要按質數減小如去掉17
ll n;
bool np[N];
int prime[N], pi[N];
int phi[PM + 1][M + 1], sz[M + 1];
int getprime() {
int cnt = 0;
np[0] = np[1] = true;
pi[0] = pi[1] = 0;
for (int i = 2; i < N; ++i) {
if (!np[i]) prime[++cnt] = i;
pi[i] = cnt;
for (int j = 1; j <= cnt && i * prime[j] < N; ++j) {
np[i * prime[j]] = true;
if (i % prime[j] == 0) break;
}
}
return cnt;
}
void init() {
getprime();
sz[0] = 1;
for (int i = 0; i <= PM; ++i) phi[i][0] = i;
for (int i = 1; i <= M; ++i) {
sz[i] = prime[i] * sz[i - 1];
for (int j = 1; j <= PM; ++j) phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1];
}
}
int sqrt2(ll x) {
ll r = (ll)sqrt(x - 0.1);
while (r * r <= x) ++r;
return int(r - 1);
}
int sqrt3(ll x) {
ll r = (ll)cbrt(x - 0.1);
while (r * r * r <= x) ++r;
return int(r - 1);
}
ll getphi(ll x, int s) {
if (s == 0) return x;
if (s <= M) return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];
if (x <= prime[s] * prime[s]) return pi[x] - s + 1;
if (x <= prime[s] * prime[s] * prime[s] && x < N) {
int s2x = pi[sqrt2(x)];
ll ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;
for (int i = s + 1; i <= s2x; ++i) ans += pi[x / prime[i]];
return ans;
}
return getphi(x, s - 1) - getphi(x / prime[s], s - 1);
}
ll getpi(ll x) {
if (x < N) return pi[x];
ll ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;
for (int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i) ans -= getpi(x / prime[i]) - i + 1;
return ans;
}
ll lehmer_pi(ll x) { //小於等於n的素數有多少個
if (x < N) return pi[x];
int a = (int)lehmer_pi(sqrt2(sqrt2(x)));
int b = (int)lehmer_pi(sqrt2(x));
int c = (int)lehmer_pi(sqrt3(x));
ll sum = getphi(x, a) + (ll)(b + a - 2) * (b - a + 1) / 2;
for (int i = a + 1; i <= b; i++) {
ll w = x / prime[i];
sum -= lehmer_pi(w);
if (i > c) continue;
ll lim = lehmer_pi(sqrt2(w));
for (int j = i; j <= lim; j++) sum -= lehmer_pi(w / prime[j]) - (j - 1);
}
return sum;
}
int main() {
ios_base::sync_with_stdio(false), cin.tie(0);
init();
while (cin >> n && n) cout << lehmer_pi(n) << "\n";
return 0;
}
OJ 上的幾道相關的題：Here
記號規定
\(\left[x\right]\) 表示對 \(x\) 下取整得到的結果。
\(p_k\) 表示第 \(k\) 個質數，\(p_1=2\)。
\(\pi\left(x\right)\) 表示 \(1\sim x\) 範圍內素數的個數。
\(\mu\left(x\right)\) 表示莫比烏斯函式。
對於集合 \(S\)，\(\# S\) 表示集合 \(S\) 的大小。
\(\delta\left(x\right)\) 表示 \(x\) 最小的質因子。
\(P^+\left(n\right)\) 表示 \(x\) 最大的質因子。
Meissel-Lehmer 演算法求 π(x)
定義 \(\phi\left(x,a\right)\) 為所有小於 \(x\) 的正整數中滿足其所有質因子都大於 \(p_a\) 的數的個數，即：
再定義 \(P_k\left(x,a\right)\) 表示為所有小於 \(x\) 的正整數中滿足可重質因子恰好有 \(k\) 個且所有質因子都大於 \(p_a\) 的數的個數，即：
特殊的，我們定義：\(P_0\left(x,a\right)=1\)，如此便有：
這個無限和式實際上是可以表示為有限和式的，因為在 \(p_a^k>x\) 時，有 \(P_k\left(x,a\right)=0\)。
設 \(y\) 為滿足 \(x^{1/3}\le y\le x^{1/2}\) 的整數，再記 \(a=\pi\left(y\right)\)。
在 \(k\ge 3\) 時，有 \(P_1\left(x,a\right)=\pi\left(x\right)-a\) 與 \(P_k\left(x,a\right)=0\)，由此我們可以推出：
這樣，計算 \(\pi\left(x\right)\) 便可以轉化為計算 \(\phi\left(x,a\right)\) 與 \(P_2\left(x,a\right)\)。
計算 P₂(x,a)
由等式 \(\left(2\right)\) 我們可以得出 \(P_2\left(x,a\right)\) 等於滿足 \(y<p\le q\) 且 \(pq\le x\) 的質數對 \(\left(p,q\right)\) 的個數。
首先我們注意到 \(p\in \left[y+1,\sqrt{x}\right]\)。此外，對於每個 \(p\)，我們都有 \(q\in\left[p,x/p\right]\)。因此：
當 \(p\in \left[y+1,\sqrt{x}\right]\) 時，我們有 \(\dfrac{x}{p}\in \left[1,\dfrac{x}{y}\right]\)。因此，我們可以篩區間 \(\left[1,\dfrac{x}{y}\right]\)，然後對於所有的的質數 \(p\in \left[y+1,\sqrt{x}\right]\) 計算 \(\pi\left(\dfrac{x}{p}\right)-\pi\left(p\right)+1\)。為了減少上述演算法的空間複雜度，我們可以考慮分塊，塊長為 \(L\)。若塊長 \(L=y\)，則我們可以在 \(O\left(\dfrac{x}{y}\log{\log{x}}\right)\) 的時間複雜度，\(O\left(y\right)\) 的空間複雜度內計算 \(P_2\left(x,a\right)\)。
計算 ϕ(x,a)
對於 \(b\le a\)，考慮所有不超過 \(x\) 的正整數，滿足它的所有質因子都大於 \(p_{b-1}\)。這些數可以被分為兩類：
- 可以被 \(p_b\) 整除的；
- 不可以被 \(p_b\) 整除的。
屬於第 \(1\) 類的數有 \(\phi\left(\dfrac{x}{p_b},b-1\right)\) 個，屬於第二類的數有 \(\phi\left(x,b\right)\) 個。
因此我們得出結論：
定理 \(5.1\)： 函式 \(\phi\) 滿足下列性質
\[\phi\left(u,0\right)=\left[u\right]\tag{5} \]\[\phi\left(x,b\right)=\phi\left(x,b-1\right)-\phi\left(\dfrac{x}{p_b},b-1\right)\tag{6} \]
計算 \(\phi\left(x,a\right)\) 的簡單方法可以從這個定理推匯出來：我們重複使用等式 \(\left(7\right)\)，知道最後得到 \(\phi\left(u,0\right)\)。這個過程可以看作從根節點 \(\phi\left(x,a\right)\) 開始建立有根二叉樹，圖 \(1\) 畫出了這一過程。通過這種方法，我們得到如下公式：
上圖表示計算 \(\phi\left(x,a\right)\) 過程的二叉樹：葉子節點權值之和就是 \(\phi\left(x,a\right)\)。
但是，這樣需要計算太多東西。因為 \(y\geq x^{1/3}\)，僅僅計算為 \(3\) 個 不超過 \(y\) 質數的乘積的數，如果按照這個方法計算，會有至少 \(\dfrac{x}{\log^3 x}\) 個項，沒有辦法我們對複雜度的需求。
為了限制這個二叉樹的“生長”，我們要改變原來的終止條件。這是原來的終止條件。
終止條件 \(1\)： 如果 \(b=0\)，則不要再對節點 \(\mu\left(n\right)\phi\left(\dfrac xn,b\right)\) 呼叫等式 \(\left(6\right)\)。
我們把它改成更強的終止條件：
終止條件 \(2\)： 如果滿足下面 \(2\) 個條件中的一個，不要再對節點 \(\mu\left(n\right)\phi\left(\dfrac xn,b\right)\) 呼叫等式 \(\left(6\right)\):
- \(b=0\) 且 \(n\le y\)；
- \(n>y\)。
我們根據 終止條件 \(2\) 將原二叉樹上的葉子分成兩種：
- 如果葉子節點 \(\mu\left(n\right)\phi\left(\dfrac xn,b\right)\) 滿足 \(n\le y\)，則稱這種葉子節點為 普通葉子；
- 如果葉子節點 \(\mu\left(n\right)\phi\left(\dfrac xn,b\right)\) 滿足 \(n>y\) 且 \(n=mp_b\left(m\le y\right)\)，則稱這種節點為 特殊葉子。
由此我們得出：
定理 \(5.2\)： 我們有：
\[\phi\left(x,a\right)=S_0+S\tag{7} \]其中 \(S_0\) 表示 普通葉子 的貢獻：
\[S_0=\sum\limits_{n\le y}{\mu\left(n\right)\left[\dfrac xn\right]}\tag{8} \]\(S\) 表示 特殊葉子 的貢獻：
\[S=\sum\limits_{n/\delta\left(n\right)\le y\le n}{\mu\left(n\right)\phi\left(\dfrac{x}{n},\pi\left(\delta\left(n\right)\right)-1 \right)}\tag{9} \]
計算 \(S_0\) 顯然是可以在 \(O\left(y\log{\log x}\right)\) 的時間複雜度內解決的，現在我們要考慮如何計算 \(S\)。
計算 S
我們有：
我們將這個等式改寫為：
其中：
注意到計算 \(S_1,S_2\) 的和式中涉及到的 \(m\) 都是質數，證明如下：
如果不是這樣，因為有 \(\delta\left(m\right)>p>x^{1/4}\)，所以有 \(m>p^2>\sqrt{x}\)，這與 \(m\le y\) 矛盾，所以原命題成立。
更多的，當 \(mp>x^{1/2}\ge y\) 時，有 \(y\le mp\)。因此我們有：
計算 S₁
因為：
所以：
所以計算 \(S_1\) 的和式中的項都是 \(1\)。所以我們實際上要計算質數對 \(\left(p,q\right)\) 的個數，滿足：\(x^{1/3}<p<q\le y\)。
因此：
有了這個等式我們便可以在 \(O\left(1\right)\) 的時間內計算 \(S_1\)。
計算 S₂
我們有：
我們將 \(S_2\) 分成 \(q>\dfrac x{p^2}\) 與 \(q\le \dfrac x{p^2}\) 兩部分：
其中：
計算 U
由 \(q>\dfrac x{p^2}\) 可得 \(p^2>\dfrac xq\le \dfrac xy,p>\sqrt{\dfrac xy}\)，因此：
因此：
因此：
因為有 \(\dfrac x{p^2}<y\)，所以我們可以預處理出所有的 \(\pi\left(t\right)\left(t\le y\right)\)，這樣我們就可以在 \(O\left(y\right)\) 的時間複雜度內計算出 \(U\)。
計算 V
對於計算 \(V\) 的和式中的每一項，我們都有 \(p\le \dfrac{x}{pq}<x^{1/2}<p^2\)。因此：
所以 \(V\) 可以被表示為：
其中：
預處理出 \(\pi\left(t\right)\left(t\le y\right)\) 我們就可以在 \(O\left(x^{1/3}\right)\) 的時間複雜度內計算出 \(V_1\)。
考慮我們如何加速計算 \(V_2\) 的過程。我們可以把 \(q\) 的貢獻拆分成若干個 \(\pi\left(\dfrac{x}{pq} \right)\) 為定值的區間上，這樣我就只需要計算出每一個區間的長度和從一個區間到下一個區間的 \(\pi\left(\dfrac{x}{pq} \right)\) 的改變數。
更準確的說，我們首先將 \(V_2\) 分成兩個部分，將 \(q\le \min\left(\dfrac x{p^2},y\right)\) 這個複雜的條件簡化：
接著我們把這個式子改寫為：
其中：
計算 W₁ 與 W₂
計算這兩個值需要計算滿足 \(y<\dfrac{x}{pq}<x^{1/2}\) 的 \(\pi\left(\dfrac{x}{pq} \right)\) 的值。可以在區間 \([1,\sqrt x]\) 分塊篩出。在每個塊中我們對於所有滿足條件的 \((p,q)\) 都累加 \(\pi\left(\dfrac x{pq}\right)\)。
計算 W₃
對於每個 \(p\)，我們把 \(q\) 分成若干個區間，每個區間都滿足它們的 \(\pi\left(\dfrac x{pq}\right)\) 是定值，每個區間我們都可以 \(O(1)\) 計算它的貢獻。當我們獲得一個新的 \(q\) 時，我們用 \(\pi(t)\)（\(t\leq y\)）的值表計算 \(\pi\left(\dfrac x{pq}\right)\)。\(y\) 以內的質數表可以給出使得 \(\pi(t)<\pi(t+1)=\pi\left(\dfrac x{pq}\right)\) 成立的 \(t\)。以此類推使得 \(\pi\left(\dfrac x{pq}\right)\) 變化的下一個 \(q\) 的值。
計算 W₄
相比於 \(W_3\)，\(W_4\) 中 \(q\) 更小，所以 \(\pi\left(\dfrac x{pq}\right)\) 改變得更快。這時候再按照計算 \(W_3\) 的方法計算 \(W_4\) 就顯得沒有任何優勢。於是我們直接暴力列舉數對 \((p,q)\) 來計算 \(W_4\)。
計算 W₅
我們像計算 \(W_3\) 那樣來計算 \(W_5\)。
計算 S₃
我們使用所有小於 \(x^{1/4}\) 的素數一次篩出區間 \(\left[1,\dfrac xy\right]\)。當我們的篩法進行到 \(p_k\) 的時候，我們算出了所有 \(m\) 滿足沒有平方因子並且 \(\delta(m)>p_k\) 的 \(-\mu(m)\phi\left(\dfrac{x}{mp_k},k-1 \right)\) 值。這個篩法是分塊進行的，我們在篩選間隔中維護一個二叉樹，以實時維護所有素數篩選到給定素數後的中間結果。這樣我們就可以只用 \(O(\log x)\) 的時間複雜度求出在篩法進行到某一個值的時候沒有被篩到的數的數量。
演算法的時空複雜度
時空複雜度被如下 \(3\) 個過程影響：
- 計算 \(P_2\left(x,a\right)\)；
- 計算 \(W_1,W_2,W_3,W_4,W_5\)；
- 計算 \(S_3\)。
計算 P₂(x,y) 的複雜度
我們已經知道了這個過程的時間複雜度為 \(O\left(\dfrac{x}{y}\log{\log x}\right)\)，空間複雜度為 \(O\left(y\right)\)。
計算 W₁,W₂,W₃,W₄,W₅ 的複雜度
計算 \(W_1,W_2\) 所進行的塊長度為 \(y\) 的篩的時間按複雜度為 \(O\left(\sqrt{x}\log{\log x}\right)\)，空間複雜度為 \(O\left(y\right)\)。
計算 \(W_1\) 所需的時間複雜度為：
計算 \(W_2\) 的時間複雜度為：
因此，計算 \(W_3\) 的時間複雜度為：
計算 \(W_4\) 的時間複雜度為：
計算 \(W_5\) 的時間複雜度為：
計算 S₃ 的複雜度
對於預處理：由於要快速查詢 \(\phi(u,b)\) 的值，我們沒辦法用普通的篩法 \(O(1)\) 求出，而是要維護一個資料結構使得每次查詢的時間複雜度是 \(O(\log x)\)，因此時間複雜度為 \(O\left(\dfrac{x}{y}\log x\log\log x\right)\)。
對於求和：對於計算 \(S_3\) 和式中的每一項，我們查詢上述資料結構，一共 \(O\left(\log x\right)\) 次查詢。我們還需要計算和式的項數，即二叉樹中葉子的個數。所有葉子的形式均為 \(\pm\phi\left(\dfrac{x}{mp_b},b-1\right)\)，其中 \(m\le y,b<\pi(x^{1/4})\)。因此，葉子的數目是 \(O\left(y\pi\left(x^{1/4}\right)\right)\) 級別的。所以計算 \(S_3\) 的總時間複雜度為：
總複雜度
這個演算法的空間複雜度為 \(O\left(y\right)\)，時間複雜度為：
我們取 \(y=x^{1/3}\log^3{x}\log{\log x}\)，就有最優時間複雜度為 \(O\left(\dfrac{x^{2/3}}{\log^2 x}\right)\)，空間複雜度為 \(O\left(x^{1/3}\log^3{x}\log{\log x}\right)\)。
一些改進
我們在這裡給出改進方法，以減少演算法的常數，提高它的實際效率。
-
在 終止條件 \(2\) 中，我們可以用一個 \(z\) 來代替 \(y\)，其中 \(z\) 滿足 \(z>y\)。我們可以證明這樣子計算 \(S_3\) 的時間複雜度可以優化到：
\[O\left(\dfrac{x}{z}\log x\log{\log x}+\dfrac{yx^{1/4}}{\log x}+z^{3/2} \right) \]這也為通過改變 \(z\) 的值來檢查計算提供了一個很好的方法。
-
為了清楚起見，我們在闡述演算法的時候選擇在 \(x^{1/4}\) 處拆分來計算總和 \(S\)，但實際上我們只需要有 \(p\le \dfrac{x}{pq}<p^2\) 就可以計算。我們可以利用這一點，漸近複雜性保持不變。
-
用前幾個素數 \(2,3,5\) 預處理計算可以節省更多的時間。
References
本文翻譯自：Computing \(\pi(x)\): the Meissel, Lehmer, Lagarias, Miller, Odlyzko method