min25篩學習總結

前言

杜教篩學了，順便把min25篩也學了吧= =剛好多校也有一道題需要補。
下面推薦幾篇博客，我之后寫一點自己的理解就是了。
傳送門1
傳送門2
傳送門3

這幾篇寫得都還是挺好的，接下來我就寫下自己對min25篩的理解吧 。

正文

簡介：

min25篩同杜教篩類似，是用來解決一類積性函數的前綴和，即\(\sum_{i=1}^nF(i)\)，並且這里的\(n\)可以達到\(10^{10}\)的規模。
但所求積性函數要求滿足以下條件：

• \(F(p)\)可以表示為簡單多項式的形式，比如\(p_1^{k_1}+p_2^{k_2}+p_3^{k_3}\)，\(p\)為質數；
• \(F(p^e)\)的值要較容易求得。

為什么是這樣？往下看就知道啦。

求解：

我們可以將上面說的多項式拆成的每項單獨來算前綴和，最后加起來即可。
設\(f(p)=p^k\)，我們現在就來求\(f(i)\)的前綴和。注意這里\(f\)函數是完全積性函數，后面我們需要用到這一性質。

Part 1:
首先設\(g(n,j)\)表示：\(1,2...n\)中，滿足條件的\(f\)之和，滿足條件是指要么為質數，要么其最小質因子大於第\(j\)個質數。
為什么要這么設？后面就知道啦。
之后我們考慮遞推求解：

• 若\(p_j^2>n\)，根據定義，\(g(n,j)=g(n,j-1)\)；
• 否則，我們要減去最小質因子為\(p_j\)的數，那么就有遞推式：\(g(n,j)=g(n,j-1)-f(p_j)*(g(\frac{n}{p_j},j-1)-sum_{j-1})\)。
• 解釋一下上式，就相當於首先提取一個\(p_j\)出來，那么只要減去剩下的質因子大於等於\(p_j\)的就行啦；但根據定義，\(g\)中還包含了\(\sum_{i=1}^{j-1}f(p_i)\)，我們減去就行了。另外上式還利用了完全積性函數的性質。

這里是min25篩的關鍵。

我們還可以用篩法的思想去考慮，我們從\(p_{j-1}\)遞推到\(p_j\)的過程，其實就是在埃氏篩中，把\(p_j\)的倍數篩去。每次得到的\(g(n,j)\)，其實就是利用前\(j\)個素數來進行埃氏篩，剩下的數以及所有質數的\(f\)之和。
這也是為啥min25后面有個“篩”的原因~

從篩法的角度來考慮，那么初始化時我們將所有的數都看作素數，然后一個一個來篩求解即可（1除外，我們先不考慮1，最后考慮即可）。

Q:那為什么要求出這個呢？
A:求出\(g\)之后，我們可以方便地求出\(\sum_pf(p)\)的值，即\(g(n,|P|)\)，其中\(|P|\)為素數集大小。

那合數怎么辦？

Part 2:
上面我們把質數的情況考慮了，利用了類似埃氏篩的思想遞推得到了\(g\)函數。
接下來我們考慮合數的情況，因為\(F\)為積性函數，那么我們直接枚舉最小質因子以及其次數來求解即可。
類似地，定義\(S(n,j)\)表示\(1,2,\cdots,n\)中，滿足條件的\(F\)之和，滿足條件是指最小質因子大於等於(這里有個等於，其實也可以沒有~)第\(j\)個質數。
那么可以直接暴力求解\(S\):

\[S(i,j)=g(i,|P|)-sum_{j-1}+\sum_{k\geq j}\sum_{e}F(p_k^e)(S(\frac{i}{p_k^e},k+1)+[e\not ={1}]) \]

稍微解釋一下：前面部分就相當於我們求出了滿足條件的\(F\)之和，這里的條件指的是大於等於\(p_j\)的質數。
接下來我們就相當於暴力枚舉最小的質因子以及其次方，並且根據其積性函數的性質，將枚舉值提出來，然后遞歸求解子問題。
但我們對於\(p_k^e\)這種數沒有統計到，所以后面加上就行了。
那么最終答案就為\(S(n,1)+F(1)\)。

所以...min25篩就這么完了。似乎也不是很難。
但我還想加個

Part 3:
雖然min25篩的思想說得差不多了，但我還想說一下其實現的一些細節。畢竟代碼也是很重要的。
下面以模板題為例：

Code

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 2e6 + 5, MOD = 1e9 + 7, inv3 = 333333336;
ll n;
ll sum1[N], sum2[N], prime[N];
ll w[N], ind1[N], ind2[N];
ll g1[N], g2[N];
bool chk[N];
int tot, cnt;
void pre(int n) { //  \sqrt
    chk[1] = 1;
    for(int i = 1; i <= n; i++) {
        if(!chk[i]) {
            prime[++tot] = i;
            sum1[tot] = (sum1[tot - 1] + i) % MOD;
            sum2[tot] = (sum2[tot - 1] + 1ll * i * i % MOD) % MOD;
        }
        for(int j = 1; j <= tot && prime[j] * i <= n; j++) {
            chk[i * prime[j]] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}
void calc_g() {
    int z = sqrt(n);
    for(ll i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i);
        w[++cnt] = n / i;
        g1[cnt] = w[cnt] % MOD;
        g2[cnt] = g1[cnt] * (g1[cnt] + 1) / 2 % MOD * (2 * g1[cnt] + 1) % MOD * inv3 % MOD - 1;
        g1[cnt] = g1[cnt] * (g1[cnt] + 1) / 2 % MOD - 1;
        if(n / i <= z) ind1[n / i] = cnt;
        else ind2[n / (n / i)] = cnt;
    }
    for(int i = 1; i <= tot; i++) {
        for(int j = 1; j <= cnt && prime[i] * prime[i] <= w[j]; j++) {
            ll tmp = w[j] / prime[i], k;
            if(tmp <= z) k = ind1[tmp]; else k = ind2[n / tmp];
            (g1[j] -= prime[i] * (g1[k] - sum1[i - 1] + MOD) % MOD) %= MOD;
            (g2[j] -= prime[i] * prime[i] % MOD * (g2[k] - sum2[i - 1] + MOD) % MOD) %= MOD;
            if(g1[j] < 0) g1[j] += MOD;
            if(g2[j] < 0) g2[j] += MOD;
        }
    }
}
ll S(ll x, int y) { // 2~x >= P_y
    if(x <= 1 || prime[y] > x) return 0;
    ll z = sqrt(n);
    ll k = x <= z ? ind1[x] : ind2[n / x];
    ll ans = (g2[k] - g1[k] + MOD - (sum2[y - 1] - sum1[y - 1]) + MOD) % MOD;
    for(int i = y; i <= tot && prime[i] * prime[i] <= x ; i++) {
        ll pe = prime[i], pe2 = prime[i] * prime[i];
        for(int e = 1; pe2 <= x; ++e, pe = pe2, pe2 *= prime[i]) {
            ans = (ans + pe % MOD * ((pe - 1) % MOD) % MOD * S(x / pe, i + 1) + pe2 % MOD * ((pe2 - 1) % MOD) % MOD) % MOD;
        }
    }
    return ans % MOD;
}
int main() {
    cin >> n;
    int tmp = sqrt(n);
    pre(tmp);
    calc_g();
    cout << (S(n, 1) + 1) % MOD ;
    return 0;
}

接下來我會說說代碼里面的一些要點：

• 預處理部分

沒什么好說的，不同的問題預處理也不同。

• 求\(g\)部分

\[g(n,j)=g(n,j-1)-f(p_j)(g(\frac{n}{p_j},j-1)-sum_{j-1}) \]

我先把式子抄下來...懶得翻上去了。
觀察這個式子可以發現，第二維每次都是從\(j-1\)優化過來，所以我們可以直接滾動掉一維。並且，每次遞推時，都是從\(\frac{n}{p_j}\)轉移過來，我們聯想到了\(\lfloor\frac{n}{i}\rfloor\)這個形式。
顯然，因為\(n\)可能會很大，我們不能直接將\(n\)求出來，質數也是，篩那么大的數，線性篩也會超時。
下面有兩個觀察：

• 觀察1：遞推中，有用的素數只有不超過\(\sqrt{n}\)的部分；
• 觀察2：因為\(\lfloor\frac{n}{i}\rfloor\)只有\(\sqrt{n}\)個不同的值，所以我們只需要預處理這\(\sqrt{n}\)個值就行了。

因為有了觀察1，我們遞推到\(\sqrt{n}\)就相當於算出了當\(j=|P|\)的情況。
但可能某些值還是很大，怎么辦？
首先我們肯定需要一個數組\(w\)來儲存所有的\(\lfloor\frac{n}{i}\rfloor\)，之后用了兩個\(ind\)數組來存儲下標，如果\(\lfloor\frac{n}{i}\rfloor<\sqrt{n}\)，那么直接存儲；否則就在另外一個數組存儲\(\lfloor\frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor\)的下標。
那么之后我們就可以通過\(\lfloor\frac{n}{i}\rfloor\)的形式來訪問相關值了。
這里類似於新加了一個映射關系。\(w\)數組映射到\(\lfloor\frac{n}{i}\rfloor\)，而\(ind\)數組用了點trick把\(\lfloor\frac{n}{i}\rfloor\)映射回\(w\)數組。

• 求\(S\)部分

\[S(i,j)=g(i,|P|)-sum_{j-1}+\sum_{k\geq j}\sum_{e}f(p_k^e)(S(\frac{i}{p_k^e},k+1)+[e\not ={1}]) \]

有一個問題是，這里怎么得到\(|P|\)。
因為我們把\(g\)滾動了的，最后雖然求出的是為\(\sqrt{n}\)的情況，但稍加思考就會發現，其值等於\(|P|\)下的值。
怎么得到\(i\)的下標？注意我們是從\(n\)開始遞推，每次會除以一個數，並且有這樣一個形式：\(\lfloor\frac{\lfloor\frac{n}{a}\rfloor}{b}\rfloor=\lfloor\frac{n}{ab}\rfloor\)，那么每次的\(i\)都一定可以表示為\(\lfloor\frac{n}{k}\rfloor\)的形式。所以直接根據剛才的映射關系來找就行啦。這也是為什么需要映射關系的原因

感覺代碼部分也說得差不多了，可能有些地方我理解的有問題或者有點復雜，各位看官也不要太過於糾結字眼...自己理解也是不錯的辦法。
反正我感覺這里利用\(\lfloor\frac{n}{i}\rfloor\)的性質來求解很關鍵，優化了代碼時間復雜度以及編寫難度（求\(g\)和\(S\)都用到了它相關性質）。十分巧妙。

時間復雜度：(我也不會)O(感覺能過)\(O(\frac{n^\frac{3}{4}}{log_n})\)。

例題：

LOJ#6053. 簡單的函數

思路：
因為有\(f(p^c)=p^c\)，注意到除開\(2\)以外的素數都為奇數，那么其余\(f(p)=p-1\)。我們處理的時候將\(f(2)\)的值也當作\(1\)來算，最后加回來就行。
基本上也是一道模板題吧~

Code

#include <bits/stdc++.h>
#define heyuhhh ok
using namespace std;
typedef long long ll;
const int N = 1e6 + 5, MOD = 1e9 + 7, inv2 = (MOD + 1) / 2;
ll n;
ll sum1[N], sum2[N], prime[N];
ll w[N], ind1[N], ind2[N];
ll g1[N], g2[N];
bool chk[N];
int tot, cnt;
void pre(int n) { //  \sqrt
    chk[1] = 1;
    for(int i = 1; i <= n; i++) {
        if(!chk[i]) {
            prime[++tot] = i;
            sum1[tot] = (sum1[tot - 1] + i) % MOD;
        }
        for(int j = 1; j <= tot && prime[j] * i <= n; j++) {
            chk[i * prime[j]] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}
void calc_g() {
    int z = sqrt(n);
    for(ll i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i);
        w[++cnt] = n / i;
        g2[cnt] = (w[cnt] - 1) % MOD;
        g1[cnt] = (w[cnt] % MOD) * ((w[cnt] + 1) % MOD) % MOD * inv2 % MOD; --g1[cnt];
        if(n / i <= z) ind1[n / i] = cnt;
        else ind2[n / (n / i)] = cnt;
    }
    for(int i = 1; i <= tot; i++) {
        for(int j = 1; j <= cnt && prime[i] * prime[i] <= w[j]; j++) {
            ll tmp = w[j] / prime[i], k;
            if(tmp <= z) k = ind1[tmp]; else k = ind2[n / tmp];
            (g1[j] -= prime[i] * (g1[k] - sum1[i - 1] + MOD) % MOD) %= MOD;
            (g2[j] -= (g2[k] - i + 1 + MOD) % MOD) %= MOD;
            if(g1[j] < 0) g1[j] += MOD;
            if(g2[j] < 0) g2[j] += MOD;
        }
    }
}
ll S(ll x, int y) { // 2~x >= P_y
    if(x <= 1 || prime[y] > x) return 0;
    ll z = sqrt(n);
    ll k = x <= z ? ind1[x] : ind2[n / x];
    ll ans = (g1[k] - sum1[y - 1] + MOD - g2[k] + y - 1 + MOD) % MOD;
    if(y == 1) ans += 2;
    for(int i = y; i <= tot && prime[i] * prime[i] <= x ; i++) {
        ll pe = prime[i];
        for(int e = 1; pe * prime[i] <= x; ++e, pe = pe * prime[i]) {
            ans = (ans + (prime[i] ^ e) * S(x / pe, i + 1) % MOD + (prime[i] ^ (e + 1))) % MOD;
        }
    }
    return ans;
}
int main() {
#ifdef heyuhhh
    freopen("input.in", "r", stdin);
#else
    ios::sync_with_stdio(false); cin.tie(0);
#endif
    cin >> n;
    int tmp = sqrt(n);
    pre(tmp);
    calc_g();
    cout << (S(n, 1) + 1) % MOD ;
    return 0;
}

牛客多校第七場K.Function
題意：
求\(\sum_{i=1}^nf(i)\)，其中，\(f(i)\)為：

\[\left\{ \begin{aligned} &3e+1,&i= p^e \ and\ p\%4=1\\ &1,&else \end{aligned} \right. \]

思路：
考慮\(min25\)篩求解。
我們不管1，2的存在，最后加上其貢獻即可。
首先求出\(g_1,g_2\)，分別質數表示\(\% 4\)為1和3的情況總數，然后一個一個來篩。篩的時候要同時考慮當前質數模4的情況以及對\(g_1,g_2\)的影響，它們都會減去某個值。
然后就直接求和就行了，基本都是套用板子。
最主要的還是求\(g\)的思路。

Code

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 2e6 + 5;
int T;
ll n;
ll sum1[N], sum2[N], prime[N];
ll w[N], ind1[N], ind2[N];
ll g1[N], g2[N];
bool chk[N];
int tot, cnt;
void pre(int n) { //  \sqrt
    chk[1] = 1;
    for(int i = 1; i <= n; i++) {
        if(!chk[i]) {
            prime[++tot] = i;
            sum1[tot] = sum1[tot - 1];
            sum2[tot] = sum2[tot - 1];
            if(i % 4 == 1) ++sum1[tot];
            if(i % 4 == 3) ++sum2[tot];
        }
        for(int j = 1; j <= tot && prime[j] * i <= n; j++) {
            chk[i * prime[j]] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}
void calc_g() {
    int z = sqrt(n);
    for(ll i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i);
        w[++cnt] = n / i;
        g1[cnt] = w[cnt] / 4 + (w[cnt] % 4 >= 1) - 1;
        g2[cnt] = w[cnt] / 4 + (w[cnt] % 4 >= 3);
        if(n / i <= z) ind1[n / i] = cnt;
        else ind2[n / (n / i)] = cnt;
    }
    for(int i = 1; i <= tot; i++) {
        for(int j = 1; j <= cnt && prime[i] * prime[i] <= w[j]; j++) {
            ll tmp = w[j] / prime[i], k;
            if(tmp <= z) k = ind1[tmp]; else k = ind2[n / tmp];
            if(prime[i] % 4 == 1) {
                g1[j] -= (g1[k] - sum1[i - 1]);
                g2[j] -= (g2[k] - sum2[i - 1]);
            } else if(prime[i] % 4 == 3) {
                g1[j] -= (g2[k] - sum2[i - 1]);
                g2[j] -= (g1[k] - sum1[i - 1]);
            }
        }
    }
}
int f(int p, int q) {
    if(p % 4 == 1) return 3 * q + 1;
    return 1;
}
ll S(ll x, int y) { // 2~x >= P_y
    if(x <= 1 || prime[y] > x) return 0;
    ll z = sqrt(n);
    ll k = x <= z ? ind1[x] : ind2[n / x];
    ll ans = 4 * g1[k] - 4 * sum1[y - 1] + g2[k] - sum2[y - 1];
    if(y == 1) ans++;
    for(int i = y; i <= tot && prime[i] * prime[i] <= x ; i++) {
        ll pe = prime[i], pe2 = prime[i] * prime[i];
        for(int e = 1; pe2 <= x; ++e, pe = pe2, pe2 *= prime[i]) {
            ans += f(prime[i], e) * S(x / pe, i + 1) + f(prime[i], e + 1);
        }
    }
    return ans;
}
int main() {
    freopen("input.in", "r", stdin);
    cin >> T;
    while(T--) {
        memset(chk, 0, sizeof(chk)); tot = cnt = 0;
        cin >> n;
        int tmp = sqrt(n);
        pre(tmp);
        calc_g();
        cout << S(n, 1) + 1 << '\n';
    }
    return 0;
}

后記

寫得還是比較匆忙，感覺也沒有很好地把心靜下來，可能有些地方寫得冗余、累贅或者有錯誤，請諒解。
總得來說，min25篩的思想以及代碼實現部分都是很巧妙的~為什么會想到搞個\(g\)來啊？為什么這樣復雜度就比較低啊？這些都是問題...
另外加一些我學習時的草稿：

因為轉移的原因，所以有用的素數只有小於等於\(\sqrt{n}\)的。最后算出來的值就是\(g(n,|P|)\)啦。

處理\(g\)時，離散化\(\lfloor\frac{n}{i}\rfloor\)，因為第一維為\(n\)不會變，並且有個這樣性質：\(\lfloor\frac{\lfloor\frac{n}{a}\rfloor}{b}\rfloor=\lfloor\frac{n}{ab}\rfloor\)，所以有用的值也就只有\(\lfloor\frac{n}{i}\rfloor\)這些。

離散化要點：利用\(w\)遞減記錄離散化的值，並且利用\(ind1,ind2\)來記錄相關數的下標，不超過\(\sqrt{n}\)。涉及整除的性質。巧妙！