bzoj 3944: Sum 杜教筛

       本来以为这种东西只能O(N)线性筛，但是大千世界，无(sang)奇(xin)不(bing)有(kuang)，确实存在更快的算法。

       省选的时候rzz讲这种东西在国内OI称为杜教筛，用来求数论函数的前缀和，课件中有一般形式。

       首先考虑对μ函数的前缀和（欧拉函数同理）：

       令f(x)=Σ(d|x) μ(d)，μ的前缀和记为s(x)，令g(x)为f的前缀和，那么有

       g(x)=Σ(i=1,x) Σ(d|i) μ(d)=Σ(d=1,x) μ(d)[x/d]=Σ(i=1,x) s(x/i)，显然有g(x)=1，那么可以得到：

       s(x)=1-Σ(i=2,x) s(x/i)，然后直接暴力hash+记忆化可以做到O(N^(3/4))，预处理前N^(2/3)就可以做到O(N^(2/3))，具体好像要用积分算我不会QAQ。。。注意后面的O(N^(2/3))不需要hash，可以发现每一项一定是x/i，而i<N^(1/3)，那么可以存储在除x之后的结果的下标处。

       小心N<2^31爆int。

AC代码如下：

#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;

int cas,n,m,cnt,c[1000005]; ll phi[2000005],mu[2000005],p[100005],q[100005]; bool vis[100005];
ll get_p(int x){
	return (x<=m)?phi[x]:p[n/x];
}
ll get_q(int x){
	return (x<=m)?mu[x]:q[n/x];
}
void solve(int x){
	if (x<=m) return; int i,j=1,t=n/x;
	if (vis[t]) return; vis[t]=1;
	p[t]=((ll)x+1)*x>>1; q[t]=1;
	while (j<x){
		i=j+1; j=x/(x/i); solve(x/i);
		p[t]-=get_p(x/i)*(j-i+1); q[t]-=get_q(x/i)*(j-i+1);
	}
}
int main(){
	scanf("%d",&cas); m=2000000;
	int i,j; phi[1]=mu[1]=1;
	for (i=2; i<=m; i++){
		if (!phi[i]){
			phi[i]=i-1; mu[i]=-1; c[++cnt]=i;
		}
		for (j=1; j<=cnt && i*c[j]<=m; j++)
			if (i%c[j]){
				phi[i*c[j]]=phi[i]*(c[j]-1); mu[i*c[j]]=-mu[i];
			} else{
				phi[i*c[j]]=phi[i]*c[j]; mu[i*c[j]]=0; break;
			}
	}
	for (i=2; i<=m; i++){ phi[i]+=phi[i-1]; mu[i]+=mu[i-1]; }
	while (cas--){
		scanf("%d",&n); memset(vis,0,sizeof(vis));
		if (n<=m) printf("%lld %lld\n",phi[n],mu[n]); else{
			solve((ll)n); printf("%lld %lld\n",p[1],q[1]);
		}
	}
}

by lych

2016.3.24