BZOJ 2154 Crash的数字表格莫比乌斯反演-优快云博客

本文链接：https://blog.youkuaiyun.com/qq_40859782/article/details/100806952

本文探讨了在n,m≤1e7的约束下，如何高效计算∑i=1n∑j=1mlcm(i,j)的问题。通过数学转换，将问题转化为更易于处理的形式，并利用线性时间预处理μ(d)d2的前缀和，实现O(n)的时间复杂度。文章详细介绍了算法步骤，包括分块技巧和具体实现代码。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

$n,m≤1e7n,m\leq1e7$ ，求 $∑i=1n∑j=1mlcm(i,j)\sum\limits_{i=1}^n\sum\limits_{j=1}^{m}lcm(i,j)$
不妨设 $n≤mn\leq m$ ，则原式
$=∑i=1n∑j=1mi∗jgcd(i,j)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{i*j}{gcd(i,j)}$
$=∑g=1n1g∑i=1n∑j=1mij[gcd(i,j)=g]=\sum\limits_{g=1}^{n}\frac{1}{g}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[gcd(i,j)=g]$
$=∑g=1ng∑i=1n/g∑j=1m/gij[gcd(i,j)=1]=\sum\limits_{g=1}^{n}g\sum\limits_{i=1}^{n/g}\sum\limits_{j=1}^{m/g}ij[gcd(i,j)=1]$
$=∑g=1ng∑i=1n/g∑j=1m/gij∑d∣gcd(i,j)μ(d)=\sum\limits_{g=1}^{n}g\sum\limits_{i=1}^{n/g}\sum\limits_{j=1}^{m/g}ij\sum\limits_{d|gcd(i,j)}\mu(d)$
$=∑g=1ng∑d=1n/gμ(d)∑i=1n/gd∑j=1m/gdijd2=\sum\limits_{g=1}^{n}g\sum\limits_{d=1}^{n/g}\mu(d)\sum\limits_{i=1}^{n/gd}\sum\limits_{j=1}^{m/gd}ijd^2$
$=∑g=1ng∑d=1n/gμ(d)d2ngd∗(ngd+1)2∗mgd∗(mgd+1)2=\sum\limits_{g=1}^{n}g\sum\limits_{d=1}^{n/g}\mu(d)d^2\frac{\frac{n}{gd}*(\frac{n}{gd}+1)}{2}*\frac{\frac{m}{gd}*(\frac{m}{gd}+1)}{2}$
$μ(d)d2\mu(d)d^2$ 可以线性时间预处理前缀和，两个分式 $O (1)$ 计算得到，对于内外两层各做一次分块。时间复杂度为 $O(n∗n)=O(n)O(\sqrt n *\sqrt n)=O(n)$ 。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int inf=0x3f3f3f3f;
const ll INF=LONG_LONG_MAX;
const int N=1e7+7;
int pri[N],tot=0;
bool ok[N];
int mu[N];
ll sum[N];
ll n,m;
const int mod=20101009; 
void getmu() {
	mu[1]=1;
	for(int i=2;i<=n;i++) {
		if(!ok[i]) pri[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*pri[j]<=n;j++) {
			ok[i*pri[j]]=1;
			if(i%pri[j]) mu[i*pri[j]]=-mu[i];
			else { mu[i*pri[j]]=0;break; }
		}
	}
	sum[1]=mu[1];
	for(int i=2;i<=n;i++) {
		sum[i]=sum[i-1]+((1LL*i*i)%mod*mu[i])%mod;
		sum[i]%=mod;	
	}
} 
ll cal(ll x) { return (x*(x+1)>>1)%mod; }
ll sub(ll x,ll y) {
	if(x-y>=0) return x-y;
	else return x-y+mod; 
}
ll solve(ll n,ll m) {
	ll ans=0;
	for(int l=1,r=0;l<=n;l=r+1) {
		r=min(n/(n/l),m/(m/l));
		ans+=sub(sum[r],sum[l-1])*cal(n/l)%mod*cal(m/l)%mod; 
		ans%=mod;
	}
	return ans;
}
int main() { 
	scanf("%lld%lld",&n,&m);
	if(n>m) swap(n,m);
	getmu();  
	ll ans=0;
	for(int l=1,r=0;l<=n;l=r+1) {
		r=min(n/(n/l),m/(m/l));
		ans+=(1LL*(l+r)*(r-l+1)/2)%mod*solve(n/l,m/l)%mod;
		ans%=mod;
	}
	printf("%lld\n",ans);
	return 0;
}