题意简述：给定序列\(a\)，\(q\)次询问在\([l1, r1],[l2, r2]\)中分别选择一个数，满足两个数互质的方案数。\(1\leq n\leq 2\times 10^4, 1\leq q\leq 2\times 10^5, 1\leq a_1, a_2, \ldots, a_n\leq 5\times 10^4\)，保证数据随机。
每个询问可以拆成四个前缀和容斥的形式，现在只需考虑在\([1, x], [1, y]\)中分别选数的方案数。
设\([1, x]\)中\(i\)出现了\(b_i\)次，\([1, y]\)中\(i\)出现了\(c_i\)次，则答案为：\(\sum_{i=1}^{\max a} \sum_{j=1}^{\max a} [\gcd(i, j)=1] b_i\cdot c_j\)。
经过一些简单的莫比乌斯反演，得到\(\sum_{d=1}^{\max a}\mu(d)\sum_{i=1}^{\lfloor\frac{\max a}{d}\rfloor}b_i \sum_{j=1}^{\lfloor\frac{\max a}{d}\rfloor}c_j\)。
设\(E_x=\sum_{i=1}^{\lfloor\frac{\max a}{i}\rfloor} b_{i\cdot x}, F_x=\sum_{i=1}^{\lfloor\frac{\max a}{i}\rfloor} c_{i\cdot x}\)，则答案为\(\sum_{d=1}^{\max a}\mu(d)\cdot E_d\cdot F_d\)。
注意到\(b_i\)或\(c_i\)增加时，会产生变化的\(E_j\)或\(F_j\)满足\(j|i\)，即只有\(\omega(i)\)项会改变，莫队维护即可。对于随机数据，复杂度期望为\(O(n\sqrt{m}\log \max a)\)，注意调整块大小为\(\frac{n}{\sqrt{m}}\)。
#include <bits/stdc++.h>
#define R register
#define ll long long
using namespace std;
const int N = 21000, M = 810000, A = 51000;
int n, q, a[N], e1[A], e2[A], bel[N], mu[A], ispr[A], prime[N], nopr, num;
ll ans[M], sum;
struct node {
int x, y, ind;
node (int x = 0, int y = 0, int i = 0) : x(x), y(y), ind(i) {}
inline bool operator < (const node &a) const {
if (bel[x] == bel[a.x])
return (bel[x] & 1) ? y < a.y : y > a.y;
return bel[x] < bel[a.x];
}
}que[M];
vector<int> fac[A];
template <class T> inline void read(T &x) {
x = 0;
char ch = getchar(), w = 0;
while (!isdigit(ch)) w = (ch == '-'), ch = getchar();
while (isdigit(ch)) x = (x << 1) + (x << 3) + (ch ^ 48), ch = getchar();
x = w ? -x : x;
return;
}
void preCalc(int n) {
mu[1] = 1;
for (R int i = 2; i <= n; ++i) {
if (!ispr[i])
prime[++nopr] = i, mu[i] = -1;
for (R int j = 1, k; j <= nopr && (k = i * prime[j]) <= n; ++j) {
ispr[k] = 1;
if (i % prime[j] == 0) break;
mu[k] = mu[i] * -1;
}
}
for (R int i = 1; i <= n; ++i)
for (R int j = 1, k; (k = i * j) <= n; ++j)
fac[k].push_back(i);
return;
}
void add(int pos, int *mod, int *que) {
int c = a[pos];
for (R int i = 0, sz = fac[c].size(); i < sz; ++i) {
sum += mu[fac[c][i]] * que[fac[c][i]];
++mod[fac[c][i]];
}
return;
}
void del(int pos, int *mod, int *que) {
int c = a[pos];
for (R int i = 0, sz = fac[c].size(); i < sz; ++i) {
sum -= mu[fac[c][i]] * que[fac[c][i]];
--mod[fac[c][i]];
}
return;
}
int main() {
int l1, r1, l2, r2;
read(n);
int maxA = 0;
for (R int i = 1; i <= n; ++i)
read(a[i]), maxA = max(maxA, a[i]);
preCalc(maxA), read(q);
int B = max(1, (int) (n / sqrt(q)));
for (R int i = 0; i < n; ++i)
bel[i + 1] = i / B + 1;
for (R int i = 1; i <= q; ++i) {
read(l1), read(r1), read(l2), read(r2);
que[++num] = node(r1, r2, i);
que[++num] = node(r1, l2 - 1, -i);
que[++num] = node(l1 - 1, r2, -i);
que[++num] = node(l1 - 1, l2 - 1, i);
}
sort(que + 1, que + 1 + num);
int p1 = 0, p2 = 0;
for (R int i = 1; i <= num; ++i) {
while (p1 < que[i].x) add(++p1, e1, e2);
while (p1 > que[i].x) del(p1--, e1, e2);
while (p2 < que[i].y) add(++p2, e2, e1);
while (p2 > que[i].y) del(p2--, e2, e1);
if (que[i].ind > 0) ans[que[i].ind] += sum;
else ans[-que[i].ind] -= sum;
}
for (R int i = 1; i <= q; ++i)
printf("%lld\n", ans[i]);
return 0;
}