杜教筛小结

首先的首先,%%%%%%杜教

……咳咳,言归正传,现在来总结一下杜教筛……

杜教筛大概是一种用于快速处理积性函数前缀和的有力工具,主要利用了狄利克雷卷积的知识

我们用高中均值不等式的复杂度分析可以得到其复杂度是$O(n^{\frac{2}{3}})$的

具体的讲解我就不写了……留几篇比较好的博客。

https://www.cnblogs.com/abclzr/p/6242020.html

这里总结一下我做过的题目吧……

1.bzoj3944

由于太水以至于没打……

可以当做板子练习题的题目。

2.bzoj4176

深刻的记住了一个式子……之前一直没有这种意识,有一些转换的确很显然但是就是没想过

\[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 }
bzoj4176

3.bzoj3512

然后这题我们发现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 }
bzoj3512

4.bzoj4916

 首先你发现$\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 }
bzoj4916

5.bzoj3930

我们先反演一波……

定义$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 }
bzoj3930

6.bzoj4652

…………不得不说是个好题

被拿去做考试题了

结果只打了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 }
bzoj4652

大概我最近做的题就这些……虽然还有另外一道模拟赛的题,式子还不错,先不放出来了……

杜教筛……怎么说呢,算是对狄利克雷卷积的灵活应用吧

作为优化的工具,不管是复杂度还是思想都很优秀

逐渐有一点点数学题的手感了知道上来无脑反演了,加油咯……

posted @ 2018-01-22 19:01 LadyLex 阅读(...) 评论(...) 编辑 收藏