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

题意

求这个(sumlimits_{i = 1}^{n}sumlimits_{j = 1}^{m}lcm(i,j) mod (20101009))

想法

以下(n > m)
在物理课上推出了柿子（逃
回家了拍草稿纸上来
(sumlimits_{i = 1}^{n}sumlimits_{j = 1}^{m}lcm(i,j))
不想写一大堆的(LaTex)下面写出的等式都是连等的

(sumlimits_{i = 1}^{n}sumlimits_{j = 1}^{m}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}sumlimits_{i = 1}^{lfloorfrac{n}{d} floor}sumlimits_{j = 1}^{lfloorfrac{m}{d} floor}i * j * d * [gcd(i,j) == 1])
再用莫反的套路
(sumlimits_{d = 1}^{n}sumlimits_{i = 1}^{lfloorfrac{n}{d} floor}sumlimits_{j = 1}^{lfloorfrac{m}{d} floor}sumlimits_{s|gcd(i,j)}i * j * d*mu(s))
分类一下
(sumlimits_{d = 1}^{n}dsumlimits_{i = 1}^{lfloorfrac{n}{d} floor}sumlimits_{j = 1}^{lfloorfrac{m}{d} floor}sumlimits_{s|gcd(i,j)}i * j * mu(s))
再把后三个求和改为枚举(s)
(sumlimits_{d = 1}^{n}dsumlimits_{s}^{lfloorfrac{n}{d} floor}sumlimits_{i = 1}^{lfloorfrac{n}{ds} floor}sumlimits_{j = 1}^{lfloorfrac{m}{ds} floor}i * j * mu(s) * s^2)
(sumlimits_{d = 1}^{n}dsumlimits_{s}^{lfloorfrac{n}{d} floor}mu(s) * s^2sumlimits_{i = 1}^{lfloorfrac{n}{ds} floor}isumlimits_{j = 1}^{lfloorfrac{m}{ds} floor}j)
后面两个求和用等差公式
此时这个柿子是(O(n^2))的
接下来的步骤为了降低复杂度
(sumlimits_{d = 1}^{n}d * ssumlimits_{s}^{lfloorfrac{n}{d} floor}mu(s) * ssumlimits_{i = 1}^{lfloorfrac{n}{ds} floor}isumlimits_{j = 1}^{lfloorfrac{m}{ds} floor}j)
改枚举(T = d * s)
(sumlimits_{T = 1}^{n}Tsumlimits_{i = 1}^{lfloorfrac{n}{T} floor}isumlimits_{j = 1}^{lfloorfrac{m}{T} floor}jsumlimits_{s|T}s*mu(s))
这样中间两个求和是(O(1))最后那个求和我们可以在筛(mu)时一并筛出(mu(s)*s)然后用狄利克雷前缀和
复杂度(O(n + N))

代码

#include<iostream>
#include<cstdio>
#define ll long long

ll n,m,mod = 20101009;
ll mu[10000010],p[10000010],f[10000010],tot;
bool v[10000010];

int N = 1e7,cnt;

void pre(){
	f[1] = mu[1] = 1;
	for(int i = 2;i <= 1e7;++i){
		if(!v[i]) mu[i] = -1,p[++tot] = i;
		f[i] = (mu[i] * i + mod) % mod;
		for(int j = 1;j <= tot && i * p[j] <= 1e7;++j){
			v[p[j] * i] = 1;
			if(i % p[j] == 0){
				mu[i * p[j]] = 0;
				break;
			}
			mu[i * p[j]] = -mu[i];
		}
	}
	f[1] = 1;
	for(int j = 1;j <= tot;++j)
	for(int i = 1;i * p[j] <= 1e7;++i)
	(f[i * p[j]] += f[i]) %= mod;
}

ll ans = 0;

ll sum(int x){
	ll ans = 1ll * x * f[x] % mod;
	ans = ans * 1ll * ((n / x + 1) * (n / x) / 2 % mod) % mod;
	ans = ans * 1ll * ((m / x + 1) * (m / x) / 2 % mod) % mod;
	return ans % mod;
}

int main(){
	scanf("%lld%lld",&n,&m);
	if(m > n)
	std::swap(m,n);
	pre();
	for(int i = 1;i <= n;++i)
	ans = (ans + sum(i)) % mod;
	std::cout<<ans<<std::endl;
}