Luogu 3911 最小公倍数之和

啊，如果你愿意的话可以去 hydro 的每日一题查看这篇文章 QWQ。

#038.2022.7.17

Luogu P3911 最小公倍数之和

标签：莫比乌斯反演，数论

Statement

给定 \(N\) 个数 \(A_1,A_2,\cdots,A_N\)，求 \(\sum_{i=1}^N\sum_{j=1}^N \mathrm{lcm}(A_i,A_j)\)
的值。

\(\textbf{Data Range :} 1 \le N \le 5\times 10^4\)，\(1 \le A_i \le 5\times 10^4\)。

Solution

solution

肯定还是要考虑莫反，但是题目中给定了 \(A_i\)，这就不能常规的进行反演。

但是其实没有影响的，我们记录每个数 \(x\) 的出现次数 \(c_x\)，然后原式子等价于下面的式子：

\[\sum_{i=1}^N\sum_{j=1}^N \text{lcm}(i,j) \times c_i \times c_j \]

于是可以进行反演了。

\[\sum_{i=1}^N\sum_{j=1}^N \text{lcm}(i,j) \times c_i \times c_j \]

\[\sum_{i=1}^N\sum_{j=1}^N \dfrac{i \times j}{\gcd(i,j)} \times c_i \times c_j \]

然后我们套路的改为枚举 \(\gcd\)。

\[\sum_{d=1}^N \sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{\frac{N}{d}} \dfrac{id \times jd}{d} \times c_{id} \times c_{jd} [(i,j)= 1] \]

\[\sum_{d=1}^N \sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{\frac{N}{d}} ijd \times c_{id} \times c_{jd} [(i,j)= 1] \]

\[\sum_{d=1}^N \sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{\frac{N}{d}} \sum_{x|(i,j)} \mu(x) \times ijd \times c_{id} \times c_{jd} \]

\[\sum_{d=1}^N \sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{\frac{N}{d}} \sum_{x|i,x | j} \mu(x) \times ijd \times c_{id} \times c_{jd} \]

之后枚举 \(x\)。

\[\sum_{d=1}^N \sum_{x=1}^{\frac{N}{d}}\mu(x) \sum_{i=1}^{\frac{N}{dx}} \sum_{j=1}^{\frac{N}{dx}} ix \times jx \times d \times c_{idx} \times c_{jdx} \]

\[\sum_{d=1}^N \sum_{x=1}^{\frac{N}{d}}\mu(x) x^2 \sum_{i=1}^{\frac{N}{dx}} \sum_{j=1}^{\frac{N}{dx}} ijd\times c_{idx} \times c_{jdx} \]

\[\sum_{d=1}^N d\sum_{x=1}^{\frac{N}{d}}\mu(x) x^2 \sum_{i=1}^{\frac{N}{dx}} \sum_{j=1}^{\frac{N}{dx}} ij\times c_{idx} \times c_{jdx} \]

令 \(T = xd\)，代入进行化简。

\[\sum_{T=1}^N \sum_{d|T} d \times \mu(\dfrac{T}{d}) \times (\dfrac{T}{d})^2 \sum_{i=1}^{\frac{N}{T}} \sum_{j=1}^{\frac{N}{T}} ij \times c_{\frac{idT}{d}} \times c_{\frac{jdT}{d}} \]

\[\sum_{T=1}^N (\sum_{d|T} d \times \mu(\dfrac{T}{d}) \times (\dfrac{T}{d})^2) \times (\sum_{i=1}^{\frac{N}{T}} \sum_{j=1}^{\frac{N}{T}} ij \times c_{\frac{idT}{d}} \times c_{\frac{jdT}{d}}) \]

对于 \(\sum_{d|T} d \times \mu(\dfrac{T}{d}) \times (\dfrac{T}{d})^2\) 这一部分，他等价于 \(T \sum_{d|T} \times \mu(\dfrac{T}{d}) \times (\dfrac{T}{d})\)，于是我们可以在进行线性筛之后暴力枚举倍数预处理就好了，复杂度是对的。

所以我们只用考虑后面一部分怎么做，也就是 \(\sum_{i=1}^{\frac{N}{T}} \sum_{j=1}^{\frac{N}{T}} ij \times c_{\frac{idT}{d}} \times c_{\frac{jdT}{d}}\)。

我们发现这东西其实不用写成 \(ij\) 两层循环，完全可以写成 \((\sum_{i=1}^{\frac{N}{T}} i \times c_{\frac{idT}{d}})^2 = (\sum_{i=1}^{\frac{N}{T}} i \times c_{iT})^2\)。

和上面一样的，莫反题到最后优化也就是整除分块和预处理了。

我们对每个 \(T\) 暴力预处理 \(\sum_{i=1}^{\frac{N}{T}} i \times c_{iT}\)，根据经典结论 \(\sum_{i=1}^N \dfrac{N}{i} = n \ln n\)。

然后枚举倍数的复杂度也是这样的，所以之前的预处理复杂度也是 \(n\ln n\) 的。

于是总预处理复杂度为 \(n \ln n\)。

最后我们直接暴力回答答案就好了，统计答案的复杂度为 \(O(n)\)。

// 德丽莎你好可爱德丽莎你好可爱德丽莎你好可爱德丽莎你好可爱德丽莎你好可爱
// 德丽莎的可爱在于德丽莎很可爱，德丽莎为什么很可爱呢，这是因为德丽莎很可爱！
// 没有力量的理想是戏言，没有理想的力量是空虚
#include <bits/stdc++.h>
#define LL long long
#define int long long
using namespace std;
char ibuf[1 << 15], *p1, *p2;
#define getchar() (p1 == p2 && (p2 = (p1 = ibuf) + fread(ibuf, 1, 1 << 15, stdin), p1==p2) ? EOF : *p1++)
inline int read() {
  char ch = getchar();  int x = 0, f = 1;
  while (ch < '0' || ch > '9')  {  if (ch == '-')  f = -1;  ch = getchar();  }
  while (ch >= '0' && ch <= '9')  x = (x << 1) + (x << 3) + (ch ^ 48), ch = getchar();
  return x * f;
}
void print(LL x) {
  if (x > 9)  print(x / 10);
  putchar(x % 10 + '0');
}
template<class T> bool chkmin(T &a, T b) { return a > b ? (a = b, true) : false; }
template<class T> bool chkmax(T &a, T b) { return a < b ? (a = b, true) : false; }
#define rep(i, l, r) for (int i = (l); i <= (r); i++)
#define repd(i, l, r) for (int i = (l); i >= (r); i--)
#define REP(i, l, r)  for (int i = (l); i < (r); i++)
const int N = 5e4 + 6, M = 5e4 + 5;
int n, vis[N], c[N], pr[N], tot, mu[N];
void prime(int n) {
  mu[1] = 1;
  rep (i, 2, n) {
    if (!vis[i]) {  pr[++tot] = i;  mu[i] = -1;  }
    rep (j, 1, tot) {
      if (pr[j] * i > n)  break;
      vis[i * pr[j]] = 1;
      if (i % pr[j] == 0) {
        mu[i * pr[j]] = 0;  break;
      }  else {
        mu[i * pr[j]] = -mu[i];
      }
    }
  }
}
int f[N], g[N];
void solve() {
  n = read();
  rep (i, 1, n) {  int x = read();  c[x] ++;  }
  prime(M);
  rep (i, 1, M) 
    for (int j = i; j <= M; j += i)  f[j] += mu[i] * i;
  rep (i, 1, M)  f[i] *= i;
  rep (T, 1, M)
    for (int i = 1; i <= M / T; i++)  g[T] += i * c[i * T];
  int ans = 0;
  rep (T, 1, M)  ans += f[T] * g[T] * g[T];
  cout << ans;
}
signed main () {
#ifdef LOCAL_DEFINE
  freopen("1.in", "r", stdin);
  freopen("1.ans", "w", stdout);
#endif
  int T = 1;  while (T--)  solve();
#ifdef LOCAL_DEFINE
  cerr << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC << " s.\n";
#endif
  return 0;
}