min_25 筛略解
min_25 筛因其发明者为 min_25 而得名。该算法脱胎自埃氏筛。
对于一类积性函数 \(f\)，它可以在低于线性的时间复杂度下计算：
-
\[\sum_{i=1}^n[i\text{ is prime}]f(i) \]
-
\[\sum_{i=1}^nf(i) \]
- （即 \(f\) 的前缀和）
并得到上式在 \(\left\{\lfloor \dfrac{n}{i}\rfloor\mid i\in \mathbf{N},i\in [1,n]\right\}\) 这些地方的值。
记号
- \(P\) 为质数集合；
- 下文默认 \(p\in P\)；
- \(p_i\) 为第 \(i\) 个质数；
- 特别地，\(p_0=1\)；
- \(\operatorname{min}(x)\) 为 \(x\) 的最小质因子；
- \(f(x)\) 为希望求出前缀和的积性函数；
- 通常要求 \(f(p)\) 为关于 \(p\) 的低次多项式，且能较快求出 \(f(p^q)\)；
- \(f'(x)\) 为 \(f'(x)=f(x)(x\in P)\) （即与 \(f\) 在质数处值相等）的一个完全积性函数；
- 通常将 \(f(p)\) 拆成多个单项式，分别作为 \(f'\)，这样 \(f'\) 可以快速求前缀和。
筛 f 在质数处的前缀和
记
先考虑我们是怎样用最原始的方法（埃氏筛）筛质数的：从小到大枚举每个质数 \(p_i\)，将 \(p_i\) 的倍数都标记为合数。第 \(i\) 轮新筛掉的数满足 \(x\neq p_i,\min(x)=p_i\)。
故记 \(g(n,i)\) 为第 \(i\) 轮筛完之后剩下的数的 \(f\) 的前缀和，即：
设 \(k\) 满足 \(p_k^2\leq n< p_{k+1}^2\)，显然 \(k\) 轮筛完后只剩质数，即 \(g(n)=g(n,k)\)。
我们如果可以通过 \(g(n,i-1)\) 求出 \(g(n,i)\)，就能通过 \(g(n,0)\) 求出 \(g(n,k)\)（即 \(g(n)\)）。
- 注意：\(g(n,i)\) 一律不包括 \(f'(1)\)。
考虑一轮中新筛去的数的贡献：新筛掉的数都满足 \(x\not=p_i\ \text{and}\ \min(x)=p_i\)，即它们有质因数 \(p_i\) ，且去掉 \(p_i\) 后其质因数仍然不小于 \(p_i\)（否则它本身就是 \(p_i\) 了，不应该筛掉）。
\(\lfloor \dfrac{n}{i}\rfloor\) 只有 \(O(\sqrt n)\) 个点值，故只存储这些地方的 \(g\)，并预处理其 \(g(i,0)\)。
\(\sum_{i=1}^{?}f'(p_i)\) 要预处理到 \(k\)（其实就是提前线筛前 \(\sqrt n\) 个数）。
求 \(f\) 的前缀和
记 \(S(n)=\sum_{i=1}^nf(i)\)。与 \(g\) 类似，定义
（注意定义略有差别，\(S(n,i)\) 并不保证包括所有质数）
则答案为 \(S(n)=S(n,1)+1\)。（\(1\) 仍然被 \(S\) 排除在外）
根据质数、合数分类计算。对于合数 \(m\)，枚举 \(m\) 的最小质因数 \(j\) 和它的次数 \(c\)，则
直接按式子递归计算，复杂度为 \(O(n^{1-\epsilon})\)。若采用和洲阁筛相同的的方式（不会），复杂度为 \(O(\dfrac{n^{0.75}}{\log n})\)。
例题
LOJ #6053 简单的函数
求 \(f(p^k)=p\oplus k\) 的前缀和。
\(f(p)=p-1(p\neq 2)\)。故第一部分筛 \(f'_1(p)=1,f'_2(p)=p\) 在质数处的前缀和。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
ll getint(){
ll ans=0,f=1;
char c=getchar();
while(c<'0'||c>'9'){
if(c=='-')f=-1;
c=getchar();
}
while(c>='0'&&c<='9'){
ans=ans*10+c-'0';
c=getchar();
}
return ans*f;
}
const int N=1e6+10,mod=1e9+7;
ll n,sq;
ll pri[N],spri[N],cnt=0;
bool boo[N];
ll w[N],g[N],h[N],m=0;
ll id1[N],id2[N];
void init(int n){
for(int i=2;i<=n;i++){
if(!boo[i]){
pri[++cnt]=i;
spri[cnt]=(spri[cnt-1]+pri[cnt])%mod;
}
for(int j=1;j<=cnt&&i*pri[j]<=n;j++){
boo[i*pri[j]]=1;
if(i%pri[j]==0)break;
}
}
}
ll S(ll x,ll y){
if(x<=1||pri[y]>x)return 0;
int k=(x<=sq?id1[x]:id2[n/x]);
ll ans=(g[k]-spri[y-1]-h[k]+y-1)%mod;
if(y==1)ans+=2;
for(int i=y;i<=cnt&&pri[i]*pri[i]<=x;i++){
ll t=pri[i];
for(int j=1;t*pri[i]<=x;j++,t*=pri[i]){
ans=(ans+(S(x/t,i+1)*(pri[i]^j)%mod)+(pri[i]^(j+1)))%mod;
}
}
return ans;
}
int main(){
n=getint(),sq=sqrt(n);
init(sq);
for(ll l=1,r=0;l<=n;l=r+1){
r=n/(n/l);
w[++m]=n/l;
h[m]=(w[m]-1)%mod; //i-1
g[m]=(w[m]+1ll)%mod*(w[m]%mod)%mod;
if(g[m]&1)g[m]+=mod;g[m]>>=1;g[m]--;//i*(i-1)/2-1
if(w[m]<=sq)id1[w[m]]=m;
else id2[r]=m;
}
for(int j=1;j<=cnt;j++){
for(int i=1;i<=m&&pri[j]*1ll*pri[j]<=w[i];i++){
ll t=w[i]/pri[j],k=(t<=sq?id1[t]:id2[n/t]);
g[i]=(g[i]-pri[j]*1ll*(g[k]-spri[j-1]))%mod;
h[i]=(h[i]-(h[k]-j+1))%mod;
}
}
ll ans=S(n,1)+1;
cout<<(ans%mod+mod)%mod;
return 0;
}
51Nod 1847 奇怪的数学题
设 \(f(n)\) 为 \(i\) 的第二大因数（\(f(n)=\dfrac{n}{\min(n)}\)），求 \(\sum_{i=1}^n \sum_{j=1}^n f(\gcd(i,j))^k\)。
简单莫反可得：
整除分块显然，\(\varphi\) 的前缀和可以杜教筛求出，问题在于 \(f(i)^k\) 的前缀和。
注意到 \(f(i)^k\) 和 \(i^k\) 相差的是 \(\min(i)^k\)，而 min_25 筛便将 \(\min(i)\) 不同的数分开来。观察 min_25 筛第一部分的式子（节选）：
红色部分便是 \(\min(j)=p_i\) 的合数 \(j\) 的 \(f(j)^k\)。
质数要单独拿出来统计。
注意到模数十分阴间，预处理自然数幂和时最好用斯特林数+下降幂。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int M=7e4+10,L=2e6+10,K=57;
const ll N=1e9+10;
#define uint unsigned int
uint n,k,sqr,lim;
uint pri[L+10],sp1[M],sp2[M],boo[L+10],ppow[M],cnt=0;
uint id1[M],id2[M],g1[M],g2[M],g[M],w[M];
uint fs[M];
uint phi[L+10],phi2[L+10];
int _(ll x){ return x<sqr?id1[x]:id2[n/x]; }
uint qpow(uint x,uint y){
uint ans=1;
while(y){
if(y&1)ans=ans*x;
x=x*x;
y>>=1;
}
return ans;
}
uint s2[K][K];
void init(int k){
s2[0][0]=1;
for(int i=1;i<=k;i++)
for(int j=1;j<=i;j++)
s2[i][j]=s2[i-1][j]*j+s2[i-1][j-1];
}
uint s(uint x){
// sum_{i=0}^{x-1}
uint ans=0;
for(int i=0;i<=k;i++){
uint p=1;
for(int j=0;j<=i;j++){
if((x-j)%(i+1)==0)p*=(x-j)/(i+1);
else p*=x-j;
}
ans+=p*s2[k][i];
}
return ans;
}
uint sphi(uint x){
if(x<=lim)return phi[x];
if(phi2[n/x])return phi2[n/x];
uint ans=x*1ll*(x+1)/2;
for(int l=2,r;l<=x;l=r+1){
r=x/(x/l);
ans-=(r-l+1)*sphi(x/l);
}
// cerr<<"! "<<x<<" "<<ans<<endl;
return phi2[n/x]=ans;
}
int main(){
cin>>n>>k;
init(k);
sqr=sqrt(n+1);
for(int i=2;i<=sqr;i++){
if(!boo[i]){
pri[++cnt]=i;
ppow[cnt]=qpow(i,k);
sp1[cnt]=(sp1[cnt-1]+ppow[cnt]);
sp2[cnt]=(sp2[cnt-1]+1);
phi[i]=i-1;
}
for(int j=1;j<=cnt&&i*pri[j]<=sqr;j++){
boo[i*pri[j]]=1;
if(i%pri[j]==0)break;
}
}
int cnt_=cnt;
phi[1]=1;
lim=pow(n+1,0.68);
cnt=0;
for(int i=2;i<=lim;i++){
if(!boo[i]){
pri[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt&&i*pri[j]<=lim;j++){
boo[i*pri[j]]=1;
if(i%pri[j])phi[i*pri[j]]=phi[i]*(pri[j]-1);
else{
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
}
}
for(int i=1;i<=lim;i++)phi[i]+=phi[i-1];
int m=1;
for(ll l=1,r;l<=n;l=r+1,m++){
r=n/(n/l);
uint t=n/l;
w[m]=t;
g1[m]=s(t+1)-1;
g2[m]=t-1;
if(t<sqr)id1[t]=m;
else id2[n/t]=m;
}
--m;
for(int i=1;i<=cnt_;i++){
int k=1;
for(int j=1;j<=m&&w[j]>=pri[i]*1ll*pri[i];j++){
ll r=_(w[j]/pri[i]);
fs[j]+=(g1[r]-sp1[i-1]);
g1[j]-=ppow[i]*(g1[r]-sp1[i-1]);
g2[j]-=(g2[r]-sp2[i-1]);
}
}
for(int i=1;i<=m;i++)fs[i]+=g2[i];
uint ans=0;
for(int i=1;i<=m;i++)
ans+=(fs[i]-fs[i+1])*(sphi(n/w[i])*2-1);
cout<<ans;
}
