$$\sum_{L\leq i<j\leq R}gcd(a[i],a[j])$$

$$\sum_{L\leq i<j\leq R}\sum_{d|gcd(a[i],a[j])}\phi(d)$$

$$\sum_{d=1}^n \sum_{L\leq i<j\leq R} \phi(d)C_{cnt[d]}^2$$

#include<bits/stdc++.h>
#define fo(i, a, b) for (Re int i = (a); i <= (b); ++i)
#define fd(i, a, b) for (Re int i = (a); i >= (b); --i)
#define edge(i, u) for (Re int i = head[u], v = e[i].v; i; i = e[i].nxt, v = e[i].v)
#define Re register
#define pb push_back
#define F first
#define S second
#define ll unsigned long long
#define lowbit(x) (x & -x)
#define mod 1000000007
#define N 20005
int T, n, a[N], m, block[N], cnt[N], l, r;
ll anss[N], ans;
std::vector<int> d[N];
struct ques{
int x, y, id;
friend bool operator < (ques x, ques y)
{
return block[x.x] == block[y.x] ? x.y < y.y: block[x.x] < block[y.x];
}
}q[N];
std::bitset<N> vis;
int p[N], tot, phi[N];
inline void init ()
{
int n = N - 3;
phi[1] = 1;
fo (i, 2, n)
{
if (!vis[i])
{
p[++tot] = i;
phi[i] = i - 1;
}
for (int j = 1; j <= tot && i * p[j] <= n; ++j)
{
vis[i * p[j]] = 1;
if (!(i % p[j]))
{
phi[i * p[j]] = phi[i] * p[j];
break;
}
phi[i * p[j]] = phi[i] * phi[p[j]];
}
}
}
inline void add (int x, bool opt)
{
if (opt)
{
int sz = d[x].size();
for (int i = 0; i < sz; ++i)
{
int now = d[x][i];
ans += cnt[now] * phi[now];
++cnt[now];
}
}
else
{
int sz = d[x].size();
for (int i = 0; i < sz; ++i)
{
int now = d[x][i];
--cnt[now];
ans -= cnt[now] * phi[now];
}
}
}
int main ()
{
init();
int up = 20000;
fo (i, 1, up)
{
int sq = (int) sqrt(i - 0.5);
fo (j, 1, sq)
if (!(i % j))
{
d[i].pb(j);
d[i].pb(i / j);
}
++sq;
if (sq * sq == i) d[i].pb(sq);
}
scanf("%d", &T);
fo (qwq, 1, T)
{
memset(cnt, 0, sizeof cnt);
scanf("%d", &n);
fo (i, 1, n) scanf("%d", &a[i]);
int sq = (int) sqrt(n);
fo (i, 1, n) block[i] = i / sq;
scanf("%d", &m);
fo (i, 1, m)
{
scanf("%d %d", &q[i].x, &q[i].y);
q[i].id = i;
}
std::sort(q + 1, q + m + 1);
l = 1; r = 0; ans = 0;
fo (i, 1, m)
{
while (l > q[i].x) add(a[--l], 1);
while (r < q[i].y) add(a[++r], 1);
while (l < q[i].x) add(a[l++], 0);
while (r > q[i].y) add(a[r--], 0);
anss[q[i].id] = ans;
}
printf("Case #%d:\n", qwq);
fo (i, 1, m)
printf("%lld\n", anss[i]);
}
return 0;
}


1. 在求最大公约数 $n$对答案的贡献的时候，我们可以将它转化为 $\sum_{d|n}\phi(d)$，这样就把贡献转移到了因子的上面。

2. 每个数的因子个数似乎也是很小的呢，$2\times 10^6$内的数平均下来才 $14$个因子。然而它的个数似乎还是找不到公式呢，反正最大的应该是小于 $sqrtn$的 (臆断)，可以现场打表检验呢反正也不麻烦。

（感觉排列随机一下更好