BZOJ.2301.[HAOI2011]Problem B(莫比乌斯反演 容斥)

[Update] 我好像现在都看不懂我当时在写什么了=-=

\(Description\)

　　求\(\sum_{i=a}^b\sum_{j=c}^d[(i,j)=k]\)

\(Solution\)

　　首先是把下界作为1.可以化为求

\[\sum_{i=1}^{\lfloor\frac{N}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{k}\rfloor}[(i,j)=1] \]

说明：大概就我不能直接看出来了。。
　　首先要求\([1,N]\)中有多少\(i，i|k\)，再求[1,j]中有多少\(j，j|k且(i,j)=1\)，显然这个\(i,j\)的上界就分别是\(\lfloor\frac{N}{k}\rfloor,\lfloor\frac{M}{k}\rfloor\)，答案就是\((i,j)=1\)的\((i,j)\)数对个数。
　　现在考虑如何求上面的式子。
　　由莫比乌斯反演，有

\[F(d)=\sum_{d|n}f(n)\Leftrightarrow f(d)=\sum_{d|n}\mu(\frac{n}{d})F(n) \]

　　设\(f[i]\)为满足\(gcd(x,y)=i\)的\((x,y)\)对数，其中\(1\leq x\leq b,1\leq y\leq d\)；
　　\(F[i]\)为满足\(i|gcd(x,y)\)的\((x,y)\)对数，其中\(1\leq x\leq b,1\leq y\leq d\)。
　　显然有\(F[i]=\sum_{i|n}f[n]\Rightarrow f[i]=\sum_{i|n}\mu(\frac{n}{i})F[n]\)
　　又显然有\(F[i]=\lfloor\frac{b}{i}\rfloor\lfloor\frac{d}{i}\rfloor\)，那么

\[f[i]=\sum_{i|n,1\leq n\leq min(b,d)}\mu(\frac{n}{i})\lfloor\frac{b}{n}\rfloor\lfloor\frac{d}{n}\rfloor \]

　　令\(k=\frac{n}{i}\)，即\(n=ki\)，令\(b'=\frac{b}{i},d'=\frac{d}{i}\)，则

\[f[i]=\sum_{k=1}^{min(b',d')}\mu(k)\lfloor\frac{b'}{k}\rfloor\lfloor\frac{d'}{k}\rfloor \]

　　(本题i就是1。)
　　上面这个式子还是\(O(n^2)\)的。。还是要分块计算。

/*
最后求解要容斥:(a,c为开区间，另外其实并不分左右)
ans=f[a~b,c~d]=f[b,d]-f[a,d]-f[b,c]+f[a,c]
*/
#include<cstdio>
#include<cctype>
#include<algorithm>
#define gc() (SS==TT&&(TT=(SS=IN)+fread(IN,1,MAXIN,stdin),SS==TT)?EOF:*SS++)
const int N=5e4+3,MAXIN=1<<17;

int cnt,P[N+2],mu[N+2],sum[N+2];
bool Not_P[N+2];
char IN[MAXIN],*SS=IN,*TT=IN;

inline int read()
{
	int now=0,f=1;register char c=gc();
	for(;!isdigit(c);c=gc()) if(c=='-') f=-1;
	for(;isdigit(c);now=now*10+c-'0',c=gc());
	return now*f;
}
void Init()
{
	mu[1]=1;
	for(int i=2;i<N;++i)
	{
		if(!Not_P[i]) P[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*P[j]<N;++j)
		{
			Not_P[i*P[j]]=1;
			if(!(i%P[j])) {mu[i*P[j]]=0; break;}
			mu[i*P[j]]=-mu[i];
		}
	}
	for(int i=1;i<N;++i) sum[i]=sum[i-1]+mu[i];
}
int calc(int n,int m)
{
	int t=std::min(n,m),ans=0;//应该不需要longlong 
//	for(int k=1;k<=t;++k) ans+=mu[k]*(n/k)*(m/k);//TLE:O(n^2)
	for(int las,i=1;i<=t;i=las+1)
	{
		las=std::min(n/(n/i),m/(m/i));
		ans+=(sum[las]-sum[i-1])*(n/i)*(m/i);
	}
	return ans;
}

int main()
{
	Init();
	int t=read(),a,b,c,d,k;
	while(t--)
		a=read()-1,b=read(),c=read()-1,d=read(),k=read(),
		a/=k,b/=k,c/=k,d/=k,printf("%d\n",calc(b,d)-calc(a,d)-calc(b,c)+calc(a,c));
	return 0;
}