Min_25筛简介
( ext{min_25})筛是一种处理一类积性函数前缀和的算法。
其中这类函数(f(x))要满足(sum_{i=1}^{n}[iin prime]cdot f(i))可以被(sum_{i=1}^{n}[iin 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}[iin prime]cdot f(i))，由于函数特性转化为了求(sum_{i=1}^{n}[iin prime]cdot i^k)。
设(g(n,i)=sum_{i=1}^{n}[iin prime or t(i)>p_i]cdot i^k)，其中(t(x))表示(x)的最小质因子，(p_i)表示第(i)个质数。
直观理解就是埃拉托斯特尼筛法进行了(i)轮之后(1sim n)剩下的数的(k)次方之和。
显然当(p_i^2>n)时(g(n,i)=g(n,i-1))，下面只考虑(p_i^2leqslant n)。
我们考虑(g(n,i))如何从(g(*,i-1))转移过来，也就是考虑第(i)轮筛掉了那些数。
那么可以得到式子：
[g(n,i)=g(n,i-1)-p_i^kcdot left(g(lfloorfrac{n}{p_i}
floor,i-1)-s_{i-1}
ight)
]
其中(s_i=sum_{j=1}^{i}p_i^k)。
解释一下这个式子，我们考虑筛掉了哪些数，然后减掉就好了，筛掉的那些数最小的质因子显然是(p_i)，而正好(g(lfloorfrac{n}{p_i} floor,i-1))代表最小质因子(geqslant p_i)的数以及质数，我们把质数减掉，拼起来就是(1sim n)。
那么写在一起就是：
[g(n,i)=
egin{cases}
g(n,i-1)&p_i^2>n\
g(n,i-1)-p_i^kcdot left(g(lfloordfrac{n}{p_i}
floor,i-1)-s_{i-1}
ight)&p_i^2leqslant n
end{cases}
]
根据常识可以知道第一维只有(2sqrt n)个取值，离散化一下然后递推就好了。
注意到这里我们顺便求出了对于每个(x),(sum_{i=1}^{lfloor n/x floor}[iin prime]cdot i^k)的值。
part 2
我们想算出(f)的前缀和，同样的我们设(h(n,i))如下：
[h(n,i)=sum_{j=1}^{n}[t(j)geqslant p_i]cdot f(j)
]
那么(h(n,1)+f(1))就是答案，注意我们(h)不会把(1)算进去。
首先还是很显然的，当(p_i>n)时，(h(n,i)=0)，同样下面只考虑(p_ileqslant n)的情况。
我们可以通过(g)很简单的算出质数的贡献，即：
[sum_{j=1}^{n}[jin prime]cdot f(j)-sum_{j=1}^{i-1}f(p_j)
]
这些都可以预处理。
然后我们考虑暴力地枚举最小的质因子是多少以及用了多少个，即：
[sum_{j=i}^{p_j^2leqslant n}sum_{k=1}^{p_j^{k+1}leqslant n}left(f(p_j^k)cdot h(lfloorfrac{n}{p_j^k}
floor,j+1)+f(p_j^{k+1})
ight)
]
解释下后面的(f(p_j^{k+1}))，这是因为(h)没有把(1)算进去，所以会漏算质数的次幂，直接加上就行了。
综合起来就是：
[h(n,i)=egin{cases}
0&p_i>n\
sum_{j=1}^{n}[jin prime]cdot f(j)-sum_{j=1}^{i-1}f(p_j)+sum_{j=i}^{p_j^2leqslant n}sum_{k=1}^{p_j^{k+1}leqslant n}left(f(p_j^k)cdot h(lfloorfrac{n}{p_j^k}
floor,j+1)+f(p_j^{k+1})
ight)&p_ileqslant n
end{cases}
]
写成这样求和符号好像有点丑...不管它
直接暴力递归算就好了，递推好像还慢些。
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('
');}
#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 [iin pri] i，其中第二位滚动掉了
h[j]=del(h[j],del(h[get_id(w[j]/pri[i])],i-1)); //h记录的是sum [iin 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;
}
这和上题也是一样的，注意别爆( ext{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('
');}
#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;
}