【BZOJ 2154】Crash的数字表格

#include<cstdio>
#include<cstring>
#include<algorithm>
#define read(x) x=getint()
using namespace std;
typedef long long LL;
const int p = 20101009;
const int N = 1E7 + 3;
int getint() {
	int k = 0, fh = 1; char c = getchar();
	for(; c < ‘0‘ || c > ‘9‘; c = getchar())
		if (c == ‘-‘) fh = -1;
	for(; c >= ‘0‘ && c <= ‘9‘; c = getchar())
		k = k * 10 + c - ‘0‘;
	return k * fh;
}
bool np[N];
int prime[N], mu[N], s[N];
void shai(int n) {
	mu[1] = 1; s[1] = 1; int num = 0;
	for(int i = 2; i <= n; ++i) {
		if (!np[i]) {
			mu[i] = -1;
			prime[++num] = i;
		}
		for(int j = 1; j <= num; ++j) {
			int t = prime[j] * i;
			if (t > n) break;
			np[t] = 1;
			if (i % prime[j] == 0) {mu[t] = 0; break;}
			mu[t] = -mu[i];
		}
		s[i] = (s[i - 1] + (1LL * i * i * mu[i]) % p) % p;
	}
}
LL sum(LL x, LL y) {return (((x * (x + 1) / 2) % p) * ((y * (y + 1) / 2) % p)) % p;}
LL F(LL x, LL y) {
	LL re = 0;
	for(int i = 1, la = 1; i <= x; i = la + 1) {
		la = min(x / (x / i), y / (y / i));
		re = (re + (s[la] - s[i - 1]) * sum(x / i, y / i) % p) % p;
	}
	return re;
}
int main() {
	int n, m;
	read(n); read(m);
	if (n > m) swap(n, m);
	shai(n);
	
	LL ret = 0;
	for(LL d = 1, la = 1; d <= n; d = la + 1) {
		la = min(n / (n / d), m / (m / d));
		ret = (ret + (d + la) * (la - d + 1) / 2 % p * F(n / d, m / d) % p) % p;
	}
	printf("%lld", (ret + p) % p);
	
	return 0;
}