luoguP1829 [国家集训队]Crash的数字表格 / JZPTAB(莫比乌斯反演)

题意

注：默认(nleqslant m)。

所求即为：(sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}lcm(i,j))
因为(i*j=gcd(i,j)*lcm(i,j)),因此原式为：
(sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}frac{i*j}{gcd(i,j)})
枚举(gcd(i,j)=d)：
(sumlimits_{d=1}^{n}frac{1}{d}*sumlimits_{i=1}^{n}sumlimits_{j=1}^mi*j*[gcd(i,j)=d])
套路地提出(d)：
(sumlimits_{d=1}^{n}frac{1}{d}sumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{m}{d}}i*d*j*d*[gcd(i,j)=1])
即：
(sumlimits_{d=1}^{n}d*sumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{m}{d}}i*j*[gcd(i,j)=1])
用莫比乌斯函数性质替换([gcd(i,j)=1])：
(sumlimits_{d=1}^{n}d*sumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{m}{d}}i*j*sumlimits_{x|gcd(i,j)}mu(x))
转而枚举(x)：
(sumlimits_{d=1}^{n}d*sumlimits_{x=1}^{frac{n}{d}}mu(x)sumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{m}{d}}i*j*[x|gcd(i,j)])
把([x|gcd(i,j)])去掉：
(sumlimits_{d=1}^{n}d*sumlimits_{x=1}^{frac{n}{d}}x^2*mu(x)sumlimits_{i=1}^{frac{n}{d*x}}sumlimits_{j=1}^{frac{m}{d*x}}i*j)
即：
(sumlimits_{d=1}^{n}d*sumlimits_{x=1}^{frac{n}{d}}x^2*mu(x)(sumlimits_{i=1}^{frac{n}{d*x}}i)(sumlimits_{j=1}^{frac{m}{d*x}}j))

求几个前缀和，第二个和第三四个都可以除法分块即可。

另一种做法

code:

#include<bits/stdc++.h>
using namespace std;
#define re register
const int maxn=10000010;
const int mod=20101009;
const int inv2=10050505;
int n,m,ans;
int mu[maxn],sum[maxn];
bool vis[maxn];
vector<int>prime;
inline void shai(int n)
{
	vis[1]=1;mu[1]=1;
	for(re int i=2;i<=n;i++)
	{
		if(!vis[i])prime.push_back(i),mu[i]=-1;
		for(re unsigned int j=0;j<prime.size()&&i*prime[j]<=n;j++)
		{
			vis[i*prime[j]]=1;
			if(i%prime[j]==0)break;
			mu[i*prime[j]]=-mu[i];
		}
	}
	for(re int i=1;i<=n;i++)sum[i]=(sum[i-1]+1ll*mu[i]*i%mod*i%mod)%mod;
}
inline int calc(int l,int r){return 1ll*(r-l+1)*(l+r)%mod*inv2%mod;}
int main()
{
	shai(10000000);
	scanf("%d%d",&n,&m);
	if(n>m)swap(n,m);
	for(re int ld=1,rd;ld<=n;ld=rd+1)
	{
		rd=min(n/(n/ld),m/(m/ld));
		int res=0;
		for(re int l=1,r;l<=n/ld;l=r+1)
		{
			r=min((n/ld)/((n/ld)/l),(m/ld)/((m/ld)/l));
			res=(res+1ll*(sum[r]-sum[l-1])*calc(1,(n/ld)/l)%mod*calc(1,(m/ld)/l)%mod)%mod;
		}
		ans=(ans+1ll*res*calc(ld,rd)%mod)%mod;
	}
	printf("%d",(ans+mod)%mod);
	return 0;
}