P1829 学习笔记

求

\[\sum_{i=1}^n \sum_{j=1}^m \operatorname{lcm}(i,j) \bmod 20101009 \]

写一个质数判断的程序不难发现，模数是一个质数。

我们把 \(\operatorname{lcm}\) 改成 \(\gcd\)

\[\sum_{i=1}^n \sum_{j=1}^m \frac{ij}{\gcd(i,j)} \]

令 \(g=\gcd(i,j)\) 并交换求和顺序

\[\sum_{g=1}^{n} g\sum_{i=1}^{\lfloor n/d \rfloor} \sum_{j=1}^{\lfloor m/d \rfloor} ij \cdot [\gcd(i,j)=1] \]

令

\[f(n,m)=\sum_{i=1}^{\lfloor n/d \rfloor} \sum_{j=1}^{\lfloor m/d \rfloor} ij \cdot [\gcd(i,j)=1] \]

把 \([\gcd(i,j)=1]\) 替换为莫比乌斯函数，即

\[f(n,m)=\sum_{i=1}^n \sum_{j=1}^m ij \sum_{d \mid i,d \mid j} \mu(d) \]

先枚举 \(d\)

\[f(n,m)=\sum_{d=1}^n \mu(d) d^2 \sum_{i=1}^{\lfloor n/d \rfloor} \sum_{j=1}^{\lfloor m/d \rfloor} ij \]

后面的

\[\sum_{i=1}^{\lfloor n/d \rfloor} \sum_{j=1}^{\lfloor m/d \rfloor} ij=\frac 14 \lfloor \frac nd \rfloor \lfloor \frac md \rfloor (\lfloor \frac nd \rfloor+1)(\lfloor \frac md \rfloor+1) \]

整除分块。

code

#include <bits/stdc++.h>
#define fast std::ios::sync_with_stdio(false); std::cin.tie(NULL); std::cout.tie(NULL);
#define Debug(x) std::cout << #x << ':' << x << std::endl;
constexpr std::int32_t mod = 20101009;
constexpr std::int32_t maxn = 1e7 + 5;

std::int32_t n, m, mu[maxn], prime[maxn], pre_sum[maxn], tot;
bool vis[maxn];
void sieve(){
	mu[1] = 1;
	for (std::int32_t i = 2; i <= std::min(n, m); i++){
		if (!vis[i]) prime[++tot] = i, mu[i] = -1;
		for (std::int32_t j = 1; j <= tot && i * prime[j] <= std::min(n, m); j++){
			vis[i * prime[j]] = true;
			if (i % prime[j] == 0){
				mu[i * prime[j]] = 0;
				break;
			} else mu[i * prime[j]] = -mu[i];
		}
	}
	for (std::int32_t i = 1; i <= std::min(n, m); i++) pre_sum[i] = (pre_sum[i - 1] + 1ll * i * i % mod * (mu[i] + mod) % mod) % mod;
}
std::int32_t f(std::int32_t x, std::int32_t y){
	return (1ll * x * (x + 1) / 2 % mod) * (1ll * y * (y + 1) / 2 % mod) % mod;
}
std::int32_t calc(std::int32_t x, std::int32_t y){
	std::int32_t res = 0;
	for (std::int32_t pl = 1, pr; pl <= std::min(x, y); pl = pr + 1){
		pr = std::min(x / (x / pl), y / (y / pl));
		res = (res + 1ll * (pre_sum[pr] - pre_sum[pl - 1] + mod) % mod * f(x / pl, y / pl) % mod) % mod;
	}
	return res;
}
std::int32_t solve(std::int32_t x, std::int32_t y){
	std::int32_t res = 0;
	for (std::int32_t pl = 1, pr; pl <= std::min(x, y); pl = pr + 1){
		pr = std::min(x / (x / pl), y / (y / pl));
		res = (res + 1ll * (pr - pl + 1) * (pl + pr) / 2 % mod * calc(x / pl, y / pl) % mod) % mod;
	}
	return res;
}

int main() {
	std::freopen("contest.in", "r", stdin);
	std::freopen("contest.out", "w", stdout);
	fast;
	std::cin >> n >> m;
	sieve();
	std::cout << solve(n, m);
	return 0;
}