Min_25 筛

不会真的有人在考了 $\sout{n}$ 次后还不会一个算法吧 不会吧不会吧

老早就学了这玩意，但是几次考试板子也写不对。于是近期重新梳理了一下算法。

算法流程

假设现在有一个函数 $f(i)$ ，我们要求下式的值：

$\sum_{i=1}^{n}f(i)$

构造函数 $F(i)$，满足如下要求：

• 方便求前缀和
• 在质数处的取值与 $f(i)$ 相同
• 完全积性函数

设 $p_i$ 表示第 $i$ 个质数。

设 $G(n,i)$ 表示「前 $n$ 个数中，质数 或者 最小质因子大于 $p_i$ 的数的 $F(i)$ 之和」。注意 $G$ 的求和不包含 $1$ 。

换一个更加直观的说法：在埃式筛中，第 $i$ 轮后仍然未被筛去的数的 $F(i)$ 之和。

可以写出转移方程：

$G(n,i)= \begin{cases} G(n,i-1)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ p_i^2>n\\ G(n,i-1)-F(p_i)(G(\lfloor\frac{n}{p_i}\rfloor,i-1)-\sum\limits_{j=1}^{i-1}F(p_j))\ \ \ \ p_i^2\leq n \end{cases}$

首先，第 $i$ 轮筛后仍未被筛去的数的 $F(i)$ 之和一定是第 $i-1$ 轮后的一部分，我们先继承过来再考虑减去多算的。可以发现，第 $i$ 轮筛一定至少从 $p_i^2$ 开始，所以若 $p_i^2>n$ 那么什么都不会发生；否则，我们从所有应该被筛去的数的函数值中提出 $F(p_i)$ （注意：此时用到了「完全积性函数」的性质），然后剩下的数的最小质因子都必须大于 $p_{i-1}$ ，所以它们的函数值和在 $G(\lfloor\frac{n}{p_i}\rfloor,i-1)$ 中。但是 $G(\lfloor\frac{n}{p_i}\rfloor,i-1)$ 中还多算了比 $p_i$ 小的质数的函数值之和，所以这部分我们也得减去。

考虑完转移，我们来考虑如何设置初值。显然，$G(n,0)$ 就相当于什么都没有筛，也就是 $F(i)$ 的前缀和。所以我们需要这个函数「方便求前缀和」。

设 $\text{limit}$ 表示 $\sqrt{n}$ 以内的所有质数， 那么由于保证了「$F(i)$ 在质数处的取值与 $f(i)$ 相同」， $G(n,\text{limit})$ 就是 $n$ 以内的所有质数的 $f(i)$ 的和。上述一大堆的工作，目的皆在于此——对于任意 $n$ ，求出 $n$ 以内的质数的 $f(i)$ 之和。

下面我们来考虑求所有合数的值。设 $S(n,i)$ 表示「前 $n$ 个数中，最小质因子大于等于 $p_i$ 的合数的 $f(i)$ 之和」。同样的，$S$ 的求和中是不包含 $1$ 的。

可以写出转移方程：

$S(n,i)= S(n,i+1)+\sum\limits_{e=1}^{p_i^{e+1}\leq n}f(p_i^e)(S(\lfloor\frac{n}{p_i^e}\rfloor,i+1)+G(\lfloor\frac{n}{p_i^e}\rfloor,\text{limit})-\sum_{j=1}^{i}f(p_j))+f(p_i^{e+1})$

首先，由于 $f(i)$ 仅仅满足是积性函数，所以我们还需要钦定 $p_i$ 的次数 $e$ ，把 $p_i$ 这个质因子完全消去（此时我们需要 $f(i)$ 满足是积性函数）。考虑一个合数除掉 $p_i^e$ 还剩什么：质数；合数；$1$ 。

我们先不管上式中 $p_i^{e+1}\leq n$ 的上界限制，考虑该怎么简便准确地求解。

首先，对于合数，新增的值就是 $f(p_i^e)S(\lfloor\frac{n}{p_i^e}\rfloor,i+1)$ 。因为此处有 $\lfloor\frac{n}{p_i^e}\rfloor$ ，所以上界可以设为 $p_{i}^{e+1}\leq n$ 。

然后，对于质数，新增的值就是 $f(p_i^e)\sum\limits_{j=i+1}^{\text{limit}}f(p_j)[p_jp_i^{e}\leq n]$ （注意这个方括号是不可或缺的！）。我们转化为 $G(\lfloor\frac{n}{p_i^e}\rfloor,\text{limit})-\sum\limits_{j=1}^{i}f(p_j)$ 。此时的细节为每一次要保证 $\lfloor\frac{n}{p_i^e}\rfloor\geq p_i$ （因为后面要减去 $f(p_j)$）。于是可以得到上界为 $p_i^{e+1}\leq n$ 。

上界统一自然是最好。所以对于 $1$ ，如果我们加在里面，我们不仅上界要改成 $p^{e}\leq n$ ，而且还要减去 $f(p)$ 。所以我们干脆提出来，变成 $f(p_i^{e+1})$ 。于是上界都统一了。

考虑完转移，我们来考虑如何设置初值。显然，$S(n,\text{limit}+1)$ 相当于什么也没有。于是我们不用赋初值。

最后 $S(n,1)+G(n,\text{limit})+f(1)$ 就是所求。

一般而言，$\rm Min\_25$ 筛的题目都是跑 $10^{10}$ 的，多一个 $10$ 会相当吃力。

实现细节

1.可以发现，如果我们直接实现，$S$ 和 $G$ 的第一维的下标都达到了 $n$ ，显然不行。我们考虑我们在整个算法流程中只用到了 $\forall i,\lfloor\frac{n}{i}\rfloor$ 的值，而这只有 $\sqrt{n}$ 个，所以我们提前预处理这 $\sqrt{n}$ 个位置，下标的范围就控制在了 $\sqrt{n}$ 级别。

2.求 $G$ 和 $S$ 中都对上界有明确要求。首先，这个上限是必须要满足的，否则可能我们在 $G(\lfloor\frac{n}{p_i}\rfloor,i-1)$ 或者 $G(\lfloor\frac{n}{p_i^e}\rfloor,\text{limit})$ 中并没有加上算某个 $p_j$ ，而之后去除质数的贡献的时候又把 $p_j$ 去掉了。除了正确性的保证，我们还可以发现，当 $p_i^2>n$ 时，两者都可以直接 break ，而这一步也是保证 $\rm Min\_25$ 筛的复杂度的关键。

3.检验时可以带入 $1$ 试试。

代码实现

例题是 $\rm LOJ6053$ 简单的函数

const int N = 1e5 + 23, M = 2e5 + 23, mod = 1e9 + 7, inv2 = (mod + 1) / 2;
using namespace std;
int m, num, valcnt;
long long n;
int ok[N], pr[N], g0[M], g1[M], s[M], id0[M], id1[M];
long long val[M];
inline void Plus(int &a, int b){
    a += b - mod, a += (a >> 31) & mod;
}
void Sieve(int n){
    for(int i = 2; i <= n; ++i){
        if(!ok[i]) pr[++num] = i;
        for(int j = 1; j <= num and pr[j] * i <= n; ++j){
            ok[i * pr[j]] = 1;
            if(i % pr[j] == 0) break;
        }
    }
}
inline int &pos(long long v){
    return v <= m ? id0[v] : id1[n / v];
}
int main(){

    fin >> n;

    if(n == 1){
        puts("1");
        return 0;
    }

    Sieve(m = sqrt(n));
    
    for(long long i = 1, now, las; i <= n; i = las + 1){
        now = n / i, las = n / now;
        val[++valcnt] = now;
        pos(now) = valcnt;
        g0[valcnt] = (now - 1) % mod;  
        g1[valcnt] = 1ll * (now % mod + 2) % mod * (now % mod - 1) % mod * inv2 % mod;
    }

    int sp0 = 0, sp1 = 0;
    for(int i = 1; i <= num; ++i){
        int p = pr[i];
        for(int j = 1; j <= valcnt; ++j){
            long long vl = val[j] / p;
            if(vl < p) break;

            int po = pos(vl);
            Plus(g0[j], (mod - g0[po] + sp0) % mod);
            Plus(g1[j], 1ll * p * (mod - g1[po] + sp1) % mod);
        }
        sp0 += 1, Plus(sp1, p);
    }

    for(int i = 1; i <= valcnt; ++i)
        Plus(g1[i], mod - g0[i] + 2);
    
    Plus(sp1, 2);
    for(int i = num; i >= 1; --i){
        int p = pr[i]; 
        for(int j = 1; j <= valcnt; ++j){
            long long vl = val[j] / p;
            if(vl < p) break;

            for(int e = 1; vl >= p; vl /= p, ++e){
                int po = pos(vl);
                Plus(s[j], 1ll * (p ^ e) * (s[po] + ((g1[po] + mod - sp1) % mod + sp0) % mod) % mod);
                Plus(s[j], p ^ (e + 1));
            }
        }
        sp0 -= 1, Plus(sp1, mod - p);
    }

    cout << (g1[1] + s[1] + 1) % mod << endl;
    return 0;
}