Min_25筛

听说这个东西能给予人力量

那就来学一学吧

功能就是筛一个积性函数\(f(i)\)的前缀和

Min_25筛好像是最近才流行起来的筛法,复杂度是非常神奇的\(O(\frac{n^{\frac{3}{4}}}{logn})\)

和杜教筛一样,使用这个筛法的也有一定要求,就是\(f(p^c)\)需要在\(O(1)\)求出

来看看这个非常力量的筛法

我们要求的东西是

\[\sum_{i=1}^nf(i) \]

我们先定义一个二元函数\(g(n,j)\)

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

\(P\)表示质数集,\(minp(i)\)表示\(i\)的质因子,\(P_j\)表示第\(j\)个质数

看起来可能没什么含义,但是我们可以理解为正在进行一个埃筛,\(P_j\)这个质数已经筛完了,这个时候所有被判定为质数的位置和所有没被筛到的位置的\(f\)函数的和,就表示\(g(n,j)\)

考虑如何推\(g(n,j)\)

如果\(P_j^2>n\),那么\(P_j\)这个数非常可怜,一个数都没有被它筛到,因为\(P_j\)筛到的数最小也得是\(P_j^2\),所以就有

\[g(n,j)=g(n,j-1) \]

如果\(P_j^2<=n\),那么我们就需要考虑\(P_j\)\(g(n,j-1)\)筛到了什么,之后减掉这些贡献

筛掉的都是那些\(minp(i)=P_j\)的数吧,根据直觉

\[g(n,j)=g(n,j-1)-f(P_j)(g(\left \lfloor \frac{n}{P_j} \right \rfloor,j-1)-\sum_{i=1}^{j-1}f(P_i)) \]

来解释一下为什么吧

首先那个\(\left \lfloor \frac{n}{P_j} \right \rfloor\)能保证这个范围里的数乘上\(P_j\)都不大于\(n\),之后这个范围里的数都\(minp\)小于等于\(P_{j-1}\)的都被干掉了,所以这个范围内剩下的数的\(minp\)都大于等于\(P_j\),之后乘上\(P_j\),最小质因子自然就是\(P_j\)

但是我们这样把质数也给算了上去,所以并不是很可行,所以后面还得把\(1\)\(j-1\)的质数加回来

同时这里有一个问题,我们将\(\left \lfloor \frac{n}{P_j} \right \rfloor\)范围的数还原的时候直接乘上\(f(P_j)\)并不科学,因为这里显然会有不互质的情况出现,所以这里是将\(f\)当做完全积性函数来做的

综上

\[g(n,j)=\begin{cases} g(n,j-1)&P_j^2\gt n\\ g(n,j-1)-f(P_j)[g(\left \lfloor \frac{n}{P_j} \right \rfloor,j-1)-\sum_{i=1}^{j-1}f(P_i)]&P_j^2\le n\end{cases} \]

求出这个\(g\)什么用都没有,我们还需要再设

\[S(n,j)=\sum_{i=1}^n[minp(i)>=P_j]f(i) \]

这里就不在把\(f(i)\)当成完全积性了

这个时候我们的答案自然就是\(S(n,1)+1\)

我们把\(S(n,j)\)分成质数贡献的一部分和合数贡献的一部分

质数哪一部分的贡献自然就是

\[g(n,|P|)-\sum_{i=1}^{j-1}f(P_j) \]

后面那个是减掉小于\(j\)的质数的贡献

而合数的部分我们通过枚举最小质因子和最小质因子出现的次数

\[\sum_{k=j}^{P_k^2<=n}\sum_{e=1}^{P_k^e<=n}f(P_k^e)(S(\left \lfloor \frac{n}{P_k^e} \right \rfloor,k+1)+[e>1]) \]

这里跟上面类似,由于后面是\(j+1\)所以会导致\(S(\left \lfloor \frac{n}{P_k^e} \right \rfloor,k+1)\)内算贡献的数都最小质因子都大于等于\(P_{k+1}\),直接乘上\(P_k^e\)不可能产生不互质的情况

由于\(1\)的贡献是最后算的,\(S(\left \lfloor \frac{n}{P_k^e} \right \rfloor,k+1)\)并没有算到\(1\)也就是还原回来的\(P_k^e\)所以在最后需要加\(1\)

特殊的是\(e=1\)本身就是质数,之前算过贡献这里就不需要加\(1\)

代码实现这个筛法就非常复杂了

来一道简单的例题

loj#6235. 区间素数个数

显然这个东西不是积性函数,但是考虑一下Min_25筛的第一步,就是只算质数部分的贡献

我们这里把\(f(i)=[i>1]\)就好了

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define maxn 1000005
#define re register
#define LL long long
const LL mod=1000000007;
LL n,w[maxn],Sqr,tot;
LL p[maxn],pre[maxn],id1[maxn],id2[maxn],h[maxn],g[maxn],m;
int f[maxn];
int main()
{
	scanf("%lld",&n);Sqr=std::sqrt(n);
	f[1]=1,pre[1]=1;
	for(re int i=2;i<=Sqr;i++) {
		if(!f[i]) p[++tot]=i,pre[tot]=(pre[tot-1]+(i^1))%mod;
		for(re int j=1;j<=tot&&p[j]*i<=Sqr;j++) {
			f[p[j]*i]=1;if(i%p[j]==0) break;
		}
	}//预处理根号范围内的质数 
	for(re LL l=1,r;l<=n;l=r+1) {
		r=n/(n/l);w[++m]=n/l;
		if(w[m]<=Sqr) id1[w[m]]=m;
			else id2[n/w[m]]=m;
		g[m]=w[m]-1;//这里的g[m]=g(w[m],0),所以初值是除去1以外的数的个数 
	}
	for(re int j=1;j<=tot;j++)
		for(re int i=1;i<=m&&p[j]*p[j]<=w[i];i++) {
			int k=(w[i]/p[j]<=Sqr)?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])];
			//求w[i]/p[j]在哪一个块里 
			g[i]=g[i]-g[k]+j-1;//这里实际上就是从g(i,j-1)推向g(i,j) 
		}
	printf("%lld\n",g[1]);
	return 0;
}

loj#6053. 简单的函数

首先注意到除去唯一的偶质数\(2\),所有的\(p\)都满足\(p⊕1=p-1\)

\(2⊕1=3\)

先设\(f(p)=p\),我们再求一个区间内素数个数求筛的时候减掉这个区间素数个数

如果有\(2\)这个质数记得把答案加\(2\)

之后也没什么了

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define maxn 300005
#define re register
#define LL long long
const LL mod=1000000007;
const LL inv2=500000004;
LL n,w[maxn],Sqr,tot;
LL p[maxn],pre[maxn],id1[maxn],id2[maxn],h[maxn],g[maxn],m;
int f[maxn];
inline LL S(LL x,int y) {
	if(x<=1||p[y]>x) return 0;
	int k=(x<=Sqr)?id1[x]:id2[n/x];
	LL ans=(h[k]-g[k]-pre[y-1]+y-1+((y==1)?2ll:0))%mod;ans=(ans+mod)%mod;
	for(re int k=y;k<=tot&&p[k]*p[k]<=x;k++) {
		LL p1=p[k];
		for(re int e=1;p1<=x;e++,p1*=p[k]) 
			ans+=(1ll*(S(x/p1,k+1)+((e>1)?1ll:0))*(p[k]^e)%mod),ans%=mod;
	}
	return ans;
}
int main()
{
	scanf("%lld",&n);Sqr=std::sqrt(n)+1;
	f[1]=1;
	for(re int i=2;i<=Sqr;i++) {
		if(!f[i]) p[++tot]=i,pre[tot]=(pre[tot-1]+i)%mod;
		for(re int j=1;j<=tot&&p[j]*i<=Sqr;j++) {
			f[p[j]*i]=1;if(i%p[j]==0) break;
		}
	}
	for(re LL l=1,r;l<=n;l=r+1) {
		r=n/(n/l);w[++m]=n/l;
		if(w[m]<=Sqr) id1[w[m]]=m;
			else id2[n/w[m]]=m;
		g[m]=(w[m]-1)%mod;
		h[m]=((w[m]+2ll)%mod)*((w[m]-1ll)%mod)%mod;h[m]=(h[m]*inv2)%mod;
	}
	for(re int j=1;j<=tot;j++)
		for(re int i=1;i<=m&&p[j]*p[j]<=w[i];i++) {
			int k=(w[i]/p[j]<=Sqr)?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])]; 
			g[i]=g[i]-g[k]+j-1ll; g[i]=(g[i]+mod)%mod;
			h[i]=(h[i]-(h[k]-pre[j-1]+mod)%mod*p[j]%mod+mod)%mod;
		}
	printf("%lld\n",1+S(n,1));
	return 0;
}
posted @ 2019-02-13 21:56  asuldb  阅读(178)  评论(0编辑  收藏  举报