杜教筛

一．杜教筛

这有一些资料
2016年国家候选队论文任之洲积性函数求和的几种方法
还有几篇博客
博客1
博客2

１．一般形式

设(f(n))为一个数论函数，求

[S(n)=sum_{i=1}^{n}f(i) ]

其中(nleq10^{10})

分析：
直接求很困难，可以采用的一个方法是找出另外一个函数(g(x))

[f*g=sum_{i=1}^{n}sum_{j|i}f(i)g(frac{i}{j})=sum_{ij leq n}f(i)g(j)=sum_{i=1}^{n}g(i)sum_{j=1}^{leftlfloor frac{n}{i} ight floor}f(j) ]

[=sum_{i=1}^{n}g(i)S(leftlfloor frac{n}{i} ight floor) ]

然后就可以得到一个关于(S(n))的递推式

[g(1)*S(n)=sum_{i=1}^{n}(f*g)(i)-sum_{i=2}^{n}g(i)S(leftlfloor frac{n}{i} ight floor) ]

如果的((f*g)(i))前缀和以及(g(i))的前缀和都很快计算，那么就可以很快地计算(S(n))，直接计算时间复杂度(O(n^{frac{3}{4}}))，线性筛预处理前(n^{frac{2}{3}})项之后，时间复杂度为(O(n^{frac{2}{3}}))

２．(mu)

求

[S(n)=sum_{i=1}^{n}mu(i) ]

分析：

设(g=1)因为 (mu*l=varepsilon)，所以，就可以得到

[S(n)=varepsilon-sum_{i=2}^{n}S(leftlfloor frac{n}{i} ight floor) ]

[S(n)=1-sum_{i=2}^{n}S(leftlfloor frac{n}{i} ight floor) ]

时间复杂度(O(n^{frac{2}{3}}))

３．(varphi)

求

[S(n)=sum_{i=1}^{n}varphi(i) ]

设(g=1)因为(varphi*1=id)，所以

[S(n)=frac{n*(n+1)}{2}-sum_{i=2}^{n}S(leftlfloor frac{n}{i} ight floor) ]

（二）习题

1.(bzoj3944)

解析上面知识部分有

所以直接上代码，一定要注意，杜教筛的时候循环从2开始

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<map>
typedef long long ll;
typedef double dd;
#define For(i,j,k) for (ll i=j;i<=k;++i)
#define Forr(i,j,k) for (ll i=j;i>=k;--i)
#define Set(a,p) memset(a,p,sizeof(a))
using namespace std;

template<typename T>bool chkmax(T& a,T b) {return a<b?a=b,1:0;}
template<typename T>bool chkmin(T& a,T b) {return a>b?a=b,1:0;}

const int maxn=2500000+100;
const ll maxx=2500000;
struct node {
	ll phi,miu;
};
ll N;
ll p[maxn],vis[maxn];
ll phi[maxn],miu[maxn];
map<ll,ll>Phi,Miu;

inline void read(ll &x) {
	x=0;
	ll p=1;
	char c=getchar();
	while (!isdigit(c)) {if (c=='-') p=-1; c=getchar();}
	while (isdigit(c)) {x=(x<<1)+(x<<3)+(c-'0'); c=getchar();}
	x*=p;
}

inline void init() {
	miu[1]=1; phi[1]=1;
	For (i,2,maxx) {
		if (!vis[i]) {
			p[++N]=i; miu[i]=-1; phi[i]=i-1;
		}
		For (j,1,N) {
			ll k=i*p[j];
			if (k>maxx) break;
			vis[k]=1;
			if (i%p[j]==0) {
				miu[k]=0; phi[k]=phi[i]*p[j];
				break;
			}
			else {
				miu[k]=-miu[i]; phi[k]=phi[i]*(p[j]-1);
			}
		}
	}
	For (i,1,maxx) {
		phi[i]=phi[i-1]+phi[i];
		miu[i]=miu[i-1]+miu[i];
	}
}

node solve(ll n) {
	if (n<=maxx) return (node){phi[n],miu[n]};
	if (Phi[n]) return (node){Phi[n],Miu[n]};
	ll ans1=n*(n+1)/2,ans2=1;
	for (ll i=2,nxt=1;i<=n;i=nxt+1) {
		node anss=solve(n/i);
		nxt=(n/(n/i));
		ans1-=(nxt-i+1)*anss.phi; ans2-=(nxt-i+1)*anss.miu;
	}
	Phi[n]=ans1; Miu[n]=ans2;
	return (node){ans1,ans2};
}

int main() {
	init();
	ll tt; read(tt);
	while (tt--) {
		ll n; read(n);
		node ans=solve(n);
		printf("%lld %lld
",ans.phi,ans.miu);
	}
	return 0;
}