【NOI2010】能量采集

题意简述

求：

\[2\sum_{i=1}^n\sum_{j=1}^m(gcd(i,j)-1)=-nm+2\sum_{i=1}^n\sum_{j=1}^mgcd(i,j) \]

题解

先说一个有点巧妙的变化,等会要用到(下面的k是给出的定值)。
首先有一个简单的结论:

\[\epsilon(n)=\sum_{d|n}\mu(d) \]

也可表示为狄利克雷卷积的形式:

\[\epsilon=\mu*1 \]

由此得出这个变化:

\[\begin{align*} \sum_{i=1}^n\sum_{j=1}^m [gcd(i, j)=1]&=\sum_{i=1}^n\sum_{j=1}^m \sum_{d|gcd(i,j)}\mu(d)\\ &=\sum_{i=1}^n\sum_{j=1}^m \sum_{d|i,d|j}\mu(d)\\ &=\sum_{d=1}^{min(n, m)}\mu(d)\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor\\ \end{align*}\]

上面的最后一步是将枚举\(i,j\)转换为枚举其公因数\(d\)，再由一个数\(d\)在\(1\sim n\)中的倍数有\(\lfloor\frac{n}{d}\rfloor\)个转化得。

还有一个关于取整函数的性质:

\[\lfloor\lfloor\frac{x}{a}\rfloor/b\rfloor=\lfloor\frac{x}{ab}\rfloor \]

所以接下来就可以这样推式子了:

\[\begin{align*} \sum_{i=1}^n\sum_{j=1}^mgcd(i,j)& =\sum_{d=1}^{min(n,m)}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i, j)=1]\\ &=\sum_{d=1}^{min(n,m)}d\sum_{k=1}^{min(\lfloor\frac{n}{d}\rfloor, \lfloor\frac{m}{d}\rfloor)}\mu(k)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor\\ \end{align*}\]

枚举\(d\)，然后后面依据\(\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor\)的取值，相同取值的区间合并，答案加上相应取值与区间和的乘积，即做一个数论分块，本题就可到此结束。

代码

#include <cstdio>
#include <cctype>
typedef long long ll;
const int maxn=1e5+10;
int prim[maxn],mu[maxn],s[maxn];
int tot;
bool vis[maxn];

int min(int x,int y) {return x<y?x:y;}
void swap(int &x,int &y) {x^=y^=x^=y;}
void prework(int n)
{
	vis[0]=vis[1]=1;
	mu[1]=s[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!vis[i])
			prim[++tot]=i,mu[i]=-1;
		s[i]=s[i-1]+mu[i];
		for (int j=1;j<=tot&&(ll)i*prim[j]<=n;j++)
		{
			vis[i*prim[j]]=1;
			if (i%prim[j]==0)
			{
				mu[i*prim[j]]=0;
				break;
			}
			mu[i*prim[j]]=-mu[i];
		}
	}
}
int read()
{
	int res=0;
	char ch=getchar();
	while(!isdigit(ch))
		ch=getchar();
	while(isdigit(ch))
		res=res*10+ch-'0',ch=getchar();
	return res;
}
int main()
{
	int n=read(),m=read();
	ll ans=0;
	if (n>m)
		swap(n, m);
	prework(n);
	for (int d=1;d<=n;d++)
	{
		for (int l=1,r=1;l<=n/d;l=r+1)
		{
			r=min(n/d/(n/d/l), m/d/(m/d/l));
			ans+=(ll)d*(s[r]-s[l-1])*(n/d/l)*(m/d/l);
		}
	}
	printf("%lld\n",ans*2-(ll)n*m);
	return 0;
}