思路
考虑没有 $a$ 限制怎么办。
\[\text{Let }n = \max(n, m) \\ \begin{align*} \sum_{i=1}^n \sum_{j=1}^m \sigma(\gcd(i, j)) &= \sum_{g=1}^n \sigma(g)\sum_{g \mid d} \mu(\frac{d}{g}) \left\lfloor\frac{n}{d}\right\rfloor \left\lfloor\frac{m}{d}\right\rfloor \\ &= \sum_{d=1}^n \sum_{g | d} \sigma(g)\mu(\frac{d}{g}) \left\lfloor\frac{n}{d}\right\rfloor \left\lfloor\frac{m}{d}\right\rfloor \\ &= \sum_{d=1}^n \left\lfloor\frac{n}{d}\right\rfloor \left\lfloor\frac{m}{d}\right\rfloor \sum_{g | d} \sigma(g) \mu(\frac{d}{g})\\ \end{align*}\]其中 $\sum\limits_{g \mid d} \sigma(g) \mu(\frac{d}{g})$ 可以先线性筛出 $\sigma, \mu$ ，然后通过枚举倍数以 $\frac{n}{1} + \frac{n}{2} + \dots + \frac{n}{n} = O(n \lg n)$ 的复杂度内计算这个和式。
再考虑加入 $a$ 限制后如何处理。我们需要动态维护上面的和式。离线处理每个询问，按照 $a$ 递增排序，每次把新的符合条件的 $\sigma$ 加入树状数组即可。
复杂度 $O(n \lg n + n \sqrt{n})$.
代码
#include <bits/stdc++.h>
#define fst first
#define snd second
using namespace std;
struct Query {
int id, n, m, a;
bool operator<(const Query &rhs) const { return a < rhs.a; }
};
const int MAXN = 1e5, MAXQ = 2e4;
int q, primecnt, primelst[MAXN+10], notprime[MAXN+10], mu[MAXN+10],
sigma[MAXN+10], smin[MAXN+10];
int pos = -1, g[MAXN+10], ans[MAXN+10]; // BIT
vector<pair<int, int> > sigmas;
Query qry[MAXQ+10];
inline int lowbit(int x) { return x & (-x); }
inline int sum(int p) {
int ret = 0;
while(p) {
ret += g[p];
p -= lowbit(p);
}
return ret;
}
inline void add(int p, int val) {
while(p <= MAXN) {
g[p] += val;
p += lowbit(p);
}
}
inline int readint() {
int f=1, r=0; char c=getchar();
while(!isdigit(c)) { if(c=='-')f=-1; c=getchar(); }
while(isdigit(c)) { r=r*10+c-'0'; c=getchar(); }
return f*r;
}
void euler_sieve() {
notprime[1] = true; mu[1] = 1;
smin[1] = sigma[1] = 1;
for(int i = 2; i <= MAXN; i++) {
if(!notprime[i]) {
primelst[++primecnt] = i;
mu[i] = -1;
smin[i] = sigma[i] = i+1;
}
for(int j = 1; j <= primecnt; j++) {
if(i * primelst[j] > MAXN) break;
notprime[i * primelst[j]] = true;
if(i % primelst[j] == 0) {
mu[i * primelst[j]] = 0;
smin[i * primelst[j]] = smin[i] * primelst[j] + 1;
sigma[i * primelst[j]] = sigma[i] / smin[i] * smin[i * primelst[j]];
break;
} else {
mu[i * primelst[j]] = -mu[i];
smin[i * primelst[j]] = smin[primelst[j]];
sigma[i * primelst[j]] = sigma[i] * smin[i * primelst[j]];
}
}
}
for(int i = 1; i <= MAXN; i++)
sigmas.push_back(make_pair(sigma[i], i));
sort(sigmas.begin(), sigmas.end());
}
int calc(int n, int m, int a) {
if(n > m) swap(n, m);
while(pos + 1 < MAXN && sigmas[pos + 1].fst <= a) {
++pos;
for(int j = 1; sigmas[pos].snd * j <= MAXN; j++) {
add(sigmas[pos].snd * j, sigmas[pos].fst * mu[j]);
}
}
int ret = 0, last = 0;
for(int i = 1; i <= n; i = last + 1) {
last = min(n/(n/i), m/(m/i));
ret += (n/i) * (m/i) * (sum(last)-sum(i-1));
}
return ret;
}
int main() {
euler_sieve();
q = readint();
for(int i = 1; i <= q; i++) {
qry[i].id = i;
qry[i].n = readint(), qry[i].m = readint(), qry[i].a = readint();
}
sort(qry + 1, qry + q + 1);
for(int i = 1; i <= q; i++)
ans[qry[i].id] = calc(qry[i].n, qry[i].m, qry[i].a) & 0x7fffffff;
for(int i = 1; i <= q; i++) printf("%d\n", ans[i]);
return 0;
}
