min_25 筛略解

min_25 筛因其发明者为 min_25 而得名。该算法脱胎自埃氏筛。

对于一类积性函数 (f)，它可以在低于线性的时间复杂度下计算：

• [sum_{i=1}^n[i ext{ is prime}]f(i) ]

• [sum_{i=1}^nf(i) ]
  • （即 (f) 的前缀和）

并得到上式在 (left{lfloor dfrac{n}{i} floormid iin mathbf{N},iin [1,n] ight}) 这些地方的值。

记号

• (P) 为质数集合；
• 下文默认 (pin 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)(xin P)) （即与 (f) 在质数处值相等）的一个完全积性函数；
  • 通常将 (f(p)) 拆成多个单项式，分别作为 (f')，这样 (f') 可以快速求前缀和。

筛 f 在质数处的前缀和

记

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

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

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

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

设 (k) 满足 (p_k^2leq 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 ot=p_i ext{and} min(x)=p_i)，即它们有质因数 (p_i) ，且去掉 (p_i) 后其质因数仍然不小于 (p_i)（否则它本身就是 (p_i) 了，不应该筛掉）。

[egin{aligned} g(n,i)=&g(n,i-1)-sum_{j=1}^n[j ot=p_i ext{and} min(j)=p_i]f'(j) quad & ext{去掉被筛的数}\ =&g(n,i-1)-f'(p_i)sum_{j=2}^{lfloor {nover p_i} floor}[min(j)>p_{i-1}]f'(j) & ext{把 $p_i$ 提出来} \=&g(n,i-1)-f'(p_i) left(sum_{j=2}^{lfloor {nover p_i} floor}[jin P ext{or} min(j)>p_{i-1}]f'(j)-sum_{j=1}^{i-1}f'(p_j) ight) \=&g(n,i-1)-f'(p_i) left(g(lfloor {nover p_i} floor,i-1)-sum_{j=1}^{i-1}f'(p_j) ight) end{aligned} ]

(lfloor dfrac{n}{i} floor) 只有 (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)=sumlimits_{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)，则

[egin{aligned} S(n,i)=&sumlimits_{j=1}^n[min(j)geq p_i]f(i)\ =&color{blue}{sumlimits_{j=i}^{p_j^2leq n} sumlimits_{c=1}^{p_j^cleq n} ([p_j^c otin P]f(p_j^c)+sumlimits_{k=2}^{lfloor {nover p_j^c} floor} [min(k)geq p_{i+1}]f(kcdot p_j^c))}+color{red}{sumlimits_{i=k}^{p_ileq n}f(p_i)} \& ext{蓝色部分为合数，红色部分为质数} \=&sumlimits_{j=i}^{p_j^2leq n} sumlimits_{c=1}^{p_j^cleq n} f(p_j^c)(color{orange}{[c>1]}+color{purple}{S({lfloor {nover p_j^c} floor},i+1)}+color{red}{g(n)-g(p_{k-1})} \=&sumlimits_{j=i}^{p_j^2leq n} sumlimits_{c=1}^{p_j^{c+1}leq n} (color{orange}{f(p_j^{c+1})}+color{purple}{f(p_j^c)S({lfloor {nover p_j^c} floor},i+1)})+color{red}{g(n)-g(p_{k-1})} \& ext{橙色部分容易理解，} \& ext{紫色部分成立的原因在于 $p_j^cleq n<p_j^{c+1}$ 时，} \& ext{$lfloor {nover p_j^c} floor<p_i<p_{i+1}$，故此时 $S({lfloor {nover p_j^c} floor},i+1)=0，$} \& ext{不必枚举 $p_j^cleq n<p_j^{c+1}$ 的情况。} end{aligned}]

直接按式子递归计算，复杂度为 (O(n^{1-epsilon}))。若采用和洲阁筛相同的的方式（不会），复杂度为 (O(dfrac{n^{0.75}}{log n}))。

例题

LOJ #6053 简单的函数

题目链接

求 (f(p^k)=poplus k) 的前缀和。

(f(p)=p-1(p eq 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+2sum_{j=1}^{lfloor frac{n}{i} floor}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 {nover p_i} floor,i-1)-sum_{j=1}^{i-1}f'(p_j) ight)} ]

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