BZOJ 2154: Crash的数字表格 莫比乌斯反演

Description

今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张N*M的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个4*5的表格如下： 1 2 3 4 5 2 2 6 4 10 3 6 3 12 15 4 4 12 4 20 看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod 20101009的值。

Input

输入的第一行包含两个正整数，分别表示N和M。

Output

输出一个正整数，表示表格中所有数的和mod 20101009的值。 

题解：

求 $\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)$

首先，定义 $n\leqslant m$

开始推公式

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

$=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{i\times j}{gcd(i,j)}$

$=\sum_{d=1}^{n}d\times \sum_{i=1}^{\left \lfloor \frac{n}{d} \right\rfloor}\sum_{j=1}^{\left \lfloor \frac{m}{d} \right\rfloor}i\times j[gcd(i,j)==1]$

定义 $calc(n,m)=\sum_{i=1}^{n}\sum_{j=1}^{m}i\times j[gcd(i,j)==1]$

$calc(n,m)=\sum_{i=1}^{n}\sum_{j=1}^{m}i\times j\sum_{d|i,d|j}\mu(d)$

$=\sum_{d=1}^{n}\mu(d)\sum_{d|i}i\sum_{d|j}j$

$=\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{\frac{n}{d}}id\sum_{j=1}^{\frac{m}{d}}jd$

$=\sum_{d=1}^{n}\mu(d)d^{2}\sum_{i=1}^{\frac{n}{d}}d\sum_{j=1}^{\frac{m}{d}}j$

发现后面那个$\sum_{i=1}^{\frac{n}{d}}d\sum_{j=1}^{\frac{m}{d}}j$可以 $O(1)$ 求出.

$=\frac{(1+\frac{n}{d})\frac{n}{d}}{2}\times\frac{(1+\frac{m}{d})\frac{m}{d}}{2}$

我们把它定义为 $Sum(\frac{n}{d},\frac{m}{d})$

可知 $calc(n,m)=\sum_{d=1}^{n}\mu(d)d^{2}Sum(\frac{n}{d},\frac{m}{d})$

可以预处理 $\mu(d)d^2$，所以这一部分的时间复杂度是 $n^{\frac{1}{2}}$ 的.

回到原式$=\sum_{d=1}^{n}d\times calc(\frac{n}{d},\frac{m}{d})$

总时间复杂度是 $O(n)$ 的. 

#include<bits/stdc++.h>
#define maxn 10010006
#define ll long long  
using namespace std;
const long long mod = 1ll*20101009;  
void setIO(string s)
{
	string in=s+".in"; 
	freopen(in.c_str(),"r",stdin); 
}
int cnt; 
int prime[maxn], vis[maxn], mu[maxn]; 
ll sumv[maxn]; 
inline void init(int hh) 
{
	int i,j;
	mu[1]=sumv[1]=1; 
	for(i=2;i<hh;++i)
	{
		if(!vis[i]) prime[++cnt]=i, mu[i]=-1; 
		for(j=1;j<=cnt&&1ll*prime[j]*i<hh;++j) 
		{
			vis[prime[j]*i]=1; 
			if(i%prime[j]!=0) mu[i*prime[j]]=-mu[i]; 
			else 
			{
				mu[i*prime[j]]=0; 
				break; 
			}
		}
		sumv[i]=(sumv[i-1]+1ll*i*i%mod*(mu[i]+mod))%mod;    
		//  sumv2[i]=(i+sumv2[i-1])%mod; 
	} 
}
inline ll getarr(ll n,ll m)
{
	ll tmp=(1ll*n*(n+1)/2%mod) * (1ll*m*(m+1)/2%mod)%mod;  
    return tmp; 
}
inline ll calc(ll n,ll m)
{
	ll i,j; 
	ll ans=0; 
	for(i=1;i<=n;i=j+1)
	{
		j=min(n/(n/i), m/(m/i));  
		ans=(ans+1ll*(sumv[j]-sumv[i-1]+mod)*getarr(n/i,m/i)%mod)%mod;  
	}
	return ans; 
}
int main()
{
// 	setIO("input");  
	ll n,m,i,j; 
	ll ans=0; 
	scanf("%lld%lld",&n,&m); 
	if(n>m)swap(n,m);        
	init(m + 1); 
	for(i=1;i<=n;i=j+1)
	{
		j=min(n/(n/i), m/(m/i));  
		ans=(ans+1ll*((j-i+1)*(i+j)/2)%mod*calc(n/i,m/i)%mod)%mod;   
		ans=(ans+mod)%mod;  
	}  
	printf("%lld\n",ans); 
	return 0; 
}