1. 程式人生 > 其它 >【瞎口胡】數學 6 杜教篩

【瞎口胡】數學 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}^{\frac nj} f(j) \]

\(S(n) = \sum \limits_{i=1}^{n} f(i)\),則原式可寫作

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

考慮 \(g(1)S(n)=\sum \limits_{i=1}^{n}g(i)S(\dfrac ni)-\sum \limits_{i=2}^{n}g(i)S(\dfrac ni)\),又 \(\sum \limits_{i=1}^{n}g(i)S(\dfrac ni)=\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(\dfrac ni) \]

這個式子有什麼用呢?如果 \(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(\dfrac ni) \]

整除分塊計算即可。

再來算 \(\mu\)。取 \(f=\mu,g=I\)。根據上一篇中講的,有 \(f*g=\epsilon\)

我們發現上面的式子變成了

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

整除分塊計算即可。

考慮上述演算法的時間複雜度。我們知道,可以對 \(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;
}