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

( ext{Problem})

求

[left( sum_{i=1}^n sum_{j=1}^m lcm(i，j) ight) mod{20101009} ]

(n le 10^7)

( ext{Analysis})

套路地推式子

[egin{aligned} sum_{i=1}^n sum_{j=1}^m lcm(i，j) &= sum_{i=1}^n sum_{j=1}^m frac{ij}{gcd(i,j)} \ &= sum_{d=1}^{min(n,m)} sum_{d|i} sum_{d|j} frac{ij}{d} [gcd(i,j)=d] \ &= sum_{d=1}^{min(n,m)} sum_{i=1}^{lfloor frac n d floor} sum_{j=1}^{lfloor frac m d floor} ijd[gcd(i,j)=1] \ &= sum_{d=1}^{min(n,m)} d sum_{i=1}^{lfloor frac n d floor} i sum_{j=1}^{lfloor frac m d floor} j sum_{g|gcd(i,j)} mu(g) \ &= sum_{d=1}^{min(n,m)} d sum_{g=1} mu(g) sum_{g|i}^{lfloor frac n d floor} sum_{g|j}^{lfloor frac m d floor} ij \ &= sum_{d=1}^{min(n,m)} d sum_{g=1} mu(g) g^2 sum_{i=1}^{lfloor frac{n}{dg} floor} i sum_{j=1}^{lfloor frac{m}{dg} floor} j \ &= sum_{d=1}^{min(n,m)} d sum_{g=1} mu(g) g^2 S(lfloor frac{n}{dg} floor,lfloor frac{m}{dg} floor) end{aligned} ]

其中 (S(n,m)=frac{n(n+1)m(m+1)}{4})
然后观察式子，显然可以先求出 (mu(g) g^2) 的前缀和
然后数论分块套数论分快，(O(n+n^{frac 3 4})) 完结

( ext{Code})

#include<cstdio>
#include<iostream>
#define re register 
using namespace std;
typedef long long LL;

const int N = 1e7, P = 20101009;
int n, m, totp, pr[N], vis[N + 5], sum[N + 5];

inline void Euler()
{
	vis[1] = sum[1] = 1;
	for(re int i = 2; i <= N; i++)
	{
		if (!vis[i]) pr[++totp] = i, sum[i] = -1;
		for(re int j = 1; j <= totp && i * pr[j] <= N; j++)
		{
			vis[i * pr[j]] = 1;
			if (!(i % pr[j])) break;
			sum[i * pr[j]] = -sum[i];
		}
	}
	for(re int i = 1; i <= N; i++) sum[i] = ((LL)sum[i - 1] + (LL)i * i % P * (sum[i] + P)) % P;
}

inline int S(int n, int m){return ((LL)n * (n + 1) / 2 % P) * ((LL)m * (m + 1) / 2 % P) % P;}

inline int F(int n, int m)
{
	LL res = 0;
	for(re int l = 1, r; l <= min(n, m); l = r + 1)
	{
		r = min(n / (n / l), m / (m / l));
		res = (res + (LL)(sum[r] - sum[l - 1] + P) * S(n / l, m / l) % P) % P;
	}
	return res;
}

int main()
{
	Euler();
	scanf("%d%d", &n, &m);
	LL ans = 0;
	for(re int l = 1, r; l <= min(n, m); l = r + 1)
	{
		r = min(n / (n / l), m / (m / l));
		ans = (ans + (LL)(l + r) * (r - l + 1) / 2 % P * F(n / l, m / l) % P) % P;
	}
	printf("%lld
", ans);
}