Min_25筛学习笔记

Min_25筛可以用来求解一类 (f(n)) 的前缀和，其中 (f) 是积性函数， (f(P)) 是关于质数 (P) 的低阶多项式，且 (forall kin N^+,f(P^k)) 可以通过某些方式快速求解。​

理论部分

首先我们求解 (sum_{i=1}^n [i in mathbb{P}]f(i)) 。

首先由于 (f(P)) 是多项式，所以我们可以把 (f(P)) 的每一项拆开来。

我们令 (g_{n,j} = sum_{i=1}^n [i in mathbb{P} space or space MinFactor(i) > P_j ]i^k) ，其中 (MinFactor(i)) 为 (i) 的最小质因子（后简写成 (MF)），(P) 为质数集，(P_j) 为第 (j) 个质数。由定义可以得到 (g_{n,|P|}) 即为我们所求, (g_{n,0} = sum_{i=2}^{n}i^k)。

那么我们该怎么转移呢？

我们分类考虑一下。如果 (P_j^2 > n) ，那么就不存在 (MF(i) > P_j) 这条限制了，毕竟不存在一个合数有大于 (sqrt n) 的质因子。

如果 (P_j^2 leq n) ，则由于 (g_{n,j}) 所囊括的被计算 (i^k) 的 (i) 必定是 (g_{n,j-1}) 所囊括的的一个子集，所以 (g_{n,j} = g_{n,j-1}-Delta)。

具体考虑 (Delta) 所包含的东西。容易发现 (g_{n,j}) 所没有包含的，而 (g_{n,j-1}) 却包含的就是最小质因子恰好等于 (P_j) 的。

[∴Delta = sum_{i=1}^n [MF(i)=p_j]i^k \ =sum_{i=1}^{lfloorfrac{n}{P_j} floor}[MF(i)geq p_j] i^k cdot P_j^k \=P_j^kcdotsum_{i=1}^{lfloorfrac{n}{P_j} floor}[MF(i)=p_j]i^k\=P_j^kcdot(g_{lfloorfrac{n}{P_j} floor,j-1}-sum_{i=1}^{j-1}[iin mathbb{P}]i^k) ]

把 (f) 的那些各个项合在一起，得到转移方程为：

[g_{n,j}=egin{cases} g_{n,j-1}&,P_j>sqrt n\ g_{n,j-1}-f(P_j) imes(g_{lfloor frac{n}{P_j} floor,j-1}-sum_{i=1}^{j-1} f(P_i))&,P_jle sqrt n end{cases} ]

其中最后的质数前缀和可以在筛质数的时候顺带处理，复杂度不计。

在整完了质数部分的答案后我们就可以合起来，考虑整体的答案了。

我们令(S_{n,j}=sum_{i=1}^n[MF(i)>P_j]f(i))，则由定义有(ans = S_{n,0}+f(1))，则

[S_{n,j} = (g_{n,|P|}-sum_{i=1}^j P_i) + (sum_{k=j+1}^{P_kleq sqrt n} sum_{e=1}^{P_k^e leq n}S_{n,k+1} + f(p_k^{e+1})) \ = (g_{n,|P|}-sum_{i=1}^j P_i) + (sum_{k=j+1}^{P_kleq sqrt n} sum_{e=1}^{P_k^e leq n}f(p_k^e)cdot S_{lfloorfrac{n}{p_k^e} floor,k+1} + f(p_k^{e+1})) ]

（第一个式子能倒腾到第二个是因为 (f(p_k^e)) 与 (S_{lfloorfrac{n}{p_k^e} floor,k+1}) 互质）

第二个式子第一个括号内的是前文中的质数部分的和，第二个括号内的则是和数部分的和。

第二个括号内的第一个 (Sigma) 相当于在枚举最小质因子，第二个则是在枚举最小质因子的指数。这个 (S_{lfloorfrac{n}{p_k^e} floor,k+1}) 是个递归，由于我们不可能从后面的东西推前面，所以我们枚举他的最小质因子和指数，然后只需要考虑除完之后剩下部分的答案就好了，即 (lfloorfrac{n}{p_k^e} floor) 。因为最小质因数已经被除完，所以剩下部分中不能再含有最小质因数，所以我们递归下去的部分要考虑的最小质因子门槛要提高，即 (k+1) 。

同时，为了补上所有被筛掉的 (p_k^e) 类的合数，我们要把他们的值加回去，即式子最后的 (f(p_k^{e+1}))) 。

然而，我们看到，光是推这个 (g) 的时空复杂度，貌似是 (O(n|P|)) 的，这很明显无法被接受 比暴力还慢的屑算法有什么存在的必要吗 ，所以说我们在实际的代码中肯定需要一些优化。

代码实现部分

首先我们看到 (g_{n,j}) 在 (jleq sqrt n) 的递推那里，即 (g_{n,j}=g_{n,j-1}-f(P_j) imes(g_{lfloor frac{n}{P_j} floor,j-1}-sum_{i=1}^{j-1} f(P_i))) 。我们可以看出，

1. 从 (j) 的角度看， (g_{n,j}) 都是从 (g_{x,j-1}) 转移来的；

2. 从 (n) 的角度看， (g_{n,j}) 都是从 (g_{lfloor frac{n}{P_j} floor, y}) 转移来的；

   这意味着什么呢？

从第一点看，我们可以压掉第二维数组。

从第二点看，我们不需要求出所有的 (g_{x}) ，只需要求 (g_{lfloor frac{n}{x} floor}) 。由整除分块的知识得到，这个东西只有可能有 (sqrt n) 种取值。于是时间复杂度降到了 (O(sqrt n imes |P|) < O(n)) （因为预处理的质数只用到 (sqrt n) 就行）。

但是由于 (n) 很大，所以我们光把下标变成 (lfloor frac{n}{x} floor) 不够，还要再离散化。具体的，我们把 (lfloor frac{n}{x} floor > sqrt n) 和 (lfloor frac{n}{x} floor leq sqrt n) 两种分成两类编号来“离散化”，这样子就把空间复杂度压下去了。

由于在 (g) 的转移方程中我们是把各个质因子拆开来想的，所以在代码里面也得这么做。具体到下面的代码就是(g1)和(g2)，分别代表一次项和二次项。

尽管洛谷的模板题不需要特殊考虑2这个唯一的偶质数，但是在有些题比如LOJ6053 简单的函数就要考虑到。

时间复杂度为 (O(n^{1-epsilon})) ，一说 (O(frac{n^{frac{3}{4}}}{log n})) ，反正我也不会证QAQ

代码（ 洛谷模板（P5325） ）

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
namespace ztd{
    using namespace std;
    typedef long long ll;
    template<typename T> inline T read(T& t) {//fast read
        t=0;short f=1;char ch=getchar();double d = 0.1;
        while (ch<'0'||ch>'9') {if (ch=='-') f=-f;ch=getchar();}
        while (ch>='0'&&ch<='9') t=t*10+ch-'0',ch=getchar();
        if(ch=='.'){ch=getchar();while(ch<='9'&&ch>='0') t+=d*(ch^48),d*=0.1,ch=getchar();}
        t*=f;
        return t;
    }
}
using namespace ztd;
const int maxn = 1e5+7;
const int mod = 1e9+7;
const ll INV6 = 166666668;//预处理6的逆元
ll n, sqrn;
inline ll f(ll x){//题面的f函数，其定义域只有prime^k
x %= mod;
    return x*(x-1)%mod;
}
ll prime[maxn], sump1[maxn], sump2[maxn], pcnt;
//sump1和sump2分别是质数的前缀和以及二次前缀和(2^2+3^2+5^2……)
bool vis[maxn];
inline void preprocess_prime(int uplimit){//线性筛筛质数
    for(int i = 2; i <= uplimit; ++i){
        if(!vis[i]){
            prime[++pcnt] = i;
            sump1[pcnt] = (sump1[pcnt-1]+i) % mod;
            sump2[pcnt] = (sump2[pcnt-1]+1ll*i*i%mod) % mod;
        } 
        for(int j = 1; j <= pcnt && i*prime[j] <= uplimit; ++j){
            vis[i*prime[j]] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}

ll g1[maxn<<1], g2[maxn<<1], match1[maxn<<1], match2[maxn<<1], w[maxn<<1], tot;
//g1,g2即上文所说， match1和match2和w就是上面所说的离散化了。
inline void preprocess_g(){//预处理g和离散化
    ll tmp;
    for(ll l = 1, r = 0; l <= n; l = r+1){
        r = n/(n/l); ++tot;
        w[tot] = n/l;
        if(w[tot] < sqrn) match1[n/l] = tot;
        else match2[n/(n/l)] = tot;
        tmp = w[tot] % mod;
        g1[tot] = tmp * (tmp+1) / 2 % mod-1;
        g2[tot] = tmp * (tmp+1) % mod * (tmp*2+1) % mod * INV6 % mod - 1;
    }
}
inline void process_g(){//g的递推
    for(int i = 1; i <= pcnt && prime[i]*prime[i] <= n; ++i){
        for(int j = 1; j <= tot && 1ll*prime[i]*prime[i] <= w[j]; ++j){
            ll pmt = w[j]/prime[i];
            ll tmp = (pmt >= sqrn) ? (match2[n/pmt]) : (match1[pmt]);
            g1[j] = (g1[j] - prime[i] * (g1[tmp]-sump1[i-1]+mod) % mod  % mod + mod) % mod;
            g2[j] = (g2[j] - prime[i]*prime[i]%mod * (g2[tmp]-sump2[i-1]+mod) % mod + mod) % mod;
        }
    }
}
ll S(ll N, ll y){//最后求答案的
    if(prime[y] > N) return 0;
    ll k = (N >= sqrn) ? (match2[n/N]) : (match1[N]);
    ll prans = (g2[k]-g1[k]+mod - (sump2[y-1]-sump1[y-1])+mod) % mod;//指质数方面的结果
    ll cpans = 0;//合数方面的结果
    for(int i = y; i <= pcnt && 1ll * prime[i]*prime[i] <= N; ++i){
        for(ll pk = prime[i]; prime[i]*pk <= N; pk *= prime[i]){
            cpans = (cpans + (f(pk) * S(N/pk, i+1) % mod + f(prime[i]*pk)) % mod) %mod;
        }
    }
    return (prans+cpans) % mod;
}
signed main(){
    read(n);  sqrn = sqrt(n);
    preprocess_prime(maxn-5);
    preprocess_g();
    process_g();
    cout << (S(n,1) + 1) % mod << '
';//这道题定义了f(1)=1
    return 0;
}