(color{#0066ff}{ 题目描述 })
对于(A_1,A_2,cdots,A_N)，求
(sum_{i=1}^Nsum_{j=1}^N lcm(A_i,A_j))的值。
(lcm(a,b)) 表示a 和b 的最小公倍数
(color{#0066ff}{输入格式})
第1 行，1 个整数N。
第2 行，N 个整数(A_1,A_2,cdots,A_N)。
(color{#0066ff}{输出格式})
1 个整数，表示所求的值。
(color{#0066ff}{输入样例})
2
2 3
(color{#0066ff}{输出样例})
17
(color{#0066ff}{数据范围与提示})
• 对于30% 的数据，(1 le N le 1000; 1 le A_i le 50000)；
• 对于另外30% 的数据，(1 le N le 50000; 1 le A_i le 1000)；
• 对于100% 的数据，(1 le N le 50000; 1 le A_i le 50000)。
(color{#0066ff}{ 题解 })
对于序列， 我们其实不太好处理
但是发现，值域只有50000
令(a_i)表示序列中等于i的数的个数
然后发现，把(a_i)弄上去后，式子变成了这样（n=50000）
[sum_{i=1}^nsum_{j=1}^n a_i*a_j*lcm(i,j)
]
这下顺眼多了，可以开始化简了
先转为gcd
[sum_{i=1}^nsum_{j=1}^n frac{a_i*a_j*i*j}{gcd(i,j)}
]
枚举gcd
[sum_{d=1}^nsum_{i=1}^nsum_{j=1}^n [gcd(i,j)==d] frac{a_i*a_j*i*j}{d}
]
把d除上去
[sum_{d=1}^nsum_{i=1}^{lfloorfrac n d
floor}sum_{j=1}^{lfloorfrac n d
floor} [gcd(i,j)==1] a_{id}*a_{jd}*i*j*d
]
把d放前面
[sum_{d=1}^n dsum_{i=1}^{lfloorfrac n d
floor}sum_{j=1}^{lfloorfrac n d
floor} [gcd(i,j)==1] a_{id}*a_{jd}*i*j
]
把(mu * 1 = e)套进去
[sum_{d=1}^n dsum_{i=1}^{lfloorfrac n d
floor}sum_{j=1}^{lfloorfrac n d
floor} sum_{k|gcd(i,j)} mu(k)* a_{id}*a_{jd}*i*j
]
枚举k
[sum_{d=1}^n dsum_{k=1}^n mu(k) sum_{i=1}^{lfloorfrac{n}{kd}
floor}sum_{j=1}^{lfloorfrac{n}{kd}
floor} a_{ikd}*a_{jkd}*i*k*j*k
]
看起来好像很复杂？ 先把k弄出来
[sum_{d=1}^n dsum_{k=1}^n mu(k)*k^2 sum_{i=1}^{lfloorfrac{n}{kd}
floor}sum_{j=1}^{lfloorfrac{n}{kd}
floor} a_{ikd}*a_{jkd}*i*j
]
然后用经典方法，kd换q
[sum_{q=1}^nsum_{k|q} mu(k)*k*qsum_{i=1}^{lfloorfrac{n}{q}
floor}sum_{j=1}^{lfloorfrac{n}{q}
floor} a_{iq}*a_{jq}*i*j
]
这。。。先把能往前提的都往前弄
[sum_{q=1}^n qsum_{k|q} mu(k)*ksum_{i=1}^{lfloorfrac{n}{q}
floor}a_{iq}*isum_{j=1}^{lfloorfrac{n}{q}
floor} a_{jq}*j
]
我去，后面两项一样诶
[sum_{q=1}^n qsum_{k|q} mu(k)*k(sum_{i=1}^{lfloorfrac{n}{q}
floor}a_{iq}*i)^2
]
然后你发现， 最后一个(sum)只跟q有关，可以枚倍数(O(nlogn))预处理
然后你又发现，中间那个也只跟q有关，也可以枚举倍数(O(nlogn))预处理
然后。。。这题还TM良心的一组数据，连数列分块都不用，直接(O(n))收集答案就行了
#include<bits/stdc++.h>
#define LL long long
LL in() {
char ch; LL x = 0, f = 1;
while(!isdigit(ch = getchar()))(ch == '-') && (f = -f);
for(x = ch ^ 48; isdigit(ch = getchar()); x = (x << 1) + (x << 3) + (ch ^ 48));
return x * f;
}
const int maxn = 50505;
const int mod = 50000;
LL a[maxn], mu[maxn], ans[maxn], tot[maxn], pri[maxn], cnt;
bool vis[maxn];
LL work() {
mu[1] = 1;
for(int i = 2; i <= mod; i++) {
if(!vis[i]) pri[++cnt] = i, mu[i] = -1;
for(int j = 1; j <= cnt&& (LL)i * pri[j] <= mod; j++) {
vis[i * pri[j]] = true;
if(i % pri[j] == 0) break;
else mu[i * pri[j]] = -mu[i];
}
}
for(int i = 1; i <= mod; i++)
for(int j = i; j <= mod; j += i) {
ans[j] += mu[i] * i;
tot[i] += a[j] * (j / i);
}
LL res = 0;
for(int i = 1; i <= mod; i++) res += 1LL * i * ans[i] * tot[i] * tot[i];
return res;
}
int main() {
int n = in();
for(int i = 1; i <= n; i++) a[in()]++;
printf("%lld
", work());
return 0;
}