●杜教筛入门(BZOJ 3944 Sum)

入门杜教筛啦。

http://blog.csdn.net/skywalkert/article/details/50500009（好文！）

可以在$O(N^{\frac{2}{3}})或O(N^{\frac{3}{4}})$的复杂度内解决求某些数论函数f(n)（或f的前缀和S(n)$）的值。

先来看看原理是什么。(接下来推导如何求数论函数f(n)的前缀和S(n)) 

---

现在有两个数论函数$f( )和g( )$

(同时定义f的前缀和函数$S(n)=\sum_{i=1}^{n}f(i)$)

有狄利克雷乘积可知：

$$f*g(n)=\sum_{i|n}f(\frac{n}{i})g(i)\quad(=\sum_{i|n}f(i)g(\frac{n}{i}))$$

那么，则有如下结论：

$$\sum_{n=1}^{N}f*g(n)=\sum_{i=1}^{N}g(i)S(\lfloor \frac{N}{i} \rfloor)$$

证明如下：

$$\begin{align*}
\sum_{n=1}^{N}f*g(n)&=\sum_{n=1}^{N}\sum_{i|n}f(\frac{n}{i})g(i)\\
&=\sum_{i=1}^{N}g(i)\sum_{i=1}^{\lfloor \frac{N}{i} \rfloor}f(i)\\
&=\sum_{i=1}^{N}g(i)S(\lfloor \frac{N}{i} \rfloor)
\end{align*}$$

然后把右边和式里的$g(1)S(N)$那一项提出来得到：

$$g(1)S(N)=\sum_{n=1}^{N}f*g(n)-\sum_{i=2}^{N}g(i)S(\lfloor \frac{N}{i} \rfloor)$$

通常令数论函数$g(n)=I(n)=1$(恒等函数$l(n)=1$，完全积性)

到目前为止，上式就是我们进行杜教筛的基础了。

因为左边的S(N)就是答案，而右边同时又可以用分块的方式计算。

不少刚刚入门的同学会有一个疑问，等式右边的后半部分确实可以分块计算，但是前半部分怎么办呢？

其实，一般前面的$\sum_{n=1}^{N}f*g(n)$都是可以O(1)计算出来的。

下面来看两个例子：

---

(一)、求莫比乌斯函数$\mu(n)$的前缀和函数$S(n)，n \leq 10^9$

首先添加一个辅助函数g(n)=l(n)=1，

然后重复上面的过程，可以得到

$$g(1)S(N)=\sum_{n=1}^{N}\mu*g(n)-\sum_{i=2}^{N}g(i)S(\lfloor \frac{N}{i} \rfloor)$$

$$S(N)=\sum_{n=1}^{N}\mu*g(n)-\sum_{i=2}^{N}S(\lfloor \frac{N}{i} \rfloor)$$

现在来看看$\sum_{n=1}^{N}\mu*g(n)$怎么求：

$$\begin{aligned}
\sum_{n=1}^{N}\mu*g(n)&=\sum_{n=1}^{N}\sum_{i|n}\mu(i)g(\frac{n}{i})\\
&=\sum_{n=1}^{N}\sum_{i|n}\mu(i)\\
&=\sum_{n=1}^{N}[n==1]\\
&=1
\end{aligned}$$

上面的化简用到了刚刚学莫比乌斯函数时的一个结论：

$$\sum_{i|n}\mu(i)=[n==1]$$

到此，我们得到：

$$S(N)=1-\sum_{i=2}^{N}S(\lfloor \frac{N}{i} \rfloor)$$

实现方式是分块计算+记忆化递归处理(用map或者hash表记忆化)

---

（二）. 求欧拉函数$\phi(n)$的前缀和函数$S(n)，n \leq 10^9$

同样地，添加一个辅助函数g(n)=l(n)=1，

然后重复上面的过程，可以得到

$$g(1)S(N)=\sum_{n=1}^{N}\phi*g(n)-\sum_{i=2}^{N}g(i)S(\lfloor \frac{N}{i} \rfloor)$$

$$S(N)=\sum_{n=1}^{N}\phi*g(n)-\sum_{i=2}^{N}S(\lfloor \frac{N}{i} \rfloor)$$

$\sum_{n=1}^{N}\phi*g(n)$又怎样求呢：

$$\begin{aligned}
\sum_{n=1}^{N}\phi*g(n)&=\sum_{n=1}^{N}\sum_{i|n}\phi(i)g(\frac{n}{i})\\
&=\sum_{n=1}^{N}\sum_{i|n}\phi(i)\\
&=\sum_{n=1}^{N}n\\
&=\frac{(1+n)n}{2}
\end{aligned}$$

上面的化简用到了这样一个结论：

$$\sum_{i|n}\phi(i)=n$$

所以我们得到：

$$S(N)=\frac{(1+n)n}{2}-\sum_{i=2}^{N}S(\lfloor \frac{N}{i} \rfloor)$$

这个同样是实现方式是分块计算+记忆化递归处理(用map或者hash表记忆化)

---

下面是代码具体实现：

关于时间复杂度的分析不太会。记了一下结论。

通常不做任何处理，就直接杜教筛的话（分块计算+记忆化递归处理），复杂度是$O(N^{\frac{3}{4}})$

但是如果预处理出前$N^{\frac{2}{3}}$个前缀和，那么总的复杂度就可以降到$O(N^{\frac{2}{3}})$

BZOJ 3944: Sum，杜教筛入门题。

多个询问，给出N，求$\sum_{n=1}^{N}\mu(n)$和$\sum_{n=1}^{N}\phi(n)$

也就是求上面的两个例子。

这里直接给出代码，用的是预处理前$N^{\frac{2}{3}}$个前缀和+hash表进行记忆化。

复杂度$O(N^{\frac{2}{3}})$

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define DJM 1664510
//#define DJM 10
#define ll long long
using namespace std;
ll phi[DJM+50],mu[DJM+50];
struct Pii{
	int x; ll a,b;
	Pii(int _x=0,ll _a=0,ll _b=0):x(_x),a(_a),b(_b){}
}nl;
struct Hash_Table{//
	#define hmod 1425367
	int nxt[hmod],head[hmod],hnt;
	Pii info[hmod];
	Hash_Table(){hnt=2;}
	void Push(Pii rtm){
		static int u; u=rtm.x%hmod;
		info[hnt]=rtm; nxt[hnt]=head[u]; head[u]=hnt++;
	}
	Pii Find(int x){
		static int u; u=x%hmod;
		for(int i=head[u];i;i=nxt[i]) if(info[i].x==x) return info[i];
		return nl;
	}
}H;
void Sieve(){
	static bool np[DJM+50];
	static int prime[DJM+50],pnt;
	phi[1]=mu[1]=1;
	for(int i=2;i<=DJM;i++){
		if(!np[i]) prime[++pnt]=i,mu[i]=-1,phi[i]=i-1;
		for(int j=1;j<=pnt&&i<=DJM/prime[j];j++){
			np[i*prime[j]]=1;
			if(i%prime[j]){mu[i*prime[j]]=-mu[i]; phi[i*prime[j]]=phi[i]*phi[prime[j]];}
			else{phi[i*prime[j]]=phi[i]*prime[j]; break;}
		}
	}
	for(int i=2;i<=DJM;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
Pii DJ_Sieve(int x){
	if(x<=DJM) return Pii(x,mu[x],phi[x]);
	if(H.Find(x).x) return H.Find(x);
	Pii tmp,now=Pii(x,1,(1ll+x)*x/2);
	for(ll i=2,last;i<=x;i=last+1){
		last=x/(x/i); tmp=DJ_Sieve(x/i);
		now.a-=tmp.a*(last-i+1); now.b-=tmp.b*(last-i+1);
	}
	H.Push(now); return now;
}
int main(){
	Sieve();
	int Case,n; Pii ans;
	scanf("%d",&Case);
	for(int i=1;i<=Case;i++){
		scanf("%d",&n);
		if(n==0) {printf("0 0\n"); continue;}
		ans=DJ_Sieve(n);
		printf("%lld %lld\n",ans.b,ans.a);
	}
	return 0;
}