简单数论(六)----------杜教筛
简单数论(六)----------杜教筛
Part I 用途
有的时候,我们需要求一些积性函数的前缀和,而且需要在低于线性时间复杂度内完成,这时就可以用到杜教筛啦
其实不是积性函数也可以但是要可以在接近线性的时间内求出来
Part II 杜教筛基础
假设我们要求积性函数f(n)的前缀和 $s(n)=\sum_{i=1}^{n}\ f(n)$
如果n够小我们可以考虑线性筛,但当n的范围为1e9呢
那杜教筛就光荣登场了
我们考虑在找到一个积性函数g
把f*g的前缀和算出来有 $ \sum _{i=1}^{n}(f*g)(i)= \sum _{i=1}^{n} \sum_{d|i}\ f(\frac{i}{d})g(d) $
可以发现这个式子等于 $ \sum_{d=1}^{n} g(d)s(\left \lfloor \frac{n}{d} \right \rfloor) $
其实上面的看不懂无所谓
其实只需要记住把我们要求的式子提出后的式子
即提出g(1)s(n)
此时有 $ g(1)s(n)= \sum _{i=1}^{n}(f*g)(i) - \sum_{d=2}^{n} g(d)s(\left \lfloor \frac{n}{d} \right \rfloor) $
当 $ m=n^{\frac{2}{3}} $ 时最优 复杂度为 O($n^{\frac{2}{3}} $)
当然如果是多组数据还是能多处理就多处理吧
一般情况下会用到的辅助的积性函数
1. g(n)=1
2.g(n) = n
Part III 例子
求 $ \sum _{i=1}^{n} \mu(i) $
构造g(n) = 1, 设f*g=h
有h(n)=[n=1]
故 s(n) = $ \sum _{i=1}^{n}(f*g)(i) - \sum_{d=2}^{n} s(\left \lfloor \frac{n}{d} \right \rfloor) $
显然 $ \sum _{i=1}^{n}(f*g)(i) $ =1
故s(n) =$ 1 - \sum_{d=2}^{n} s(\left \lfloor \frac{n}{d} \right \rfloor) $
整除分块+记忆化搜索即可
代码是洛谷模板的代码多了一个求欧拉函数的代码
#include<bits/stdc++.h> using namespace std; #define R register const int MAXN = 7000000 + 10; inline int read() { int f=1,x=0; char ch; do { ch=getchar(); if(ch=='-') f=-1; }while(ch<'0'||ch>'9'); do { x=(x<<3)+(x<<1)+ch-'0'; ch=getchar(); }while(ch>='0'&&ch<='9'); return f*x; } long long t; int n; unsigned long long phi[MAXN]; int mu[MAXN]; int prime[MAXN],prime_cnt; bool notprime[MAXN]; map<long long ,long long > w2; map<unsigned long long ,long long >w1; inline void init() { phi[1]=1;notprime[1]=1;mu[1]=1; for(R int i=2;i<=MAXN;i++) { if(notprime[i]==0){ ++prime_cnt; prime[prime_cnt]=i; phi[i]=i-1;mu[i]=-1; } for(R int j=1;j<=prime_cnt;j++) { if(i*prime[j]>MAXN) { break; } notprime[i*prime[j]] = 1; if(i%prime[j]==0) { phi[i*prime[j]]=phi[i]*prime[j]; mu[i*prime[j]]=0; break; } else { mu[i*prime[j]]=-mu[i]; phi[i*prime[j]]=phi[i]*(prime[j]-1); } } } for(R int i=2;i<=MAXN;i++) mu[i]+=mu[i-1]; for(R int i=2;i<=MAXN;i++) phi[i]+=phi[i-1]; } inline unsigned long long djs_phi(long long x) { if(x<=MAXN) return phi[x]; if(w1[x]) return w1[x]; long long sum=x*(x+1)/2; for(R int l=2,r;l<=x;l=r+1) { r=(x/(x/l)); sum-=(r-l+1)*djs_phi(x/l); } return w1[x]=sum; } inline long long djs_mu(long long x) { if(x<=MAXN) return mu[x]; if(w2[x]) return w2[x]; unsigned long long sum=1; for(R unsigned long long l=2,r;l<=x;l=r+1) { r=(x/(x/l)); sum-=(r-l+1)*djs_mu(x/l); } return w2[x]=sum; } int main() { init(); int t=read(); while(t--) { int n=read(); printf("%lld %lld\n",djs_phi(n),djs_mu(n)); } }