《積性函數求和的幾種方法》這篇paper大概就是講了杜教篩和任之州一種神奇的自創做法。%%%IOI爺
分別復雜度是O(n^(2/3))和O(n^(3/4)/logn)的。
在一般情況下，后者的常數和復雜度都更加優秀。
這篇就先講杜教篩好了
①杜教篩
運用Dircichlet卷積來完成復雜度的化簡。
可以參考唐教的介紹 http://blog.csdn.net/skywalkert/article/details/50500009
還是上幾個例題。
A. 求n個正整數的約數之和，其中n<=10^12。
運用莫比烏斯反演的那個技巧即可做到O(sqrt(n))的復雜度。
其實這跟杜教篩並沒有什么關系…
#include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; long long n; int main() { cin>>n; long long ed=1,ans=0; for(long long a=1;a<=n;a=ed+1) { ed=(n/(n/a)); ans+=(a+ed)*(ed-a+1)/2*(n/a); } cout<<ans; }
B. 求前n個正整數的歐拉函數之和，其中n<=10^11。
因為
，所以
。
設
，則
那么我們只要算出
個
的值就可以計算出歐拉函數的前綴和。
我們運用主定理進行分析：
根據唐教說只要展開一層就行了.
那么
最后一步可以直接根據積分得到(symbolab)。
如果我們用篩法處理前k個前綴和，那么復雜度就變成了O($\frac{n}{\sqrt{k}}$)，當k=$n^\frac{2}{3}$時可以取到極值。
為什么我們會這么考慮呢？
因為Dircichlet卷積：
例如我們知道莫比烏斯反演
則
也就是說
。
設$id(n)=n$，我們知道$\varphi(i)*I=id$。
所以
。
那總結一下如果Dircichlet卷積的左邊其中一個函數是待求前綴和的，而且卷積的另外兩個函數比較好計算前綴和，那么就可以簡化計算的過程。
51nod 1239
實現的時候要注意一些細節
由於常數問題比較嚴重，還是建議用map稍微記憶化意思意思
upd on 2017-02-01：
這篇文章有一些年頭了...這里說一下這個map其實是可以不用的...
注意到每個參數都是n的約數，有一個儲存上的trick，大致如下：
![KVR$M]RNG`~DNWXE@YXPALS](/image/aHR0cHM6Ly9pbWFnZXMyMDE1LmNuYmxvZ3MuY29tL2Jsb2cvOTAzNjUyLzIwMTcwMi85MDM2NTItMjAxNzAyMDExNDAyMjMzMzktMjE0NTM2NzAyLnBuZw==.png "KVR$M]RNG`~DNWXE@YXPALS")
其中s大約為$\sqrt{n}$。這樣儲存就少了一個log了（洲閣篩也需要這個trick，但是實際運行速度好像相差不大）
#include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; typedef long long ll; ll MOD=1000000007; #define MM 5000000 int ps[MM+5],pn=0,mu[MM+5],eul[MM+5]; bool np[MM+5]; ll qzheul[MM+5]; void shai() { np[1]=mu[1]=eul[1]=1; for(int i=2;i<=MM;i++) { if(!np[i]) {mu[i]=-1; ps[++pn]=i; eul[i]=i-1;} for(int j=1;j<=pn&&i*ps[j]<=MM;j++) { np[i*ps[j]]=1; if(i%ps[j]) { mu[i*ps[j]]=-mu[i]; eul[i*ps[j]]=eul[i]*(ps[j]-1); } else { mu[i*ps[j]]=0; eul[i*ps[j]]=eul[i]*ps[j]; break; } } } for(int i=1;i<=MM;i++) qzheul[i]=(qzheul[i-1]+eul[i])%MOD; } ll qp(ll a,ll b) { if(b==0) return 1; ll hf=qp(a,b>>1); hf=hf*hf%MOD; if(b&1) hf=hf*(a%MOD)%MOD; return hf; } ll rv2=qp(2,MOD-2),n; #define G 233333 ll p1[G],p2[G]; ll &gm(ll x) { if(x<G) return p1[x]; return p2[n/x]; } ll eulsum(ll x) { if(x<=MM) return qzheul[x]; ll &ps=gm(x); if(ps!=-1) return ps; ll ans,lst; if(x%MOD!=0) ans=(x%MOD)*((x%MOD)+1)%MOD*rv2%MOD; else ans=0; for(ll p=2;p<=x;p=lst+1) { lst=x/(x/p); ans-=((lst-p+1)%MOD)*eulsum(x/p)%MOD; ans%=MOD; } ans=(ans%MOD+MOD)%MOD; return ps=ans; } int main() { memset(p1,-1,sizeof(p1)); memset(p2,-1,sizeof(p2)); shai(); cin>>n; cout<<eulsum(n)<<"\n"; }
C. 求前n個正整數的mu函數（莫比烏斯函數）之和，n<=10^11。
咦好像$\sum_{d|n}\mu(d)=[n=1]$ 
設單位函數
也就是說$\mu*1=dw$。
我們來手推一發
$1={\sum_{i=1}^n\sum_{d|i}\mu(d)}={\sum_{d=1}^n\mu(d)*\lfloor\frac{n}{d}\rfloor}={\sum_{d=1}^n\mu(d)*\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}1}=\sum_{i=1}^{n}\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}\mu(d)$
（latex打得好辛苦啊QAQ
所以也可以用同樣的方法搞搞
51nod 1244
#include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; typedef long long ll; ll MOD=1000000007; #define MM 6000000 int ps[MM+5],pn=0,mu[MM+5]; bool np[MM+5]; ll qzhmu[MM+5]; void shai() { np[1]=mu[1]=1; for(int i=2;i<=MM;i++) { if(!np[i]) {mu[i]=-1; ps[++pn]=i;} for(int j=1;j<=pn&&i*ps[j]<=MM;j++) { np[i*ps[j]]=1; if(i%ps[j]) mu[i*ps[j]]=-mu[i]; else{mu[i*ps[j]]=0; break;} } } for(int i=1;i<=MM;i++) qzhmu[i]=qzhmu[i-1]+mu[i]; } #define G 233333 ll p1[G],p2[G],n; ll &gm(ll x) { if(x<G) return p1[x]; return p2[n/x]; } ll musum_(ll x) { if(x<=MM) return qzhmu[x]; ll &ps=gm(x); if(ps!=p2[0]) return ps; ll ans=1,lst; for(ll p=2;p<=x;p=lst+1) { lst=x/(x/p); ans-=(lst-p+1)*musum_(x/p); } return ps=ans; } ll musum(ll x) { memset(p1,-2333,sizeof(p1)); memset(p2,-2333,sizeof(p2)); n=x; return musum_(x); } int main() { shai(); ll a,b; cin>>a>>b; cout<<musum(b)-musum(a-1)<<"\n"; }
還有唐教博客上的一些例題也順帶做了一下
bzoj3944 雙倍經驗，略
hdu5608
還是杜教篩搞搞。
設g(n)=n^2-3n+2
容易得到f*1=g。
則$\sum_{i=1}^ng(i)=\sum_{i=1}^n\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}f(d)$。（推導和上面一樣）
而$\sum_{i=1}^ng(i)=\frac{n(n+1)(2n+1)}{6}-\frac{3n(n+1)}{2}+2n=\frac{n(n-1)(n-2)}{3}$。
初始化的時候我們可以先線篩，然后莫比烏斯反演一下暴力更新。
$f(n)=\sum_{d|n}{\mu(\frac{n}{d})g(d)}$
對於每一個g暴力更新它的倍數，由調和級數顯然是O(plogp)的（假設你預處理到p
這種題出到bc div2簡直了
#include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; typedef long long ll; ll MOD=1000000007; #define MM 600000 int ps[MM/5],pn=0,mu[MM+5]; bool np[MM+5]; ll smf[MM+5],fqzh[MM+5]; void shai() { np[1]=mu[1]=1; for(int i=2;i<=MM;i++) { if(!np[i]) {mu[i]=-1; ps[++pn]=i;} for(int j=1;j<=pn&&i*ps[j]<=MM;j++) { np[i*ps[j]]=1; if(i%ps[j]) mu[i*ps[j]]=-mu[i]; else {mu[i*ps[j]]=0; break;} } } for(int d=1;d<=MM;d++) { ll g=((ll)d*d%MOD-3*d%MOD+2)%MOD; for(int n=d;n<=MM;n+=d) smf[n]=(smf[n]+mu[n/d]*g%MOD)%MOD; } for(int i=1;i<=MM;i++) fqzh[i]=(fqzh[i-1]+smf[i])%MOD; } ll qp(ll a,ll b) { if(b==0) return 1; ll hf=qp(a,b>>1); hf=hf*hf%MOD; if(b&1) hf=hf*(a%MOD)%MOD; return hf; } ll rv3=qp(3,MOD-2); #define G 23333 ll p1[G],p2[G],n; ll &gm(ll x) { if(x<G) return p1[x]; return p2[n/x]; } ll gf(ll x) { if(x<=MM) return fqzh[x]; ll &ps=gm(x); if(ps!=p2[0]) return ps; ll lst,ans=x*(x-1)%MOD*(x-2)%MOD*rv3%MOD; for(ll p=2;p<=x;p=lst+1) { lst=x/(x/p); ans=ans-(lst-p+1)*gf(x/p)%MOD; ans%=MOD; } ans=(ans%MOD+MOD)%MOD; return ps=ans; } ll fs(ll x) { memset(p1,-2333,sizeof(p1)); memset(p2,-2333,sizeof(p2)); n=x; return gf(x); } int main() { shai(); int T,n; scanf("%d",&T); while(T--) { scanf("%d",&n); ll ans=fs(n); printf("%d\n",int((ans%MOD+MOD)%MOD)); } }
2017-02-01：所有代碼已經去掉了map= =