【瞎口胡】數學 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) \]調換求和順序
記 \(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\) 和 \(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;
}