令 \(s_n\) 表示最大公因数恰好为 \(n\) 的集合个数，\(f_n\) 表示与 \(n\) 互质的数的个数。
那么我们要求的就是 \(\sum_i s_if_i\)。
考虑求出 \(s\) 和 \(f\)。
考虑计算 \(f_i\)。令 \(c_i\) 表示数 \(i\) 的出现次数。
则直接莫比乌斯反演：
\[f_n=\sum_{i}[\gcd(i,n)=1]c_i \]\[=\sum_{i}c_i\sum_{d\mid i,d\mid n}\mu(d) \]\[=\sum_{d\mid n}\mu(d)\sum_{d\mid i}c_i \]令 \(g_k=\sum\limits_{k\mid i}c_i,g_{k}^{′}=\mu(k)g_k\)，那么 \(f_n=\sum\limits_{d\mid n}g_{d}^{′}\)
发现这个 \(g_i\) 就表示 \(i\) 的倍数的出现次数。
再考虑求 \(s_i\)。我们令 \(s_i^′\) 表示最大公因数为 \(i\) 的倍数的方案数。要使最大公因数为 \(i\) 的倍数，那么每个选的数必须是 \(i\) 的倍数。所以容易得到 \(s_i^′=2^{g_i}-1\)。
又有 \(s_i^{′}=\sum\limits_{i\mid d}s_d\)。
观察上面的式子的形式，它们都是这样的形式：
\[b_k=\sum_{i\mid k}a_i \]\[b_k=\sum_{k\mid i}a_i \]其中，上面的式子为狄利克雷前缀和，下面的式子为狄利克雷后缀和。
这都可以在 \(\mathcal O(w\log\log w)\) 的时间内求出，可以先去做一下【模板】Dirichlet 前缀和（这个实际上就是一个关于质因子分解后的指数的高维前缀和，使用类似 FMT 的方法就可以了）。
计算 \(g\) 和 \(f\) 都是标准的狄利克雷前/后缀和形式。
但还有个问题，已知 \(a\) 求 \(b\) 能做了，但是那个 \(s\) 的计算，是已知 \(b\) 求 \(a\) 的。
前面是加，那现在就把它改成减就可以了。当然也可以用高维前缀和的知识解释。
总时间复杂度 \(\mathcal O(n+w\log\log w)\)。
具体细节看代码。
Code：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 500005, M = 10000000, mod = 1e9 + 7;
int n, ans;
int a[N], cnt[M+5];
int mu[M+5];
int pri[M / 10], tot;
bool vis[M+5];
int pw[N];
int f[M], g[M];
inline void add(int &a, int b) {
a += b;
if (a >= mod) a -= mod;
if (a < 0) a += mod;
}
void init(int n) {
pw[0] = 1;
for (int i = 1; i <= 500000; ++i) pw[i] = 1ll * pw[i - 1] * 2 % mod;
mu[1] = 1;
for (int i = 2; i <= n; ++i) {
if (!vis[i]) pri[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot && i * pri[j] <= n; ++j) {
vis[i * pri[j]] = 1;
if (i % pri[j] == 0) {
mu[i * pri[j]] = 0;
break;
}
mu[i * pri[j]] = -mu[i];
}
}
}
int main() {
init(M);
scanf("%d", &n);
for (int i = 1; i <= n; ++i) scanf("%d", &a[i]), ++cnt[a[i]];
for (int i = 1; i <= tot; ++i)
for (int j = M / pri[i]; j; --j)
add(cnt[j], cnt[j * pri[i]]);
for (int i = 1; i <= M; ++i) g[i] = cnt[i] * mu[i];
for (int i = 1; i <= tot; ++i)
for (int j = 1; j * pri[i] <= M; ++j)
add(g[j * pri[i]], g[j]);
for (int i = 1; i <= M; ++i) f[i] = pw[cnt[i]] - 1;
for (int i = 1; i <= tot; ++i)
for (int j = 1; j * pri[i] <= M; ++j)
add(f[j], -f[j * pri[i]]);
for (int i = 2; i <= M; ++i) add(ans, 1ll * f[i] * g[i] % mod);
printf("%d", ans);
return 0;
}