杜教筛
比更快还更快
杜教筛
字符串好难啊,还是数论比较有趣.
从洛谷试炼场找了一个莫反的题做,非常套路,熟练的做完之后发现线性复杂度过不了,所以不得不学一下这种奇妙的筛法了.
设要求前缀和的函数$f$,是一个积性函数.
凭借灵感或者是经验找出另一个积性函数$g$来,两个函数卷一下.
$h=f\times g$
$\sum_{i=1}^nh(i)=\sum_{i=1}^n\sum_{d|i}f(\frac{i}{d})g(d)$
莫反常见套路:改变枚举顺序;
$\sum_{i=1}^nh(i)=\sum_{d=1}^ng(d)\sum_{k=1}^{\frac{n}{d}}f(k)$
$\sum_{i=1}^nh(i)=\sum_{d=1}^ng(d)S(\frac{n}{d})$
$\sum_{i=1}^nh(i)=g(1)S(n)+\sum_{d=2}^ng(d)S(\frac{n}{d})$
$S(n)g(1)=\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)S(\frac{n}{d})$
那么$\sum_{i=1}^nh(i)$怎么求呢?积性函数的卷积还是积性函数,显然可以线性筛!
不过要是这样我们为什么不直接筛$f$…
这就是考察技巧的时候了,需要选择一个巧妙的函数$g$使得$h$的前缀和非常好算,举几个例子:
$$\sum_{i=1}^n\varepsilon(n)=[n>=1]$$
$$\sum_{i=1}^n1=n$$
$$\sum_{i=1}^ni=\frac{n(n+1)}{2}$$
$$\sum_{i=1}^ni^2=\frac{n(n+1)(2n+1)}{6}$$
$$\sum_{i=1}^n i^3=\frac{n^2(n+1)^2}{4}$$
自然数的$k$次幂求和是一个关于$n$的$k+1$次多项式.可以用高斯消元暴力求每一项的系数,也可以用拉格朗日插值法做…好像说远了,回到杜教筛吧.
再说一点:今天烜神仙问我为什么一定就是$k+1$次多项式?目瞪口呆.jpg,因为之前学的课件上也只是写了”通过观察可以发现”,但是经过不懈的yy,终于想到一个靠谱一点的证明:
首先明确一点:N维空间中的直线解析式是一个N次多项式,没有什么疑问吧.
对于k次方求和,我们将它表示到k+1维空间中,坐标为(n,n,…,n^k),构成了一个$n+1$维”正方体”,那么它的体对角线长度就是
$$\sqrt[n+1]{(n+1)x^{n+1}}=\sqrt[n+1]{n+1}x$$
是一个关于$x$的线性函数,而且每个立方体体对角线的角度都是相同的,沿着体对角线的方向这些点就构成了一条直线,它们的和就是在这条直线上一些等距离的点的求和,是一个k+1次多项式.
这么说可能不是很好理解,举两个简单的例子:
$i^1$: $i^2$:用电脑不大好画…
$$S(n)g(1)=\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)S(\frac{n}{d})$$
前半部分$O(1)$求,后半部分除法分块,递归调用这个式子,一定要记忆化！
单纯考察杜教筛的题目不是很多，但是部分莫反的最后20-40分必须要低于线性才能做，每当看到$n<=10^{9}$或者$n<=10^{10}$的时候… …
复杂度分析:不会积分.jpg，所以直接抄结论了,首先线性筛预处理出$n^{\frac{2}{3}}$的答案,小于时直接返回,记忆化时不用$map$,而是用$f[x]$表示$f[n/x]$,可以做到$O(n^{\frac{2}{3}}$的复杂度.出于这样的分析,我们可以发现需要记忆化的东西并不是很多,即使是多组询问也很快,复杂度几乎不变.实际运用中为了更快,可以在在不影响复杂度的情况下多线筛一些.注意…除法分块应从2开始分，至于为什么…看定义…
莫比乌斯函数求和：https://www.51nod.com/Challenge/Problem.html#!#problemId=1244
题意概述:求$\sum_{i=a}^{b}\mu(i),a,b<=10^{10}$
首先给莫比乌斯函数卷上一个别的函数,使得前缀和好求.
这个函数还是比较好找的,毕竟 $\mu$ 的某个定义式就已经给了我们一些启发,
$$\sum_{d|n}\mu(d)=\varepsilon(n)$$
你可能说这哪里卷积了? 我们可以假装它乘了一个一嘛.
$$\sum_{d|n}\mu(d)1(\frac{n}{d})=\varepsilon(n)$$
然后套模板就好了,不过这道题存在的意义就是为了讲解模板写法.
1 # include <cstdio> 2 # include <iostream> 3 # include <cstring> 4 # include <string> 5 # define ll long long 6 # define R register int 7 8 using namespace std; 9 10 const int maxn=5000000; 11 int N=5000000,pri[maxn],h,vis[maxn+4],dvis[100005]; 12 ll m[maxn+3],f[100005]; 13 ll a,b,n; 14 15 void init() 16 { 17 m[1]=1; 18 for (R i=2;i<=N;++i) 19 { 20 if(!vis[i]) pri[++h]=i,m[i]=-1; 21 for (R j=1;j<=h&&i*pri[j]<=N;++j) 22 { 23 vis[ i*pri[j] ]=1; 24 if(i%pri[j]==0) break; 25 m[ i*pri[j] ]=-m[i]; 26 } 27 } 28 for (R i=2;i<=N;++i) 29 m[i]+=m[i-1]; 30 } 31 32 ll get_mu (ll x) { return (x<=N)?m[x]:f[n/x]; } 33 34 void solve (ll x) 35 { 36 if(x<=N) return; 37 ll i=2,l,t=n/x; 38 if(dvis[t]) return; 39 dvis[t]=1; 40 f[t]=1; 41 while(i<=x) 42 { 43 l=x/(x/i); 44 solve(x/i); 45 f[t]-=(l-i+1)*get_mu(x/i); 46 i=l+1; 47 } 48 } 49 50 int main() 51 { 52 scanf("%lld%lld",&a,&b); 53 if(a>b) swap(a,b); 54 a--; 55 if(a<0) a=0; 56 init(); 57 ll ans=0; 58 if(b<=N) 59 { 60 ans+=m[b]; 61 n=a; 62 if(a<=N) ans-=m[a]; else solve(a),ans-=f[1]; 63 } 64 else 65 { 66 n=b; 67 solve(b); 68 ans+=f[1]; 69 n=a; 70 if(a<=N) ans-=m[a]; 71 else 72 { 73 memset(dvis,0,sizeof(dvis)); 74 solve(a); 75 ans-=f[1]; 76 } 77 } 78 printf("%lld",ans); 79 return 0; 80 }
1244
欧拉函数求和:https://www.51nod.com/Challenge/Problem.html#!#problemId=1239
题意概述:求$\sum_{i=1}^{n}\varphi(i),n<=10^{10}$
考虑之前证明过的一个性质,可以快速的想出合理的卷积函数.
$$id=\varphi \times 1$$
然后就是套板子啦.注意因为欧拉函数的和的范围会比较大,要开$ll$,有的时候要用慢速乘.
1 # include <cstdio> 2 # include <iostream> 3 # include <cstring> 4 # include <map> 5 # define R register int 6 # define ll long long 7 8 using namespace std; 9 10 const int maxn=5000006; 11 const ll mod=1000000007; 12 int N=5000000,pri[maxn],h,phi[maxn]; 13 ll n,vis[100005],f[100005]; 14 map <ll,ll> m; 15 16 void init() 17 { 18 phi[1]=1; 19 for (R i=2;i<=N;++i) 20 { 21 if(!phi[i]) pri[++h]=i,phi[i]=i-1; 22 for (R j=1;j<=h&&i*pri[j]<=N;++j) 23 { 24 if(i%pri[j]) phi[ i*pri[j] ]=1LL*phi[i]*(pri[j]-1)%mod; 25 else 26 { 27 phi[ i*pri[j] ]=1LL*phi[i]*pri[j]%mod; 28 break; 29 } 30 } 31 } 32 for (R i=1;i<=N;++i) 33 phi[i]=(phi[i]+phi[i-1])%mod; 34 } 35 36 ll get_phi (ll x) 37 { 38 return (x<=N)?phi[x]:f[n/x]; 39 } 40 41 ll mul (ll a,ll b) 42 { 43 ll s=0; 44 while(b) 45 { 46 if(b&1LL) s=(s+a)%mod; 47 a=(a+a)%mod; 48 b>>=1; 49 } 50 return s; 51 } 52 53 void solve (ll x) 54 { 55 if(x<=N) return; 56 ll i=2,l,t=n/x; 57 if(vis[t]) return; 58 f[t]=0; 59 vis[t]=1; 60 if(m[x]) return; 61 while(i<=x) 62 { 63 l=x/(x/i); 64 solve(x/i); 65 f[t]=(f[t]+(l-i+1)*get_phi(x/i)%mod)%mod; 66 i=l+1; 67 } 68 ll y=x+1; 69 if(y%2LL) x>>=1; else y>>=1; 70 f[t]=((mul(x,y)-f[t])%mod+mod)%mod; 71 } 72 73 int main() 74 { 75 scanf("%lld",&n); 76 if(n<=N) N=n; 77 init(); 78 if(n<=N) printf("%d",phi[n]); 79 else 80 { 81 solve(n); 82 printf("%lld",f[1]); 83 } 84 return 0; 85 }
1239
Sum:https://www.lydsy.com/JudgeOnline/problem.php?id=3944
题意概述:上述两个题的二合一.
1 # include <cstdio> 2 # include <iostream> 3 # include <cstring> 4 # define R register int 5 # define ll long long 6 7 using namespace std; 8 9 const int maxn=4000000; 10 int t,N=4000000,pri[maxn+2],h,vis[100005]; 11 ll n,mu[maxn+2],phi[maxn+2],m[100005],p[100005]; 12 13 void init() 14 { 15 mu[1]=1,phi[1]=1; 16 for (R i=2;i<=N;++i) 17 { 18 if(!phi[i]) phi[i]=i-1,pri[++h]=i,mu[i]=-1; 19 for (R j=1;j<=h&&i*pri[j]<=N;++j) 20 { 21 if(i%pri[j]==0) 22 { 23 phi[ i*pri[j] ]=phi[i]*pri[j]; 24 break; 25 } 26 phi[ i*pri[j] ]=phi[i]*phi[ pri[j] ]; 27 mu[ i*pri[j] ]=-mu[i]; 28 } 29 } 30 for (R i=1;i<=N;++i) 31 mu[i]+=mu[i-1],phi[i]+=phi[i-1]; 32 } 33 34 ll smu (ll x) { return x<=N?mu[x]:m[n/x]; } 35 ll sphi (ll x) { return x<=N?phi[x]:p[n/x]; } 36 37 void solve (ll x) 38 { 39 if(x<=N) return; 40 ll i=2,l,t=n/x; 41 if(vis[t]) return; 42 vis[t]=true; 43 p[t]=(x*(x+1))/2LL,m[t]=1; 44 while(i<=x) 45 { 46 l=x/(x/i); 47 solve(x/i); 48 p[t]-=sphi(x/i)*(l-i+1); 49 m[t]-=smu(x/i)*(l-i+1); 50 i=l+1; 51 } 52 } 53 54 int main() 55 { 56 scanf("%d",&t); 57 init(); 58 while(t--) 59 { 60 scanf("%lld",&n); 61 if(n<=N) 62 printf("%lld %lld\n",phi[n],mu[n]); 63 else 64 { 65 memset(vis,0,sizeof(vis)); 66 solve(n); 67 printf("%lld %lld\n",p[1],m[1]); 68 } 69 } 70 return 0; 71 }
3944
许多莫反的题都可以用杜教筛来优化,举个例子:
GCDSUM：https://www.51nod.com/Challenge/Problem.html#!#problemId=1237
题意概述:求$\sum_{i=1}^n\sum_{j=1}^n(i,j),n<=10^{10}$.
套路题，枚举答案，前面除法分块，后面杜教筛求和。注意，因为这道题其实是变相的多组询问，所以不要用之前的技巧(要多次清空数组)，反而是unordered_map快一些.
$$\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}} \varphi(i)$$
$$\sum_{d=1}^nS(\frac{n}{d})d$$
这样化出来的式子会少算一半,答案乘二后,如果$i,j$相同又会多算一次,减掉就好了.
1 # include <cstdio> 2 # include <iostream> 3 # include<tr1/unordered_map> 4 # define ll long long 5 # define R register int 6 7 using namespace std; 8 9 const ll mod=1000000007; 10 const int maxn=5000006; 11 int N=5000000,h,cnt,pri[maxn],phi[maxn]; 12 ll n; 13 tr1::unordered_map <ll,ll> m; 14 15 void init() 16 { 17 phi[1]=1; 18 for (R i=2;i<=N;++i) 19 { 20 if(!phi[i]) pri[++h]=i,phi[i]=i-1; 21 for (R j=1;j<=h&&i*pri[j]<=N;++j) 22 { 23 if(i%pri[j]) 24 phi[ i*pri[j] ]=1LL*phi[i]*(pri[j]-1)%mod; 25 else 26 { 27 phi[ i*pri[j] ]=1LL*phi[i]*pri[j]%mod; 28 break; 29 } 30 } 31 } 32 for (R i=1;i<=N;++i) 33 phi[i]=(phi[i]+phi[i-1])%mod; 34 } 35 36 ll mul (ll a,ll b) 37 { 38 ll s=0; 39 while(b) 40 { 41 if(b&1LL) s=(s+a)%mod; 42 a=(a+a)%mod; 43 b>>=1LL; 44 } 45 return s; 46 } 47 48 ll get_phi (ll x) { return (x<=N)?phi[x]:m[x]; } 49 50 ll S (ll x) 51 { 52 ll y=x+1; 53 if(y%2) x=x/2; else y=y/2; 54 return mul(x,y); 55 } 56 57 ll solve (ll x) 58 { 59 if(x<=N) return phi[x]; 60 ll i=2,l,ans=0; 61 if(m[x]) return m[x]; 62 while(i<=x) 63 { 64 l=x/(x/i); 65 solve(x/i); 66 ans=(ans+(l-i+1)*get_phi(x/i)); 67 i=l+1; 68 } 69 ans=((S(x)-ans)%mod+mod)%mod; 70 return m[x]=ans; 71 } 72 73 int main() 74 { 75 scanf("%lld",&n); 76 init(); 77 ll i=1,l,ans=0,a; 78 while(i<=n) 79 { 80 l=n/(n/i); 81 a=solve(n/i); 82 ans=(ans+mul((S(l)-S(i-1)+mod)%mod,a))%mod; 83 i=l+1; 84 } 85 ans=ans*2%mod; 86 ans=(ans-S(n)+mod)%mod; 87 printf("%lld",ans); 88 return 0; 89 }
1237
—shzr