AT_agc038_c 做题笔记

题目链接

莫反好题，不仅仅是莫反，还有很多思维含量。

由于推式子过程太过于漫长了，所以我仅仅讲下大概。

题目是给你一个长度为 $n$ 的数组，请求出 $\sum\limits_{i=1}^n\sum\limits_{j=i+1}^n \operatorname{lcm}(A_i, A_j)$

莫反通常是对于值域考虑，直接推是不可行的，所以开一个桶 $b$，$b_k$ 存储 $A_i=k$ 的个数。

对值域分类的话就只能把忽略大小关系的 $i$ $j$ 对算出来后再容斥得到答案，即：

$\sum\limits_{i=1}^V\sum\limits_{j=1}^V \operatorname{lcm}(i,j)\cdot b_i \cdot b_j$。

多说一句，我刚开始没忽略大小关系，是 $\sum\limits_{i=1}^V\sum\limits_{j=i+1}^V \operatorname{lcm}(i,j)\cdot b_i \cdot b_j$，但这个实际上是错的，我当时没注意，推到最后发现没办法推了才想到这个办法的。。。

推式子，大概就是先把 $\operatorname{lcm}$ 变成相乘除以最大公约数之后枚举最大公约数，然后再同时除以最大公约数，然后再用莫比乌斯反演之后再枚举。

此时有一个整除的条件，枚举的时候枚举倍数就行了，最后把所有 $\sum$ 内的能拿出来的都拿出来，用乘法分配律拿。

推完后是：$\sum\limits_{d=1}^V (d\cdot \sum\limits_{x=1}^{\lfloor\frac{V}{d}\rfloor}\cdot (\mu(x)\cdot x^2\cdot\sum\limits_{i=1}^{\lfloor\frac{V}{dx}\rfloor}i\cdot b_{i\cdot d\cdot x}\sum\limits_{j=1}^{\lfloor\frac{V}{dx}\rfloor}j\cdot b_{j\cdot d\cdot x})$。

此时令 $f(x)=\sum\limits_{i=1}^{\lfloor\frac{V}{x}\rfloor}i\cdot b_{i\cdot x}$，然后就可以化为 $\sum\limits_{d=1}^V (d\cdot \sum\limits_{x=1}^{\lfloor\frac{V}{d}\rfloor} (\mu(x)\cdot x^2\cdot f^2(dx)))$。

对 $f$ 进行预处理，然后调和级数 $n\log n$ 就可求得，注意最后要除以二，用逆元，还有减去 $\sum_{i=1}^n a_i$，因为我们多算了要容斥掉。

代码：

#include <bits/stdc++.h>
#define int long long
#define For(i, a, b) for (int i = (a); i <= (b); i ++)
#define foR(i, a, b) for (int i = (a); i >= (b); i --)
using namespace std;
int n, x, ans, cnt;
int b[1000005], f[1000005];
int mu[1000005], primes[200005];
bool prime[1000005];
const int mod = 998244353;
void init () {
    mu[1] = 1;
    For (i, 2, 1000000) {
        if (!prime[i]) {
            primes[++ cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * primes[j] <= 1000000; j ++) {
            prime[i * primes[j] ] = 1;
            if (i % primes[j] == 0) break;
            mu[i * primes[j] ] = -mu[i];
        }
    }
}
void solve () {
    init ();
    cin >> n;
    For (i, 1, n) {
        cin >> x;
        ans -= x;
        ++ b[x];
    }
    ans = (ans % mod + mod) % mod;
    For (i, 1, 1000000) {
        For (j, 1, 1000000 / i) f[i] += j * b[i * j];
        f[i] %= mod;
    }
    For (i, 1, 1000000) {
        int x = 1000000 / i, res = 0;
        For (j, 1, x) {
            res += mu[j] * j * j % mod * f[i * j] % mod * f[i * j] % mod;
            if (res >= mod) res -= mod;
            if (res <= mod) res += mod;
        }
        ans += i * res % mod;
        if (ans >= mod) ans -= mod;
    }
    cout << ans * 499122177 % mod;
}
signed main () {
    int _ = 1;
//    cin >> _;
    while (_ --) {
        solve ();
        cout << '\n';
    }
    return 0;
}