【BZOJ 1101】[POI2007] Zap

题目传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=1101
离线版题目：http://paste.ubuntu.com/20151271/

莫比乌斯反演第一题！ 1A！ *★,°*:.☆\(￣▽￣)/$:*.°★*
然而，之前翻来覆去看题解不下5遍QAQ

一开始想拿phi函数来搞，结果发现搞不了。
于是还是老老实实搞莫比乌斯：
不难发现，我们实际是要求这个东西：\(\sum\limits_{i = 1}^{\frac{a}{d}} {\sum\limits_{j = 1}^{\frac{b}{d}} {[gcd(i,j) = 1]} } \)
然后利用这个式子：\(\sum\limits_{d|n} {\mu (d) = \left\{ \begin{array}{l}1\left( {n = 1} \right)\\0(n > 1)\end{array} \right.} \)
就可以化成下面这个式子：\(\sum\limits_{i = 1}^{\frac{a}{d}} {\sum\limits_{j = 1}^{\frac{b}{d}} {\sum\limits_{k|\gcd (i,j)} {\mu (k)} } } \)
更进一步，因为 k|gcd(i,j)\( \Leftrightarrow \)k|i && k|j ,所以就可以化成这个式子：\(\sum\limits_{i = 1}^{\frac{a}{d}} {\sum\limits_{j = 1}^{\frac{b}{d}} {\sum\limits_{k|i} {\sum\limits_{k|j} {\mu (k)} } } } \)
然后把莫比乌斯函数前提，得到：\(\sum\limits_{k = 1}^{\min (\frac{a}{d},\frac{b}{d})} {\mu (k)\sum\limits_{i = 1}^{\frac{a}{d}} {\sum\limits_{j = 1}^{\frac{b}{d}} {\sum\limits_{k|i} {\sum\limits_{k|j} 1 } } } } \)
最后化简得到：\(\sum\limits_{k = 1}^{\min (\frac{a}{d},\frac{b}{d})} {\mu (k)\left\lfloor {\frac{a}{{d \cdot k}}} \right\rfloor \left\lfloor {\frac{b}{{d \cdot k}}} \right\rfloor } \)
然后就是分块了，因为\({\left\lfloor {\frac{a}{{d \cdot k}}} \right\rfloor }\)只会有\(\sqrt a \)种取值

#include<iostream>
#include<cstdio>
#include<bitset>
using namespace std;

const int MAXN = 50000+9;
const int MX = 50000;

int mu[MAXN],pri[MAXN],tot,sum[MAXN];
bitset<MAXN> tag;

inline int read(){
	char c=getchar(); int buf=0,f=1;
	while (c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
	while (c<='9'&&c>='0'){buf=buf*10+c-'0';c=getchar();}
	return buf*f;
}

inline void Get_mu(){
	mu[1] = 1;
	for (int i=2;i<=MX;i++){
		if (!tag[i]) pri[++tot] = i, mu[i] = -1;
		for (int j=1;j<=tot && pri[j]*i <= MX;j++) {
			tag[pri[j]*i] = 1;
			if (i % pri[j]) mu[pri[j]*i] = -mu[i];
			else {mu[pri[j]*i] = 0; break;}
		}
	}
	for (int i=1;i<=MX;i++) sum[i] = sum[i-1] + mu[i];
}

int main(){
	Get_mu();
	int T = read(),a,b,k,ret; 
	while (T--) {
		a = read(); b = read(); k = read();
		a = a/k; b = b/k; ret = 0;
		for (int i=1,lim=min(a,b),tmp;i<=lim;i=tmp+1){
			tmp = min(a/(a/i),b/(b/i));
			ret += (sum[tmp]-sum[i-1])*(a/i)*(b/i);
		}
		printf("%d\n",ret);
	}
	return 0;
}