【学习笔记】min25筛

话说我们现在要求一个函数\(f\)的前缀和。即求\(F(n)=\sum_{i=1}^nf(i)\)。

min25筛这个算法的主要思想是把\(1...n\)这些数按质数和合数分类，然后分别考虑质数和合数的贡献。

STEP1 质数贡献

我们尝试先解决一个小问题：求\(G(m)=\sum_{i=1}^m[i\in prime]f(i)\)，即\(m\)以下的质数的\(f\)值之和。其中\(m\in\{\lfloor\frac{n}{1}\rfloor,\lfloor\frac{n}{2}\rfloor,...,\lfloor\frac{n}{n}\rfloor\}\)，共有\(O(\sqrt{n})\)种取值。

这玩意很难直接求，于是考虑DP。设\(g(j,m)=\sum_{i=1}^m[i\in prime \; or \; minp(i)>p[j]]f'(i)\)，其中\(minp(i)\)表示\(i\)的最小质因子。这个式子的含义是“\(i\)为质数或\(i\)的最小质因子\(>p[j]\)的函数值\(f'(i)\)”之和。其中\(f'(i)\)并不一定是我们要求的\(f(i)\)，它可以是我们构造出来的另一个函数，但它必须满足以下条件：

\(f'\)在质数处的取值与\(f\)相同。即\(\forall p\in prime\)有\(f'(p)=f(p)\)。
\(f'\)是完全积性函数。
\(f'\)可以快速求前缀和。

显然DP的边界是\(g(0,m)=\sum_{i=2}^mf'(i)\)。

考虑转移，求\(g(j,m)\)。当\(p[j]^2>m\)时，\(m\)以内不可能有任何最小质因子是\(p[j]\)的合数，即\(p[j]\)无法筛掉任何数，因此\(g(j,m)=g(j-1,m)\quad (p[j]^2>m)\)。

否则，我们只需要在\(g(j-1,m)\)里“筛掉”\(\leq m\)的最小质因子为\(p[j]\)的合数即可。而它们都可以表示为“\(p[j]\times\)一个\(\lfloor\frac{m}{p[j]}\rfloor\)以内的最小质因子大于等于\(p[j]\)的数”。于是有：

\[g(j,m)=g(j-1,m)-f'(p[j])\times \big(g(j-1,\lfloor\frac{m}{p[j]}\rfloor)-\sum_{i=1}^{j-1}f'(p[i])\big)\qquad (p[j]^2\leq m) \]

这时\(f'\)作为一个完全积性函数的好处就体现出来了，它可以直接乘以后面的一坨东西而不用担心是否互质。请读者认真理解这个式子的含义。

令\(|P|\)表示\(\leq m\)的质数个数，则\(g(|P|,m)\)就是我们前面要求的\(G(m)\)的值了。显然，我们可以把DP的第一维滚掉。

在具体实现时，我们先预处理出\(\sqrt{n}\)以内的质数（线性筛），以及所有\(m\)的值（用数论分块）。暴力的做法我们可以用STL-map来记录所有\(m\)的值的编号，但这里可以更巧妙一些。我们开两个数组id1和id2，大小都是\(\sqrt{n}\)，巧妙地利用“对所有\(>\sqrt{n}\)的数\(x\)，\(n/x\)一定\(<\sqrt{n}\)”，就可以记下所有\(m\)的编号了。

int n,val[N*2],id1[N],id2[N];
//主函数（预处理）
int sqrt_n=sqrt(n),tot=0;
for(int i=1,j;i<=n;i=j+1) {
    j=n/(n/i);
    int w=n/i;val[++tot]=w;
    if(w<=sqrt_n) id1[w]=tot;
    else id2[n/w]=tot;
}
//查询某个m的编号
inline int get_id(int m) {
    if(m<=sqrt_n) return id1[m];
    else return id2[n/m];
}

这一部分的时间复杂度被证明是\(O(\frac{n^{\frac{3}{4}}}{\log n})\)的。

STEP2 贡献加和

现在我们用刚刚搞出来的\(G(m)=\sum_{i=1}^m[i\in prime]f(i)\)，求出所有\(f\)的和。

设\(S(n,j)=\sum_{i=1}^n[minp(i)\geq p[j]]f(i)\)。则显然最终答案是\(f(1)+S(n,1)\)。

正如一开始所说的，我们分别考虑质数和合数对\(S(n,j)\)的贡献。显然质数的贡献是\(G(n)-\sum_{i=1}^{j-1}f(p[i])\)。对于合数的贡献，我们可以枚举这个合数的“最小质因子及其次数”，然后进行递推。容易列出递推式：

\[S(n,j)=G(n)-\sum_{i=1}^{j-1}f(p[i])+\sum_{i=j}^{p[i]^2\leq n}\sum_{e=1}^{p[i]^{e+1}\leq n}\Big(f(p[i]^e)S(\lfloor\frac{n}{p[i]^e}\rfloor,i+1)+f(p[i]^{e+1})\Big) \]

式子里的二重循环就是在枚举合数的“最小质因子及其次数”，注意合数可以有多个不同的质因子（继续递归）或只有一个质因子的若干次方，这些都在式子里得以体现，请读者仔细理解。

在求\(S\)时我们直接递归，不用记忆化。这一部分的时间复杂度，被证明也是\(O(\frac{n^{\frac{3}{4}}}{\log n})\)。反正你只要知道它能跑很大的数据（\(10^{10}\)的数据大约1秒搞定）就可以了。

一些栗子

wqy大神的良心模板题

预处理\(f_1(i)=i^2\)和\(f_2(i)=i\)在质数处的前缀和，在需要用到时减一下即可。

参考代码：

//P5325
#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define mk make_pair
#define pii pair<int,int>
#define fst first
#define scd second
using namespace std;
inline ll read() {
	ll f=1,x=0;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
inline void write(ll x) {
	if(!x)putchar('0');if(x<0)x=-x,putchar('-');
	static int sta[20];register int tot=0;
	while(x)sta[tot++]=x%10,x/=10;
	while(tot)putchar(sta[--tot]+48);
}
const int MAXN=1e6+5;
const int MOD=1e9+7,INV6=166666668;
ll n,sqrtn,prm[MAXN],sump[MAXN],sum2p[MAXN],cnt,g1[MAXN],g2[MAXN];
bool v[MAXN];
void sieve(int n) {
	cnt=0;
	for(int i=2;i<=n;++i) {
		if(!v[i]) {
			prm[++cnt]=i;
			sump[cnt]=(sump[cnt-1]+i)%MOD;
			sum2p[cnt]=(sum2p[cnt-1]+(ll)i*i%MOD)%MOD;
		}
		for(int j=1;j<=cnt && prm[j]*i<=n;++j) {
			v[prm[j]*i]=1;
			if(i%prm[j]==0) break;
		}
	}
}
ll val[MAXN];
int id1[MAXN],id2[MAXN];
inline int get_id(ll x){return ((x<=sqrtn)?id1[x]:id2[n/x]);}
ll getS(ll x,ll y) {
	if(x<prm[y] || x<=1) return 0;
	int k=get_id(x);
	ll res=((g2[k]-g1[k]+MOD)%MOD-(sum2p[y-1]-sump[y-1]+MOD)%MOD+MOD)%MOD;
	for(int i=y;i<=cnt && prm[i]*prm[i]<=x;++i) {
		ll t1=prm[i],t2=prm[i]*prm[i];
		for(int j=1;t2<=x;++j,t1=t2,t2*=prm[i]) {
			ll tt1=t1%MOD,tt2=t2%MOD;
			res=(res+getS(x/t1,i+1)*tt1%MOD*(tt1-1)%MOD+tt2*(tt2-1)%MOD)%MOD;
		}
	}
	return res%MOD;
}
int main() {
	ios::sync_with_stdio(0);/*syn加速*/
	n=read();
	sieve(sqrtn=sqrt(n));
	int tot=0;/*O(sqrt(n))级别*/
	for(ll i=1,j;i<=n;i=j+1) {
		j=n/(n/i);
		ll w=n/i;
		val[++tot]=w;
		if(w<=sqrtn) id1[w]=tot;
		else id2[n/w]=tot;
		
		w%=MOD;
		g1[tot]=w*(w+1)/2%MOD;
		g2[tot]=w*(w+1)%MOD*(2LL*w+1)%MOD*INV6%MOD;
		
		g1[tot]=(g1[tot]-1+MOD)%MOD;
		g2[tot]=(g2[tot]-1+MOD)%MOD;
	}
	for(int j=1;j<=cnt;++j) {
		for(int i=1;i<=tot && prm[j]*prm[j]<=val[i];++i) {
			int k=get_id(val[i]/prm[j]);
			g1[i]=(g1[i]-prm[j]*(g1[k]-sump[j-1]+MOD)%MOD+MOD)%MOD;
			g2[i]=(g2[i]-prm[j]*prm[j]%MOD*(g2[k]-sum2p[j-1]+MOD)%MOD+MOD)%MOD;
		}
	}
	//for(int i=1;i<=tot;++i)cout<<val[i]<<" "<<g1[i]<<" "<<g2[i]<<endl;
	write((getS(n,1)+1)%MOD);puts("");
	return 0;
}

意外搞定了杜教筛的模板

先看\(\sum_{i=1}^n\varphi(i)\)。首先\(\varphi(p)=p-1\quad (p\in prime)\)。于是我们预处理\(f_1(i)=i\)和\(f_2(i)=1\)的前缀和，按与上一题相同的方法相减即可求出\(G(m)\)。

另外，\(\varphi(p^c)=p^{c-1}(p-1)\)（自己yy可得），这样在求\(S\)时，质数的整次幂的函数值也很好求。

然后是\(\mu(i)\)，显然\(\mu(p)=-1\quad (p\in prime)\)。于是我们只要对前面已经求好的\(f_2(i)=1\)求个相反数即可。

另外，\(\mu(p^k)=0\quad (k>1)\)，因此求\(S\)时不需要枚举最小质因子的次数。

参考代码：

//P4213 min_25 
#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define mk make_pair
#define pii pair<int,int>
#define fst first
#define scd second
using namespace std;
/*  -----  by:duyi  -----  */
const int N=5e4;
const int MAXN=N*2+5;
int n,sn,cnt,tot,p[MAXN],sum[MAXN],id1[MAXN],id2[MAXN],val[MAXN];
ll g1[MAXN],g2[MAXN];
bool v[MAXN];
void sieve() {
	v[1]=1;
	for(int i=2;i<=N;++i) {
		if(!v[i]) p[++cnt]=i;
		for(int j=1;j<=cnt && (ll)i*p[j]<=N;++j) {
			v[i*p[j]]=1;
			if(i%p[j]==0) break;
		}
	}
	for(int i=1;i<=cnt;++i) sum[i]=sum[i-1]+p[i];
}
inline int get_id(int x) {
	if(x<=sn) return id1[x];
	else return id2[n/x];
}
ll S_phi(int x,int y) {
	if(x<=1 || p[y]>x) return 0;
	ll res=g1[get_id(x)]-g2[get_id(x)]-(sum[y-1]-(y-1));
	for(int i=y;i<=cnt && (ll)p[i]*p[i]<=x;++i) {
		ll pre=1,cur=p[i];
		for(int j=1;cur*p[i]<=x;++j) {
			res+=pre*(p[i]-1LL)*S_phi(x/cur,i+1)+cur*(p[i]-1LL);
			pre=cur;cur=cur*p[i];
		}
	}
	return res;
}
ll S_mu(int x,int y) {
	if(x<=1 || p[y]>x) return 0;
	ll res=-g2[get_id(x)]+y-1;
	for(int i=y;i<=cnt && (ll)p[i]*p[i]<=x;++i) {
		res+=(-S_mu(x/p[i],i+1));
	}
	return res;
}
void solve() {
	cin>>n;
	sn=sqrt(n);tot=0;
	for(int i=1,j;i<=n;i=j+1) {
		j=n/(n/i);int w=n/i;
		val[++tot]=w;
		if(w<=sn) id1[w]=tot;
		else id2[n/w]=tot;
		g1[tot]=(ll)w*(w+1LL)/2LL-1LL;
		g2[tot]=w-1;
	}
	for(int i=1;i<=cnt;++i) {
		for(int j=1;j<=tot && (ll)p[i]*p[i]<=val[j];++j) {
			int t=get_id(val[j]/p[i]);
			g1[j]-=(ll)p[i]*(g1[t]-sum[i-1]);
			g2[j]-=(g2[t]-(i-1));
		}
	}
	cout<<S_phi(n,1)+1LL<<" "<<S_mu(n,1)+1LL<<endl;
}
int main() {
	ios::sync_with_stdio(0);/*syn加速*/
	sieve();
	int T;cin>>T;while(T--)solve();
	return 0;
}

总结一下

用min_25做题主要是要想清楚两个问题：

该函数在质数处的取值怎么求和？（找到合适的\(f'\)）
质数的次幂处的取值\(f(p^c)\)能否快速求？

本文写得比较粗浅，欢迎大家补充。如有纰漏欢迎指正。

posted @ 2020-08-06 21:19 duyiblue 阅读(1815) 评论(3) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

望舒草

就让失散的誓言飞舞吧，随西风飘荡

【学习笔记】min25筛

STEP1 质数贡献

STEP2 贡献加和

一些栗子

总结一下

公告