杜教篩小結

時間 2019-12-14

標籤小結简体版

原文原文鏈接

首先的首先，%%%%%%杜教php

……咳咳，言歸正傳，如今來總結一下杜教篩……html

杜教篩大概是一種用於快速處理積性函數前綴和的有力工具，主要利用了狄利克雷卷積的知識算法

咱們用高中均值不等式的複雜度分析能夠獲得其複雜度是$O(n^{\frac{2}{3}})$的ide

具體的講解我就不寫了……留幾篇比較好的博客。函數

http://www.javashuo.com/article/p-gslvwtiw-dw.html工具

這裏總結一下我作過的題目吧……優化

1.bzoj3944ui

因爲太水以致於沒打……spa

能夠當作板子練習題的題目。.net

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 }