●BZOJ 2154 Crash的数字表格

题链：

http://www.lydsy.com/JudgeOnline/problem.php?id=2154

题解：

莫比乌斯反演。

题意还是很清楚的，就不赘述了。

显然有

$ANS=sum_{i=1;j=1}^{n;m} lcm(i,j)$

化为较为熟悉的gcd形式：

$quadquad=sum_{i=1;j=1}^{n;m} frac{i imes j}{gcd(i,j)}$

令$g$为gcd的值，$F(n,m)=sum i imes j，满足1leq i leq lfloor frac{n}{g} floor，1leq j leq lfloor frac{m}{g} floor，且gcd(i,j)=1$

那么，$ANS=sum_{g=1}^{min(n,m)} g imes F(n,m)$

如果已知$F(n,m)$的值，那么ANS就可以在$Osqrt N$的复杂度内求出。

---

接下来看看对于一个确定的g，$F(n,m)$怎么求

即我们要求满足$1leq i leq lfloor frac{n}{g} floor，1leq j leq lfloor frac{m}{g} floor$，且$gcd(i,j)=1$的所有$i imes j$的和

还记得这个题么？●HDU 1695 GCD，让求的是满足上述条件的$(i，j)$的对数。

我们仍可以类似地去求：（令$x=lfloor frac{n}{g} floor，y=lfloor frac{m}{g} floor$）

令$t(k)为gcd(i,j)=k的i imes j的和$

$T(k)为gcd(x,y)=lambda k的i imes j的和$

显然$T(k)=sum_{k|d}{f(d)}$，即T为t的倍数关系和函数

考虑一下如何计算T(k)。

首先i的取值有：$k，2k，3k，cdots，lfloor frac{x}{k} floor k$

首先j的取值有：$k，2k，3k，cdots，lfloor frac{y}{k} floor k$

任意的两两组合都是T(k)的一部分。

令$sum(a，b)=frac{(1+a)a}{2} imesfrac{(1+b)b}{2}$

所以$T(k)=k^2 imes sum(lfloor frac{x}{k} floor，lfloor frac{y}{k} floor)$

那么由莫比乌斯反演公式的形式二(倍数关系那个式子)：

$t(k)=sum_{k|d}mu(frac{d}{k})T(d)$

$quadquad=sum_{k|d}mu(frac{d}{k})d^2 imes sum(lfloor frac{x}{d} floor，lfloor frac{y}{d} floor)$

而我们要求的是t(1)，所以

$t(1)=sum_{d=1}^{min(lfloor frac{n}{g} floor，lfloor frac{m}{g} floor)}mu(d)d^2 imes sum(lfloor frac{x}{d} floor，lfloor frac{y}{d} floor)$

显然，这个式子也可以在$Osqrt N$的复杂度内求出。

所以综上，时间复杂度为$O(N)$

代码：

#include<cstdio>
#include<cstring>
#include<iostream>
#define MAXN 10000050
using namespace std;
const int mod=20101009;
int mu[MAXN],pmui2[MAXN];
void Sieve(int n){
	static bool np[MAXN];
	static int prime[MAXN],pnt;
	mu[1]=pmui2[1]=1;
	for(int i=2;i<=n;i++){
		if(!np[i]) prime[++pnt]=i,mu[i]=-1;
		for(int j=1;j<=pnt&&i<=n/prime[j];j++){
			np[i*prime[j]]=1;
			if(i%prime[j]) mu[i*prime[j]]=-mu[i];
			else{mu[i*prime[j]]=0; break;}
		}
		pmui2[i]=(pmui2[i-1]+1ll*mu[i]*i%mod*i%mod)%mod;
	}
}
int sum(int n,int m){
	return (1ll*(1+n)*n/2%mod)*(1ll*(1+m)*m/2%mod)%mod;
}
int F(int n,int m){
	int mini=min(n,m),ret=0;
	for(int d=1,last;d<=mini;d=last+1){
		last=min(n/(n/d),m/(m/d));
		ret=(1ll*ret+1ll*(pmui2[last]-pmui2[d-1]+mod)%mod*sum(n/d,m/d)%mod)%mod;
	}
	return ret;
}
int main(){
	int n,m,mini,ans=0;
	scanf("%d%d",&n,&m); mini=min(n,m);
	Sieve(mini);
	for(int g=1,last;g<=mini;g=last+1){
		last=min(n/(n/g),m/(m/g));
		ans=(1ll*ans+1ll*(g+last)*(last-g+1)/2%mod*F(n/g,m/g)%mod)%mod;
	}
	printf("%d
",(ans+mod)%mod);
	return 0;
}