2024 年春节集训 _ 第三课 - 莫比乌斯反演

练习 \(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