[算法]Min_25篩

前言

本篇文章中使用的字母\(p\)，指\(\text{任意的} p \in \text{素數集合}\)

應用場景

若函數\(f(x)\)滿足，

1. \(f(x)\)是積性函數
2. \(f(p)\)可以使用多項式表示。
3. 已知\(f(p)\)，要能在常數級的時間內計算\(f(p^x),x \in N^+\)。

Min_25篩可以在\(\Theta(\frac{n^{\frac{3}{4}}}{log_2n})\)的時間復雜度內計算\(f(x)\)的前綴和
或者說\(\Theta(n ^ {1 - \epsilon})\)?個人傾向上面那種。

算法

分類

\[\sum_{i=1}^n f(i)=\sum_{p=1}^nf(p),p \in prime\ + \sum_{q=1}^nf(q) ,q \notin prime \]

prime指素數集合

輔助函數\(g(x,y)\)的構造與遞推

首先線性篩出質數，線性篩就不多說了（如果您線性篩都不會建議先去學習基礎算法）
我們從小到大設\(p_i\)為從小到大排列的第i個質數
例如，\(p_1=2,p_2=3,\dots\)
設素數集合為prime，i的最小質因子為\(MPF_i\) (minimal prime fact of i)。
\(||\)為或者，即\(or\); \([a]\)表示a成立是為\(1\)，否則為\(0\)
設函數\(g(x,y)\)，使得

\[g(x,y)=\sum_{i=1}^x[i\in prime\ ||\ MPF_i>p_y]f(i) \]

考慮埃拉托斯特尼篩法，每次選出一個質數，篩掉它的倍數。
注:埃拉托斯特尼篩法的百度百科
我們發現\(g(x,y)\)有一個神奇的性質，我們的\(g(x,y)\)就是\([1,x]\)運行y次篩法以后剩余所有數之和加上所有的\(f(p),p<x\)的和。
我們所求的(prime為質數集合)

\[\sum_{i=1}^x[i\in prime] \]

實際上等同於\(g(x,|prime|)\)，|prime|表示\([1,x]\)之內的質數集合的大小。
首先我們要了解\(g(x,y)\)的初值
\(g(x,0)\)表示所有數的和，也就是把所有數帶入\(f(x)\)計算出的結果。
那么接下來就是\(g(x,y)\)的遞推了。

1. \(p_y^2 > x\)。 此時埃篩的第\(j\)次沒有篩去任何質數(理論上第\(y\)次應該篩掉的最小質數為\(p_y^2\))，所以\(g(x,y)=g(x,y-1)\)
2. \(p_y^2 \leq x\)。此時埃篩篩去了所有大於\(p_y\)倍的\(p_y\)的倍數，若我們從\(g(x,y-1)\)遞推，顯然有多計算的部分，該部分就是

\[f(p_y)\times(g(\frac{x}{p_y},y-1)-\sum_{i=1}^{y-1}f(i)) \]

表示成公式就是:

\[g(x,y) = \left\{ \begin{matrix} g(x,y-1) & p_y^2>x \\ \\ g(x,y-1)-f(p_y)\times(g(\frac{x}{p_y},y-1)-\sum_{i=1}^{y-1}f(i)) & p_y^2 \leq x\end{matrix} \right. \]

前綴函數\(S(x,y)\)的構造與遞推

好了現在我們已經有了一個輔助函數\(g(x,y)\)，但這玩意貌似一點用都沒有....
我們設

\[S(x,y)=\sum_{i=1}^x[MPF_p \geq p_y]f(i) \]

講人話就是所有最小質因子大於等於\(p_j\)的\(f\)值之和。
如果我們要求\(\sum f(i),i \in [1,n]\)，我們要求的東西是

\[S(n,1)+f(1) \]

鑒於所有質數對答案的貢獻我們已經完成計算，它的貢獻是

\[g(n,j)-\sum_{i=1}^{j-1}f(p_i) \]

(我們要保證最小質因子大於等於\(p_j\)所以要把小於它的質數減去)
考慮合數對答案的貢獻，枚舉合數的最小質因子和它的出現次數，然后直接計算。
合數對答案的貢獻是:

\[\sum_{i=y}^{p_i^2 \leq x} \sum_{j = 1}^{p_i^{j + 1} \leq x}S(\frac{x}{p_i^j},i+1)+f(p_{i}^{j+1}) \]

總結起來就是

\[S(x,y) = g(x,|prime|) - \sum_{i=1}^{y-1}f(p_i) + \sum_{j=y}^{p_j^2 \leq x} \sum_{k = 1}^{p_j^{k + 1} \leq x}S(\frac{x}{p_j^k},j+1)+f(p_{j}^{k+1}) \]

\(prime\)指素數集合，\(||\)表示取當前的最大質因子。

復雜度比較

借用一下WC2019課件中的一張圖片

例題

[LOJ6053]簡單的函數

代碼

#include <cstdio>
#include <cmath> 
#define ll long long
#define MAXN 1000005
#define MOD 1000000007

long long sqrtn;

ll read(){
    ll x = 0; int zf = 1; char ch = ' ';
    while (ch != '-' && (ch < '0' || ch > '9')) ch = getchar();
    if (ch == '-') zf = -1, ch = getchar();
    while (ch >= '0' && ch <= '9') x = x * 10 + ch - '0', ch = getchar(); return x * zf;
}

// the number n, cac the sum of f(i) i belongs to [1, n]
ll n; 
// primes   number of primes in total
ll p[MAXN], pcnt;
// g is the function g, the sum of the f(prime).
ll g[MAXN];
// the number of the prime less the pn_i
ll pn[MAXN];
// the id of every helpful number
ll id[2][MAXN];
// prefix of p
ll prep[MAXN];

// the different value of zcfk, or (n/i) 
// zcfk -> https://www.cnblogs.com/linzhengmin/p/11061244.html
ll w[MAXN], wcnt;
// is not a prime -> array in the Euler Sieve
bool np[MAXN];

// Euler Sieve
void getPrime(int n){
    np[0] = np[1] = 0;
    for (int i = 2; i <= n; ++i){
        if (!np[i]) p[++pcnt] = i, prep[pcnt] = (prep[pcnt - 1] + i) % MOD;
        for (int j = 1; j <= pcnt && i * p[j] <= n; ++j){
            np[p[j] * i] = 1;
            if (!(i % p[j])) break;
        }
    }
}

// cac function s 
int S(ll x, int y){
    if (x <= 1 || p[y] > x) return 0;
    int k = (x <= sqrtn) ? id[0][x] : id[1][n / x];
    int res = (((ll)g[k] - pn[k] - prep[y - 1] + y - 1) % MOD + MOD) % MOD;
    if (y == 1) res += 2;
    for (int i = y; i <= pcnt && (ll)p[i] * p[i] <= x; ++i){
        ll p1 = p[i], p2 = (ll)p[i] * p[i];
        for (int j = 1; p2 <= x; ++j, p1 = p2, p2 *= p[i])
            (res += (1ll * S(x / p1, i + 1) * (p[i] ^ j) % MOD + (p[i] ^ (j + 1))) % MOD) %= MOD;
    }
    return res;
}

int main(){
    n = read(); sqrtn = sqrt(n); getPrime(sqrtn);
    // zcfk -> caculate function g, init array w, id
    for (ll i = 1, j; i <= n; i = j + 1){
        w[++wcnt] = n / i, j = n / w[wcnt]; 
        (w[wcnt] <= sqrtn) ? id[0][w[wcnt]] = wcnt : id[1][j] = wcnt;
        pn[wcnt] = (w[wcnt] - 1) % MOD, g[wcnt] = (w[wcnt] % MOD) * ((w[wcnt] + 1) % MOD) % MOD;
        if (g[wcnt] & 1) g[wcnt] = g[wcnt] + MOD;
        g[wcnt] >>= 1; g[wcnt]--;
    }
    for (int j = 1; j <= pcnt; ++j)  
        for (int i = 1; i <= wcnt && p[j] * p[j] <= w[i]; ++i){
            int k = (w[i] / p[j] <= sqrtn) ? id[0][w[i] / p[j]] : id[1][n/(w[i] / p[j])];
            // minus the extra part
            (g[i] -= (ll)p[j] * (g[k] - prep[j - 1]) % MOD) %= MOD;
            // refresh pn
            (pn[i] -= pn[k] - j + 1) %= MOD;
        }
    // answer is S(n,1) plus f(1) = S(n,1) + 1
    int ans = S(n, 1) + 1;
    printf("%d\n", (ans + MOD) % MOD);
    return 0;
}