【BZOJ4176】Lucas的数论（杜教筛）

点此看题面

大致题意： 令\(f(i)\)表示\(i\)的约数个数，求\(\sum_{i=1}^n\sum_{j=1}^nf(ij)\)。

前言

刚写完【BZOJ2627】JZPKIL再来写这道题，发现其他数论题都好简单啊。。。

话说上面这道题是数论大杂烩，却缺少了此题的线性筛+杜教筛，这么看来我是不是一个下午把数论算法基本复习了一遍。。。

\(f(ij)\)的性质

考虑我们枚举\(x|i,y|j\)以求\(f(ij)\)，显然不可能直接对于每一对\(x,y\)计算一次贡献，这样一定会出现重复的因数。

因此，我们规定\(\frac ix\)与\(y\)互质时才能计算贡献，容易发现，这样一来贡献的计算就不重不漏了。

而实际上，由于我们枚举的是\(x|i\)，因此直接看作\(x\)与\(y\)互质时算贡献也未尝不可，而且还方便了许多。

也就是说：

\[f(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)=1] \]

推式子

把上面的式子代入原式：

\[\sum_{i=1}^n\sum_{j=1}^n\sum_{x|i}\sum_{y|j}[gcd(x,y)=1] \]

喜闻乐见枚举\(gcd\)，得到：

\[\sum_{d=1}^n\mu(d)\sum_{x=1}^{\lfloor\frac nd\rfloor}\lfloor\frac n{dx}\rfloor\sum_{y=1}^{\lfloor\frac nd\rfloor}\lfloor\frac n{dy}\rfloor=\sum_{d=1}^n\mu(d)(\sum_{x=1}^{\lfloor\frac nd\rfloor}\lfloor\frac n{dx}\rfloor)^2 \]

于是我们对于\(d\)除法分块，然后杜教筛筛出\(\mu\)，再套个除法分块处理后面的和。

总复杂度就感性理解一下吧，反正推完式子感觉应该能过就写了一发，结果跑的还挺快（单个点最多也就\(0.7s\)）。

一开始线性筛范围手贱多摁了个0结果MLE了，我说怎么样例都要跑一秒多。。。

代码

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define X 1000000007
#define Inc(x,y) ((x+=(y))>=X&&(x-=X))
using namespace std;
int n,sn;
class DuSieve
{
	private:
		#define SZ 5000000
		int Pt,P[SZ+5],mu[SZ+5],f1[SZ+5],f2[SZ+5];
	public:
		I DuSieve()//线性筛预处理
		{
			mu[1]=1;for(RI i=2,j;i<=SZ;++i)
				for(!P[i]&&(mu[P[++Pt]=i]=X-1),j=1;j<=Pt&&i*P[j]<=SZ;++j)
					if(P[i*P[j]]=1,i%P[j]) mu[i*P[j]]=X-mu[i];else break;
			for(RI i=2;i<=SZ;++i) Inc(mu[i],mu[i-1]);
		}
		I int SMu(CI x)//杜教筛筛μ
		{
			if(x<=SZ) return mu[x];//较小数据直接返回
			if(x<=sn) {if(f1[x]) return f1[x];}else {if(f2[n/x]) return f2[n/x];}//记忆化
			RI l,r,t=1;for(l=2;l<=x;l=r+1) r=x/(x/l),t=(t-1LL*(r-l+1)*SMu(x/l)%X+X)%X;
			return (x<=sn?f1[x]:f2[n/x])=t;
		}
}D;
I int Sum(CI x)//除法分块求和
{
	RI l,r,t=0;for(l=1;l<=x;l=r+1) r=x/(x/l),t=(1LL*(r-l+1)*(x/l)+t)%X;return t;
}
int main()
{
	RI l,r,k,t=0;for(scanf("%d",&n),sn=sqrt(n),l=1;l<=n;l=r+1)//除法分块
		r=n/(n/l),k=Sum(n/l),t=(1LL*k*k%X*(D.SMu(r)-D.SMu(l-1)+X)+t)%X;
	return printf("%d",t),0;
}