Min 25 筛

Min 25 筛是用来求出积性函数前缀和的方法。

如果你不知道什么是积性函数：

  • 积性函数：∀a⊥b:f(ab)=f(a)f(b)\forall a\perp b:f(ab)=f(a)f(b)；
  • 完全积性函数：∀a,b:f(ab)=f(a)f(b)\forall a,b:f(ab)=f(a)f(b)。

本文中，若无特殊说明，pp 表示质数。

第一步 - 求出函数质数部分的前缀和

做法

此步我们要求 ff 是完全积性函数，在题目中 f(p)f(p) 往往是关于 pp 的多项式，这样就可以拆成若干个完全积性函数，也就是说，我们期望 f(p)f(p) 的表达式较简单。

即求 g(n)=∑p≤nf(p)g(n) = \sum_{p\le n} f(p)。
设 pjp_j 为第 jj 个素数（即p1=2p_1=2），特别地，令 p0=1p_0=1。

设 g(n,j)g(n, j) = ∑i=1n[i是质数∨i 的最小质因子>pj]f(i)\sum_{i=1}^n [i\text{是质数}\vee i\text{ 的最小质因子}>p_j]f(i)。

考虑从 g(n,j−1)g(n, j-1) 递推到 g(n,j)g(n,j)，就是要减去最小质因子为 pjp_j 的合数的贡献。

考虑所有 pjp_j 的倍数，f(pj)⋅g(⌊npj⌋,j−1)f(p_j)\cdot g\left(\left \lfloor \frac{n}{p_j}\right \rfloor, j-1 \right)，表示两个部分：

  1. pjp_j 的倍数，且剩下的质因子都 ≥pj\ge p_j 的贡献；
  2. pjp_j 的倍数，且剩下的部分是一个 <pj<p_j 的质数的贡献。

需要把第二部分的再减掉，这一项就是

f(pj)⋅g(pj−1,j−1)f(p_j)\cdot g\left(p_{j-1}, j-1\right)

因此，能得到递推式：

g(n,j)=g(n,j−1)−f(pj)(g(⌊npj⌋,j−1)−g(pj−1,j−1))g(n,j) = g(n, j-1) - f(p_j)\left(g\left(\left \lfloor \dfrac{n}{p_j}\right \rfloor, j-1 \right) - g\left(p_{j-1}, j-1\right)\right)

可以发现，上述式子只会用到在 ⌊nx⌋\left \lfloor \frac{n}{x} \right \rfloor 处的 dp 值，那么我们知道，这样的数个数是 O(n)O(\sqrt n) 的（具体来说，小于等于 2n2\sqrt n）。

注意到我们还用到了小于等于 n\sqrt{n} 的质数处的值，但事实上，可以证明 ∀1≤i≤n\forall 1\le i\le \sqrt n，都存在 xx，使得 ⌊nx⌋=i\left\lfloor \frac n x \right\rfloor = i。

证明：
令 x=⌊ni⌋x = \left\lfloor \frac n i \right\rfloor，下证 ⌊nx⌋=i\left\lfloor \frac n x \right\rfloor = i。

由取整函数定义，有 x≤ni<x+1x\le \frac n i <x+1，可得 i≤nxi\le \frac n x。
只需证明 nx<i+1\frac n x < i+1。

若 nx≥i+1\frac n x \ge i+1，则 x≤ni+1<ni<x+1x \le \frac n {i+1} < \frac n i < x+1。
得 x(i+1)≤n<i(x+1)x(i+1)\le n < i(x+1)，得到 x<ix<i，即 ⌊ni⌋<i\left\lfloor \frac n i \right\rfloor < i。
则 ni<i\frac n i<i，i2>ni^2 >n 矛盾！□\square

因此不需要考虑求出这些位置的 gg。

我们求出这些数并离散化，顺带求出 g(n,0)g(n,0)，

复杂度为我们用到的所有 g(x,⋅)g(x,\cdot) 中的 π(x)\pi(\sqrt x) 之和，即 O(n34log⁡n)O\left(\dfrac{n^{\frac 3 4}}{\log n}\right)。

例题

LOJ#6235. 区间素数个数

求 1∼n1\sim n 之间素数的个数，2≤n≤10112\le n\le 10^{11}。

由上文我们只要求出 f(x)=1f(x)=1 函数质数部分的前缀和就能得到答案。

在具体实现方面注意几点：

  1. 在离散化的时候，我们不希望引入 log⁡\log，可以用两个数组 id1 和 id2 来储存，若 x≤nx\le \sqrt n，那么就存在 id1[x] 中，否则存在 id2[n/x] 中。
  2. 在递推这个 g(n,j)g(n,j) 的过程中，直接滚动数组，按照 jj 从小到大枚举。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include<bits/stdc++.h>
using namespace std;

using ll=long long;
const int N = 4e5+5;
int id1[N],id2[N];

int p[N],tot;
bool vis[N];
void init(){ // 筛出前 sqrt(n) 的质数
    for(int i=2;i<N;i++){
        if(!vis[i]){
            p[++tot] = i;
        }

        for(int j=1;j<=tot && p[j]*i<N;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0)break;
        }
    }
};

ll v[N*2], g[N*2];
int main(){
    ll n; cin>>n;
    int m=0;

    init();
    for(ll l=1,r;l<=n;l=r+1){ // 数论分块
        r = n/(n/l); 
        v[++m] = n/l; 
        if(v[m]<N)id1[v[m]]=m;
        else id2[n/v[m]]=m;

        g[m] = v[m]-1;
    }
    
    auto get = [&](ll x){
        return x<N?id1[x]:id2[n/x];
    };

    for(int j=1;j<=tot;j++){ // 递推求出所有的 g
        for(int i=1;i<=m&&p[j]<=v[i]/p[j];i++){
            g[i] -= g[get(v[i]/p[j])] - g[get(p[j-1])];
        }
    }

    cout << g[get(n)] << '\n';
    return 0;
}

练习：https://qoj.ac/contest/1871/problem/9867。

第二步 - 求出函数的前缀和

此步我们就不要求 ff 是完全积性函数了，但是我们仍需能够快速计算出 f(pk)f(p^k)。

做法一

设 S(n,j)S(n, j) = ∑i=2n[i 的最小质因子>pj]f(i)\sum_{i=2}^n [i\text{ 的最小质因子}>p_j]f(i)。

最后答案就是 S(n,0)+f(1)S(n,0) + f(1)。

把 S(n,j)S(n,j) 分成质数部分和合数部分进行计算。

S(n,j)=g(n)−g(pj)+∑j<k,pk≤n,1≤e,pke≤nf(pke)(S(⌊npke⌋,k)+[e≠1])S(n,j) = g(n) - g(p_j) + \sum_{j<k,p_k\le\sqrt n, 1\le e,p_k^e\le n }f(p_k^e)\left ( S\left ( \left \lfloor \dfrac{n}{p_k^e} \right \rfloor,k \right ) + [e\ne 1] \right )

直接暴力计算即可。

这个式子可以修改成

S(n,j)=g(n)−g(pj)+∑j<k,pk≤n,1≤e,pke+1≤n(f(pke)S(⌊npke⌋,k)+f(pke+1))S(n,j) = g(n) - g(p_j) + \sum_{j<k,p_k\le\sqrt n, 1\le e,p_k^{e+1}\le n }\left( f(p_k^e) S\left ( \left \lfloor \dfrac{n}{p_k^e} \right \rfloor,k \right ) + f(p_k^{e+1}) \right)

正确性显然。

例题：
P5325 【模板】Min_25 筛
求积性函数 f(pk)=pk(pk−1)f(p^k) = p^k(p^k-1) 的前缀和。n≤1010n\le 10^{10}。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#include<bits/stdc++.h>
using namespace std;

using ll=long long;

const ll mod = 1e9+7;
constexpr ll ksm(ll a,ll b){
    ll r = 1;
    while(b){
        if(b&1)r=r*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return r;
}
constexpr ll inv2 = ksm(2, mod-2);
constexpr ll inv6 = ksm(6, mod-2);

constexpr ll sq(ll n){n%=mod;return n*n%mod;}
constexpr ll sum(ll n){n%=mod;return n*(n+1)%mod*inv2%mod;}
constexpr ll sum_sq(ll n){n%=mod;return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;}

const int N = 2e5+5;


ll p[N],tot; bool vis[N];
void init_sieve(){
    for(int i=2;i<N;i++){
        if(!vis[i])p[++tot]=i;
        for(int j=1;j<=tot && i*p[j]<N;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0)break;
        }
    }
}

ll v[N],m;
int id1[N],id2[N];
ll g1[N],g2[N];
ll n;
ll get(ll x){
    return x<N?id1[x]:id2[n/x];
}
void init_sqrt(){
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        v[++m]=n/l;
        if(v[m]<N)id1[v[m]]=m;
        else id2[n/v[m]]=m;
        g1[m]=(sum_sq(v[m])+mod-1)%mod, g2[m]=(sum(v[m])+mod-1)%mod;
    }
}

void calc_g(){
    for(int j=1;j<=tot;j++){
        for(int i=1;i<=m && p[j]<=v[i]/p[j];i++){
            g1[i] = (g1[i] - sq(p[j])*g1[get(v[i]/p[j])]%mod + sq(p[j])*g1[get(p[j-1])] + mod) % mod;
            g2[i] = (g2[i] - p[j]*g2[get(v[i]/p[j])]%mod + p[j]*g2[get(p[j-1])] + mod) % mod;
        }
    }
}

ll f(ll x){x%=mod;return (x*x%mod-x+mod)%mod;}
ll S(ll x,int j){
    if(p[j] >= x) return 0;
    ll res = (g1[get(x)]-g2[get(x)]-g1[get(p[j])]+g2[get(p[j])]+mod+mod)%mod;
    for(int k=j+1;k<=tot && p[k]<=x/p[k];k++){
        ll w = p[k];
        for(;w<=x/p[k];w*=p[k]){
            res = (res + f(w) * S(x/w, k)%mod + f(w*p[k]))%mod;
        }
    }
    return res;
}   

int main(){
    ios::sync_with_stdio(0);cin.tie(0);
    cin>>n;
    init_sieve();
    init_sqrt();
    calc_g();

    cout << ((S(n, 0) + 1)%mod+mod)%mod << '\n';
    return 0;
}

做法二

考虑继承第一步的状态设计。

设 S(n,j)S(n, j) = ∑i=2n[i是质数∨i 的最小质因子>pj]f(i)\sum_{i=2}^n [i\text{是质数}\vee i\text{ 的最小质因子}>p_j]f(i)。

写出递推式：

S(n,j)=S(n,j+1)+∑pj+1≤n,1≤e,pj+1e+1≤nf(pj+1e)(S(⌊npj+1e⌋,j+1)−g(pj+1))+f(pj+1e+1)S(n,j) = S(n,j+1) + \sum_{p_{j+1}\le\sqrt n,1\le e,p_{j+1}^{e+1}\le n}f(p_{j+1}^e)\left(S\left ( \left \lfloor \dfrac{n}{p_{j+1}^e} \right \rfloor ,j+1 \right ) -g(p_{j+1})\right) + f(p_{j+1}^{e+1})

和第一步中一样，滚动数组 dp 即可，边界条件 S(x,⋅)=g(x)S(x,\cdot)=g(x)。

注意到做法二可以求出所有 S(⌊nx⌋)S(\left\lfloor\frac n x\right\rfloor)，这是做法一无法完成的。

例题：
P5325 【模板】Min_25 筛
求积性函数 f(pk)=pk(pk−1)f(p^k) = p^k(p^k-1) 的前缀和。n≤1010n\le 10^{10}。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include<bits/stdc++.h>
using namespace std;

using ll=long long;

const ll mod = 1e9+7;
constexpr ll ksm(ll a,ll b){
    ll r = 1;
    while(b){
        if(b&1)r=r*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return r;
}
constexpr ll inv2 = ksm(2, mod-2);
constexpr ll inv6 = ksm(6, mod-2);

constexpr ll sq(ll n){n%=mod;return n*n%mod;}
constexpr ll sum(ll n){n%=mod;return n*(n+1)%mod*inv2%mod;}
constexpr ll sum_sq(ll n){n%=mod;return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;}

const int N = 2e5+5;


ll p[N],tot; bool vis[N];
void init_sieve(){
    for(int i=2;i<N;i++){
        if(!vis[i])p[++tot]=i;
        for(int j=1;j<=tot && i*p[j]<N;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0)break;
        }
    }
}

ll v[N],m;
int id1[N],id2[N];
ll g1[N],g2[N];
ll n;
ll get(ll x){
    return x<N?id1[x]:id2[n/x];
}
void init_sqrt(){
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        v[++m]=n/l;
        if(v[m]<N)id1[v[m]]=m;
        else id2[n/v[m]]=m;
        g1[m]=(sum_sq(v[m])+mod-1)%mod, g2[m]=(sum(v[m])+mod-1)%mod;
    }   
}

void calc_g(){
    for(int j=1;j<=tot;j++){
        for(int i=1;i<=m && p[j]<=v[i]/p[j];i++){
            g1[i] = (g1[i] - sq(p[j])*g1[get(v[i]/p[j])]%mod + sq(p[j])*g1[get(p[j-1])] + mod) % mod;
            g2[i] = (g2[i] - p[j]*g2[get(v[i]/p[j])]%mod + p[j]*g2[get(p[j-1])] + mod) % mod;
        }
    }
}

ll f(ll x){x%=mod;return (x*x%mod-x+mod)%mod;}
ll g[N],S[N];

int main(){
    ios::sync_with_stdio(0);cin.tie(0);
    cin>>n;
    init_sieve();
    init_sqrt();
    calc_g();

    for(int i=1;i<=m;i++)S[i]=g[i]=(g1[i]-g2[i]+mod)%mod;
    for(int j=tot-1;j>=0;j--){
        for(int i=1;i<=m && p[j+1]<=v[i]/p[j+1];i++){
            for(ll w=p[j+1];w<=v[i]/p[j+1];w*=p[j+1]){
                S[i] = (S[i] + f(w)*(S[get(v[i]/w)]-g[get(p[j+1])]+mod)%mod + f(w*p[j+1]))%mod;
            }
        }
    }
    cout << S[get(n)]+1 << '\n';
    return 0;
}

练习

LOJ#6784. 简单的函数 10^12 版

给定积性函数 f(pk)=p⊕kf(p^k)=p\oplus k（f(1)=1f(1)=1），S(m)∑i=1mf(i) mod (109+7)S(m)\sum_{i=1}^m f(i) \bmod (10^9 + 7)。
给定 nn，求 ⨁{S(⌊n/i⌋)∣i∈[n]}\bigoplus \{S(\lfloor n/i \rfloor)\mid i\in[n]\}。（注意：）

n≤1012n\le 10^{12}。

因为这个函数是积性函数，且 f(pk)f(p^k) 是很容易计算的，第二步照抄做法二即可。

第一步在于快速求出质数函数值的前缀和，可以发现

f(p)={p+1p=2p−1elsef(p)=\begin{cases} p+1 &p=2 \\ p-1 &\text{else} \end{cases}

那么，把 f(p)f(p) 当成 p−1p-1，然后拆成 pp 和 11 两个完全积性函数去计算即可，最后求出的 gg 函数在处理一下 f(2)f(2) 少算了 22 的问题。

于是就做完了！

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include<bits/stdc++.h>
using namespace std;

using ll=long long;

const int N = 2e6+10;


const ll mod = 1e9+7;
constexpr ll ksm(ll a,ll b){
    ll r = 1;
    while(b){
        if(b&1)r=r*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return r;
}
constexpr ll inv2 = ksm(2, mod-2);
constexpr ll sum(ll n){n%=mod;return n*(n+1)%mod*inv2%mod;}

ll n;

ll p[N],tot;
bool vis[N];

ll id1[N],id2[N],v[N],m;
ll get(ll x){return x<N?id1[x]:id2[n/x];}

ll g1[N],g2[N],g[N];

ll f(ll P,ll K){return (P^K)%mod;}
ll S[N];
int main(){
    for(int i=2;i<N;i++){
        if(!vis[i])p[++tot]=i;
        for(int j=1;j<=tot && i*p[j]<N;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0)break;
        }
    }

    cin>>n;

    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        v[++m]=n/l;
        if(v[m]<N)id1[v[m]]=m;
        else id2[n/v[m]]=m;

        g1[m]=(sum(v[m])+mod-1)%mod, g2[m]=(v[m]-1+mod)%mod;
    }
    for(int j=1;j<=tot;j++){
        for(int i=1;i<=m && p[j]<=v[i]/p[j];i++){
            g1[i] = (g1[i] - p[j]*(g1[get(v[i]/p[j])]-g1[get(p[j-1])]+mod)%mod + mod) % mod;
            g2[i] = (g2[i] - (g2[get(v[i]/p[j])] - g2[get(p[j-1])])%mod + mod) % mod;
        }

    }

    for(int i=1;i<=m;i++)g[i]=(g1[i]-g2[i]+mod + (v[i]>=2?2:0))%mod;

    for(int i=1;i<=m;i++)S[i]=g[i];
    for(int j=tot-1;j>=0;j--){
        for(int i=1;i<=m && p[j+1]<=v[i]/p[j+1];i++){
            for(ll w=p[j+1],k=1;w<=v[i]/p[j+1];w*=p[j+1],++k){
                S[i] = (S[i] + f(p[j+1],k)*(S[get(v[i]/w)]-g[get(p[j+1])]+mod)%mod + f(p[j+1],k+1))%mod;
            }
        }
    }

    for(int i=1;i<=m;i++)S[i]=(S[i]+1)%mod;
    sort(S+1,S+1+m);
    m=unique(S+1,S+1+m)-S-1;  

    ll ans = 0;
    for(int i=1;i<=m;i++)ans^=S[i];
    cout << ans << '\n';
    return 0;
}

参考资料