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 在质数处的前缀和

记

\[g(n)=\sum_{i=1}^n[i\in P]f(i)=\sum_{i=1}^n[i\in P]f'(i) \]

先考虑我们是怎样用最原始的方法（埃氏筛）筛质数的：从小到大枚举每个质数 \(p_i\)，将 \(p_i\) 的倍数都标记为合数。第 \(i\) 轮新筛掉的数满足 \(x\neq p_i,\min(x)=p_i\)。

故记 \(g(n,i)\) 为第 \(i\) 轮筛完之后剩下的数的 \(f\) 的前缀和，即：

\[g(n,i)=\sum_{j=1}^n[j\in P\ \text{or}\ \min(j)>p_i]f(i) \]

设 \(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\) 了，不应该筛掉）。

\[\begin{aligned} g(n,i)=&g(n,i-1)-\sum_{j=1}^n[j\not=p_i\ \text{and}\ \min(j)=p_i]f'(j) \quad &\text{去掉被筛的数}\\ =&g(n,i-1)-f'(p_i)\sum_{j=2}^{\lfloor {n\over p_i}\rfloor}[\min(j)>p_{i-1}]f'(j) &\text{把 $p_i$ 提出来} \\=&g(n,i-1)-f'(p_i) \left(\sum_{j=2}^{\lfloor {n\over p_i}\rfloor}[j\in P\ \text{or}\ \min(j)>p_{i-1}]f'(j)-\sum_{j=1}^{i-1}f'(p_j) \right) \\=&g(n,i-1)-f'(p_i) \left(g(\lfloor {n\over p_i}\rfloor,i-1)-\sum_{j=1}^{i-1}f'(p_j) \right) \end{aligned} \]

\(\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)=\sum\limits_{j=2}^n[\min(j)\geq p_i]f(i) \]

（注意定义略有差别，\(S(n,i)\) 并不保证包括所有质数）

则答案为 \(S(n)=S(n,1)+1\)。（\(1\) 仍然被 \(S\) 排除在外）

根据质数、合数分类计算。对于合数 \(m\)，枚举 \(m\) 的最小质因数 \(j\) 和它的次数 \(c\)，则

\[\begin{aligned} S(n,i)=&\sum\limits_{j=1}^n[\min(j)\geq p_i]f(i)\\ =&\color{blue}\sum\limits_{j=i}^{p_j^2\leq n} \sum\limits_{c=1}^{p_j^c\leq n} ([p_j^c\not\in P]f(p_j^c)+\sum\limits_{k=2}^{\lfloor {n\over p_j^c}\rfloor} [\min(k)\geq p_{i+1}]f(k\cdot p_j^c))\color{black}+\color{red}\sum\limits_{i=k}^{p_i\leq n}f(p_i) \\&\text{蓝色部分为合数，红色部分为质数} \\=&\sum\limits_{j=i}^{p_j^2\leq n} \sum\limits_{c=1}^{p_j^c\leq n} f(p_j^c)(\color{orange}[c>1]\color{black}+\color{purple}S({\lfloor {n\over p_j^c}\rfloor},i+1)\color{black})+\color{red}g(n)-g(p_{k-1}) \\=&\sum\limits_{j=i}^{p_j^2\leq n} \sum\limits_{c=1}^{p_j^{c+1}\leq n} (\color{orange}f(p_j^{c+1})\color{black}+\color{purple}f(p_j^c)S({\lfloor {n\over p_j^c}\rfloor},i+1)\color{black})+\color{red}g(n)-g(p_{k-1}) \\&\text{橙色部分容易理解，} \\&\text{紫色部分成立的原因在于 $p_j^c\leq n<p_j^{c+1}$ 时，} \\&\text{$\lfloor {n\over p_j^c}\rfloor<p_i<p_{i+1}$，故此时 $S({\lfloor {n\over p_j^c}\rfloor},i+1)=0，$} \\&\text{不必枚举 $p_j^c\leq n<p_j^{c+1}$ 的情况。} \end{aligned}\]

直接按式子递归计算，复杂度为 \(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\)。

简单莫反可得：

\[ans=\sum_{i=2}^{n}f(i)^k(-1+2\sum_{j=1}^{\lfloor \frac{n}{i}\rfloor}\varphi(i)) \]

整除分块显然，\(\varphi\) 的前缀和可以杜教筛求出，问题在于 \(f(i)^k\) 的前缀和。

注意到 \(f(i)^k\) 和 \(i^k\) 相差的是 \(\min(i)^k\)，而 min_25 筛便将 \(\min(i)\) 不同的数分开来。观察 min_25 筛第一部分的式子（节选）：

\[g(n,i-1)-f'(p_i) \color{red}\left(g(\lfloor {n\over p_i}\rfloor,i-1)-\sum_{j=1}^{i-1}f'(p_j) \right) \]

红色部分便是 \(\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;
}