[国家集训队] Crash的数字表格

考虑到(lcm(i,j)=frac{ij}{gcd(i,j)})

(sum_{i=1}^nsum_{j=1}^mfrac{ij}{gcd(i,j)})

(sum_{d=1}^{n}sum_{i=1}^nsum_{j=1}^m[gcd(i,j)==d]frac{ij}{d})

(sum_{d=1}^{n}sum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]{ijd})

(=sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]{ij})

看后面
(sum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]{ij})
(=sum_{i=1}^{x}sum_{j=1}^{y}[gcd(i,j)==1]{ij})

考虑反演
(f(d)=sum_{i=1}^{x}sum_{j=1}^{y}[gcd(i,j)==d]{ij})

(G(d)=sum_{i=1}^{x}sum_{j=1}^{y}[d|gcd(i,j)]{ij})

(G(d)=d^2sum_{i=1}^{x/d}sum_{j=1}^{y/d}[1|gcd(i,j)]{ij})
(G(d)=d^2sum_{i=1}^{x/d}sum_{j=1}^{y/d}{ij})

于是
(f(1)=sum_{i=1}^xmu(i)G(i))

(ans=sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]{ij})

在这里做第一次分块，而后
(f(1)=sum_{i=1}^xmu(i)G(i))

(=sum_{i=1}^xmu(i)i^2sum_{i=1}^{x/d}sum_{j=1}^{y/d}{ij})

然后内外各做一次整除分块即可

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 12000005;
const int mod = 20101009;
int mu[N+5],sum[N+5],pr[N+5],is[N+5],cnt;

int h(int n,int m) {
    int l=1,r,ans=0;
    while(l<=n) {
        r=min(n/(n/l),m/(m/l));
        ans += (sum[r]-sum[l-1]+mod)%mod
            * ((n/l)*(n/l+1)/2%mod)%mod
            * ((m/l)*(m/l+1)/2%mod)%mod;
        ans %= mod;
        l=r+1;
    }
    return ans;
}

signed main() {
    mu[0]=mu[1]=1; is[1]=1;
	for(int i=2;i<=N;i++) {
		if(is[i]==0) {
			pr[++cnt]=i;
			mu[i]=-1;
		}
		for(int j=1; j<=cnt&&pr[j]*i<N; ++j) {
			is[pr[j]*i]=1;
			if(i%pr[j]==0) {
				mu[pr[j]*i]=0;
				break;
			}
			else {
				mu[pr[j]*i]=-mu[i];
			}
		}
	}
	for(int i=1;i<=N;i++) sum[i]=(sum[i-1]+mu[i]*i*i%mod)%mod;
    int n,m;
    cin>>n>>m;
    if(n>m) swap(n,m);
    int l=1,r,ans=0;
    while(l<=n) {
        r=min(n/(n/l),m/(m/l));
        ans+=(r-l+1)*(l+r)/2%mod*h(n/l,m/l)%mod;
        ans%=mod;
        l=r+1;
    }
    cout<<ans<<endl;
}