练习 \(5\) \(\color{orange}{\texttt{E -> link}}\)
求
\[\sum _{i=1} ^n \sum _{j=1} ^m lcm(i,j)
\]
\(n,m \leq 10^7 ,\ T\leq 10^4\)
贴个照片。
及其丑陋的照片（我的草稿）
如上化简最后可以得到
\[\sum _{d=1} ^n d\sum_{k=1} ^{\left[\dfrac{n}{d}\right]} \mu(k)k^2 \sum_{i=1} ^{\left[\dfrac{n}{dk}\right]} i \sum_{j=1} ^{\left[\dfrac{m}{dk}\right]} j
\]
好了然后注意到后面的是等差序列，可以丢尽分块处理，所以考虑 设 \(T=dk, f(x)=\dfrac{x(x+1)}{2}\)
代入进去得到：
\[\sum _{d=1} ^n d\sum_{k=1} ^{\left[\dfrac{n}{d}\right]} \mu(k)k^2 f(\left[ \dfrac{n}{T}\right])f(\left[ \dfrac{m}{T}\right]) \
\]
把后面两个 \(f\) 往前面扔。枚举 \(T.\) 就像之前 \(cbh\) 同学往前枚举 \(pd\) 是同样的道理。得到
\[\sum_{T=1}^{n} f(\left[ \dfrac{n}{T}\right])f(\left[ \dfrac{m}{T}\right]) \sum_{d|T} d\mu(d) T
\]
设 \(F(T) = \sum_{d|T} d\mu(d) T.\)
那么原式就可以直接暴力二维整除分块了。
至于 \(F\) 函数在 \(\texttt{Euler}\) 筛法中处理。详见程式。
程式
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const ll N = 1e7 + 5, mod = 1e8 + 9;
int f[N], prime[N / 10];
ll s[N];
bitset<N> vis;
int top, t, n, m;
inline void sieve(void) {
vis[1] = s[1] = f[1] = 1;
for (int i = 2; i < N; ++i) {
if (!vis[i])
vis[i] = true, prime[++top] = i, f[i] = 1 - i;
for (int j = 1; j <= top; ++j) {
if (1ll * i * prime[j] > N)
break;
vis[1ll * i * prime[j]] = 1;
if (!(1ll * i % prime[j])) {
f[1ll * i * prime[j]] = f[i];
break;
}
f[1ll * i * prime[j]] = 1ll * f[i] * f[prime[j]] % mod;
}
s[i] = s[i - 1] + 1ll * f[i] * i, f[i] %= mod, s[i] %= mod;
}
return;
}
inline ll func(int x) { return 1ll * x * (x + 1) / 2 % mod; }
inline ll calc(void) {
ll lt = 1, rt = 1, ans = 0;
for (lt = 1; lt <= min(n, m); lt = rt + 1) {
rt = min(n / (n / lt), m / (m / lt));
(ans += 1ll * func(n / lt) % mod * func(m / lt) % mod * (s[rt] - s[lt - 1]) % mod + mod) % mod;
}
return ans % mod;
}
int main() {
scanf("%d", &t);
sieve();
while (t--) {
scanf("%d%d", &n, &m);
if (n > m)
swap(n, m);
printf("%lld\n", calc());
}
return 0;
}
// 1 29 47
