BZOJ 3944: Sum(杜教筛)

传送门

解题思路

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

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

　　把第一项提出来：

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

　　移项

\[g(1)S(n)=\sum\limits_{i=1}^nh(i)-\sum\limits_{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('\n');
    }
    return 0;   
}