[BZOJ2154]Crash的数字表格

题目求“是质数”，我们可以枚举每个质数然后统计答案，假设$n\leq m$

$$\begin{align*}\sum\limits_{\substack{p是质数\\p\leq n}}\ \sum\limits_{i=1}^n\sum\limits_{j=1}^m\left[(i,j)=p\right]&=\sum\limits_{\substack{p是质数\\p\leq n}}\ \sum\limits_{i=1}^{\left\lfloor\frac np\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac mp\right\rfloor}\left[(i,j)=1\right]\\&=\sum\limits_{\substack{p是质数\\p\leq n}}\ \sum\limits_{i=1}^{\left\lfloor\frac np\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac mp\right\rfloor}\sum\limits_{d|(i,j)}\mu(d)\\&=\sum\limits_{\substack{p是质数\\p\leq n}}\ \sum\limits_{i=1}^{\left\lfloor\frac np\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac mp\right\rfloor}\sum\limits_{\substack{d|i\\d|j}}\mu(d)\\&=\sum\limits_{\substack{p是质数\\p\leq n}}\ \sum\limits_{d=1}^{\left\lfloor\frac np\right\rfloor}\mu(d)\left\lfloor\dfrac n{dp}\right\rfloor\left\lfloor\dfrac m{dp}\right\rfloor\end{align*}$$

“$p$是质数”这个限制条件很烦人，我们要把它放到内层的sigma里，令$k=dp$，注意到原来内层对$d$的限制是$d\leq\left\lfloor\dfrac np\right\rfloor$，即$k\leq n$，所以我们可以把它改写为$\sum\limits_{k=1}^n\sum\limits_{\substack{p是质数\\p|k}}\mu\left(\dfrac kp\right)\left\lfloor\dfrac nk\right\rfloor\left\lfloor\dfrac mk\right\rfloor$

求和符号的形式决定了它可以快速求和，问题在于中间那个$\sum\limits_{\substack{p是质数\\p|k}}\mu\left(\dfrac kp\right)$，显然我们是要预处理它的

用比较暴力的方式就可以了，即枚举$p$和它的所有倍数，因为一个数$\leq n$的倍数数量是均摊$O(\ln n)$的（$\sum\limits_{i=2}^n\left\lfloor\dfrac ni\right\rfloor\lt n\int_1^n\dfrac 1xdx=n\ln n-n$），而且$\pi(x)=\Theta\left(\dfrac x{\ln x}\right)$（一个上下界很松的估计是$\left(\dfrac{\ln2}3\right)\dfrac x{\ln x}\lt\pi(x)\lt(6\ln2)\dfrac x{\ln x}$，证明可参考《初等数论》中“$\pi(x)$的上、下界估计”），所以枚举所有$\leq n$的质数的倍数是均摊$O(n)$的，这种看似暴力的做法其实并不会超时

于是就做完了，学到一个套路：出现$[x=1]$这种东西的时候可以把它转化为莫比乌斯函数看是否能简化计算

这是...反演套反演么==

#include<stdio.h>
#define mod 20101009
#define ll long long
int mu[10000010],pr[10000010];
bool np[10000010];
int ad(int a,int b){return(a+b)%mod;}
int mul(int a,int b){return a*(ll)b%mod;}
void sieve(int n){
	int i,j,m=0;
	np[1]=1;
	mu[1]=1;
	for(i=2;i<=n;i++){
		if(!np[i]){
			m++;
			pr[m]=i;
			mu[i]=-1;
		}
		for(j=1;j<=m;j++){
			if(pr[j]*(ll)i>n)break;
			np[i*pr[j]]=1;
			if(i%pr[j]==0)break;
			mu[i*pr[j]]=-mu[i];
		}
	}
	for(i=1;i<=n;i++)mu[i]=mul(mu[i],mul(i,i));
	for(i=2;i<=n;i++)mu[i]=ad(mu[i],mu[i-1]);
}
int min(int a,int b){return a<b?a:b;}
void swap(int&a,int&b){a^=b^=a^=b;}
int S(int n){return mul(mul(n,n+1),10050505);}
int f(int n,int m){
	int i,nex,s=0;
	for(i=1;i<=n;i=nex+1){
		nex=min(n/(n/i),m/(m/i));
		s=ad(s,mul((mu[nex]-mu[i-1])%mod,mul(S(n/i),S(m/i))));
	}
	return s;
}
int mob(int n,int m){
	int i,nex,s=0;
	for(i=1;i<=n;i=nex+1){
		nex=min(n/(n/i),m/(m/i));
		s=ad(s,mul((S(nex)-S(i-1))%mod,f(n/i,m/i)));
	}
	return s;
}
int main(){
	int n,m;
	scanf("%d%d",&n,&m);
	if(n>m)swap(n,m);
	sieve(n);
	printf("%d",(mob(n,m)+mod)%mod);
}