筛

前言

好像没什么要说的。

线性筛

用来线性时间内筛质数，也可以求积性函数的值。

每个数会被最小质因子筛掉。

void init(){
	for(int i=2;i<=n;++i){
		if(!notp[i]) pri[++prs]=i;
		for(int j=1;j<=prs && i*pri[j]<=n;++j){
			notp[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
}

一下几个筛都可以做到在亚线性时间内求积性函数前缀和。

杜教筛

设需要求 \(S(n)=\sum_{i=1}^n f(i)\)，找出一个积性函数 \(g\)，考虑 \(f*g\) 的前缀和：

\[\begin{aligned} \sum_{i=1}^n (f*g)(i)&=\sum_{i=1}^n\sum_{d|i}g(d)f(\frac{i}{d})\\ \sum_{i=1}^n (f*g)(i)&=\sum_{i=1}^ng(i)S(\lfloor\frac{n}{i}\rfloor)\\ S(n)&=\sum_{i=1}^n (f*g)(i)-\sum_{i=2}^ng(i)S(\lfloor\frac{n}{i}\rfloor) \end{aligned} \]

所以如果 \(g\) 和 \(f*g\) 的前缀和都好求，那么我们可以用整除分块处理这个问题。

这个做法的时间是 \(\mathcal{O}(n^{\frac{3}{4}})\)，如果可以把 \(m\le n^{\frac{2}{3}}\) 的 \(S(i)\) 都先用线性筛算出来，则这个做法的时间可以做到 \(\mathcal{O}(n^{\frac{2}{3}})\)。

构造 \(g\) 的话可以利用 \(\varphi * 1=id,\mu *1=\varepsilon\)，比如 \(f(x)=x\varphi(x)\) 可以构造 \(g(x)=x\)，则 \((f*g)(x)=x^2\)。

以下是求 \(\mu\) 和 \(\varphi\) 的代码。

int GetSum_phi(int lim){
	if(lim<N) return phi[lim];
	if(sphi[lim]) return sphi[lim];
	int res=(1+lim)*lim/2;
	for(int l=2,r;l<=lim;l=r+1){
		r=lim/(lim/l);
		res-=(r-l+1)*GetSum_phi(lim/l);
	}
	return sphi[lim]=res;
}
int GetSum_mu(int lim){
	if(lim<N) return mu[lim];
	if(smu[lim]) return smu[lim];
	int res=1;
	for(int l=2,r;l<=lim;l=r+1){
		r=lim/(lim/l);
		res-=(r-l+1)*GetSum_mu(lim/l);
	}
	return smu[lim]=res;
}

Min_25 筛

我们以 \(f(p^e)=(p^e)^k-1\) 为例。

Min_25 筛的核心思想就是所有合数的最小质因子小于等于根号。

我们小于等于 \(n\) 的数分成质数和合数，对于合数，我们枚举它的最小质因子以及幂指数：

\[\begin{aligned} \sum_{i=1}^nf(i)=\sum_{p_i\le n}f(p_i)+\sum_{p_i\le \sqrt n,p_i^e\le n} f(p_i^e)\sum_{j\le n,minp(j)>p_i}f(j) \end{aligned} \]

Step 1

第一个问题是 \(n\) 很大，如何求 \(\sum_{p_i\le n}f(p_i)\)。再次运用合数的最小质因子小于等于根号，我们可以用小于等于根号的质数筛掉所有合数。

因为 \(f(p)=p^k-1\)，\(p^k\) 和 \(1\) 可以分别求和。

现在计算 \(p^k\)，注意到我们可以去筛完全积性函数 \(f'(x)=x^k\)，而不用管 \(f\)，因为我们需要求的只是质数位置的和。

设 \(g(i,j)\) 表示 \([1,i]\) 中的 \(f'(x)\) 之和，满足 \(x\) 是质数或者最小质因子大于等于 \(p_j\)（既目前我们已经筛完了前 \(j-1\) 个质数）。

\[\begin{aligned} g(i,j)=g(i,j-1)-p_{j-1}^k(g(\lfloor\frac{i}{p_{j-1}}\rfloor,j-1)-sp(j-2)) \end{aligned} \]

其中 \(sp(j)\) 代表小于等于 \(p_j\) 的质数的 \(p^k\) 之和，这个显然可以预处理。

注意到第一维都是 \(\lfloor\frac{n}{i}\rfloor\) 的形式，可以离散化成一个 \(2\sqrt n\) 大小的数组，具体见代码。

空间可以滚动数组，这个直接直接做时间复杂度是 \(\mathcal{O}(\frac{n}{\log n})\) 的，但是上面那个转移有意义需要 \(p_{j-1}^2\le i\)，所以复杂度是 \(\mathcal{O}(\frac{n^{\frac{3}{4}}}{\log n})\) 的。

Step 2

接下来就是求答案，设 \(S(n,x)\) 表示 \([1,n]\) 中最小质因子大于等于 \(p_x\) 的 \(f\) 之和，则：

\[S(n,x)=\sum_{p_x\le p_i\le n}f(p_i)+\sum_{i\ge x,p_i\le \sqrt n,p_i^e\le n}f(p_i^e)S(\lfloor\frac{n}{p_i^e}\rfloor,i+1)+f(p_i^e)[e\not=1] \]

需要注意一下 1 的问题，上面这些都是没有算 1 的，所以答案要加 1。

这是 lg 的模板题的代码。

#include<iostream>
#include<stdio.h>
#include<ctype.h>
#include<math.h>
#define int long long
#define N 200005
using namespace std;
inline int read(){
	int x=0,f=0; char ch=getchar();
	while(!isdigit(ch)) f|=(ch==45),ch=getchar();
	while(isdigit(ch)) x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return f?-x:x;
}
const int mo=1e9+7;
inline int qpow(int x,int b){
	int res=1;
	for(;b;x=x*x%mo,b>>=1) if(b&1) res=res*x%mo;
	return res;
}
int inv2=qpow(2,mo-2),inv6=qpow(6,mo-2);
inline void red(int &x){x>=mo?x-=mo:0;}
int gn,lim,pri[N],prs,notp[N],lis[N],cnt,spk[N],spk2[N],G[N][2],pk[N],le[N],ge[N];
inline void init(){
	for(int i=2;i<=lim;++i){
		if(!notp[i]){
			pri[++prs]=i;
			red(spk[prs]=spk[prs-1]+i);
			red(spk2[prs]=spk2[prs-1]+i*i%mo);
		}
		for(int j=1;j<=prs && i*pri[j]<=lim;++j){
			notp[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
	for(int l=1,r;l<=gn;l=r+1){
		r=gn/(gn/l);
		lis[++cnt]=gn/l;
		if(gn/l<=lim) le[gn/l]=cnt;
		else ge[r]=cnt;
		int v=lis[cnt]%mo;
		red(G[cnt][0]=v*(1+v)%mo*inv2%mo-1+mo);
		red(G[cnt][1]=v*(v+1)%mo*(2*v+1)%mo*inv6%mo-1+mo);
	}
}
inline int id(int x){return (x<=lim?le[x]:ge[gn/x]);}
int Fprime[N];
inline void calcFprime(){
	for(int k=1;k<=prs;++k){
		int T=pri[k]*pri[k];
		for(int i=1;lis[i]>=T;++i){
			int I=id(lis[i]/pri[k]);
			red(G[i][0]+=mo-pri[k]*(G[I][0]-spk[k-1]+mo)%mo);
			red(G[i][1]+=mo-pri[k]*pri[k]%mo*(G[I][1]-spk2[k-1]+mo)%mo); 
		}
	}
	for(int i=1;i<=cnt;++i) red(Fprime[i]=G[i][1]-G[i][0]+mo);
}
inline int Fk(int p){p%=mo;return p*(p-1)%mo;}
int F(int k,int n){
	if(pri[k]>n) return 0;
	int I=id(n);
	int ans=(Fprime[I]-spk2[k-1]+spk[k-1]+mo)%mo;
	for(int i=k;i<=prs && pri[i]*pri[i]<=n;++i){
		int pc=pri[i];
		for(int c=1;pc*pri[i]<=n;++c,pc*=pri[i]){
			red(ans+=Fk(pc)*F(i+1,n/pc)%mo+Fk(pc*pri[i])),red(ans);
		}
	}
	return ans;
}
signed main(){
	gn=read();
	lim=sqrt(gn)+10;
	init();
	calcFprime();
	cout<<(F(1,gn)%mo+mo+1)%mo;
	return 0;
}

Powerful Number 筛

如果一个数的幂指数都大于等于 2，那么这就是个 Powerful Number。显然每个 Powerful Number 都与 \(a^2b^3\) 相对应，积分一下可以证明 Powerful Number 的个数是 \(\mathcal{O}(\sqrt n)\) 的。

我们要求 \(S(n)=\sum_{i=1}^n f(i)\)。我们找到一个在质数处取值和 \(f\) 相等的积性函数 \(g\)。一般来说 \(g\) 的前缀和要好求，比如用可以杜教筛求或者可以 \(\mathcal{O}(1)\) 求之类的。

考虑 \(h=\frac{f}{g}\)，对于质数 \(p\)，则有：

\[\begin{aligned} f(p)&=g(1)h(p)+h(1)g(p)\\ h(p)&=f(p)-g(p)\\ h(p)&=0 \end{aligned} \]

又因为 \(h\) 是积性函数，所以 \(h\) 只有在 Powerful Number 处取值才不为 0。

回到 \(f\) 求和：

\[\begin{aligned} \sum_{i=1}^n f(i)&=\sum_{i=1}^n \sum_{d|i}g(d)h(\frac{i}{d})\\ &=\sum_{i}h(i)\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}g(j) \end{aligned} \]

因为 \(g\) 的前缀和可以快速求，所以我们可以爆搜所有的 Powerful Number（显然只用小于等于根号的质数即可）。

接下来是如何求 \(h(x)\) 的问题，显然我们只需要求出 \(h(p^k)\) 即可：

\[\begin{aligned} f(p^k)&=\sum_{i=0}^k h(p^i)g(p^{k-i})\\ h(p^k)&=f(p^k)-\sum_{i=0}^{k-1}h(p^i)g(p^{k-i}) \end{aligned} \]

暴力计算这个的复杂度是 \(\mathcal{O}(\sum_{p\in prime}\lfloor \log^2_p \sqrt n\rfloor)\)，推一下：

\[\begin{aligned} \mathcal{O}(\sum_{p\in prime}\lfloor \log^2_p \sqrt n\rfloor)&=\mathcal{O}(\sum_{p\in prime,p\le n^{\frac{1}{4}}}\lfloor \log^2_p \sqrt n\rfloor+\sum_{p\in prime,p>n^{\frac{1}{4}}} 1)\\ &\le\mathcal{O}(\frac{n^{\frac{1}{4}}}{\ln n^\frac{1}{4}}\cdot \log_2^2 \sqrt n+\frac{\sqrt n}{\ln \sqrt n})\\ &=\mathcal{O}(n^{\frac{1}{4}}\log n+\frac{\sqrt n}{\ln n}) \end{aligned} \]

所以这个时间应该是比 \(\mathcal{O}(\sqrt n)\) 还要快的。

空间复杂度瓶颈在于存 \(h(p^k)\)，空间复杂度为 \(\mathcal{O}(\frac{\sqrt n}{\log n}\cdot \log n)=\mathcal{O}(\sqrt n)\)。

如果 \(g\) 的前缀和可以 \(\mathcal{O}(1)\) 求的话，整个算法的复杂度可以做到 \(\mathcal{O}(\sqrt n)\)。如果 \(g\) 用的是杜教筛求的话，复杂度瓶颈在于杜教筛，即 \(\mathcal{O}(n^{\frac{2}{3}})\)。