首先设 \(lim = \max\limits_{i = 1}^N a_i,c_k = \sum\limits_{i = 1}^N [a_i = k]\)。
于是问题转化为求解 \(\sum\limits_{i = 1}^{lim} \sum\limits_{j = 1}^{lim} \text{lcm}(i,j) c_i c_j\)。
\[\begin{align*} & \sum\limits_{i = 1}^{lim} \sum\limits_{j = 1}^{lim} \text{lcm}(i,j) c_i c_j \\ = & \sum\limits_{i = 1}^{lim} \sum\limits_{j = 1}^{lim} \dfrac{ijc_ic_j}{\gcd(i,j)} \\ = & \sum\limits_{d = 1}^{lim} d \sum\limits_{id \le lim} \sum\limits_{jd \le lim} [\gcd(i,j) = 1]ijc_{id}c_{jd} \end{align*}\]
设 \(f(x) = \sum\limits_{id \le lim} \sum\limits_{jd \le lim} [\gcd(i,j) = x]ijc_{id}c_{jd},F(x) = \sum\limits_{x | d} f(d) = \sum\limits_{id \le lim} \sum\limits_{jd \le lim} [x | \gcd(i,j)]ijc_{id}c_{jd}\)。
\[\begin{align*} F(x) = & \sum\limits_{id \le lim} \sum\limits_{jd \le lim} [x | \gcd(i,j)]ijc_{id}c_{jd} \\ = & x^2 \sum\limits_{xid \le lim} \sum\limits_{xjd \le lim} ijc_{ixd}c_{jxd} \\ = & x^2 \left(\sum\limits_{xid \le lim} ic_{ixd}\right)^2 \end{align*}\]
\[\begin{align*} & \sum\limits_{d = 1}^{lim} d \sum\limits_{id \le lim} \sum\limits_{jd \le lim} [\gcd(i,j) = 1]ijc_{id}c_{jd} \\ = & \sum\limits_{d = 1}^{lim} d f(1) \\ = & \sum\limits_{d = 1}^{lim} d \sum\limits_{id \le lim} i^2 \mu(i) \left(\sum\limits_{ijd \le lim} jc_{ijd}\right)^2 \end{align*}\]
于是可以预处理 \(\left(\sum\limits_{xd \le lim} jc_{xd}\right)^2\) 然后这题就没了。
代码： 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
using namespace std;
const int N = 5e4;
int n,lim,c[N + 5];
int vis[N + 5],prime[N + 5],cnt,mu[N + 5];
long long f[N + 5];
long long ans;
int main()
{
mu[1] = 1;
for(register int i = 2;i <= N;++i)
{
if(!vis[i])
mu[prime[++cnt] = i] = -1;
for(register int j = 1;j <= cnt && i * prime[j] <= N;++j)
{
vis[i * prime[j]] = 1;
if(!(i % prime[j]))
break;
else
mu[i * prime[j]] = -mu[i];
}
}
scanf("%d",&n);
int x;
for(register int i = 1;i <= n;++i)
scanf("%d",&x),++c[x],lim = max(lim,x);
for(register int i = 1;i <= lim;f[i] *= f[i],++i)
for(register int j = 1;i * j <= lim;++j)
f[i] += (long long)j * c[i * j];
for(register int i = 1;i <= lim;++i)
for(register int j = 1;i * j <= lim;++j)
ans += (long long)i * j * j * mu[j] * f[i * j];
printf("%lld\n",ans);
}