min_25筛

发表于 更新于 分类于

Min_25筛简介

$\text{min_25}$筛是一种处理一类积性函数前缀和的算法。

其中这类函数$f(x)$要满足$\sum_{i=1}^{n}[i\in prime]\cdot f(i)$可以被$\sum_{i=1}^{n}[i\in prime]\cdot i^k$简单表示或者快速计算，其中$k$为较小的常数。

时间复杂度好像是$O(\frac{n^{0.75}}{\log n})$，不过据说被证伪了…也有人说是$O(n^{1-\epsilon})$，反正可以做$n=1e10$，貌似跑的比杜教筛快。

空间复杂度$O(\sqrt{n})$。

part 1

首先我们要处理的问题是$\sum_{i=1}^{n}[i\in prime]\cdot f(i)$，由于函数特性转化为了求$\sum_{i=1}^{n}[i\in prime]\cdot i^k$。

设$g(n,i)=\sum_{i=1}^{n}[i\in prime\ or \ t(i)>p_i]\cdot i^k$，其中$t(x)$表示$x$的最小质因子，$p_i$表示第$i$个质数。

直观理解就是埃拉托斯特尼筛法进行了$i$轮之后$1\sim n$剩下的数的$k$次方之和。

显然当$p_i^2>n$时$g(n,i)=g(n,i-1)$，下面只考虑$p_i^2\leqslant n$。

我们考虑$g(n,i)$如何从$g(*,i-1)$转移过来，也就是考虑第$i$轮筛掉了那些数。

那么可以得到式子：

其中$s_i=\sum_{j=1}^{i}p_i^k$。

解释一下这个式子，我们考虑筛掉了哪些数，然后减掉就好了，筛掉的那些数最小的质因子显然是$p_i$，而正好$g(\lfloor\frac{n}{p_i}\rfloor,i-1)$代表最小质因子$\geqslant p_i$的数以及质数，我们把质数减掉，拼起来就是$1\sim n$。

那么写在一起就是：

根据常识可以知道第一维只有$2\sqrt n$个取值，离散化一下然后递推就好了。

注意到这里我们顺便求出了对于每个$x$,$\sum_{i=1}^{\lfloor n/x\rfloor}[i\in prime]\cdot i^k$的值。

part 2

我们想算出$f$的前缀和，同样的我们设$h(n,i)$如下：

那么$h(n,1)+f(1)$就是答案，注意我们$h$不会把$1$算进去。

首先还是很显然的，当$p_i>n$时，$h(n,i)=0$，同样下面只考虑$p_i\leqslant n$的情况。

我们可以通过$g$很简单的算出质数的贡献，即：

这些都可以预处理。

然后我们考虑暴力地枚举最小的质因子是多少以及用了多少个，即：

解释下后面的$f(p_j^{k+1})$，这是因为$h$没有把$1$算进去，所以会漏算质数的次幂，直接加上就行了。

综合起来就是：

写成这样求和符号好像有点丑…不管它

直接暴力递归算就好了，递推好像还慢些。

code

[LOJ6053] 简单的函数

$f(p)=p-1$，$p=2$的情况特殊处理下就好了。

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
#include<bits/stdc++.h>
using namespace std;

void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
 
void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}

#define lf double
#define ll long long 

#define pii pair<int,int >
#define vec vector<int >

#define pb push_back
#define mp make_pair
#define fr first
#define sc second

#define FOR(i,l,r) for(int i=l,i##_r=r;i<=i##_r;i++) 

const int maxn = 1e6+10;
const int inf = 1e9;
const lf eps = 1e-8;
const int mod = 1e9+7;
const int inv2 = 5e8+4;

int add(int x,int y) {return x+y>=mod?x+y-mod:x+y;}
int del(int x,int y) {return x-y<0?x-y+mod:x-y;}
int mul(int x,int y) {return 1ll*x*y-1ll*x*y/mod*mod;}

ll n,w[maxn];
int b,tot,pri[maxn],vis[maxn],id[maxn],id2[maxn],m,g[maxn],h[maxn],f[maxn],p[maxn];

void sieve() {
    for(int i=2;i<=b<<1;i++) {
        if(!vis[i]) pri[++tot]=i,p[tot]=add(p[tot-1],pri[tot]);
        for(int j=1;j<=tot&&i*pri[j]<=b<<1;j++) {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        }
    }
}

int get_id(ll x) {return x<=b?id[x]:id2[n/x];}

void get_g() {
    for(ll l=1;l<=n;) {
        ll v=n/l,r=n/v;
        if(v<=b) id[v]=++m;else id2[n/v]=++m;  //利用整除分块来离散化
        w[m]=v,v%=mod;g[m]=mul(v+2,mul(v-1,inv2)),h[m]=del(v,1);l=r+1;  
    }
    for(int i=1;pri[i]<b;i++)
        for(int j=1;j<=m;j++) {
            if(1ll*pri[i]*pri[i]>w[j]) break;   //剪枝
            g[j]=del(g[j],mul(pri[i],del(g[get_id(w[j]/pri[i])],p[i-1]))); //g函数记录的是\sum [i\in pri] i，其中第二位滚动掉了
            h[j]=del(h[j],del(h[get_id(w[j]/pri[i])],i-1));     //h记录的是\sum [i\in pri] 1
        }
    for(int i=1;i<=m;i++) f[i]=del(g[i],h[i]);
}

int calc(ll x,int i) {
    if(x<=1||pri[i]>x) return 0;
    int res=del(f[get_id(x)],del(p[i-1],i-1));if(i==1) res=add(res,2);  //这里处理了f(2)=2+1的情况
    for(int j=i;1ll*pri[j]*pri[j]<=x;j++)
        for(ll k=1,e=pri[j];e*pri[j]<=x;e*=pri[j],k++)
            res=add(res,add(mul(pri[j]^k,calc(x/e,j+1)),pri[j]^(k+1))); //直接暴力算
    return res;
}

int main() {
    scanf("%lld",&n);b=ceil(sqrt(n));
    sieve();get_g();write(add(calc(n,1),1));  //注意加 1
    return 0;
}

洛谷P5325 【模板】Min_25筛

这和上题也是一样的，注意别爆$\text{long long}$。

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;

void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
 
void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}

#define lf double
#define ll long long 

#define pii pair<int,int >
#define vec vector<int >

#define pb push_back
#define mp make_pair
#define fr first
#define sc second

#define FOR(i,l,r) for(int i=l,i##_r=r;i<=i##_r;i++) 

const int maxn = 1e6+10;
const int inf = 1e9;
const lf eps = 1e-8;
const int mod = 1e9+7;
const int inv2 = 5e8+4;
const int inv6 = 166666668;

int add(int x,int y) {return x+y>=mod?x+y-mod:x+y;}
int del(int x,int y) {return x-y<0?x-y+mod:x-y;}
int mul(int x,int y) {return 1ll*x*y-1ll*x*y/mod*mod;}

ll n,w[maxn];
int b,tot,pri[maxn],vis[maxn],id[maxn],id2[maxn],m,g[maxn],h[maxn],f[maxn],p[maxn],p2[maxn];

void sieve() {
    for(int i=2;i<=b<<1;i++) {
        if(!vis[i]) pri[++tot]=i,p[tot]=add(p[tot-1],pri[tot]),p2[tot]=add(p2[tot-1],mul(i,i));
        for(int j=1;j<=tot&&i*pri[j]<=b<<1;j++) {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        }
    }
}

int get_id(ll x) {return x<=b?id[x]:id2[n/x];}

void get_g() {
    for(ll l=1;l<=n;) {
        ll v=n/l,r=n/v;
        if(v<=b) id[v]=++m;else id2[n/v]=++m;
        w[m]=v,v%=mod;
        g[m]=del(mul(mul(v,mul(v+1,(2*v+1)%mod)),inv6),1);
        h[m]=del(mul(v,mul(v+1,inv2)),1);l=r+1;
    }
    for(int i=1;pri[i]<b;i++)
        for(int j=1;j<=m;j++) {
            if(1ll*pri[i]*pri[i]>w[j]) break;
            g[j]=del(g[j],mul(mul(pri[i],pri[i]),del(g[get_id(w[j]/pri[i])],p2[i-1])));  //sum of square of prime
            h[j]=del(h[j],mul(pri[i],del(h[get_id(w[j]/pri[i])],p[i-1]))); //sum of prime
        }
    for(int i=1;i<=m;i++) f[i]=del(g[i],h[i]); 
}

int calc(ll x,int i) {
    if(x<=1||pri[i]>x) return 0;
    int res=del(f[get_id(x)],del(p2[i-1],p[i-1]));
    for(int j=i;1ll*pri[j]*pri[j]<=x;j++)
        for(ll k=1,e=pri[j],ee=e;e*pri[j]<=x;e*=pri[j],k++,ee=e%mod) //注意e过大可能会爆，先模一下就好了
            res=add(res,add(mul(mul(ee,ee-1),calc(x/e,j+1)),mul(mul(ee,pri[j]),del(mul(ee,pri[j]),1))));
    return res;
}

int main() {
    scanf("%lld",&n);b=ceil(sqrt(n));
    sieve();get_g();write(add(calc(n,1),1));
    return 0;
}