杜教筛小结
首先的首先，%%%%%%杜教
……咳咳，言归正传，现在来总结一下杜教筛……
杜教筛大概是一种用于快速处理积性函数前缀和的有力工具，主要利用了狄利克雷卷积的知识
我们用高中均值不等式的复杂度分析可以得到其复杂度是$O(n^{\frac{2}{3}})$的
具体的讲解我就不写了……留几篇比较好的博客。
https://www.cnblogs.com/abclzr/p/6242020.html
这里总结一下我做过的题目吧……
由于太水以至于没打……
可以当做板子练习题的题目。
深刻的记住了一个式子……之前一直没有这种意识，有一些转换的确很显然但是就是没想过
\[d\mid ij\Leftrightarrow \frac{d}{(i,d)}\mid j\]
……自己打公式太麻烦了，放弃……
……其实怎么演都能演出来，关键就看第一步转化能不能出来
代码：
1 #include <cstdio> 2 #include <map> 3 #include <cstring> 4 using namespace std; 5 #define mod 1000000007 6 #define N 1000000 7 #define LL long long 8 int n,prime[N/10],tot,vis[N+10],mu[N+10]; 9 inline void intn() 10 { 11 register int i,j,tmp; 12 for(mu[1]=1,i=2;i<=N;++i) 13 { 14 if(!vis[i])prime[++tot]=i,mu[i]=-1; 15 for(j=1;j<=tot&&(tmp=i*prime[j])<=N;++j) 16 {vis[tmp]=1;if(i%prime[j]==0){mu[tmp]=0;break;}mu[tmp]=-mu[i];} 17 } 18 for(i=1;i<=N;++i)mu[i]+=mu[i-1]; 19 } 20 map<int,int> mmp; 21 inline int getmu(int x) 22 { 23 if(x<=N)return mu[x]; 24 if(mmp.count(x))return mmp[x]; 25 register int i,last,ret=1; 26 for(i=2;i<=x;i=last+1) 27 last=x/(x/i),ret=(ret-((LL)(last-i+1)*getmu(x/i)%mod))%mod; 28 return mmp[x]=ret; 29 } 30 inline int calc2(int x) 31 { 32 register int i,last,ret=0; 33 for(i=1;i<=x;i=last+1) 34 last=x/(x/i),ret=(ret+((LL)(last-i+1)*(x/i)%mod))%mod; 35 return (LL)ret*ret%mod; 36 } 37 int main() 38 { 39 // freopen("Ark.in","r",stdin); 40 register int i,last,ans=0; 41 scanf("%d",&n);intn(); 42 for(i=1;i<=n;i=last+1) 43 last=n/(n/i),ans=(ans+((LL)(getmu(last)-getmu(i-1))*calc2(n/i)%mod) )%mod; 44 printf("%d\n",(ans+mod)%mod); 45 }
然后这题我们发现n比较小，m却比较大
所以考虑与枚举n有关的算法……
设$S(n,m)=\sum _{i=1} ^{n} \varphi (i*n)$
那么我们知道$\varphi$是个积性函数……所以上来质因数分解一发
设$n=\prod p_{i} ^{a_{i}}$
那么有$S(n,m)=p_{i}^{a_{i}-1} * S(\frac{n}{p_{i}^{a_{i}-1}},m)$
然后我们把p都干下去了之后，
考虑$\varphi(p*q)=\varphi(p)*\varphi(q),(p,q)==1$
那么我们有$S(n,m)=\varphi(p_{i}) * \sum _{i=1} ^{m} \varphi(\frac{n}{p_{i}}*i)*[(i,p_{i})==1]+p_{i}* \sum _{i=1} ^{m} \varphi(n*\frac{i}{p_{i}})*[p_{i}|i]$
又因为$\varphi(p)==p-1$
所以我们可以凑个整……
$S(n,m)=\varphi(p_{i}) * \sum _{i=1} ^{m} \varphi(\frac{n}{p_{i}}*i)+\sum _{i=1} ^{ \left \lfloor \frac{m}{p_{i}} \right \rfloor } \varphi(n*i)$
$S(n,m)=S(\frac{n}{p_{i}},m)+S(n, \left \lfloor \frac{m}{p_{i}} \right \rfloor)$
这样我们就可以记忆化搞了
边界条件$n==1$就是$\varphi$前缀和，杜教筛之。
1 #include <cstdio> 2 #include <cstring> 3 #include <map> 4 using namespace std; 5 #define N 1000000 6 #define LL long long 7 #define mod 1000000007 8 int prime[N/10],tot,mind[N+10],phi[N+10],sum[N+10]; 9 inline void intn() 10 { 11 register int i,j,tmp; 12 phi[1]=1; 13 for(i=2;i<=N;++i) 14 { 15 if(!mind[i])prime[++tot]=i,mind[i]=i,phi[i]=i-1; 16 for(j=1;j<=tot&&(tmp=i*prime[j])<=N;++j) 17 { 18 mind[tmp]=prime[j]; 19 if(i%prime[j]==0){phi[tmp]=phi[i]*prime[j];break;} 20 phi[tmp]=phi[i]*phi[prime[j]]; 21 } 22 } 23 for(i=1;i<=N;++i)sum[i]=(sum[i-1]+phi[i])%mod; 24 } 25 inline int quick_mod(int di,int mi) 26 { 27 int ans=1; 28 for(di%=mod;mi;mi>>=1,di=(LL)di*di%mod) 29 if(mi&1)ans=(LL)ans*di%mod; 30 return ans; 31 } 32 inline int divide(int x) 33 { 34 register int ret=1,i; 35 for(i=1;(LL)prime[i]*prime[i]<=x;++i) 36 if(x%prime[i]==0) 37 { 38 ret*=prime[i]; 39 while(x%prime[i]==0)x/=prime[i]; 40 } 41 return ret*x; 42 } 43 #define N1 100010 44 map<int,int>mmp; 45 int mem[N1]; 46 inline int getphi(int x) 47 { 48 if(x<=N)return sum[x]; 49 else if(mmp.count(x))return mmp[x]; 50 int ret=((LL)x*(x+1)/2%mod); 51 register int i,last; 52 for(i=2;i<=x;i=last+1) 53 last=x/(x/i),ret=( ret+mod-( (LL)(last-i+1)*getphi(x/i)%mod ) )%mod; 54 return mmp[x]=ret; 55 } 56 inline int calc(int a,int b) 57 { 58 if(!b)return 0; 59 if(b==1)return phi[a]; 60 if(a==1)return getphi(b); 61 return ((LL)calc(a/mind[a],b)*(mind[a]-1)%mod+calc(a,b/mind[a]))%mod; 62 } 63 int main() 64 { 65 // freopen("Ark.in","r",stdin); 66 register int n,m,ans=0,tmp,x,i; 67 scanf("%d%d",&n,&m),intn(); 68 memset(mem,-1,sizeof(mem)); 69 for(i=1;i<=n;++i) 70 { 71 x=divide(i),tmp=i/x; 72 if(mem[x]==-1)mem[x]=calc(x,m); 73 ans=(ans+(LL)tmp*mem[x]%mod)%mod; 74 } 75 printf("%d\n",ans); 76 }
首先你发现$\sum \mu(i^{2})$是出来搞笑的，输出1即可
然后我们考虑一下这个$\sum _{i=1} ^{n} \varphi(i^{2})$
由于$\varphi(i)$是积性函数，所以我们考虑
$\varphi(i)=\prod \varphi(p_{i})*p_{i} ^{a_{i}-1}$
$\varphi(i^{2})=\prod \varphi(p_{i})*p_{i} ^{a_{i}-1+a_{i}}$
$\varphi(i^{2})=\prod \varphi(p_{i})*p_{i} ^{a_{i}-1} *p_{i}^{a_{i}}$
所以……$\varphi(i^{2})=i*\varphi(i)$
设$f(i)=\varphi(i^{2})$
那么我们不是知道$\sum _{d|n} \varphi(d) = n$ 嘛
那么转化一下 $\sum _{d|n} d*\varphi(d) * \frac {n}{d} = n * \sum _{d|n} \varphi(d) =n^{2}$
所以我们有积性函数$f(i)=i\varphi(i),g(i)=i,h(i)=i^{2}$
这样就可以杜教筛了，设大写字母为前缀和，过程略……
最后结果是$F(i)=H(i)-\sum _{j=2} ^{m} g(j)*F(\left \lfloor \frac{m}{j} \right \rfloor )$
其中，$G(i)=\frac {n*(n+1)} {2} , H(i)= \frac {n*(n+1)*(2n+1)} {6} $
代码：
1 #include <cstdio> 2 #include <map> 3 #include <cstring> 4 using namespace std; 5 #define N 1000000 6 #define mod 1000000007 7 int n,inv6,vis[N+10],prime[N/10],tot,phi[N+10],sum[N+10]; 8 map<int,int> mmp; 9 #define LL long long 10 inline void intn() 11 { 12 register int i,j,tmp; 13 for(phi[1]=1,i=2;i<=N;++i) 14 { 15 if(!vis[i])prime[++tot]=i,phi[i]=i-1; 16 for(j=1;j<=tot&&(tmp=i*prime[j])<=N;++j) 17 { 18 if(i%prime[j]==0) 19 {vis[tmp]=1,phi[tmp]=phi[i]*prime[j];break;} 20 vis[tmp]=1,phi[tmp]=phi[i]*phi[prime[j]]; 21 } 22 } 23 for(i=1;i<=N;++i) 24 sum[i]=(sum[i-1]+(LL)i*phi[i]%mod)%mod; 25 } 26 inline int quick_mod(int di,int mi) 27 { 28 int ans=1; 29 for(;mi;mi>>=1,di=(LL)di*di%mod) 30 if(mi&1)ans=(LL)ans*di%mod; 31 return ans; 32 } 33 inline int getf(int x) 34 { 35 if(x<=N)return sum[x]; 36 if(mmp.count(x))return mmp[x]; 37 register int i,last,ret=(LL)x*(x+1)%mod*(2*x+1)%mod*inv6%mod; 38 for(i=2;i<=x;i=last+1) 39 last=x/(x/i), 40 ret=( ret +mod -(LL)getf(x/i)*( (LL)(last-i+1)*(last+i)/2 %mod )%mod )%mod; 41 return mmp[x]=ret; 42 } 43 int main() 44 { 45 // freopen("Ark.in","r",stdin); 46 inv6=quick_mod(6,mod-2); 47 scanf("%d",&n),intn(); 48 printf("1\n%d\n",getf(n)); 49 }
我们先反演一波……
定义$f(i)$为$gcd(a_{1},a_{2},......,a_{n})==i$的方案数,$F(i)$为$i|gcd(a_{1},a_{2},......,a_{n})$的方案数
那么有$F(i)=\sum _{i|n} f(n)$
然后我们演他一下，可以得到：
$f(n)=\sum _{n|d} \mu(\frac{d}{n})*F(d)$
我们考虑一下……$F(d)=( \left \lfloor \frac{r}{d} \right \rfloor - \left \lfloor \frac{l-1}{d} \right \rfloor )^{n}$
然后我们考虑枚举$i=\frac{d}{n}$
$f(n)=\sum_{i=1}^{ \left \lfloor \frac{r}{n} \right \rfloor } \mu(i) * ( \left \lfloor \frac{r}{i*n} \right \rfloor - \left \lfloor \frac{l-1}{i*n} \right \rfloor )^{n}$
然后我们对$\left \lfloor \frac{r}{i*n} \right \rfloor$和$\left \lfloor \frac{l-1}{i*n} \right \rfloor$除法分块
再杜教筛一个$\mu$的前缀和就行了。
代码：
1 #include <cstdio> 2 #include <map> 3 #include <cstring> 4 using namespace std; 5 #define mod 1000000007 6 #define LL long long 7 #define inf 0x3f3f3f3f 8 #define N 1000000 9 int prime[N],tot,mu[N+10],vis[N+10]; 10 int n,d,l,r; 11 inline void intn() 12 { 13 register int i,j,tmp; 14 for(mu[1]=1,i=2;i<=N;++i) 15 { 16 if(!vis[i])prime[++tot]=i,mu[i]=-1; 17 for(j=1;j<=tot&&(tmp=i*prime[j])<=N;++j) 18 { 19 vis[tmp]=1; 20 if(i%prime[j]==0)break; 21 mu[tmp]=-mu[i]; 22 } 23 } 24 for(i=2;i<=N;++i)mu[i]+=mu[i-1]; 25 } 26 inline int min(int a,int b){return a<b?a:b;} 27 map<int,int> mmp; 28 inline LL getmu(int x) 29 { 30 if(x<=N)return mu[x]; 31 if(mmp.count(x))return mmp[x]; 32 LL ret=1; 33 for(register int i=2,last;i<=x;i=last+1) 34 last=x/(x/i),ret=( ret-getmu(x/i)*(LL)(last-i+1)%mod )%mod; 35 return mmp[x]=ret; 36 } 37 inline LL quick_mod(LL di,int mi) 38 { 39 LL ans=1; 40 for(;mi;mi>>=1,di=di*di%mod) 41 if(mi&1)ans=ans*di%mod; 42 return ans; 43 } 44 signed main() 45 { 46 register int i,last; 47 scanf("%d%d%d%d",&n,&d,&l,&r); 48 r=r/d,l=(l-1)/d,intn(); 49 LL ans=0; 50 for(i=1;i<=r;i=last+1) 51 { 52 last= (l/i) ? min( r/(r/i), l/(l/i) ) : r/(r/i) , 53 ans=(ans+ ( getmu(last)-getmu(i-1) )*quick_mod((r/i-l/i),n) )%mod; 54 } 55 printf("%lld\n",(ans%mod+mod)%mod); 56 }
…………不得不说是个好题
被拿去做考试题了
结果只打了24分，而不是常见的84分暴力……原因待会再说
介绍一下考试历程
首先我观察了一下十进制下都以谁为分母是纯循环的……
发现2，5，10都不是
然后发现1无论何时是纯循环的
然后发现这是个完全积性函数！！！哇太棒了！
于是我就把这个函数命名为$f(i)$，当$\frac{1}{i}$纯循环时$f(i)$为1，否则为0
然后我就开始反演……
$\sum _{i=1} ^{n} \sum _{j=1} ^{m} f(j)[(i,j)==1]$
$\sum _{i=1} ^{n} \sum _{j=1} ^{m} f(j) \sum _{d|(i,j)}\mu(d)$
然后我们枚举$d$得到
$\sum _{d=1} ^{n} \mu(d)\left \lfloor \frac{n}{d} \right \rfloor \sum _{j=1}^{ \left \lfloor \frac{m}{d} \right \rfloor }f(j*d)$
由于$f(i)$是完全积性函数,所以$f(j*d)=f(j)*f(d)$
所以可以写成$\sum _{d=1} ^{n} \mu(d)f(d)\left \lfloor \frac{n}{d} \right \rfloor \sum _{j=1}^{ \left \lfloor \frac{m}{d} \right \rfloor }f(j)$
然后我成功的yy出了怎么线性筛$f(i)$的值
......然后不会处理$\sum _{j=1}^{ \left \lfloor \frac{m}{d} \right \rfloor }f(j)$这一部分了……我以为可以杜教筛这部分，结果发现不行
最后在空间允许的范围内打了时间复杂度为$O(n)$的24分算法
考完试和$wq$交流之后我发现......
$f(i)\Leftrightarrow[(i,k)==1]$
……当时整个人傻了
后来的后来发现自己式子的某个主流题解是一样的……
更加心里苦
我们还是来证明一下上面那个式子为啥成立吧……
设循环节长度为$L$
那么根据题给的循环节生成方式，应该有$x \% y=x*k^{L} \% y$
由于题面要求被计算的$x$和$y$互质，所以我们可以消去$x$
从而有$k^{L} \equiv 1(mod y)$
这就意味着$k$在模$y$的意义下有逆元
如果一个数在模意义下有逆元……我们回想一下$exgcd$的过程
那么就会有$(k,y)==1$
这就好说了……真是心累
然后还有一点我考试没有想到的，就是快速计算$F(n)=\sum _{i=1} ^{n} [(i,k)==1]$
由于我们知道$(a,b)==(a-k*b,b)$
所以$\sum _{i=1} ^{k} [(i,k)==1]==\sum _{i=k+1} ^{2*k} [(i,k)==1]$
由此我们可以归纳出
$F(n)=\left \lfloor \frac{n}{k} \right \rfloor F(k)+F(n \% k)$
我们回到刚才我演的那个式子
$\sum _{d=1} ^{n} \mu(d) [(k,d)==1] \left \lfloor \frac{n}{d} \right \rfloor F( \left \lfloor \frac{m}{d} \right \rfloor ) $
这样，我们预处理1～k内的$F$值，就可以$O(n)$回答询问了，可以得到84分
然而事实证明剩下的16分分值和难度不成正比……
我们考虑上面这个东西，$\left \lfloor \frac{n}{d} \right \rfloor$和$\left \lfloor \frac{m}{d} \right \rfloor$是可以除法分块的
那么我们如果能处理出来一个$\sum _{i=1} ^{n} \mu(i) f(i)$前缀和就能$O(\sqrt{n})$回答了……
这里面的f(i)是我上面定义的$f(i)\Leftrightarrow[(i,k)==1]$
我们考虑怎么处理这个前缀和……
设其为$G(n,k)$,考虑化简为子问题递归解决
对于$k$的某一个质因子$p$，我们可以把$k$写成$k=p^{a}q((p,q)==1)$的形式
那么如果$(i,k)==1$,则$(i,q)==1$且$(i,p)==1$
如果我们处理出$G(n,q)$，再减去与$q$互质但是与$p$不互质的数的个数，就可以得到$G(n,k)$了
我们假设某个与$q$互质但是与$p$不互质的数$i$可以表示为$i=p^{c}y((y,q)==1,(y,p)==1)$
由于与$p$不互质，所以$c>0$又因为$c>=2$时$\mu(i)==0$所以$c==1$
因此我们考虑的是一系列形如$i=py$的$i$
我们写一下这个式子……
$G(n,k)=G(n,q)-\sum _{y=1} ^{\left \lfloor \frac{n}{p} \right \rfloor } \mu(py) [(y,q)==1] [(y,p)==1]$
因为$(p,y)==1,\mu(py)=\mu(p)\mu(y)$
因此就会有$G(n,k)=G(n,q)-\sum _{y=1} ^{\left \lfloor \frac{n}{p} \right \rfloor } \mu(p)\mu(y) [(y,p)==1][(y,q)==1]$
我们把$\mu(p)=-1$带入，有$G(n,k)=G(n,q)+\sum _{y=1} ^{\left \lfloor \frac{n}{p} \right \rfloor } \mu(y) [(y,p)==1][(y,q)==1]$
又由于我们一开始利用的条件，“如果$(i,k)==1$,则$(i,q)==1$且$(i,p)==1$”
所以$[(y,p)==1][(y,q)==1]\Leftrightarrow[(y,k)==1]$
所以原式可以写成$G(n,k)=G(n,q)+\sum _{y=1} ^{\left \lfloor \frac{n}{p} \right \rfloor } \mu(y) [(y,k)==1]$
即$G(n,k)=G(n,q)+G(\left \lfloor \frac{n}{p} \right \rfloor ,k)$
这样我们就可以记忆化搜索我们的$G$值，边界条件是$n==0$时$G(0,k)=0$,$n==1$时$G(1,k)=1$；
$k==1$时$G(n,1)=\sum _{i=1} ^{n} \mu(i)$，杜教筛之。
这个部分和上面的bzoj3309一样，都是对递归的巧妙应用
这种定义一个函数求递推式来帮助思考的方法很有效！
这样这题就被解决了！的确是优美的题目啊！
不知道为什么跑的贼慢的代码：
1 #include <map> 2 #include <cstdio> 3 #include <cstring> 4 using namespace std; 5 #define N 20000000 6 #define K 2010 7 #define LL long long 8 int n,m,k; 9 int prime[N/10+10],tot,mu[N+10],sum[K]; 10 bool vis[N+10]; 11 inline int min(int a,int b){return a<b?a:b;} 12 inline int gcd(int a,int b){return b==0?a:gcd(b,a%b);} 13 inline void init() 14 { 15 register int i,j,tmp; 16 for(i=1;i<=k;++i)sum[i]=sum[i-1]+(gcd(i,k)==1); 17 for(mu[1]=1,i=2;i<=N;++i) 18 { 19 if(!vis[i])prime[++tot]=i,mu[i]=-1; 20 for(j=1;j<=tot&&(tmp=i*prime[j])<=N;++j) 21 { 22 vis[tmp]=1; 23 if(i%prime[j]==0){mu[tmp]=0;break;} 24 mu[tmp]=-mu[i]; 25 } 26 } 27 for(i=1;i<=N;++i)mu[i]+=mu[i-1]; 28 } 29 int bin[100],ge; 30 inline void divide(int kkk) 31 { 32 register int i; 33 for(i=1;prime[i]*prime[i]<=kkk;++i) 34 if(kkk%prime[i]==0) 35 { 36 bin[++ge]=prime[i]; 37 while(kkk%prime[i]==0)kkk/=prime[i]; 38 } 39 if(kkk>1)bin[++ge]=kkk; 40 } 41 inline int getF(int x) 42 { 43 if(x<=k)return sum[x]; 44 return (x/k)*sum[k]+sum[x%k]; 45 } 46 map<int,int>smu; 47 inline int getmu(int x) 48 { 49 if(x<=N)return mu[x]; 50 if(smu.count(x))return smu[x]; 51 register int i,last,ret=1; 52 for(i=2;i<=x;i=last+1) 53 last=x/(x/i),ret-=(LL)getmu(x/i)*(last-i+1); 54 return smu[x]=ret; 55 } 56 map<int,int>mmp[7]; 57 inline int getsum(int x,int layer) 58 { 59 if(x<=1)return x; 60 if(layer==0)return getmu(x); 61 if(mmp[layer].count(x))return mmp[layer][x]; 62 return mmp[layer][x]=getsum(x,layer-1)+getsum(x/bin[layer],layer); 63 } 64 int main() 65 { 66 scanf("%d%d%d",&n,&m,&k); 67 init(),divide(k); 68 register int i,to=min(n,m),last; 69 LL ans=0; 70 // if(n>m)n^=m,m^=n,n^=m; 71 for(i=1;i<=to;i=last+1) 72 last=min(n/(n/i),m/(m/i)), 73 ans+=(getsum(last,ge)-getsum(i-1,ge))*(LL)getF(m/i)*(n/i); 74 printf("%lld\n",ans); 75 }
大概我最近做的题就这些……虽然还有另外一道模拟赛的题，式子还不错，先不放出来了……
杜教筛……怎么说呢，算是对狄利克雷卷积的灵活应用吧
作为优化的工具，不管是复杂度还是思想都很优秀
逐渐有一点点数学题的手感了知道上来无脑反演了，加油咯……
