Min_25 篩這個東西,完全理解花了我很長的時間,所以寫點東西來記錄一些自己的理解。
它能做什么
對於某個數論函數 \(f\),如果滿足以下幾個條件,那么它就可以用 Min_25 篩來快速求出這個函數的前綴和。
-
它是一個積性函數
-
對於一個質數 \(p\) ,\(f(p)\) 的表達式必須是一個項數比較小的多項式。即 \(\displaystyle f(p) = \sum a_ip^{b_i}\)。
-
對於一個質數 \(p\) ,\(f(p^k)\) 的表達式必須可以由 \(f(p)\) 快速得到。
看起來這些條件比較緊,其實仔細想想許多函數還是滿足這些性質的。接下來討論它的實現過程。
實現過程
Min_25 篩的核心思想就是對於 \(\le n\) 的數中,除掉質數之后,每個數的最小質因子大小不超過 \(\sqrt{n}\) 。即類似埃氏篩的方式,從每個質因子篩掉一部分的貢獻。整個過程可以分成三部分,即質數的部分的貢獻之和與合數的部分的貢獻和,以及 \(1\) 的特判。
記號約定
這里用 \(P\) 表示質數集合,\(P_i\) 表示第 \(i\) 小的質數。
質數部分
首先我們不妨設 \(f(p) = p^k\) 。對於更加一般的項數小的多項式形式,我們可以分別記每一個項的貢獻,然后分別相加。注意到這里 \(f(p)\) 是一個完全積性函數。對於 \(p\) 為合數的情況,為了方便,我們先姑且同樣定義 \(f(p) = p^k\) ,到最后你會明白這樣定義的原因。
我們再設 \(g(n,i)\) 表示 \(\le n\) 的所有數 \(x\) 中,\(x\) 的最小質因子 \(> P_i\) 或者 \(x\) 本身為質數的所有 \(f(x)\) 之和,這樣 \(g(n,|P|)\) 就是所有質數的貢獻之和了。而 \(g(n,0)\) 的答案也十分顯然,即 \(\displaystyle \sum_{i = 2} ^{n} i^k\) 。我們只需要考慮,如何從 \(g(*,i - 1)\) 如何得到 \(g(*,i)\) 。
考慮我們從 \(i - 1\) 推到 \(i\) 的過程中,會減少哪個部分的貢獻,容易發現最小質因子恰好為 \(P_i\) 的貢獻被減掉了。如果 \(P_i^2 > n\) ,那么這一部分是沒有貢獻的,即 \(g(n,i) = g(n,i - 1)\)。
對於 \(P_i^2 \le n\) 的情況,我們需要把所有 \(P_i\) 的倍數,且滿足其最小質因子除去 \(P_i\) 之外 \(> P_{i-1}\) 的貢獻全部減掉。如果直接減去 \(f(P_i) \cdot g(\lfloor\frac{n}{P_i}\rfloor,i - 1)\) 你會發現貢獻是減多了的,原因就是 \(g\) 函數的值還包含了所有質數的貢獻之和,仔細分析一下,可以發現,對於 \(x = P_iP_j(j < i)\) 的貢獻在 \(j\) 這個地方已經減過了,所以還要加回來,即 \(\displaystyle f(P_i) \cdot g(P_i - 1,i - 1) = \sum_{j = 1}^{i - 1} f(P_i)f(P_j)\) 。
式子匯總起來也就是
回頭來,我們再來看最初我們提出的一些問題。
我們最開始為什么要定義,對於一個合數 \(p\) ,\(f(p) = p^k\) 呢?無非就是方便計算,方便利用完全積性函數的性質。你發現到最后,合數的貢獻還是沒有算進去,所以對於合數的定義是無所謂的。
我們到底在哪里利用了 \(f(p)\) 是一個完全積性函數的性質呢?這個性質是在上式中第二種情況中利用到的。注意到我們在減貢獻的時候,並沒有枚舉 \(P_i\) 的次數,這也就意味着 \(g(\lfloor\frac{n}{P_i}\rfloor,i - 1)\) 中可能包含了某些數仍然是 \(P_i\) 的倍數。所以這一步需要有一個對於所有質數 \(p,f(p^2) = f^2(p)\) 的條件。由於之前我們假定了 \(f(p)\) 對於 \(p\) 為合數的情況也滿足 \(f(p) = p^k\) ,因此這個條件是得到滿足了的。
遞歸實現很顯然。考慮一下非遞歸的實現。首先對於這個步驟,我們有用的 \(n\) 只有所有不同的 \(\lfloor\frac{n}{i}\rfloor,(1 \le i \le n)\),所以可以離散掉這些值,用一個大小為 \(2\sqrt{n}\) 的數組存下來。由於這個遞推是一層一層推的,所以實現上可以不用開二維數組,可以直接離散之后用一個一維數組存下來,需要注意枚舉順序。
合數部分
考慮我們在求質數部分的時候,除了質數的答案之和之外,還求出了些什么東西。事實上,對於所有的 \(\lfloor\frac{n}{i}\rfloor,(1 \le i \le n)\) ,我們求出了 \(g(\lfloor\frac{n}{i}\rfloor,|P|)\) 的值。
我們設 \(S(n,i)\) 表示 \(\le n\) 的所有數 \(x\) 中,\(x\) 的最小質因子 \(\ge P_i\) 的 \(f(x)\) 之和。注意這里的 \(f(x)\) 在 \(x\) 為合數的時候不再沿用質數部分時的定義,也注意這里的定義與 \(g\) 的區別。
接下來我們先算上所有滿足條件的質數貢獻之和,即 \(\displaystyle g(n,|P|) - \sum_{j = 1}^{i - 1} f(P_j)\) 。對於合數,我們直接枚舉其最小質因子以及質因子的個數,直接遞歸計算。注意在這里形如 \(p^k,(p \in P)\) 的貢獻並沒有算進去,所以還要單獨加一下。式子的形式如下:
直接對着式子看,應該還是很好理解的。對這個式子,直接遞歸暴力計算即可,跑得會比第一部分快。
至此 Min_25 篩的全過程就講完了。時間復雜度是 \(\displaystyle O(\frac{n^{\frac{3}{4}}}{\log n})\) 的,不大會證……
例題
首先可以用一個小問題來感受 Min_25 篩的威力。
SPOJ DIVCNT2
題意
求
\[\sum_{i = 1}^{n} \sigma_0(i^2) \]其中 \(\sigma_0(i)\) 表示 \(i\) 的約數個數。
\(n \le 10^{12},TimeLimit = 20s\) 。
題解
傳統的杜教篩做法需要考慮將式子拆為每個質因子貢獻相乘,然后反演 + 杜教篩。如果我們有 Min_25 篩呢?
我們發現 \(f(x) = \sigma_0(x^2)\) 滿足 Min_25 篩的三個基本條件,\(f(p) = 3,f(p^k) = 2k + 1\) ,因此就是一道裸題了。
我們發現 SPOJ 還有兩道題叫 DIVCNT3 / DIVCNTK。於是此題三倍經驗……
代碼
#include<bits/stdc++.h>
#define For(i,l,r) for(int i = (l),i##end = (r);i <= i##end;i++)
#define Fordown(i,r,l) for(int i = (r),i##end = (l);i >= i##end;i--)
#define debug(x) cout << #x << " = " << x << endl
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
template <typename T> inline bool chkmin(T &x,T y) { return y < x ? x = y,1 : 0; }
template <typename T> inline bool chkmax(T &x,T y) { return x < y ? x = y,1 : 0; }
const int INF = 0x3f3f3f3f;
const int N = 2e6 + 10;
int prime[N],vis[N],id1[N],id2[N];
ull B[N],Ans[N];
ll val[N],n;
int pcnt = 0;
inline ll read() {
ll x = 0,flag = 1;
char ch = getchar();
while(!isdigit(ch) && ch != '-') ch = getchar();
if(ch == '-') flag = -1,ch = getchar();
while(isdigit(ch)) x = (x << 3) + (x << 1) + (ch - '0'),ch = getchar();
return x * flag;
}
inline void init(int Lim) {
For(i,2,Lim) {
if(!vis[i]) prime[++pcnt] = i;
for(int j = 1;j <= pcnt && i * prime[j] <= Lim;j++) {
vis[i * prime[j]] = true;
if(i % prime[j] == 0) break;
}
}
}
#define id(x) (x <= d ? id1[x] : id2[n / (x)])
inline ull S(ll Lim,int cnt) {
if(Lim <= 1 || prime[cnt] > Lim) return 0;
int d = sqrt(n);
ull ans = 3 * (B[id(Lim)] - (cnt - 1));
for(int i = cnt;i <= pcnt && 1ll * prime[i] * prime[i] <= Lim;i++) {
ll prod = prime[i];
for(int times = 1;1ll * prod * prime[i] <= Lim;times++,prod = 1ll * prod * prime[i])
ans += (2 * times + 1) * S(Lim / prod,i + 1) + (2 * times + 3);
}
return ans;
}
inline void Calc(ll Lim) {
int d = sqrt(Lim),cnt = 0;
for(ll i = 1;i <= Lim;i = Lim / (Lim / i) + 1) {
val[++cnt] = Lim / i,id(val[cnt]) = cnt;
B[cnt] = val[cnt] - 1;
}
for(int i = 1;i <= pcnt && 1ll * prime[i] * prime[i] <= Lim;i++)
for(int j = 1;j <= cnt && 1ll * prime[i] * prime[i] <= val[j];j++)
B[j] = B[j] - B[id(val[j] / prime[i])] + (i - 1);
}
int main() {
#ifndef ONLINE_JUDGE
freopen("DIVCNT2.in","r",stdin);
freopen("DIVCNT2.out","w",stdout);
#endif
int Case = read();
init(2e6 + 5);
while(Case--) {
n = read();
Calc(n);
ull res = S(n,1) + 1;
printf("%llu\n",res);
}
return 0;
}
Min_25 篩在一些質因子相關的問題上,也有很強的適用性,比如下面的這些題。
UOJ188 Sanrd
題意
設 \(f(x)\) 為 \(x\) 的次大質因子,重復質因子算多個。求
\[\sum_{i = l} ^ {r} f(i) \]\(l,r \le 10^{11},TimeLimt = 5s\) 。
題解
這個函數顯然不是積性函數了,但是這個次大質因子的形式,和 Min_25 篩枚舉最小質因子的形式比較契合,因此可以想到在 Min_25 篩過程中統計答案。
考慮我們在遞歸計算 \(S\) 的過程中,假設當前遞歸到了 \(S(a,b)\) ,那么這些 \(\le a\) 且 \(> P_b\) 的質數,在最開始的時候,一定是 \(P_{b - 1}\) 的倍數(具體原因可以仔細看看 \(S\) 的遞歸計算的式子)。所以這些數的次大質因子都是 \(P_{b - 1}\) 。因此就可以在 Min_25 篩過程中直接計算了。
LOJ572 Misaka Network 與求和
題意
設 \(f(x)\) 為 \(x\) 的次大質因子,重復質因子算多個。求
\[\sum_{i = 1}^{n}\sum_{j = 1}^{n}f(\gcd(i,j))^k \bmod 2^{32} \]\(n,k \le 2 \times 10^9,TimeLimit = 4s\) 。
題解
前面經典的反演比較簡單,式子可以轉化為
設 \(F(x) = f(x)^k\) ,則所求即為
於是我們只要能夠快速求出 \(F * \mu\) 的前綴和,就可以整除分塊了。我們設 \(\displaystyle S(n) = \sum_{i = 1}^{n} (F * \mu)(n)\),可以寫出一個杜教篩的式子。
由於 \(\mu * 1 = \epsilon\) ,容易想到令 \(g = 1\) 。那么 \((F * \mu * g)(i) = (F * (\mu * g))(i) = (F * \epsilon)(i) = F(i)\)
對於 \(F(i)\) 的前綴和,直接套用 Min_25 篩即可。