转化下题面。
求
[sum_{i=1}^n sum_{j=1}^n lcm(i,j) * cnt_i * cnt_j
]
其中， (n) 、 (cnt_k (1 leq k leq n)) 已事先给出规模均为 (5e4) 。
[sum_{i=1}^n sum_{j=1}^n lcm(i,j) * cnt_i * cnt_j
]
[= sum_{i=1}^n sum_{j=1}^n frac{i*j}{gcd(i,j)} * cnt_i * cnt_j
]
[= sum_{d=1}^n sum_{i=1}^n sum_{j=1}^n [gcd(i,j)==d] * frac{i*j}{d} * cnt_i * cnt_j
]
[= sum_{d=1}^n sum_{i=1}^{lfloor frac{n}d
floor} sum_{j=1}^{lfloor frac{n}d
floor} [gcd(i,j)==1] * i * j * d * cnt_{id} * cnt_{jd}
]
[= sum_{d=1}^n sum_{i=1}^{lfloor frac{n}d
floor} sum_{j=1}^{lfloor frac{n}d
floor} sum_{s|gcd(i,j)} mu(s) * i * j * d * cnt_{id} * cnt_{jd}
]
[= sum_{d=1}^n d*sum_{i=1}^{lfloor frac{n}d
floor} sum_{j=1}^{lfloor frac{n}d
floor} sum_{s|gcd(i,j)} mu(s) * i * j * cnt_{id} * cnt_{jd}
]
[= sum_{d=1}^n d* sum_{s=1}^{lfloor frac{n}d
floor} mu(s) sum_{i=1}^{lfloor frac{n}{sd}
floor} sum_{j=1}^{lfloor frac{n}{sd}
floor} s^2 * i * j * cnt_{ids} * cnt_{jds}
]
[= sum_{d=1}^n d* sum_{s=1}^{lfloor frac{n}d
floor} mu(s) * s^2 sum_{i=1}^{lfloor frac{n}{sd}
floor} sum_{j=1}^{lfloor frac{n}{sd}
floor} i * j * cnt_{ids} * cnt_{jds}
]
[= sum_{d=1}^n d* sum_{s=1}^{lfloor frac{n}d
floor} mu(s) * s^2 sum_{i=1}^{lfloor frac{n}{sd}
floor} i*cnt_{ids} sum_{j=1}^{lfloor frac{n}{sd}
floor} j * cnt_{jds}
]
[= sum_{d=1}^n d* sum_{s=1}^{lfloor frac{n}d
floor} mu(s) * s^2 Big( sum_{i=1}^{lfloor frac{n}{sd}
floor} i*cnt_{ids} Big)^2
]
[= sum_{sd=1}^n sd sum_{s|sd} ig( mu(s) * s ig) Big( sum_{i=1}^{lfloor frac{n}{sd}
floor} i*cnt_{ids} Big)^2
]
把 (sd) 换成 一个漂亮点的数 (T) ， 最后的式子就是
[sum_{T=1}^n T *ig( sum_{s|T} mu(s) * s ig) * Big( sum_{i=1}^{lfloor frac{n}{T}
floor} i*cnt_{iT} Big)^2
]
中间的括号可以 (O(n log n)) 预处理， 后面的括号暴力算的总复杂度是个调和级数， 复杂度也是 (O(n log n))。
luogu数据AC代码：
#include<bits/stdc++.h>
using namespace std;
const int maxn = 50005;
#define li long long
int prime[maxn], v[maxn], mu[maxn], m;
li f[maxn];
void euler(int n) {
mu[1] = 1;
for(int i=2; i<=n; ++i) {
if(!v[i]) {
v[prime[++m] = i] = i;
mu[i] = -1;
}
for(int j=1; j<=m; ++j) {
if(prime[j] > n/i || prime[j] > v[i]) break;
v[prime[j] * i] = prime[j];
mu[prime[j] * i] = mu[i] * (i%prime[j] ? -1 : 0);
}
}
for(int i=1; i<=n; ++i)
for(int j=1; j<=n/i; ++j)
f[i*j] += i * mu[i];
}
int n, a[maxn], c[maxn];
int N;
int main()
{
scanf("%d", &n);
for(int i=1;i<=n;++i) {
scanf("%d", &a[i]);
++c[a[i]];
N = max(N, a[i]);
}
euler(N);
li ans = 0ll;
for(int T=1; T<=N; ++T) {
li nowans = 0ll;
for(int i=1; i<=N/T; ++i) nowans += i*c[i*T];
nowans = nowans * nowans;
ans += T * f[T] * nowans;
}
cout << ans;
return 0;
}