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
\(f(p)=p-1\)，\(p=2\)的情况特殊处理下就好了。
#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;
}
这和上题也是一样的，注意别爆\(\text{long long}\)。
#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;
}
