题意
Sol
首先不考虑\(a\)的限制
我们要求的是
\[\sum_{i = 1}^n \sum_{j = 1}^m \sigma(gcd(i, j))\]
用常规的套路可以化到这个形式
\[\sum_{d = 1}^n \sigma (d) \sum_{k = 1}^{\frac{n}{d}} \mu(k) \frac{n}{kd} \frac{m}{kd}\]
设\(kd = T\)
那么
\(\sum_{T = 1}^n \left\lfloor \frac{n}{T} \right\rfloor \left\lfloor \frac{m}{T} \right\rfloor \sum_{d \ | T} \sigma(d) \mu(\frac{T}{d})\)
然后按照套路筛出后面的卷积。
但是现在有\(a\)的限制了。。
那么可以直接离线之后树状数组维护一下。
时间复杂度:\(O(T\sqrt{n}logn)\)
约数个数和用埃氏筛两行就筛出来了
#include<bits/stdc++.h>
#define chmax(a, b) (a = (a > b ? a : b))
#define chmin(a, b) (a = (a < b ? a : b))
//#define int long long
using namespace std;
const int MAXN = 3e5 + 10;
inline int read() {
char c = getchar(); int x = 0, f = 1;
while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
return x * f;
}
int T, mx, mu[MAXN], phi[MAXN], prime[MAXN], tot, vis[MAXN], si[MAXN], out[MAXN], p[MAXN];
struct Query{
int N, M, a, id;
bool operator < (const Query &rhs) const {
return a < rhs.a;
}
}q[MAXN];
int comp(const Query &a, const Query &b) {
return a.id < b.id;
}
int comp2(const int &a, const int &b) {
return si[a] < si[b];
}
void Get(int N) {
for(int i = 1; i <= N; i++)
for(int j = i; j <= N; j += i) si[j] += i;
vis[1] = 1; mu[1] = 1;
for(int i = 2; i <= N; i++) {
if(!vis[i]) prime[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
vis[i * prime[j]] = 1;
if(!(i % prime[j])) {mu[i * prime[j]] = 0; break;}
else mu[i * prime[j]] = -mu[i];
}
}
}
#define lb(x) (x & (-x))
int Tr[MAXN];
void Add(int x, int val) {
while(x <= mx) Tr[x] += val, x += lb(x);
}
int Sum(int x) {
int ans = 0;
while(x) ans += Tr[x], x -= lb(x);
return ans;
}
int solve(int n, int m) {
int rt = 0, las = 0, now;
for(int i = 1, nxt; i <= n; i = nxt + 1) {
nxt = min(m / (m / i), n / (n / i));
now = Sum(nxt);
rt += (n / i) * (m / i) * (now - las);
las = now;
}
return rt;
}
signed main() {
T = read();
for(int i = 1; i <= T; i++) {
q[i].N = read(), q[i].M = read(), q[i].a = read(), q[i].id = i;
if(q[i].N > q[i].M) swap(q[i].N, q[i].M);;
chmax(mx, q[i].N);
}
Get(mx);
for(int i = 1; i <= mx; i++) p[i] = i;
sort(p + 1, p + mx + 1, comp2);
sort(q + 1, q + T + 1);
for(int t = 1, d = 1; t <= T; t++) {
for(; d <= mx && si[p[d]] <= q[t].a; d++) {
for(int k = p[d]; k <= mx; k += p[d]) {
if(!mu[k / p[d]]) continue;
Add(k, si[p[d]] * mu[k / p[d]]);
}
}
out[q[t].id] = solve(q[t].N, q[t].M);
}
for(int i = 1; i <= T; i++)
printf("%d\n", out[i] < 0 ? out[i] + 2147483648 : out[i]);
//printf("%d\n", out[i] <);
return 0;
}