【瞎口胡】数学 6 杜教筛

杜教筛用于求一类积性函数的前缀和。其具体思想是考虑狄利克雷卷积的应用。

例题 Luogu P4213

求 \(\sum\limits_{i=1}^{n} \varphi(i)\) 和 \(\sum \limits_{i=1}^{n} \mu(i)\)。

多组数据，数据组数 \(T \leq 10\)，\(1 \leq n \leq 2^{31}-1\)

---

杜教筛的思想是，对于数论函数 \(f\)，找到另外一个数论函数 \(g\)。

考虑 \(f*g\) ：

\[\sum \limits_{i=1}^{n} (f*g)(i) \]

由狄利克雷卷积的定义可知原式等于

\[\sum \limits_{i=1}^{n} \sum \limits_{d\mid i} f(d)g(\dfrac nd) \]

调换求和顺序

\[\sum \limits_{i=1}^{n} g(i) \sum \limits_{j=1}^{\left \lfloor\frac nj \right \rfloor} f(j) \]

记 \(S(n) = \sum \limits_{i=1}^{n} f(i)\)，则原式可写作

\[\sum \limits_{i=1}^{n} g(i) S(\left \lfloor\dfrac ni \right \rfloor) \]

考虑 \(g(1)S(n)=\sum \limits_{i=1}^{n}g(i)S(\left \lfloor\dfrac ni\right \rfloor)-\sum \limits_{i=2}^{n}g(i)S(\left \lfloor\dfrac ni \right \rfloor)\)，又 \(\sum \limits_{i=1}^{n}g(i)S( \left \lfloor \dfrac ni \right \rfloor)=\sum \limits_{i=1}^{n}(f*g)(i)\)，那么则有

\[g(1)S(n)= \sum \limits_{i=1}^{n} (f*g)(i)-\sum \limits_{i=2}^{n}g(i)S(\left \lfloor\dfrac ni \right \rfloor) \]

这个式子有什么用呢？如果 \(g\) 和 \(f*g\) 的前缀和易求，那么我们可以通过整除分块来计算 \(g(1)S(n)\)，从而得到 \(S(n)\)，即要求的函数的前缀和。

让我们试试！

首先来算 \(\varphi\)。取 \(f = \varphi,g=I\)。根据上一篇中讲的，有 \(f * g = Id\)。

我们发现上面的式子变成了

\[S(n)= \dfrac{n(n-1)}{2}-\sum \limits_{i=2}^{n}S(\left \lfloor\dfrac ni \right \rfloor) \]

整除分块计算即可。

再来算 \(\mu\)。取 \(f=\mu,g=I\)。根据上一篇中讲的，有 \(f*g=\epsilon\)。

我们发现上面的式子变成了

\[S(n)= 1- \sum \limits_{i=2}^{n}S(\left \lfloor\dfrac ni \right \rfloor) \]

整除分块计算即可。

考虑上述算法的时间复杂度。我们知道，可以对 \(S(n)\) 记忆化求解。那么加上记忆化之后的时间复杂度为 \(O(n^{\frac 34})\)。如果我们通过线性筛求出 \(S(i)(i\leq n^{\frac 23})\) 再求解，那么时间复杂度为 \(O(n^{\frac 23})\)。证明较为复杂，感兴趣的读者可以自行查阅资料。这里给出参考。

Code

# include <bits/stdc++.h>

const int N=3000010,INF=0x3f3f3f3f,MAXN=N-10;
typedef long long ll;
int prime[N],psum;
ll phi[N],mu[N];
std::unordered_map <int,ll> caphi,camu; // 当然，也可以使用 n/x 最多有 \sqrt n 种取值的 trick
                                        // 在 min_25 筛里面可以看见这个 trick 的应用
bool vis[N]; 

inline int read(void){
	int res,f=1;
	char c;
	while((c=getchar())<'0'||c>'9')
		if(c=='-')f=-1;
	res=c-48;
	while((c=getchar())>='0'&&c<='9')
		res=res*10+c-48;
	return res*f;
}

inline void init(void){ // 预处理
	mu[1]=phi[1]=1;
	for(int i=2;i<=MAXN;++i){
		if(!vis[i]){
			prime[++psum]=i;
			phi[i]=i-1,mu[i]=-1;
		}
		for(int j=1;j<=psum&&i*prime[j]<=MAXN;++j){
			vis[i*prime[j]]=true;
			if(i%prime[j]==0){
				phi[i*prime[j]]=phi[i]*prime[j],mu[i*prime[j]]=0;
				break;
			}else{
				phi[i*prime[j]]=phi[i]*phi[prime[j]],mu[i*prime[j]]=-mu[i];
			}
		}
	}
	for(int i=2;i<=MAXN;++i){
		mu[i]+=mu[i-1],phi[i]+=phi[i-1];
	}
	return;
}
ll getmu(int x){
	if(x<=MAXN) // 预处理过了
		return mu[x];
	if(camu[x]) // 算过了
		return camu[x];
	ll res=1ll;
	for(ll l=2ll,r;l<=x;l=r+1ll){ // 整除分块
		r=x/(x/l);
		res-=(r-l+1ll)*getmu(x/l);
	}
//	if(x==2147483647){
//		printf("%lld\n",res);
//	}
	return camu[x]=res;
}
ll getphi(int x){
	if(x<=MAXN)
		return phi[x];
	if(caphi[x])
		return caphi[x];
	ll res=1ll*x*(x+1ll)/2ll;
	for(ll l=2ll,r;l<=x;l=r+1ll){
		r=x/(x/l);
		res-=(r-l+1ll)*getphi(x/l);
	}
	return caphi[x]=res;
}
int main(void){
//	freopen("4213.in","r",stdin);
	init();
	int T=read();
	while(T--){
		int n=read();
		printf("%lld %lld\n",getphi(n),getmu(n));
	}
	return 0;
}