\(\text{Problem}\)

\[\left( \sum_{i=1}^n \sum_{j=1}^m lcm(i，j) \right) \bmod{20101009} \]

\(n \le 10^7\)

\(\text{Analysis}\)

套路地推式子

\[\begin{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 \rfloor} \sum_{j=1}^{\lfloor \frac m d \rfloor} ijd[\gcd(i,j)=1] \\ &= \sum_{d=1}^{\min(n,m)} d \sum_{i=1}^{\lfloor \frac n d \rfloor} i \sum_{j=1}^{\lfloor \frac m d \rfloor} 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 \rfloor} \sum_{g|j}^{\lfloor \frac m d \rfloor} ij \\ &= \sum_{d=1}^{\min(n,m)} d \sum_{g=1} \mu(g) g^2 \sum_{i=1}^{\lfloor \frac{n}{dg} \rfloor} i \sum_{j=1}^{\lfloor \frac{m}{dg} \rfloor} j \\ &= \sum_{d=1}^{\min(n,m)} d \sum_{g=1} \mu(g) g^2 S(\lfloor \frac{n}{dg} \rfloor,\lfloor \frac{m}{dg} \rfloor) \end{aligned} \]

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

\(\text{Code}\)

#include
#include
#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\n", ans);
}