杜教筛进阶+洲阁筛讲解+SPOJ divcnt3

Part 1:杜教筛进阶
在了解了杜教筛基本应用，如$sum_{i=1}^nvarphi(i)$的求法后，我们看一些杜教筛较难的应用。
求$sum_{i=1}^nvarphi(i)*i$
考虑把它与$id$函数狄利克雷卷积后的前缀和。
$$sum_{i=1}^nsum_{d|i}varphi(d)*d*frac id=sum_{i=1}^ni^2$$枚举$T=frac id$，原式化为
$$sum_{T=1}^nTsum_{d=1}^{lfloorfrac nT floor}varphi(d)*d=sum_{i=1}^ni^2$$移项，得
$$sum_{i=1}^nvarphi(i)*i=sum_{i=1}^ni^2-sum_{T=2}^nTsum_{d=1}^{lfloorfrac nT floor}varphi(d)*d$$右边的$sum_{d=1}^{lfloorfrac nT floor}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=316241,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
",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
",a1+f2[1]);
    }
    return 0;
}