1. 程式人生 > >杜教篩進階+洲閣篩講解+SPOJ divcnt3

杜教篩進階+洲閣篩講解+SPOJ divcnt3

前綴 con bre sca else rac poj class algo

Part 1:杜教篩進階
在了解了杜教篩基本應用,如$\sum_{i=1}^n\varphi(i)$的求法後,我們看一些杜教篩較難的應用。
求$\sum_{i=1}^n\varphi(i)*i$
考慮把它與$id$函數狄利克雷卷積後的前綴和。
$$\sum_{i=1}^n\sum_{d|i}\varphi(d)*d*\frac id=\sum_{i=1}^ni^2$$枚舉$T=\frac id$,原式化為
$$\sum_{T=1}^nT\sum_{d=1}^{\lfloor\frac nT\rfloor}\varphi(d)*d=\sum_{i=1}^ni^2$$移項,得
$$\sum_{i=1}^n\varphi(i)*i=\sum_{i=1}^ni^2-\sum_{T=2}^nT\sum_{d=1}^{\lfloor\frac nT\rfloor}\varphi(d)*d$$右邊的$\sum_{d=1}^{\lfloor\frac nT\rfloor}\varphi(d)*d$遞歸求就行了。
總結:當遇到一些不好求前綴和的函數時,一般將其與一個易於求前綴和的函數進行狄利克雷卷積,得到另一個易於求前綴和的函數,然後通過簡單的數學變換,得到可以遞歸的式子。
Part 2:洲閣篩講解


有一篇博客講的挺好:
http://debug18.com/posts/calculate-the-sum-of-multiplicative-function
Part 3:SPOJ divcnt3

洲閣篩的簡單應用。

#include <cstdio>
#include <algorithm>
using namespace std;

typedef long long ll;
const int N=316228,p=1000003;
int T,e,tt,t2,pr[N],hd[p],nx[p],w[p],mx[N],ci[N],s[N],D[p];
ll n,a1,to[p],d[p],g[p],f[N],f2[p],sf[N];
void ins(int x,ll y) {int h=y%p; to[++e]=y,w[e]=x,nx[e]=hd[h],hd[h]=e;} int qr(ll x) {for(int i=hd[x%p];i;i=nx[i]) if(to[i]==x) return w[i]; return 0;} void sol() { e=t2=0; for(ll i=1;i<=n;i=n/(n/i)+1) hd[n/i%p]=0; for(ll i=1;i<=n;i=n/(n/i)+1) ins(++t2,n/i),d[t2]=g[t2]=n/i,D[t2]=0; for(int
i=1;i<=tt;i++) for(int j=1;j<=t2&&(ll)pr[i]*pr[i]<=d[j];j++) { int k=qr(d[j]/pr[i]); g[j]-=g[k]-(i-1-D[k]),D[j]=i; } } void sol2() { for(int i=1;i<=t2;i++) f2[t2]=1; for(int i=tt;i;i--) for(int j=1;j<=t2&&(ll)pr[i]*pr[i]<=d[j];j++) { if(pr[i+1]>d[j]) f2[j]=1; else if((ll)pr[i+1]*pr[i+1]>d[j]) f2[j]=(s[min(N-1LL,d[j])]-s[pr[i+1]-1])*4+1; for(ll pi=pr[i],c=1;d[j]>=pi;pi*=pr[i],c++) { ll k=d[j]/pi,k2; if(pr[i+1]>k) k2=1; else if((ll)pr[i+1]*pr[i+1]>k) k2=(s[min(N-1LL,k)]-s[pr[i+1]-1])*4+1; else k2=f2[qr(k)]; f2[j]+=k2*(3*c+1); } } } int main() { scanf("%d",&T),f[1]=sf[1]=1; for(int i=2;i<N;i++) { s[i]=s[i-1]+!f[i]; if(!f[i]) pr[++tt]=i,f[i]=4,mx[i]=i,ci[i]=1; for(int j=1,k;j<=tt&&(k=i*pr[j])<N;j++) { if(i%pr[j]) f[k]=f[i]*4,mx[k]=pr[j],ci[k]=1; else {f[k]=f[i/mx[i]]*(ci[i]*3+4),mx[k]=mx[i]*pr[j],ci[k]=ci[i]+1; break;} } sf[i]=sf[i-1]+f[i]; } pr[tt+1]=316241; while(T--) { scanf("%lld",&n); if(n<N) {printf("%lld\n",sf[n]); continue;} a1=0,sol(),sol2(); for(int i=1,r;i<N;i=r+1) { int j=qr(n/i); ll k; if(pr[tt+1]>n/i) k=1; else k=g[j]-(tt-D[j]); a1+=(sf[r=min(N-1LL,n/(n/i))]-sf[i-1])*(k-1)*4; } printf("%lld\n",a1+f2[1]); } return 0; }

杜教篩進階+洲閣篩講解+SPOJ divcnt3