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

$n,m\le 1e7$，求$\sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j)$
不妨设$n\le m$，则原式
$=\sum _{i=1}^n\sum _{j=1}^m\frac{i\ast j}{gcd(i,j)}$
$=\sum _{g=1}^n\frac{1}{g}\sum _{i=1}^n\sum _{j=1}^{m}ij[gcd(i,j)=g]$
$=\sum _{g=1}^{n}g\sum _{i=1}^{n/g}\sum _{j=1}^{m/g}ij[gcd(i,j)=1]$
$=\sum _{g=1}^{n}g\sum _{i=1}^{n/g}\sum _{j=1}^{m/g}ij\sum _{d|gcd(i,j)}\mu (d)$
$=\sum _{g=1}^{n}g\sum _{d=1}^{n/g}\mu (d)\sum _{i=1}^{n/gd}\sum _{j=1}^{m/gd}ij{d}^2$
$=\sum _{g=1}^{n}g\sum _{d=1}^{n/g}\mu (d)d^2\frac{\frac{n}{gd}\ast (\frac{n}{gd}+1)}{2}\ast \frac{\frac{m}{gd}\ast (\frac{m}{gd}+1)}{2}$
$\mu (d)d^2$可以线性时间预处理前缀和，两个分式$O(1)$计算得到，对于内外两层各做一次分块。时间复杂度为$O(\sqrt{n}\ast \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;
}