Crash的数字表格

【题目描述】
今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数(a)和(b)，(operatorname{lcm}(a,b))表示能同时被(a)和(b)整除的最小正整数。例如，(operatorname{lcm}(6,8)) = (24)。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张(N imes M)的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为(operatorname{lcm}(i,j))。一个(4 imes 5)的表格如下：
(egin{matrix}1&2&3&4&5\2&2&6&4&10\3&6&3&12&15\4&4&12&4&20\ end{matrix})
看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和 (operatorname{mod} 20101009) 的值。

【输入格式】
输入的第一行包含两个正整数，分别表示(N)和(M)。

【输出格式】
输出一个正整数，表示表格中所有数的和 (operatorname{mod} 20101009) 的值。

(n le 10^7)
万幸20101009是个质数

题解

不妨设$nle m$

那就直接开始推式子吧。。。

懵逼警告

需要求的是 (sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}operatorname{lcm}(i,j))

(=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}frac{i*j}{gcd(i,j)})

枚举(gcd(i,j))的值

(=sumlimits_{g=1}^{n}frac{1}{g}*sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=g]*i*j)

把后半部分提出来一个(g^2)

(=sumlimits_{g=1}^{n}g*sumlimits_{i=1}^{lfloorfrac{n}{g} floor}sumlimits_{j=1}^{lfloorfrac{m}{g} floor}[gcd(i,j)=1]*i*j)

套用公式([n=1]=sumlimits_{d|n}mu(d))

(=sumlimits_{g=1}^{n}g*sumlimits_{i=1}^{lfloorfrac{n}{g} floor}sumlimits_{j=1}^{lfloorfrac{m}{g} floor}i*j*sumlimits_{d|gcd(i,j)}mu(d))

如法炮制 枚举(d)的值

(=sumlimits_{g=1}^{n}g*sumlimits_{d=1}^{lfloorfrac{n}{g} floor}mu(d)*sumlimits_{i=1}^{lfloorfrac{n}{g*d} floor}sumlimits_{j=1}^{lfloorfrac{m}{g*d} floor}d^2*i*j)

提出(d^2)

(=sumlimits_{g=1}^{n}g*sumlimits_{d=1}^{lfloorfrac{n}{g} floor}mu(d)*d^2*sumlimits_{i=1}^{lfloorfrac{n}{g*d} floor}sumlimits_{j=1}^{lfloorfrac{m}{g*d} floor}i*j)

后面的两个(sum)可以用等差数列求和的方式表示

(=sumlimits_{g=1}^{n}g*sumlimits_{d=1}^{lfloorfrac{n}{g} floor}mu(d)*d^2*frac{lfloorfrac{n}{g*d} floor*(lfloorfrac{n}{g*d} floor+1)}{2}*frac{lfloorfrac{m}{g*d} floor*(lfloorfrac{m}{g*d} floor+1)}{2})

好了 这个式子就是最终的形式了

观察一下 枚举(g) 我们只需要知道 (frac{n}{g}) 及 (frac{m}{g}) 是多少， (sumlimits_{d=1}^{lfloorfrac{n}{g} floor}mu(d)*d^2*frac{lfloorfrac{n}{g*d} floor*(lfloorfrac{n}{g*d} floor+1)}{2}*frac{lfloorfrac{m}{g*d} floor*(lfloorfrac{m}{g*d} floor+1)}{2}) 这一串就能算出来

所以用数论分块来枚举(g)就可以了 记(lfloorfrac{n}{g} floor =n_0,lfloorfrac{m}{g} floor =m_0)

同理 只需要知道(frac{n_0}{d})和(frac{m_0}{d})， 里面那串 (frac{lfloorfrac{n}{g*d} floor*(lfloorfrac{n}{g*d} floor+1)}{2}*frac{lfloorfrac{m}{g*d} floor*(lfloorfrac{m}{g*d} floor+1)}{2}) 也就能算出来了

所以(d)也可以用数论分块枚举 处理一下 (mu(d)*d^2) 这个东西的前缀和就可以了

时间复杂度是(O(sqrt{n}^2)=O(n))，很完美

注意处理负数 (mu(d)*d^2)的前缀和会有负数

推导很冗长 代码很简洁

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

inline ll read() {
	ll x = 0, f = 1; char ch = getchar();
	for (; ch > '9' || ch < '0'; ch = getchar()) if (ch == '-') f = -1;
	for (; ch <= '9' && ch >= '0'; ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ '0');
	return x * f;
}

const ll mod = 20101009;
const ll div2 = mod / 2 + 1;
ll n, m, ans;

bool np[10000005];
ll prime[1000005], mob[10000005], sum[10000005], tot;

inline void getprime() {
	mob[1] = 1;
	for (ll i = 2; i <= max(n, m); i++) {
		if (!np[i]) prime[++tot] = i, mob[i] = -1;
		for (ll j = 1; j <= tot && i * prime[j] <= max(n, m); j++) {
			np[i * prime[j]] = 1;
			if (i % prime[j] == 0) {
				mob[i * prime[j]] = 0;
				break;
			} else mob[i * prime[j]] = -mob[i];
		}
	}
	for (ll i = 1; i <= max(n, m); i++) {
		sum[i] = (sum[i-1] + i * i % mod * mob[i] % mod) % mod;
	}
}

inline ll calc(ll n0, ll m0) {
	ll ret = 0;
	for (ll l = 1, r; l <= min(n0, m0); l = r + 1) {
		r = min(n0 / (n0 / l), m0 / (m0 / l));
		ll s = (sum[r] - sum[l-1] + mod) % mod;
		ll s1 = (n0 / l) * (n0 / l + 1) % mod * div2 % mod, s2 = (m0 / l) * (m0 / l + 1) % mod * div2 % mod;
		ret = ((ret + s * s1 % mod * s2 % mod) % mod + mod) % mod;
	}	
	return ret;
}

int main() {
	n = read(); m = read();
	getprime(); 
	for (ll l = 1, r; l <= min(n, m); l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		ans = (((r - l + 1) * (r + l) % mod * div2 % mod * calc(n / l, m / l) % mod + ans) % mod + mod) % mod;
	}
	printf("%lld
", ans); 
	return 0; 
}