「专题总结」杜教筛
这个专题也是第二次学了。
而且上一次做的时候还是只背了$\mu$和$\varphi$的式子,居然能做的动题。。。
所以这次还是去大概研究了一下原理,重新把之前水过的题再做一遍(发现自己不会了)
这个专题里多数还是板子,但是《选数》挺有意思有两种解法,而《DZY Loves Maths IV》是真的大神题。
还是挺有意思的。
未强调的除法均向下取整。
先假装明白的总结两句:
狄利克雷$Dirichlet$卷积是啥?就是$(f*g)(i) = \sum\limits_{d|i} f(d) g(\frac{i}{d})$
这个形式的卷积就是狄利克雷卷积。平时好像也没啥用。
再提一嘴在杜教筛部分常用的函数$\mu,\varphi,\sigma,d$这四个是积性函数,而$I,id,\epsilon$是完全积性函数。
光看符号挺懵的,它们的含义是莫比乌斯函数,欧拉函数,约数和,约数个数,恒等函数(1),单位函数$id(i)=i$,元函数$\epsilon(x)=[x]$
然后杜教筛的通常套路就是拿这些东西和目标函数凑一凑看能不能卷出好求的玩意。
$\sum\limits_{i=1}^{n} (f*g)(i)=\sum\limits_{i=1}^{n} \sum\limits_{d|i} f(d)g(\frac{i}{d}) =\sum\limits_{d=1}^{n} \sum\limits_{i=1}^{\frac{n}{d}} f(i)$
如果单纯让你求右边那个东西,你会数论分块+快速求$f$前缀和来得到。而左边也是卷积函数的前缀和。
把右边$i=1$单独提出来移项,设$S(n)=\sum\limits_{i=1}^{n} f(i)$就可以得到:$g(1)S(n) =\sum\limits_{i=1}^{n} (f*g)(i) - \sum\limits_{i=2}^{n} g(i)S(\frac{n}{i})$
然后右边的$S$就可以递归求解了。也是数论分块+求$g$的前缀和。
所以如果我们能在合法的时间复杂度内求出$g$和$f*g$的前缀和,我们就可以快速的得到$f$的前缀和。复杂度$O(n^{0.75})$
然而如果加一发$unordered \ map$记忆化并且线筛预处理一部分$f$前缀和的话,复杂度就会下降,理论最优参数是$\frac{2}{3}$,复杂度下降为$O(n^{\frac{2}{3}})$
但是因为它本身的常数不小,所以线筛预处理的范围可以稍大一些,个人认为$0.7$这个参数也不错了。
常用结论:$\mu * I = \epsilon,\varphi * I =id,\mu * id =\varphi$
选数:
$Description:$
我们知道,从区间$[L,H]$(L和H为整数)中选取N个整数,总共有$(H-L)^N$种方案。小$Z$很好奇这样选出的数的最大公约数的规律,他决定对每种方案选出的N个整数都求一次最大公约数,以便进一步研究。然而他很快发现工作量太大了,于是向你寻求帮助。你的任务很简单,小 Z 会告诉你一个整数K,你需要回答他最大公约数刚好为 的选取方案有多少个。由于方案数较大,你只需要输出其除以$10^9+7$的余数即可。
$H-L \le 10^5$,其余均为$10^9$
解法一:暴力杜教筛
作为一个思维僵化的文化课选手我自然会疯狂搬套路。
$\sum\limits_{A_i=L}^{H} [gcd(A)=K]$
$=\sum\limits_{A_i=\frac{L-1}{K}+1}^{\frac{H}{K}} [gcd(A)]$
$=\sum\limits_{A_i=\frac{L-1}{K}+1}^{\frac{H}{K}} \sum\limits_{d|A} \mu(d)$
$=\sum\limits_{d} \mu(d) \sum\limits_{A_i=\frac{L-1}{Kd}+1}^{\frac{H}{Kd}} 1$
$=\sum\limits_{d} \mu(d) (\frac{H}{Kd}-\frac{L-1}{Kd}-1)^N$
然后数论分块+杜教筛$\mu$解决问题。
1 #include<bits/stdc++.h> 2 using namespace std; 3 unordered_map<int,int>Mu; 4 int dp[100005],L,R,k,n,C,mu[10000005],M,p[3333333],pc,ans;char np[10000005]; 5 #define mod 1000000007 6 int pow(int b,int t,int a=1){for(;t;t>>=1,b=1ll*b*b%mod)if(t&1)a=1ll*a*b%mod;return a;} 7 int MU(int n){ 8 if(n<M)return mu[n];if(Mu[n])return Mu[n];int a=1; 9 for(int i=2,l,N;N=n/i,i<=n;i=l+1)l=n/N,a-=MU(N)*(l-i+1);return Mu[n]=a; 10 } 11 main(){ 12 cin>>n>>k>>L>>R;L--;R/=k;L/=k;M=pow(R,0.7);mu[1]=1; 13 for(int i=2;i<M;++i){ 14 if(!np[i])p[++pc]=i,mu[i]=-1; 15 for(int j=1,x;j<=pc&&(x=p[j]*i)<M;++j) 16 if(i%p[j])np[x]=1,mu[x]=-mu[i]; 17 else{np[x]=1;break;} 18 }for(int i=2;i<M;++i)mu[i]+=mu[i-1]; 19 for(int i=1,l,le,ri;le=L/i,ri=R/i,i<=R;i=l+1) l=min(le?L/le:mod,R/ri),ans=(ans+mod+1ll*(MU(l)-MU(i-1))*pow(ri-le,n))%mod; 20 cout<<ans<<endl; 21 }
解法二:然而更加简单优雅快速省空间代码短的方法是dp。
因为$H-L$很小,所以不难意识到,对于任意大于$H-L$的数,如果它是$[L,H]$中某个数的因子,那一定就不是其它数的因子。
也就是这些大数作为因子只能出现在特定一个数内。所以如果最终数列的$gcd$是这个大数的话我们就能确定整个数列都是同一个数。
而既然整个数列都相同了,那么它们的$gcd$就是这个数列中的任意一项。
所以,如果整个数列都相同且方案合法,当且仅当整个数列都是$K$,且$L \le K \le H$。只有这一种方案。
剩下的所有情况都保证数列并不是完全由同一个数组成的了,那么它们的$gcd$必然小于等于$H-L$。
这样的话暴力枚举$gcd$用莫比乌斯函数容斥算一下就是最终答案。
1 #include<bits/stdc++.h> 2 using namespace std; 3 int dp[100005],L,R,k,n,C; 4 #define mod 1000000007 5 int pow(int b,int t,int a=1){for(;t;t>>=1,b=1ll*b*b%mod)if(t&1)a=1ll*a*b%mod;return a;} 6 main(){ 7 cin>>n>>k>>L>>R;((--L)/=k)++;R/=k; 8 if(R<L)return puts("0"),0;C=R-L+1; 9 for(int i=1;i<=C;++i){int x=R/i-(L-1)/i;if(x>0)dp[i]=(pow(x,n)-x+mod)%mod;} 10 for(int i=C;i;--i)for(int j=i+i;j<=C;j+=i)dp[i]=(dp[i]-dp[j]+mod)%mod; 11 cout<<dp[1]+(L==1)<<endl; 12 }
神犇和蒟蒻:
$Description:$
很久很久以前,有一只神犇叫yzy;很久很久之后,有一只蒟蒻叫lty;请你输出一个整数$A=\sum_{i=1}^{N}\mu (i^2)$;请你输出一个整数$B=\sum_{i=1}^{N}\varphi (i^2)$;$N \le 10^9$
题面挺有意思?蒟蒻是不是就是指第一问的答案一定是$1$啊。
然后是神犇。考虑这个$\varphi(i^2)=\varphi(i) \times i$
用各种积性函数卷一卷。想到求普通欧拉函数时总卷$I$,带进来试一下得到$(f*id)(x)=x^2$
前缀和是自然数幂和。它的值等于$\frac{n(n+1)(2n+1)}{6}$。
两个前缀和会求了,就可以筛了。
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define mod 1000000007 4 #define S 10000007 5 unordered_map<int,int>Sum; 6 int n,M,p[S],pc,phi[S];char np[S]; 7 int cal(int n){return n*(n+1ll)%mod*(2*n+1)%mod*166666668%mod;} 8 int Cal(int n){return n*(n+1ll)/2%mod;} 9 int sum(int n){ 10 if(n<M)return phi[n];if(Sum[n])return Sum[n];int a=cal(n); 11 for(int i=2,N,l;N=n/i,i<=n;i=l+1)l=n/N,a=(a+mod-sum(N)*(Cal(l)-Cal(i-1)+mod+0ll)%mod)%mod;return Sum[n]=a; 12 } 13 int main(){ 14 cin>>n;M=pow(n,0.7);phi[1]=1; 15 for(int i=2;i<M;++i){ 16 if(!np[i])phi[p[++pc]=i]=i*(i-1ll)%mod; 17 for(int j=1,x;j<=pc&&(x=p[j]*i)<M;++j) 18 if(i%p[j])phi[x]=1ll*phi[i]*p[j]%mod*(p[j]-1)%mod,np[x]=1; 19 else{phi[x]=1ll*phi[i]*p[j]%mod*p[j]%mod;np[x]=1;break;} 20 }for(int i=2;i<M;++i)phi[i]=(phi[i-1]+phi[i])%mod; 21 puts("1");cout<<sum(n)<<endl; 22 }
DZY Loves Math IV:
$Description:$
给定$n,m$,求$\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{m} \varphi(ij)$
$n \le 10^5,m \le 10^9$
神题啊。没有用到反演,但是需要熟练掌握欧拉函数的性质
$n,m$范围不一样大,一定有问题,于是就在最外层枚举$n$,设$S(n,m)=\sum\limits_{i=1}^{m} \varphi(ni)$。总问题是$\sum\limits_{i=1}^{n} S(i,m)$
现在问题转化为求$S(n,m)$。下面的$n$的含有发生了变化(就是总问题里枚举的$i$而不是原来的$n$了)
首先对于欧拉函数它有一些熟悉但不很常用的性质。
一个因子出现多次的话可以提出去一个,所以我们只保留$n$中每个因子第一次出现的累乘,设为$x$
而对于多余的因子,我们也单独提出来,称为$r$
于是$S(n,m)=r\sum\limits_{i=1}^{m} \varphi(xi)$
沿用套路,在$x$和$i$中的重复因子也可以直接提出来,提出来之后两者就互质了,根据积性函数,可以拆开。
$=r\sum\limits_{i=1}^{m} gcd(x,i) \varphi(\frac{x}{gcd(i,x)}) \varphi(i)$
然后再套用一个我们经常见到但是不知道怎么用的式子,就是$p=\sum\limits_{d|p} \varphi(d)$.用它替换掉那个$gcd$
$=r\sum\limits_{i=1}^{m} \sum\limits_{d|i \ and \ d|x} \varphi(d) \varphi(i) \varphi(\frac{x}{gcd(i,x)})$
我们又知道因为原来$x$中每个因子只有一次,所以$\frac{x}{gcd(i,x)}$与$gcd(i,x)$是互质的,由于积性函数而合并。
$=r\sum\limits_{i=1}^{m} \varphi(i) \sum\limits_{d|i \ and \ d|x} \varphi(\frac{x}{d})$
改变枚举顺序:
$=r\sum\limits_{d=1}^{m} \varphi(\frac{x}{d}) \sum\limits_{i=1}^{\frac{m}{d}} \varphi(di)$
好像化不动了?但是突然发现最后一个和式有些熟悉?那不是和$S$的定义一样吗?递归求解。
$=r\sum\limits_{d=1}^{m} \varphi(\frac{x}{d}) S(e,\frac{m}{d})$
枚举约数+递归求解,杜教筛解决$\varphi$的前缀和以便处理$n=1$的递归边界。
因为$n$较大的$S(n,m)$会调用$n$较小的,所以复杂度只用算一次杜教筛的复杂度,即$O(m^{\frac{2}{3}})$
然后内部的枚举约数复杂度是$O(n \sqrt{m})$的。理论复杂度很大,建议线筛预处理出每个数的最小因子然后就可以$2^c$枚举约数了。常数会小一些?
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define Z 10000005 4 #define mod 1000000007 5 unordered_map<int,int>Phi,s[Z>>6]; 6 int n,m,M,phi[Z],p[Z],pc,ans,lp[Z];char np[Z]; 7 int PHI(int n){ 8 if(n<M)return phi[n];if(Phi[n])return Phi[n];int a=n*(n+1ll)/2%mod; 9 for(int i=2,l,N;N=n/i,i<=n;i=l+1)l=n/N,a=(a+mod-(l-i+1ll)*PHI(N)%mod)%mod; 10 return Phi[n]=a; 11 }int aphi(int x){return (phi[x]-phi[x-1]+mod)%mod;} 12 int S(int n,int m){ 13 if(n==1)return PHI(m);if(m==1)return aphi(n);if(!m)return 0;if(s[n][m])return s[n][m]; 14 int x=n,w=1,y=1,L,a=0;vector<int>d; 15 while(x^1){L=lp[x];w*=L;x/=L;d.push_back(L);while(x%L==0)x/=L,y*=L;} 16 for(int s=0,e;e=1,s<1<<d.size();++s){ 17 for(int i=0;i<d.size();++i)if(s&1<<i)e*=d[i]; 18 a=(a+1ll*aphi(w/e)*S(e,m/e))%mod; 19 }return s[n][m]=1ll*a*y%mod; 20 } 21 int main(){ 22 cin>>n>>m;M=max(n+1.1,pow(M,0.7));phi[1]=1; 23 for(int i=2;i<M;++i){ 24 if(!np[i])phi[p[++pc]=i]=i-1,lp[i]=i; 25 for(int j=1,x;j<=pc&&(x=p[j]*i)<M;++j) 26 if(i%p[j])phi[x]=phi[i]*(p[j]-1),np[x]=1,lp[x]=p[j]; 27 else{np[x]=1;phi[x]=phi[i]*p[j];lp[x]=p[j];break;} 28 }for(int i=2;i<M;++i)phi[i]=(phi[i-1]+phi[i])%mod; 29 for(int i=1;i<=n;++i)ans=(ans+S(i,m))%mod;cout<<ans<<endl; 30 }
Lucas的数论:
$Description:$
去年的Lucas非常喜欢数论题,但是一年以后的Lucas却不那么喜欢了。
在整理以前的试题时,发现了这样一道题目“求Sigma(f(i)),其中1<=i<=N”,其中 表示i的约数个数。他现在长大了,题目也变难了。求如下表达式的值:其中f(ij)表示ij的约数个数。他发现答案有点大,只需要输出模1000000007的值。$1 \le n \le 10^9$
《约数个数和》数据加强版。还是搬结论反演然后套个式子杜教筛莫比乌斯函数前缀和配合数论分块就可以。
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define S 2332333 4 #define mod 1000000007 5 unordered_map<int,int>Mu; 6 int mu[S],p[S],pc,n,ans,M;char np[S]; 7 int sum(int n,int a=1){ 8 if(n<M)return mu[n];if(Mu[n])return Mu[n]; 9 for(int i=2,l,N;N=n/i,i<=n;i=l+1)l=n/N,a-=sum(N)*(l-i+1);return Mu[n]=a; 10 } 11 int cal(int n,int a=0){for(int i=1,l,N;N=n/i,i<=n;i=l+1)l=n/N,a=(a+1ll*N*(l-i+1))%mod;return 1ll*a*a%mod;} 12 int main(){ 13 cin>>n;M=pow(n,0.7);mu[1]=1; 14 for(int i=2;i<M;++i){ 15 if(!np[i])mu[p[++pc]=i]=-1; 16 for(int j=1,x;j<=pc&&(x=p[j]*i)<M;++j) 17 if(i%p[j])mu[x]=-mu[i],np[x]=1; 18 else{np[x]=1;break;} 19 }for(int i=2;i<M;++i)mu[i]+=mu[i-1]; 20 for(int i=1,l,N;N=n/i,i<=n;i=l+1)l=n/N,ans=(ans+1ll*(sum(l)-sum(i-1)+mod)*cal(N))%mod; 21 cout<<ans<<endl; 22 }
Sum:
$Description:$
求$\sum\limits_{i=1}^{n} \mu(i)$和$\sum\limits_{i=1}^{n} \varphi(i)$
$n \le 10^9$
板子,不说了。
因为调用轨迹一致,所以用$pair$同时解决两个会快好多。
1 #include<bits/stdc++.h> 2 using namespace std; 3 unordered_map<int,long long>Phi; 4 unordered_map<int,int>Mu; 5 #define S 3333333 6 int mu[S],p[S],pc;char np[S];long long phi[S]; 7 int sum_mu(int n){ 8 if(n<S)return mu[n]; 9 if(Mu.find(n)!=Mu.end())return Mu[n]; 10 int ans=1; 11 for(unsigned int i=2,l,N;N=n/i,i<=n;i=l+1)l=n/N,ans-=sum_mu(N)*(l-i+1); 12 return Mu[n]=ans; 13 } 14 long long sum_phi(int n){ 15 if(n<S)return phi[n]; 16 if(Phi.find(n)!=Phi.end())return Phi[n]; 17 long long ans=n*(n+1ll)>>1; 18 for(unsigned int i=2,l,N;N=n/i,i<=n;i=l+1)l=n/N,ans-=sum_phi(N)*(l-i+1); 19 return Phi[n]=ans; 20 } 21 int main(){ 22 int t,n;mu[1]=phi[1]=1; 23 for(int i=2;i<S;++i){ 24 if(!np[i])p[++pc]=i,phi[i]=i-1,mu[i]=-1; 25 for(int j=1,x;j<=pc&&(x=p[j]*i)<S;++j) 26 if(i%p[j])np[x]=1,mu[x]=-mu[i],phi[x]=phi[i]*phi[p[j]]; 27 else{np[x]=1;phi[x]=phi[i]*p[j];break;} 28 }for(int i=2;i<S;++i)phi[i]+=phi[i-1],mu[i]+=mu[i-1]; 29 cin>>t;while(t--)cin>>n,cout<<sum_phi(n)<<' '<<sum_mu(n)<<endl; 30 }