BZOJ 3944: Sum(杜教筛)

传送门

解题思路

　　机房里的神仙们一万年前就会的东西拿出来学一学。杜教筛可以在(O(n^{2/3}))的时间内积性函数前缀和，做法如下。
　　首先设要求的是(sumlimits_{i=1}^n f(i))。设(h=f*g)，(S(x)=sumlimits_{i=1}^nf(i))，那么可以得出$$sumlimits_{i=1}^{nh(i)=sumlimits_{i=1}}nsumlimits_{d|n}g(d)f(frac{n}{d})$$
　　换一下枚举顺序：
　　$$sumlimits_{i=1}^{nh(i)=sumlimits_{d=1}}ng(d)sumlimits_{i=1}^{n/d}f(i)$$

[sumlimits_{i=1}^nh(i)=sumlimits_{d=1}^ng(d)S(n/d) ]

　　把第一项提出来：

[sumlimits_{i=1}^nh(i)=g(1)S(n)+sumlimits_{d=2}^ng(d)S(n/d) ]

　　移项

[g(1)S(n)=sumlimits_{i=1}^nh(i)-sumlimits_{d=2}^ng(d)S(n/d) ]

　　发现要求的东西已经被放在左边了，只要找出一个好求前缀和的函数(h)，后面部分可以数论分块解决，这个可以递归实现。这就是杜教筛的思路。现在要解决的问题就是如何找(f,g)，然而坑的是这个玩意并没有什么通用的，似乎只能做题时自行(yy)。
　　这道题要求(mu)和(phi)，而(mu*I=epsilon)，(epsilon)是元函数，当且仅当(i=1)时(epsilon(i)=1)，否则为(0)，而(I)是恒等函数，(I(i)=1)。而(phi*id=epsilon)，(id)是单位函数，(id(i)=i)。这样这道题就可以做了。在下人丑常数大还懒，所以这道题严重卡常，用(map)存平添(log)，交了一页终于用循环展开大法过了。。

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<map>
 
using namespace std;
const int N=3000005;
typedef long long LL;

template<class T> void rd(T &x){
	x=0; char ch=getchar();
	while(!isdigit(ch)) ch=getchar();
	while(isdigit(ch)) x=(x<<1)+(x<<3)+ch-'0',ch=getchar();	
}
template<class T> void out(T x){
	if(!x) return; out(x/10); putchar('0'+x%10);
}
 
int T,prime[N],miu[N],lim=3000000,cnt,sum2[N];
LL sum1[N],n,phi[N],ans1,ans2;
bool vis[N];
map<LL,LL> mp1;
map<int,int> mp2;
 
struct Data{
	LL tmp1,tmp2;	
	Data(LL A=0,LL B=0){
		tmp1=A; tmp2=B;	
	}
};
 
void prework(){
    miu[1]=1; phi[1]=1;
    for(register int i=2;i<=lim;++i){
        if(!vis[i]) prime[++cnt]=i,miu[i]=-1,phi[i]=i-1;
        for(register int j=1;j<=cnt && prime[j]*i<=lim;++j){
            vis[prime[j]*i]=1;
            if(!(i%prime[j])) {
                vis[i*prime[j]]=1;
                phi[i*prime[j]]=phi[i]*prime[j];    
                break;
            }
            miu[i*prime[j]]=-miu[i];
            phi[i*prime[j]]=phi[i]*phi[prime[j]];
        }
    }
    register int i=1;
    for(;i<=lim;i+=16){
    	sum1[i]=sum1[i-1]+phi[i];
    	sum1[i+1]=sum1[i]+phi[i+1];
    	sum1[i+2]=sum1[i+1]+phi[i+2];
    	sum1[i+3]=sum1[i+2]+phi[i+3];
    	sum1[i+4]=sum1[i+3]+phi[i+4];
    	sum1[i+5]=sum1[i+4]+phi[i+5];
    	sum1[i+6]=sum1[i+5]+phi[i+6];
    	sum1[i+7]=sum1[i+6]+phi[i+7];
    	sum1[i+8]=sum1[i+7]+phi[i+8];
    	sum1[i+9]=sum1[i+8]+phi[i+9];
    	sum1[i+10]=sum1[i+9]+phi[i+10];
    	sum1[i+11]=sum1[i+10]+phi[i+11];
    	sum1[i+12]=sum1[i+11]+phi[i+12];
    	sum1[i+13]=sum1[i+12]+phi[i+13];
    	sum1[i+14]=sum1[i+13]+phi[i+14];
    	sum1[i+15]=sum1[i+14]+phi[i+15];
        sum2[i]=sum2[i-1]+miu[i];
        sum2[i+1]=sum2[i]+miu[i+1];
        sum2[i+2]=sum2[i+1]+miu[i+2];
        sum2[i+3]=sum2[i+2]+miu[i+3];
        sum2[i+4]=sum2[i+3]+miu[i+4];
        sum2[i+5]=sum2[i+4]+miu[i+5];
        sum2[i+6]=sum2[i+5]+miu[i+6];
        sum2[i+7]=sum2[i+6]+miu[i+7];
        sum2[i+8]=sum2[i+7]+miu[i+8];
        sum2[i+9]=sum2[i+8]+miu[i+9];
        sum2[i+10]=sum2[i+9]+miu[i+10];
        sum2[i+11]=sum2[i+10]+miu[i+11];
        sum2[i+12]=sum2[i+11]+miu[i+12];
        sum2[i+13]=sum2[i+12]+miu[i+13];
        sum2[i+14]=sum2[i+13]+miu[i+14];
        sum2[i+15]=sum2[i+14]+miu[i+15];
    }
    for(;i<=lim;++i) 
    	sum1[i]=sum1[i-1]+phi[i],
    	sum2[i]=sum2[i-1]+miu[i];
}
 
Data solve(LL x){
    if(x<=lim) {return Data(sum1[x],sum2[x]);}
    if(mp1[x]) {return Data(mp1[x],mp2[x]);}
    Data tmp,nxt; tmp.tmp1=(LL)x*(x+1)/2;
    tmp.tmp2=1;
    for(register LL l=2,r;l<=x;l=r+1){
        r=x/(x/l); nxt=solve(x/l); 
        tmp.tmp1-=(r-l+1)*nxt.tmp1;
        tmp.tmp2-=(r-l+1)*nxt.tmp2;
    }
    mp1[x]=tmp.tmp1; mp2[x]=tmp.tmp2;
    return tmp;
}
 
signed main(){
	rd(T); prework(); Data zz;
    while(T--){
    	rd(n); zz=solve(n);
    	ans1=zz.tmp1; ans2=zz.tmp2;
    	if(!ans1) putchar('0');
    	else out(ans1); putchar(' ');
    	if(!ans2) putchar('0');
    	else {
    		if(ans2<0) putchar('-'),ans2=-ans2;
			out(ans2);	
    	} putchar('
');
    }
    return 0;   
}